OGS
RichardsComponentTransportFEM-impl.h
Go to the documentation of this file.
1
16
17namespace ProcessLib
18{
19namespace RichardsComponentTransport
20{
21template <typename ShapeFunction, int GlobalDim>
23 MeshLib::Element const& element,
24 std::size_t const local_matrix_size,
25 NumLib::GenericIntegrationMethod const& integration_method,
26 bool is_axially_symmetric,
27 RichardsComponentTransportProcessData const& process_data,
28 ProcessVariable const& transport_process_variable)
29 : _element(element),
30 _process_data(process_data),
31 _integration_method(integration_method),
32 _transport_process_variable(transport_process_variable)
33{
34 // This assertion is valid only if all nodal d.o.f. use the same shape
35 // matrices.
36 assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
37 (void)local_matrix_size;
38
39 unsigned const n_integration_points =
41 _ip_data.reserve(n_integration_points);
42
43 auto const shape_matrices =
45 element, is_axially_symmetric, _integration_method);
46
47 for (unsigned ip = 0; ip < n_integration_points; ip++)
48 {
49 auto const& sm = shape_matrices[ip];
50 double const integration_factor = sm.integralMeasure * sm.detJ;
51 _ip_data.emplace_back(
52 sm.N, sm.dNdx,
54 integration_factor,
55 sm.N.transpose() * sm.N * integration_factor *
57 }
58}
59
60template <typename ShapeFunction, int GlobalDim>
62 double const t, double const dt, std::vector<double> const& local_x,
63 std::vector<double> const& /*local_x_prev*/,
64 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
65 std::vector<double>& local_b_data)
66{
67 auto const local_matrix_size = local_x.size();
68 // This assertion is valid only if all nodal d.o.f. use the same shape
69 // matrices.
70 assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
71
73 local_M_data, local_matrix_size, local_matrix_size);
75 local_K_data, local_matrix_size, local_matrix_size);
77 local_b_data, local_matrix_size);
78
79 unsigned const n_integration_points =
80 _integration_method.getNumberOfPoints();
81
82 auto p_nodal_values = Eigen::Map<const NodalVectorType>(
83 &local_x[pressure_index], pressure_size);
84
85 auto const& b = _process_data.specific_body_force;
86
88
89 GlobalDimMatrixType const& I(
90 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
91
92 // Get material properties
93 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
94 auto const& phase = medium.phase("AqueousLiquid");
95 auto const& component =
96 phase.component(_transport_process_variable.getName());
97
98 auto KCC = local_K.template block<concentration_size, concentration_size>(
99 concentration_index, concentration_index);
100 auto MCC = local_M.template block<concentration_size, concentration_size>(
101 concentration_index, concentration_index);
102 auto Kpp = local_K.template block<pressure_size, pressure_size>(
103 pressure_index, pressure_index);
104 auto Mpp = local_M.template block<pressure_size, pressure_size>(
105 pressure_index, pressure_index);
106 auto Bp = local_b.template block<pressure_size, 1>(pressure_index, 0);
107
108 for (std::size_t ip(0); ip < n_integration_points; ++ip)
109 {
110 auto const& ip_data = _ip_data[ip];
111 auto const& N = ip_data.N;
112 auto const& dNdx = ip_data.dNdx;
113 auto const& w = ip_data.integration_weight;
114
116 std::nullopt, _element.getID(),
119 _element, N))};
120
121 double C_int_pt = 0.0;
122 double p_int_pt = 0.0;
123 // Order matters: First C, then p!
124 NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
125
126 vars.capillary_pressure = -p_int_pt;
128 .template value<double>(vars, pos, t, dt);
129
130 double const dSw_dpc =
132 .template dValue<double>(
134 pos, t, dt);
135
136 vars.concentration = C_int_pt;
137 vars.liquid_phase_pressure = p_int_pt;
138 // setting pG to 1 atm
139 // TODO : rewrite equations s.t. p_L = pG-p_cap
140 vars.gas_phase_pressure = 1.0e5;
141
142 // \todo the argument to getValue() has to be changed for non
143 // constant storage model
144 auto const specific_storage =
146 .template value<double>(vars, pos, t, dt);
147 // \todo the first argument has to be changed for non constant
148 // porosity model
149 auto const porosity =
151 .template value<double>(vars, pos, t, dt);
152
153 auto const retardation_factor =
155 .template value<double>(vars, pos, t, dt);
156
157 auto const solute_dispersivity_transverse =
159 .template value<double>(vars, pos, t, dt);
160 auto const solute_dispersivity_longitudinal =
162 .template value<double>(vars, pos, t, dt);
163
164 // Use the fluid density model to compute the density
165 auto const density = phase[MaterialPropertyLib::PropertyType::density]
166 .template value<double>(vars, pos, t, dt);
167 vars.density = density;
168
169 auto const decay_rate =
171 .template value<double>(vars, pos, t, dt);
172 auto const pore_diffusion_coefficient =
175 .value(vars, pos, t, dt));
176
179 vars, pos, t, dt));
180 vars.liquid_saturation = Sw;
181 auto const k_rel =
183 .template value<double>(vars, pos, t, dt);
184 // Use the viscosity model to compute the viscosity
186 .template value<double>(vars, pos, t, dt);
187 auto const K_times_k_rel_over_mu = K * (k_rel / mu);
188
189 GlobalDimVectorType const velocity =
190 _process_data.has_gravity
191 ? GlobalDimVectorType(-K_times_k_rel_over_mu *
192 (dNdx * p_nodal_values - density * b))
193 : GlobalDimVectorType(-K_times_k_rel_over_mu * dNdx *
194 p_nodal_values);
195
196 double const velocity_magnitude = velocity.norm();
197 GlobalDimMatrixType const hydrodynamic_dispersion =
198 velocity_magnitude != 0.0
200 porosity * pore_diffusion_coefficient +
201 solute_dispersivity_transverse * velocity_magnitude * I +
202 (solute_dispersivity_longitudinal -
203 solute_dispersivity_transverse) /
204 velocity_magnitude * velocity * velocity.transpose())
205 : GlobalDimMatrixType(porosity * pore_diffusion_coefficient +
206 solute_dispersivity_transverse *
207 velocity_magnitude * I);
208
209 // matrix assembly
210 KCC.noalias() +=
211 (dNdx.transpose() * hydrodynamic_dispersion * dNdx +
212 N.transpose() * velocity.transpose() * dNdx +
213 N.transpose() * decay_rate * porosity * retardation_factor * N) *
214 w;
215 MCC.noalias() += w * N.transpose() * porosity * retardation_factor * N;
216 Kpp.noalias() += w * dNdx.transpose() * K_times_k_rel_over_mu * dNdx;
217 // \TODO Extend to pressure dependent density.
218 double const drhow_dp(0.0);
219 Mpp.noalias() += (specific_storage * Sw + porosity * Sw * drhow_dp -
220 porosity * dSw_dpc) *
221 ip_data.mass_operator;
222
223 if (_process_data.has_gravity)
224 {
225 Bp += w * density * dNdx.transpose() * K_times_k_rel_over_mu * b;
226 }
227 /* with Oberbeck-Boussing assumption density difference only exists
228 * in buoyancy effects */
229 }
230}
231
232template <typename ShapeFunction, int GlobalDim>
233std::vector<double> const&
235 const double t,
236 std::vector<GlobalVector*> const& x,
237 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
238 std::vector<double>& cache) const
239{
240 unsigned const n_integration_points =
241 _integration_method.getNumberOfPoints();
242
243 constexpr int process_id = 0; // monolithic scheme
244 auto const indices =
245 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
246 assert(!indices.empty());
247 auto const local_x = x[process_id]->get(indices);
248
249 cache.clear();
250 auto cache_mat = MathLib::createZeroedMatrix<
251 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
252 cache, GlobalDim, n_integration_points);
253
255
256 // Get material properties
257 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
258 auto const& phase = medium.phase("AqueousLiquid");
259
260 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
261 &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
262
263 for (unsigned ip = 0; ip < n_integration_points; ++ip)
264 {
265 auto const& ip_data = _ip_data[ip];
266 auto const& N = ip_data.N;
267 auto const& dNdx = ip_data.dNdx;
268
270 std::nullopt, _element.getID(),
273 _element, N))};
274
275 auto const dt = std::numeric_limits<double>::quiet_NaN();
278 vars, pos, t, dt));
280 .template value<double>(vars, pos, t, dt);
281
282 double C_int_pt = 0.0;
283 double p_int_pt = 0.0;
284 NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
285
286 // saturation
287 vars.capillary_pressure = -p_int_pt;
289 .template value<double>(vars, pos, t, dt);
290
291 vars.liquid_saturation = Sw;
292 auto const k_rel =
294 .template value<double>(vars, pos, t, dt);
295
296 cache_mat.col(ip).noalias() = -dNdx * p_nodal_values;
297 if (_process_data.has_gravity)
298 {
299 vars.concentration = C_int_pt;
300 vars.liquid_phase_pressure = p_int_pt;
301 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
302 .template value<double>(vars, pos, t, dt);
303 auto const b = _process_data.specific_body_force;
304 // here it is assumed that the vector b is directed 'downwards'
305 cache_mat.col(ip).noalias() += rho_w * b;
306 }
307 cache_mat.col(ip).noalias() = k_rel / mu * (K * cache_mat.col(ip));
308 }
309 return cache;
310}
311
312template <typename ShapeFunction, int GlobalDim>
313Eigen::Map<const Eigen::RowVectorXd>
315 const unsigned integration_point) const
316{
317 auto const& N = _ip_data[integration_point].N;
318
319 // assumes N is stored contiguously in memory
320 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
321}
322
323template <typename ShapeFunction, int GlobalDim>
324std::vector<double> const&
326 const double t,
327 std::vector<GlobalVector*> const& x,
328 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
329 std::vector<double>& cache) const
330{
332
333 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
334
335 unsigned const n_integration_points =
336 _integration_method.getNumberOfPoints();
337
338 constexpr int process_id = 0; // monolithic scheme
339 auto const indices =
340 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
341 assert(!indices.empty());
342 auto const local_x = x[process_id]->get(indices);
343
344 cache.clear();
345 auto cache_vec = MathLib::createZeroedVector<
346 Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>>(
347 cache, n_integration_points);
348
349 for (unsigned ip = 0; ip < n_integration_points; ++ip)
350 {
351 auto const& ip_data = _ip_data[ip];
352 auto const& N = ip_data.N;
353
355 std::nullopt, _element.getID(),
358 _element, N))};
359
360 double C_int_pt = 0.0;
361 double p_int_pt = 0.0;
362 NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
363
364 // saturation
365 vars.capillary_pressure = -p_int_pt;
366 auto const dt = std::numeric_limits<double>::quiet_NaN();
368 .template value<double>(vars, pos, t, dt);
369 cache_vec[ip] = Sw;
370 }
371
372 return cache;
373}
374
375} // namespace RichardsComponentTransport
376} // namespace ProcessLib
constexpr double getWeight() const
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType > > > _ip_data
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool is_axially_symmetric, RichardsComponentTransportProcessData const &process_data, ProcessVariable const &transport_process_variable)
std::vector< double > const & getIntPtSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ retardation_factor
specify retardation factor used in component transport process.
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)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
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::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)