OGS
HeatTransportBHELocalAssemblerSoil-impl.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <vector>
14
24#include "SecondaryData.h"
25
26namespace ProcessLib
27{
28namespace HeatTransportBHE
29{
30template <typename ShapeFunction>
33 MeshLib::Element const& e,
34 NumLib::GenericIntegrationMethod const& integration_method,
35 bool const is_axially_symmetric,
36 HeatTransportBHEProcessData& process_data)
37 : _process_data(process_data),
38 _integration_method(integration_method),
39 _element_id(e.getID())
40{
41 unsigned const n_integration_points =
43
44 _ip_data.reserve(n_integration_points);
45 _secondary_data.N.resize(n_integration_points);
46
49 3 /* GlobalDim */>(e, is_axially_symmetric,
51
53 x_position.setElementID(_element_id);
54
55 // ip data initialization
56 for (unsigned ip = 0; ip < n_integration_points; ip++)
57 {
58 // create the class IntegrationPointDataBHE in place
59 auto const& sm = _shape_matrices[ip];
60 double const w = _integration_method.getWeightedPoint(ip).getWeight() *
61 sm.integralMeasure * sm.detJ;
62 _ip_data.push_back({sm.N, sm.dNdx, w});
63
64 _secondary_data.N[ip] = sm.N;
65 }
66}
67
68template <typename ShapeFunction>
70 double const t, double const dt, std::vector<double> const& local_x,
71 std::vector<double> const& /*local_x_prev*/,
72 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
73 std::vector<double>& /*local_b_data*/)
74{
75 assert(local_x.size() == ShapeFunction::NPOINTS);
76 (void)local_x; // Avoid unused arg warning.
77
79 local_M_data, ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
81 local_K_data, ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
82
84 pos.setElementID(_element_id);
85
86 auto const& medium = *_process_data.media_map.getMedium(_element_id);
87 auto const& solid_phase = medium.phase("Solid");
88 auto const& liquid_phase = medium.phase("AqueousLiquid");
89
91
92 unsigned const n_integration_points =
93 _integration_method.getNumberOfPoints();
94
95 for (unsigned ip = 0; ip < n_integration_points; ip++)
96 {
97 auto& ip_data = _ip_data[ip];
98 auto const& N = ip_data.N;
99 auto const& dNdx = ip_data.dNdx;
100 auto const& w = ip_data.integration_weight;
101
102 double T_int_pt = 0.0;
103 NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt);
104
105 vars.temperature = T_int_pt;
106
107 // for now only using the solid and liquid phase parameters
108 auto const density_s =
110 .template value<double>(vars, pos, t, dt);
111
112 auto const heat_capacity_s =
113 solid_phase
114 .property(
116 .template value<double>(vars, pos, t, dt);
117
118 auto const density_f =
119 liquid_phase.property(MaterialPropertyLib::PropertyType::density)
120 .template value<double>(vars, pos, t, dt);
121
122 auto const heat_capacity_f =
123 liquid_phase
124 .property(
126 .template value<double>(vars, pos, t, dt);
127
128 auto const porosity =
130 .template value<double>(vars, pos, t, dt);
131
132 auto const velocity =
133 liquid_phase
135 .template value<Eigen::Vector3d>(vars, pos, t, dt);
136
137 // calculate the hydrodynamic thermodispersion tensor
138 auto const thermal_conductivity =
140 medium
141 .property(
143 .value(vars, pos, t, dt));
144
145 auto thermal_conductivity_dispersivity = thermal_conductivity;
146
147 double const velocity_magnitude = velocity.norm();
148
149 if (velocity_magnitude >= std::numeric_limits<double>::epsilon())
150 {
151 auto const thermal_dispersivity_longitudinal =
152 medium
154 thermal_longitudinal_dispersivity)
155 .template value<double>();
156 auto const thermal_dispersivity_transversal =
157 medium
159 thermal_transversal_dispersivity)
160 .template value<double>();
161
162 auto const thermal_dispersivity =
163 density_f * heat_capacity_f *
164 (thermal_dispersivity_transversal * velocity_magnitude *
165 Eigen::Matrix3d::Identity() +
166 (thermal_dispersivity_longitudinal -
167 thermal_dispersivity_transversal) /
168 velocity_magnitude * velocity * velocity.transpose());
169 thermal_conductivity_dispersivity += thermal_dispersivity;
170 }
171
172 // assemble Conductance matrix
173 local_K.noalias() +=
174 (dNdx.transpose() * thermal_conductivity_dispersivity * dNdx +
175 N.transpose() * velocity.transpose() * dNdx * density_f *
176 heat_capacity_f) *
177 w;
178
179 // assemble Mass matrix
180 local_M.noalias() += N.transpose() * N * w *
181 (density_s * heat_capacity_s * (1 - porosity) +
182 density_f * heat_capacity_f * porosity);
183 }
184
185 if (_process_data._mass_lumping)
186 {
187 // only mass lumping at the BHE connected soil elements
188 if (_process_data.mass_lumping_soil_elements[_element_id])
189 {
190 local_M = local_M.colwise().sum().eval().asDiagonal();
191 }
192 }
193
194 // debugging
195 // std::string sep = "\n----------------------------------------\n";
196 // Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
197 // std::cout << local_K.format(CleanFmt) << sep;
198 // std::cout << local_M.format(CleanFmt) << sep;
199}
200
201template <typename ShapeFunction>
203 double const t, double const dt, std::vector<double> const& local_x,
204 std::vector<double> const& local_x_prev,
205 std::vector<double>& local_rhs_data, std::vector<double>& local_Jac_data)
206{
207 assert(local_x.size() == ShapeFunction::NPOINTS);
208 auto const local_matrix_size = local_x.size();
209 // initialize x and x_prev
210 auto x =
211 Eigen::Map<NodalVectorType const>(local_x.data(), local_matrix_size);
212 auto x_prev = Eigen::Map<NodalVectorType const>(local_x_prev.data(),
213 local_matrix_size);
214 // initialize local_Jac and local_rhs
216 local_Jac_data, local_matrix_size, local_matrix_size);
218 local_rhs_data, local_matrix_size);
219
220 std::vector<double> local_M_data;
221 std::vector<double> local_K_data;
222 assemble(t, dt, local_x, local_x_prev, local_M_data, local_K_data,
223 local_rhs_data /*not going to be used*/);
224
225 // convert to matrix
227 local_M_data, local_matrix_size, local_matrix_size);
229 local_K_data, local_matrix_size, local_matrix_size);
230
231 // Jac matrix and rhs vector operation
232 local_Jac.noalias() += local_K + local_M / dt;
233 local_rhs.noalias() -= local_K * x + local_M * (x - x_prev) / dt;
234
235 local_M.setZero();
236 local_K.setZero();
237}
238
239} // namespace HeatTransportBHE
240} // namespace ProcessLib
constexpr double getWeight() const
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setElementID(std::size_t element_id)
std::vector< ShapeMatrices, Eigen::aligned_allocator< ShapeMatrices > > _shape_matrices
void assemble(double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
void assembleWithJacobian(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
std::vector< IntegrationPointDataSoil< NodalRowVectorType, GlobalDimNodalMatrixType >, Eigen::aligned_allocator< IntegrationPointDataSoil< NodalRowVectorType, GlobalDimNodalMatrixType > > > _ip_data
HeatTransportBHELocalAssemblerSoil(HeatTransportBHELocalAssemblerSoil const &)=delete
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N