OGS
ThermoRichardsFlowFEM-impl.h
Go to the documentation of this file.
1
10#pragma once
11
12#include <cassert>
13
28
29namespace ProcessLib
30{
31namespace ThermoRichardsFlow
32{
33template <typename ShapeFunction, int GlobalDim>
36 MeshLib::Element const& e,
37 std::size_t const /*local_matrix_size*/,
38 NumLib::GenericIntegrationMethod const& integration_method,
39 bool const is_axially_symmetric,
41 : _process_data(process_data),
42 _integration_method(integration_method),
43 _element(e),
44 _is_axially_symmetric(is_axially_symmetric)
45{
46 unsigned const n_integration_points =
48
49 _ip_data.reserve(n_integration_points);
50
51 auto const shape_matrices =
52 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType, GlobalDim>(
53 e, is_axially_symmetric, _integration_method);
54
55 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
56
57 for (unsigned ip = 0; ip < n_integration_points; ip++)
58 {
59 auto const& sm = shape_matrices[ip];
60 _ip_data.emplace_back();
61 auto& ip_data = _ip_data[ip];
62 _ip_data[ip].integration_weight =
64 sm.integralMeasure * sm.detJ;
65
66 ip_data.N = sm.N;
67 ip_data.dNdx = sm.dNdx;
68
69 ParameterLib::SpatialPosition const x_position{
70 std::nullopt, _element.getID(), ip,
73 _element, sm.N))};
74 // Initial porosity. Could be read from integration point data or mesh.
75 ip_data.porosity = medium[MPL::porosity].template initialValue<double>(
76 x_position,
77 std::numeric_limits<double>::quiet_NaN() /* t independent */);
78 }
79}
80
81template <typename ShapeFunction, int GlobalDim>
83 setIPDataInitialConditions(std::string_view const name,
84 double const* values,
85 int const integration_order)
86{
87 if (integration_order !=
88 static_cast<int>(_integration_method.getIntegrationOrder()))
89 {
91 "Setting integration point initial conditions; The integration "
92 "order of the local assembler for element {:d} is different "
93 "from the integration order in the initial condition.",
94 _element.getID());
95 }
96
97 if (name == "saturation")
98 {
99 return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
100 &IpData::saturation);
101 }
102 if (name == "porosity")
103 {
104 return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
105 &IpData::porosity);
106 }
107 return 0;
108}
109
110template <typename ShapeFunction, int GlobalDim>
112 setInitialConditionsConcrete(Eigen::VectorXd const local_x,
113 double const t,
114 int const /*process_id*/)
115{
116 assert(local_x.size() == temperature_size + pressure_size);
117
118 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
119
120 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
121 MPL::VariableArray variables;
122
123 unsigned const n_integration_points =
124 _integration_method.getNumberOfPoints();
125 for (unsigned ip = 0; ip < n_integration_points; ip++)
126 {
127 auto const& N = _ip_data[ip].N;
128
129 ParameterLib::SpatialPosition const x_position{
130 std::nullopt, _element.getID(), ip,
133 _element, N))};
134
135 double p_cap_ip;
136 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
137
138 variables.capillary_pressure = p_cap_ip;
139 variables.liquid_phase_pressure = -p_cap_ip;
140 // setting pG to 1 atm
141 // TODO : rewrite equations s.t. p_L = pG-p_cap
142 variables.gas_phase_pressure = 1.0e5;
143
144 // Note: temperature dependent saturation model is not considered so
145 // far.
146 _ip_data[ip].saturation_prev =
147 medium[MPL::PropertyType::saturation].template value<double>(
148 variables, x_position, t,
149 std::numeric_limits<double>::quiet_NaN());
150 }
151}
152
153template <typename ShapeFunction, int GlobalDim>
155 assembleWithJacobian(double const t, double const dt,
156 std::vector<double> const& local_x,
157 std::vector<double> const& local_x_prev,
158 std::vector<double>& /*local_M_data*/,
159 std::vector<double>& /*local_K_data*/,
160 std::vector<double>& local_rhs_data,
161 std::vector<double>& local_Jac_data)
162{
163 auto const local_matrix_dim = pressure_size + temperature_size;
164 assert(local_x.size() == local_matrix_dim);
165
166 auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
167 temperature_size> const>(local_x.data() + temperature_index,
168 temperature_size);
169 auto const p_L = Eigen::Map<
170 typename ShapeMatricesType::template VectorType<pressure_size> const>(
171 local_x.data() + pressure_index, pressure_size);
172
173 auto const T_prev =
174 Eigen::Map<typename ShapeMatricesType::template VectorType<
175 temperature_size> const>(local_x_prev.data() + temperature_index,
176 temperature_size);
177 auto const p_L_prev = Eigen::Map<
178 typename ShapeMatricesType::template VectorType<pressure_size> const>(
179 local_x_prev.data() + pressure_index, pressure_size);
180
181 auto local_Jac = MathLib::createZeroedMatrix<
182 typename ShapeMatricesType::template MatrixType<local_matrix_dim,
183 local_matrix_dim>>(
184 local_Jac_data, local_matrix_dim, local_matrix_dim);
185
186 auto local_rhs = MathLib::createZeroedVector<
187 typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
188 local_rhs_data, local_matrix_dim);
189
191 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
192 temperature_size);
194 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
195 pressure_size);
197 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
198 temperature_size);
200 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
201 pressure_size);
202 typename ShapeMatricesType::NodalMatrixType dK_TT_dp =
203 ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
204 pressure_size);
206 ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
207 temperature_size);
208 typename ShapeMatricesType::NodalMatrixType laplace_p =
209 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
210 typename ShapeMatricesType::NodalMatrixType laplace_T =
211 ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
212 temperature_size);
213 typename ShapeMatricesType::NodalMatrixType storage_p_a_p =
214 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
215
216 typename ShapeMatricesType::NodalMatrixType storage_p_a_S_Jpp =
217 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
218
219 typename ShapeMatricesType::NodalMatrixType storage_p_a_S =
220 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
221
222 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
223 auto const& liquid_phase = medium.phase("AqueousLiquid");
224 auto const& solid_phase = medium.phase("Solid");
225 MPL::Phase const* gas_phase =
226 medium.hasPhase("Gas") ? &medium.phase("Gas") : nullptr;
227 MPL::VariableArray variables;
228 MPL::VariableArray variables_prev;
229
230 unsigned const n_integration_points =
231 _integration_method.getNumberOfPoints();
232 for (unsigned ip = 0; ip < n_integration_points; ip++)
233 {
234 auto const& w = _ip_data[ip].integration_weight;
235
236 auto const& N = _ip_data[ip].N;
237 auto const& dNdx = _ip_data[ip].dNdx;
238
239 ParameterLib::SpatialPosition const x_position{
240 std::nullopt, _element.getID(), ip,
243 _element, N))};
244
245 double T_ip;
247
248 double p_cap_ip;
249 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
250
251 double p_cap_prev_ip;
252 NumLib::shapeFunctionInterpolate(-p_L_prev, N, p_cap_prev_ip);
253
254 variables.capillary_pressure = p_cap_ip;
255 variables.liquid_phase_pressure = -p_cap_ip;
256 // setting pG to 1 atm
257 // TODO : rewrite equations s.t. p_L = pG-p_cap
258 variables.gas_phase_pressure = 1.0e5;
259 variables.temperature = T_ip;
260
261 auto& S_L = _ip_data[ip].saturation;
262 auto const S_L_prev = _ip_data[ip].saturation_prev;
263 auto const alpha =
264 medium[MPL::PropertyType::biot_coefficient].template value<double>(
265 variables, x_position, t, dt);
266
267 auto& solid_elasticity = *_process_data.simplified_elasticity;
268 // TODO (buchwaldj)
269 // is bulk_modulus good name for bulk modulus of solid skeleton?
270 auto const beta_S =
271 solid_elasticity.bulkCompressibilityFromYoungsModulus(
272 solid_phase, variables, x_position, t, dt);
273 auto const beta_SR = (1 - alpha) * beta_S;
274 variables.grain_compressibility = beta_SR;
275
276 auto const rho_LR =
277 liquid_phase[MPL::PropertyType::density].template value<double>(
278 variables, x_position, t, dt);
279 variables.density = rho_LR;
280 auto const& b = _process_data.specific_body_force;
281
282 double const drho_LR_dp =
283 liquid_phase[MPL::PropertyType::density].template dValue<double>(
284 variables, MPL::Variable::liquid_phase_pressure, x_position, t,
285 dt);
286 auto const beta_LR = drho_LR_dp / rho_LR;
287
288 S_L = medium[MPL::PropertyType::saturation].template value<double>(
289 variables, x_position, t, dt);
290 variables.liquid_saturation = S_L;
291 variables_prev.liquid_saturation = S_L_prev;
292
293 // tangent derivative for Jacobian
294 double const dS_L_dp_cap =
295 medium[MPL::PropertyType::saturation].template dValue<double>(
296 variables, MPL::Variable::capillary_pressure, x_position, t,
297 dt);
298 // secant derivative from time discretization for storage
299 // use tangent, if secant is not available
300 double const DeltaS_L_Deltap_cap =
301 (p_cap_ip == p_cap_prev_ip)
302 ? dS_L_dp_cap
303 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
304
305 auto chi_S_L = S_L;
306 auto chi_S_L_prev = S_L_prev;
307 auto dchi_dS_L = 1.0;
308 if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
309 {
310 auto const chi = [&medium, x_position, t, dt](double const S_L)
311 {
312 MPL::VariableArray variables;
313 variables.liquid_saturation = S_L;
314 return medium[MPL::PropertyType::bishops_effective_stress]
315 .template value<double>(variables, x_position, t, dt);
316 };
317 chi_S_L = chi(S_L);
318 chi_S_L_prev = chi(S_L_prev);
319
320 dchi_dS_L = medium[MPL::PropertyType::bishops_effective_stress]
321 .template dValue<double>(
322 variables, MPL::Variable::liquid_saturation,
323 x_position, t, dt);
324 }
325 // TODO (buchwaldj)
326 // should solid_grain_pressure or effective_pore_pressure remain?
327 // double const p_FR = -chi_S_L * p_cap_ip;
328 // variables.solid_grain_pressure = p_FR;
329
330 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
331 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
332
333 auto& phi = _ip_data[ip].porosity;
334 { // Porosity update
335
336 variables_prev.porosity = _ip_data[ip].porosity_prev;
337 phi = medium[MPL::PropertyType::porosity].template value<double>(
338 variables, variables_prev, x_position, t, dt);
339 variables.porosity = phi;
340 }
341
342 if (alpha < phi)
343 {
344 OGS_FATAL(
345 "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
346 "porosity {} in element/integration point {}/{}.",
347 alpha, phi, _element.getID(), ip);
348 }
349
350 double const k_rel =
351 medium[MPL::PropertyType::relative_permeability]
352 .template value<double>(variables, x_position, t, dt);
353 auto const mu =
354 liquid_phase[MPL::PropertyType::viscosity].template value<double>(
355 variables, x_position, t, dt);
356
357 auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
358 medium[MPL::PropertyType::permeability].value(variables, x_position,
359 t, dt));
360
361 GlobalDimMatrixType const Ki_over_mu = K_intrinsic / mu;
362 GlobalDimMatrixType const rho_Ki_over_mu = rho_LR * Ki_over_mu;
363
364 auto const K_pT_thermal_osmosis =
365 (solid_phase.hasProperty(
367 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
368 solid_phase
369 [MPL::PropertyType::thermal_osmosis_coefficient]
370 .value(variables, x_position, t, dt))
371 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
372
373 // Consider anisotropic thermal expansion.
374 // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
375 // component.
376 Eigen::Matrix<double, 3, 3> const
377 solid_linear_thermal_expansion_coefficient =
379 solid_phase
381 .value(variables, x_position, t, dt));
382
383 auto const rho_SR =
384 solid_phase[MPL::PropertyType::density].template value<double>(
385 variables, x_position, t, dt);
386
387 //
388 // pressure equation, pressure part.
389 //
390 laplace_p.noalias() +=
391 dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
392 laplace_T.noalias() +=
393 dNdx.transpose() * rho_LR * K_pT_thermal_osmosis * dNdx * w;
394 const double alphaB_minus_phi = alpha - phi;
395 double const a0 = alphaB_minus_phi * beta_SR;
396 double const specific_storage_a_p =
397 S_L * (phi * beta_LR + S_L * a0 +
398 chi_S_L * alpha * alpha *
399 solid_elasticity.storageContribution(
400 solid_phase, variables, x_position, t, dt));
401 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
402
403 double const dspecific_storage_a_p_dp_cap =
404 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0 +
405 alpha * alpha *
406 solid_elasticity.storageContribution(
407 solid_phase, variables, x_position, t, dt) *
408 (chi_S_L + dchi_dS_L * S_L));
409 double const dspecific_storage_a_S_dp_cap =
410 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
411
412 storage_p_a_p.noalias() +=
413 N.transpose() * rho_LR * specific_storage_a_p * N * w;
414
415 storage_p_a_S.noalias() -= N.transpose() * rho_LR *
416 specific_storage_a_S * DeltaS_L_Deltap_cap *
417 N * w;
418
419 local_Jac
420 .template block<pressure_size, pressure_size>(pressure_index,
421 pressure_index)
422 .noalias() += N.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
423 rho_LR * dspecific_storage_a_p_dp_cap * N * w;
424
425 storage_p_a_S_Jpp.noalias() -=
426 N.transpose() * rho_LR *
427 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
428 specific_storage_a_S * dS_L_dp_cap) /
429 dt * N * w;
430
431 double const dk_rel_dS_L =
432 medium[MPL::PropertyType::relative_permeability]
433 .template dValue<double>(variables,
434 MPL::Variable::liquid_saturation,
435 x_position, t, dt);
436 GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
437 local_Jac
438 .template block<pressure_size, pressure_size>(pressure_index,
439 pressure_index)
440 .noalias() += dNdx.transpose() * rho_Ki_over_mu * grad_p_cap *
441 dk_rel_dS_L * dS_L_dp_cap * N * w;
442
443 local_Jac
444 .template block<pressure_size, pressure_size>(pressure_index,
445 pressure_index)
446 .noalias() += dNdx.transpose() * rho_LR * rho_Ki_over_mu * b *
447 dk_rel_dS_L * dS_L_dp_cap * N * w;
448
449 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
450 dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
451
452 //
453 // pressure equation, temperature part.
454 //
455 double const fluid_volumetric_thermal_expansion_coefficient =
456 MPL::getLiquidThermalExpansivity(liquid_phase, variables, rho_LR,
457 x_position, t, dt);
458 const double eff_thermal_expansion =
459 S_L * (alphaB_minus_phi *
460 solid_linear_thermal_expansion_coefficient.trace() +
461 phi * fluid_volumetric_thermal_expansion_coefficient +
462 alpha * solid_elasticity.thermalExpansivityContribution(
463 solid_linear_thermal_expansion_coefficient,
464 solid_phase, variables, x_position, t, dt));
465 M_pT.noalias() -=
466 N.transpose() * rho_LR * eff_thermal_expansion * N * w;
467
468 //
469 // temperature equation.
470 //
471 {
472 auto const specific_heat_capacity_fluid =
474 .template value<double>(variables, x_position, t, dt);
475
476 auto const specific_heat_capacity_solid =
477 solid_phase
479 .template value<double>(variables, x_position, t, dt);
480
481 M_TT.noalias() +=
482 w *
483 (rho_SR * specific_heat_capacity_solid * (1 - phi) +
484 (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
485 N.transpose() * N;
486
487 auto const thermal_conductivity =
488 MaterialPropertyLib::formEigenTensor<GlobalDim>(
490 thermal_conductivity]
491 .value(variables, x_position, t, dt));
492
493 GlobalDimVectorType const velocity_L = GlobalDimVectorType(
494 -Ki_over_mu * k_rel * (dNdx * p_L - rho_LR * b) -
495 K_pT_thermal_osmosis * dNdx * T);
496
497 K_TT.noalias() += (dNdx.transpose() * thermal_conductivity * dNdx +
498 N.transpose() * velocity_L.transpose() * dNdx *
499 rho_LR * specific_heat_capacity_fluid) *
500 w;
501
502 //
503 // temperature equation, pressure part
504 //
505 K_Tp.noalias() +=
506 dNdx.transpose() * T_ip * K_pT_thermal_osmosis * dNdx * w;
507 dK_TT_dp.noalias() -= rho_LR * specific_heat_capacity_fluid *
508 N.transpose() * (dNdx * T).transpose() *
509 k_rel * Ki_over_mu * dNdx * w;
510
511 dK_TT_dp.noalias() -= rho_LR * specific_heat_capacity_fluid *
512 N.transpose() * velocity_L.dot(dNdx * T) /
513 k_rel * dk_rel_dS_L * dS_L_dp_cap * N * w;
514 }
515 if (gas_phase && S_L < 1.0)
516 {
517 variables.density = rho_LR;
518
519 double const rho_wv =
520 gas_phase->property(MPL::PropertyType::density)
521 .template value<double>(variables, x_position, t, dt);
522
523 double const drho_wv_dT =
524 gas_phase->property(MPL::PropertyType::density)
525 .template dValue<double>(variables,
526 MPL::Variable::temperature,
527 x_position, t, dt);
528 double const drho_wv_dp =
529 gas_phase->property(MPL::PropertyType::density)
530 .template dValue<double>(
531 variables, MPL::Variable::liquid_phase_pressure,
532 x_position, t, dt);
533 auto const f_Tv =
534 gas_phase
535 ->property(
536 MPL::PropertyType::thermal_diffusion_enhancement_factor)
537 .template value<double>(variables, x_position, t, dt);
538
539 variables.porosity = phi;
540 auto const tortuosity =
541 medium.property(MPL::PropertyType::tortuosity)
542 .template value<double>(variables, x_position, t, dt);
543 double const D_v =
544 phi * (1.0 - S_L) * tortuosity *
545 gas_phase->property(MPL::PropertyType::diffusion)
546 .template value<double>(variables, x_position, t, dt);
547
548 double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
549 double const D_pv = D_v * drho_wv_dp;
550
551 GlobalDimVectorType const grad_T = dNdx * T;
552 GlobalDimVectorType const vapour_flux =
553 -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
554 double const specific_heat_capacity_vapour =
556 .template value<double>(variables, x_position, t, dt);
557
558 M_TT.noalias() +=
559 w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
560 N.transpose() * N;
561
562 K_TT.noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
563 specific_heat_capacity_vapour * w;
564
565 double const storage_coefficient_by_water_vapor =
566 phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
567
568 storage_p_a_p.noalias() +=
569 N.transpose() * storage_coefficient_by_water_vapor * N * w;
570
571 double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
572 M_pT.noalias() += N.transpose() * vapor_expansion_factor * N * w;
573
574 local_Jac
575 .template block<pressure_size, temperature_size>(
576 pressure_index, temperature_index)
577 .noalias() += dNdx.transpose() * f_Tv_D_Tv * dNdx * w;
578
579 local_rhs.template segment<pressure_size>(pressure_index)
580 .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
581
582 laplace_p.noalias() += dNdx.transpose() * D_pv * dNdx * w;
583
584 //
585 // Latent heat term
586 //
587 if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
588 {
589 double const factor = phi * (1 - S_L) / rho_LR;
590 // The volumetric latent heat of vaporization of liquid water
591 double const L0 =
592 gas_phase->property(MPL::PropertyType::specific_latent_heat)
593 .template value<double>(variables, x_position, t, dt) *
594 rho_LR;
595
596 double const drho_LR_dT =
597 liquid_phase.property(MPL::PropertyType::density)
598 .template dValue<double>(variables,
599 MPL::Variable::temperature,
600 x_position, t, dt);
601
602 double const rho_wv_over_rho_L = rho_wv / rho_LR;
603 M_TT.noalias() +=
604 factor * L0 *
605 (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
606 N.transpose() * N * w;
607
608 M_Tp.noalias() +=
609 (factor * L0 *
610 (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
611 L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
612 N.transpose() * N * w;
613
614 // temperature equation, temperature part
615 K_TT.noalias() +=
616 L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
617 // temperature equation, pressure part
618 K_Tp.noalias() +=
619 L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
620 }
621 }
622 }
623
624 if (_process_data.apply_mass_lumping)
625 {
626 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
627 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
628 storage_p_a_S_Jpp =
629 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
630 }
631
632 //
633 // -- Jacobian
634 //
635 // temperature equation.
636 local_Jac
637 .template block<temperature_size, temperature_size>(temperature_index,
638 temperature_index)
639 .noalias() += M_TT / dt + K_TT;
640 // temperature equation, pressure part
641 local_Jac
642 .template block<temperature_size, pressure_size>(temperature_index,
643 pressure_index)
644 .noalias() += K_Tp + dK_TT_dp;
645
646 // pressure equation, pressure part.
647 local_Jac
648 .template block<pressure_size, pressure_size>(pressure_index,
649 pressure_index)
650 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
651
652 // pressure equation, temperature part (contributed by thermal expansion).
653 local_Jac
654 .template block<pressure_size, temperature_size>(pressure_index,
655 temperature_index)
656 .noalias() += M_pT / dt + laplace_T;
657
658 //
659 // -- Residual
660 //
661 // temperature equation
662 local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
663 M_TT * (T - T_prev) / dt + K_TT * T;
664 local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
665 K_Tp * p_L;
666
667 // pressure equation
668 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
669 laplace_p * p_L + laplace_T * T +
670 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
671 M_pT * (T - T_prev) / dt;
672 if (gas_phase)
673 {
674 if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
675 {
676 // Jacobian: temperature equation, pressure part
677 local_Jac
678 .template block<temperature_size, pressure_size>(
679 temperature_index, pressure_index)
680 .noalias() += M_Tp / dt;
681 // RHS: temperature part
682 local_rhs.template segment<temperature_size>(temperature_index)
683 .noalias() -= M_Tp * (p_L - p_L_prev) / dt;
684 }
685 }
686}
687
688template <typename ShapeFunction, int GlobalDim>
690 double const t, double const dt, std::vector<double> const& local_x,
691 std::vector<double> const& local_x_prev, std::vector<double>& local_M_data,
692 std::vector<double>& local_K_data, std::vector<double>& local_rhs_data)
693{
694 auto const local_matrix_dim = pressure_size + temperature_size;
695 assert(local_x.size() == local_matrix_dim);
696
697 auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
698 temperature_size> const>(local_x.data() + temperature_index,
699 temperature_size);
700 auto const p_L = Eigen::Map<
701 typename ShapeMatricesType::template VectorType<pressure_size> const>(
702 local_x.data() + pressure_index, pressure_size);
703
704 auto const p_L_prev = Eigen::Map<
705 typename ShapeMatricesType::template VectorType<pressure_size> const>(
706 local_x_prev.data() + pressure_index, pressure_size);
707
708 auto local_K = MathLib::createZeroedMatrix<
709 typename ShapeMatricesType::template MatrixType<local_matrix_dim,
710 local_matrix_dim>>(
711 local_K_data, local_matrix_dim, local_matrix_dim);
712
713 auto local_M = MathLib::createZeroedMatrix<
714 typename ShapeMatricesType::template MatrixType<local_matrix_dim,
715 local_matrix_dim>>(
716 local_M_data, local_matrix_dim, local_matrix_dim);
717
718 auto local_rhs = MathLib::createZeroedVector<
719 typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
720 local_rhs_data, local_matrix_dim);
721
722 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
723 auto const& liquid_phase = medium.phase("AqueousLiquid");
724 auto const& solid_phase = medium.phase("Solid");
725 MPL::Phase const* gas_phase =
726 medium.hasPhase("Gas") ? &medium.phase("Gas") : nullptr;
727 MPL::VariableArray variables;
728 MPL::VariableArray variables_prev;
729
730 unsigned const n_integration_points =
731 _integration_method.getNumberOfPoints();
732 for (unsigned ip = 0; ip < n_integration_points; ip++)
733 {
734 auto const& w = _ip_data[ip].integration_weight;
735
736 auto const& N = _ip_data[ip].N;
737 auto const& dNdx = _ip_data[ip].dNdx;
738
739 ParameterLib::SpatialPosition const x_position{
740 std::nullopt, _element.getID(), ip,
743 _element, N))};
744
745 double T_ip;
747
748 double p_cap_ip;
749 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
750
751 double p_cap_prev_ip;
752 NumLib::shapeFunctionInterpolate(-p_L_prev, N, p_cap_prev_ip);
753
754 variables.capillary_pressure = p_cap_ip;
755 variables.liquid_phase_pressure = -p_cap_ip;
756 // setting pG to 1 atm
757 // TODO : rewrite equations s.t. p_L = pG-p_cap
758 variables.gas_phase_pressure = 1.0e5;
759 variables.temperature = T_ip;
760
761 auto& S_L = _ip_data[ip].saturation;
762 auto const S_L_prev = _ip_data[ip].saturation_prev;
763 auto const alpha =
764 medium[MPL::PropertyType::biot_coefficient].template value<double>(
765 variables, x_position, t, dt);
766
767 auto& solid_elasticity = *_process_data.simplified_elasticity;
768 // TODO (buchwaldj)
769 // is bulk_modulus good name for bulk modulus of solid skeleton?
770 auto const beta_S =
771 solid_elasticity.bulkCompressibilityFromYoungsModulus(
772 solid_phase, variables, x_position, t, dt);
773 auto const beta_SR = (1 - alpha) * beta_S;
774 variables.grain_compressibility = beta_SR;
775
776 auto const rho_LR =
777 liquid_phase[MPL::PropertyType::density].template value<double>(
778 variables, x_position, t, dt);
779 auto const& b = _process_data.specific_body_force;
780
781 double const drho_LR_dp =
782 liquid_phase[MPL::PropertyType::density].template dValue<double>(
783 variables, MPL::Variable::liquid_phase_pressure, x_position, t,
784 dt);
785 auto const beta_LR = drho_LR_dp / rho_LR;
786
787 S_L = medium[MPL::PropertyType::saturation].template value<double>(
788 variables, x_position, t, dt);
789 variables.liquid_saturation = S_L;
790 variables_prev.liquid_saturation = S_L_prev;
791
792 // tangent derivative for Jacobian
793 double const dS_L_dp_cap =
794 medium[MPL::PropertyType::saturation].template dValue<double>(
795 variables, MPL::Variable::capillary_pressure, x_position, t,
796 dt);
797 // secant derivative from time discretization for storage
798 // use tangent, if secant is not available
799 double const DeltaS_L_Deltap_cap =
800 (p_cap_ip == p_cap_prev_ip)
801 ? dS_L_dp_cap
802 : (S_L - S_L_prev) / (p_cap_ip - p_cap_prev_ip);
803
804 auto chi_S_L = S_L;
805 auto chi_S_L_prev = S_L_prev;
806 if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
807 {
808 auto const chi = [&medium, x_position, t, dt](double const S_L)
809 {
810 MPL::VariableArray variables;
811 variables.liquid_saturation = S_L;
812 return medium[MPL::PropertyType::bishops_effective_stress]
813 .template value<double>(variables, x_position, t, dt);
814 };
815 chi_S_L = chi(S_L);
816 chi_S_L_prev = chi(S_L_prev);
817 }
818 // TODO (buchwaldj)
819 // should solid_grain_pressure or effective_pore_pressure remain?
820 // double const p_FR = -chi_S_L * p_cap_ip;
821 // variables.solid_grain_pressure = p_FR;
822
823 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
824 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
825
826 auto& phi = _ip_data[ip].porosity;
827 { // Porosity update
828
829 variables_prev.porosity = _ip_data[ip].porosity_prev;
830 phi = medium[MPL::PropertyType::porosity].template value<double>(
831 variables, variables_prev, x_position, t, dt);
832 variables.porosity = phi;
833 }
834
835 if (alpha < phi)
836 {
837 OGS_FATAL(
838 "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
839 "porosity {} in element/integration point {}/{}.",
840 alpha, phi, _element.getID(), ip);
841 }
842
843 auto const K_pT_thermal_osmosis =
844 (solid_phase.hasProperty(
846 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
847 solid_phase
848 [MPL::PropertyType::thermal_osmosis_coefficient]
849 .value(variables, x_position, t, dt))
850 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
851
852 double const k_rel =
853 medium[MPL::PropertyType::relative_permeability]
854 .template value<double>(variables, x_position, t, dt);
855 auto const mu =
856 liquid_phase[MPL::PropertyType::viscosity].template value<double>(
857 variables, x_position, t, dt);
858
859 auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
860 medium[MPL::PropertyType::permeability].value(variables, x_position,
861 t, dt));
862
863 GlobalDimMatrixType const Ki_over_mu = K_intrinsic / mu;
864 GlobalDimMatrixType const rho_Ki_over_mu = rho_LR * Ki_over_mu;
865
866 // Consider anisotropic thermal expansion.
867 // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
868 // component.
869 Eigen::Matrix<double, 3, 3> const
870 solid_linear_thermal_expansion_coefficient =
872 solid_phase
874 .value(variables, x_position, t, dt));
875
876 auto const rho_SR =
877 solid_phase[MPL::PropertyType::density].template value<double>(
878 variables, x_position, t, dt);
879
880 //
881 // pressure equation, pressure part.
882 //
883 local_K
884 .template block<pressure_size, pressure_size>(pressure_index,
885 pressure_index)
886 .noalias() += dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
887
888 const double alphaB_minus_phi = alpha - phi;
889 double const a0 = alphaB_minus_phi * beta_SR;
890 double const specific_storage_a_p =
891 S_L * (phi * beta_LR + S_L * a0 +
892 chi_S_L * alpha * alpha *
893 solid_elasticity.storageContribution(
894 solid_phase, variables, x_position, t, dt));
895 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
896
897 local_M
898 .template block<pressure_size, pressure_size>(pressure_index,
899 pressure_index)
900 .noalias() += N.transpose() * rho_LR *
901 (specific_storage_a_p -
902 specific_storage_a_S * DeltaS_L_Deltap_cap) *
903 N * w;
904
905 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
906 dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
907
908 //
909 // pressure equation, temperature part.
910 //
911 double const fluid_volumetric_thermal_expansion_coefficient =
912 MPL::getLiquidThermalExpansivity(liquid_phase, variables, rho_LR,
913 x_position, t, dt);
914 const double eff_thermal_expansion =
915 S_L * (alphaB_minus_phi *
916 solid_linear_thermal_expansion_coefficient.trace() +
917 phi * fluid_volumetric_thermal_expansion_coefficient +
918 alpha * solid_elasticity.thermalExpansivityContribution(
919 solid_linear_thermal_expansion_coefficient,
920 solid_phase, variables, x_position, t, dt));
921
922 local_K
923 .template block<pressure_size, temperature_size>(pressure_index,
924 temperature_index)
925 .noalias() +=
926 dNdx.transpose() * rho_LR * K_pT_thermal_osmosis * dNdx * w;
927
928 local_M
929 .template block<pressure_size, temperature_size>(pressure_index,
930 temperature_index)
931 .noalias() -=
932 N.transpose() * rho_LR * eff_thermal_expansion * N * w;
933
934 //
935 // temperature equation.
936 //
937 {
938 auto const specific_heat_capacity_fluid =
940 .template value<double>(variables, x_position, t, dt);
941
942 auto const specific_heat_capacity_solid =
943 solid_phase
945 .template value<double>(variables, x_position, t, dt);
946
947 local_M
948 .template block<temperature_size, temperature_size>(
949 temperature_index, temperature_index)
950 .noalias() +=
951 w *
952 (rho_SR * specific_heat_capacity_solid * (1 - phi) +
953 (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
954 N.transpose() * N;
955
956 auto const thermal_conductivity =
957 MaterialPropertyLib::formEigenTensor<GlobalDim>(
959 thermal_conductivity]
960 .value(variables, x_position, t, dt));
961
962 GlobalDimVectorType const velocity_L = GlobalDimVectorType(
963 -Ki_over_mu * k_rel * (dNdx * p_L - rho_LR * b) -
964 K_pT_thermal_osmosis * dNdx * T);
965
966 local_K
967 .template block<temperature_size, temperature_size>(
968 temperature_index, temperature_index)
969 .noalias() += (dNdx.transpose() * thermal_conductivity * dNdx +
970 N.transpose() * velocity_L.transpose() * dNdx *
971 rho_LR * specific_heat_capacity_fluid) *
972 w;
973 local_K
974 .template block<temperature_size, pressure_size>(
975 temperature_index, pressure_index)
976 .noalias() +=
977 dNdx.transpose() * T_ip * K_pT_thermal_osmosis * dNdx * w;
978 }
979 if (gas_phase && S_L < 1.0)
980 {
981 variables.density = rho_LR;
982
983 double const rho_wv =
985 .template value<double>(variables, x_position, t, dt);
986
987 double const drho_wv_dT =
989 .template dValue<double>(variables,
990 MPL::Variable::temperature,
991 x_position, t, dt);
992 double const drho_wv_dp =
994 .template dValue<double>(
995 variables, MPL::Variable::liquid_phase_pressure,
996 x_position, t, dt);
997 auto const f_Tv =
998 gas_phase
999 ->property(
1000 MPL::PropertyType::thermal_diffusion_enhancement_factor)
1001 .template value<double>(variables, x_position, t, dt);
1002
1003 variables.porosity = phi;
1004 auto const tortuosity =
1005 medium.property(MPL::PropertyType::tortuosity)
1006 .template value<double>(variables, x_position, t, dt);
1007 double const D_v =
1008 phi * (1.0 - S_L) * tortuosity *
1009 gas_phase->property(MPL::PropertyType::diffusion)
1010 .template value<double>(variables, x_position, t, dt);
1011
1012 double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
1013 double const D_pv = D_v * drho_wv_dp;
1014
1015 GlobalDimVectorType const grad_T = dNdx * T;
1016 GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
1017 GlobalDimVectorType const vapour_flux =
1018 -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
1019 double const specific_heat_capacity_vapour =
1021 .template value<double>(variables, x_position, t, dt);
1022
1023 local_M
1024 .template block<temperature_size, temperature_size>(
1025 temperature_index, temperature_index)
1026 .noalias() +=
1027 w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
1028 N.transpose() * N;
1029
1030 local_K
1031 .template block<temperature_size, temperature_size>(
1032 temperature_index, temperature_index)
1033 .noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
1034 specific_heat_capacity_vapour * w;
1035
1036 double const storage_coefficient_by_water_vapor =
1037 phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
1038 local_M
1039 .template block<pressure_size, pressure_size>(pressure_index,
1040 pressure_index)
1041 .noalias() +=
1042 N.transpose() * storage_coefficient_by_water_vapor * N * w;
1043
1044 double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
1045 local_M
1046 .template block<pressure_size, temperature_size>(
1047 pressure_index, temperature_index)
1048 .noalias() += N.transpose() * vapor_expansion_factor * N * w;
1049
1050 local_rhs.template segment<pressure_size>(pressure_index)
1051 .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
1052
1053 local_K
1054 .template block<pressure_size, pressure_size>(pressure_index,
1055 pressure_index)
1056 .noalias() += dNdx.transpose() * D_pv * dNdx * w;
1057
1058 //
1059 // Latent heat term
1060 //
1061 if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
1062 {
1063 double const factor = phi * (1 - S_L) / rho_LR;
1064 // The volumetric latent heat of vaporization of liquid water
1065 double const L0 =
1066 gas_phase->property(MPL::PropertyType::specific_latent_heat)
1067 .template value<double>(variables, x_position, t, dt) *
1068 rho_LR;
1069
1070 double const drho_LR_dT =
1071 liquid_phase.property(MPL::PropertyType::density)
1072 .template dValue<double>(variables,
1073 MPL::Variable::temperature,
1074 x_position, t, dt);
1075
1076 double const rho_wv_over_rho_L = rho_wv / rho_LR;
1077 local_M
1078 .template block<temperature_size, temperature_size>(
1079 temperature_index, temperature_index)
1080 .noalias() +=
1081 factor * L0 *
1082 (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
1083 N.transpose() * N * w;
1084
1085 local_M
1086 .template block<temperature_size, pressure_size>(
1087 temperature_index, pressure_index)
1088 .noalias() +=
1089 (factor * L0 *
1090 (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
1091 L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
1092 N.transpose() * N * w;
1093
1094 // temperature equation, temperature part
1095 local_K
1096 .template block<temperature_size, temperature_size>(
1097 temperature_index, temperature_index)
1098 .noalias() +=
1099 L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
1100 // temperature equation, pressure part
1101 local_K
1102 .template block<temperature_size, pressure_size>(
1103 temperature_index, pressure_index)
1104 .noalias() +=
1105 L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
1106 }
1107 }
1108 }
1109
1110 if (_process_data.apply_mass_lumping)
1111 {
1112 auto Mpp = local_M.template block<pressure_size, pressure_size>(
1113 pressure_index, pressure_index);
1114 Mpp = Mpp.colwise().sum().eval().asDiagonal();
1115 }
1116}
1117
1118template <typename ShapeFunction, int GlobalDim>
1119std::vector<double> const&
1122 const double /*t*/,
1123 std::vector<GlobalVector*> const& /*x*/,
1124 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1125 std::vector<double>& cache) const
1126{
1127 unsigned const n_integration_points =
1128 _integration_method.getNumberOfPoints();
1129
1130 cache.clear();
1131 auto cache_matrix = MathLib::createZeroedMatrix<
1132 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1133 cache, GlobalDim, n_integration_points);
1134
1135 for (unsigned ip = 0; ip < n_integration_points; ip++)
1136 {
1137 cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
1138 }
1139
1140 return cache;
1141}
1142
1143template <typename ShapeFunction, int GlobalDim>
1144std::vector<double> ThermoRichardsFlowLocalAssembler<
1145 ShapeFunction, GlobalDim>::getSaturation() const
1146{
1147 std::vector<double> result;
1148 getIntPtSaturation(0, {}, {}, result);
1149 return result;
1150}
1151
1152template <typename ShapeFunction, int GlobalDim>
1153std::vector<double> const&
1155 const double /*t*/,
1156 std::vector<GlobalVector*> const& /*x*/,
1157 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1158 std::vector<double>& cache) const
1159{
1161 _ip_data, &IpData::saturation, cache);
1162}
1163
1164template <typename ShapeFunction, int GlobalDim>
1165std::vector<double>
1167{
1168 std::vector<double> result;
1169 getIntPtPorosity(0, {}, {}, result);
1170 return result;
1171}
1172
1173template <typename ShapeFunction, int GlobalDim>
1174std::vector<double> const&
1176 const double /*t*/,
1177 std::vector<GlobalVector*> const& /*x*/,
1178 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1179 std::vector<double>& cache) const
1180{
1182 &IpData::porosity, cache);
1183}
1184
1185template <typename ShapeFunction, int GlobalDim>
1186std::vector<double> const&
1189 const double /*t*/,
1190 std::vector<GlobalVector*> const& /*x*/,
1191 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
1192 std::vector<double>& cache) const
1193{
1195 _ip_data, &IpData::dry_density_solid, cache);
1196}
1197
1198template <typename ShapeFunction, int GlobalDim>
1200 computeSecondaryVariableConcrete(double const t, double const dt,
1201 Eigen::VectorXd const& local_x,
1202 Eigen::VectorXd const& local_x_prev)
1203{
1204 auto const T =
1205 local_x.template segment<temperature_size>(temperature_index);
1206
1207 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
1208
1209 auto p_L_prev =
1210 local_x_prev.template segment<pressure_size>(pressure_index);
1211
1212 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
1213 auto const& liquid_phase = medium.phase("AqueousLiquid");
1214 auto const& solid_phase = medium.phase("Solid");
1215 MPL::VariableArray variables;
1216 MPL::VariableArray variables_prev;
1217
1218 unsigned const n_integration_points =
1219 _integration_method.getNumberOfPoints();
1220
1221 double saturation_avg = 0;
1222 double porosity_avg = 0;
1223
1224 for (unsigned ip = 0; ip < n_integration_points; ip++)
1225 {
1226 auto const& N = _ip_data[ip].N;
1227
1228 ParameterLib::SpatialPosition const x_position{
1229 std::nullopt, _element.getID(), ip,
1232 _element, N))};
1233
1234 double T_ip;
1236
1237 double p_cap_ip;
1238 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
1239
1240 double p_cap_prev_ip;
1241 NumLib::shapeFunctionInterpolate(-p_L_prev, N, p_cap_prev_ip);
1242
1243 variables.capillary_pressure = p_cap_ip;
1244 variables.liquid_phase_pressure = -p_cap_ip;
1245 // setting pG to 1 atm
1246 // TODO : rewrite equations s.t. p_L = pG-p_cap
1247 variables.gas_phase_pressure = 1.0e5;
1248
1249 variables.temperature = T_ip;
1250
1251 auto& S_L = _ip_data[ip].saturation;
1252 auto const S_L_prev = _ip_data[ip].saturation_prev;
1253 S_L = medium[MPL::PropertyType::saturation].template value<double>(
1254 variables, x_position, t, dt);
1255 variables.liquid_saturation = S_L;
1256 variables_prev.liquid_saturation = S_L_prev;
1257
1258 auto chi_S_L = S_L;
1259 auto chi_S_L_prev = S_L_prev;
1260 if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
1261 {
1262 auto const chi = [&medium, x_position, t, dt](double const S_L)
1263 {
1264 MPL::VariableArray variables;
1265 variables.liquid_saturation = S_L;
1266 return medium[MPL::PropertyType::bishops_effective_stress]
1267 .template value<double>(variables, x_position, t, dt);
1268 };
1269 chi_S_L = chi(S_L);
1270 chi_S_L_prev = chi(S_L_prev);
1271 }
1272 variables.effective_pore_pressure = -chi_S_L * p_cap_ip;
1273 variables_prev.effective_pore_pressure = -chi_S_L_prev * p_cap_prev_ip;
1274
1275 auto const alpha =
1276 medium[MPL::PropertyType::biot_coefficient].template value<double>(
1277 variables, x_position, t, dt);
1278
1279 auto& solid_elasticity = *_process_data.simplified_elasticity;
1280 auto const beta_S =
1281 solid_elasticity.bulkCompressibilityFromYoungsModulus(
1282 solid_phase, variables, x_position, t, dt);
1283 auto const beta_SR = (1 - alpha) * beta_S;
1284 variables.grain_compressibility = beta_SR;
1285
1286 auto& phi = _ip_data[ip].porosity;
1287 { // Porosity update
1288 variables_prev.porosity = _ip_data[ip].porosity_prev;
1289 phi = medium[MPL::PropertyType::porosity].template value<double>(
1290 variables, variables_prev, x_position, t, dt);
1291 variables.porosity = phi;
1292 }
1293
1294 auto const mu =
1295 liquid_phase[MPL::PropertyType::viscosity].template value<double>(
1296 variables, x_position, t, dt);
1297 auto const rho_LR =
1298 liquid_phase[MPL::PropertyType::density].template value<double>(
1299 variables, x_position, t, dt);
1300
1301 auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
1302 medium[MPL::PropertyType::permeability].value(variables, x_position,
1303 t, dt));
1304
1305 double const k_rel =
1306 medium[MPL::PropertyType::relative_permeability]
1307 .template value<double>(variables, x_position, t, dt);
1308
1309 GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
1310
1311 auto const rho_SR =
1312 solid_phase[MPL::PropertyType::density].template value<double>(
1313 variables, x_position, t, dt);
1314 _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
1315
1316 auto const& b = _process_data.specific_body_force;
1317
1318 auto const K_pT_thermal_osmosis =
1319 (solid_phase.hasProperty(
1321 ? MaterialPropertyLib::formEigenTensor<GlobalDim>(
1322 solid_phase
1323 [MPL::PropertyType::thermal_osmosis_coefficient]
1324 .value(variables, x_position, t, dt))
1325 : Eigen::MatrixXd::Zero(GlobalDim, GlobalDim));
1326
1327 // Compute the velocity
1328 auto const& dNdx = _ip_data[ip].dNdx;
1329 _ip_data[ip].v_darcy.noalias() = -K_over_mu * dNdx * p_L -
1330 K_pT_thermal_osmosis * dNdx * T +
1331 rho_LR * K_over_mu * b;
1332
1333 saturation_avg += S_L;
1334 porosity_avg += phi;
1335 }
1336 saturation_avg /= n_integration_points;
1337 porosity_avg /= n_integration_points;
1338
1339 (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
1340 (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
1341}
1342
1343template <typename ShapeFunction, int GlobalDim>
1345 ShapeFunction, GlobalDim>::getNumberOfIntegrationPoints() const
1346{
1347 return _integration_method.getNumberOfPoints();
1348}
1349
1350} // namespace ThermoRichardsFlow
1351} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
Property const & property(PropertyType const &p) const
Definition Phase.cpp:53
bool hasProperty(PropertyType const &p) const
Definition Phase.cpp:69
double getWeight() const
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void assembleWithJacobian(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &, std::vector< double > &, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
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::size_t setIPDataInitialConditions(std::string_view const name, double const *values, int const integration_order) override
std::vector< double > const & getIntPtDryDensitySolid(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
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
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_rhs_data) override
ThermoRichardsFlowLocalAssembler(ThermoRichardsFlowLocalAssembler const &)=delete
std::vector< double > const & getIntPtPorosity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void setInitialConditionsConcrete(Eigen::VectorXd const local_x, double const t, int const process_id) override
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
void computeSecondaryVariableConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev) override
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
double getLiquidThermalExpansivity(Phase const &phase, VariableArray const &vars, const double density, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
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::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
MaterialPropertyLib::MaterialSpatialDistributionMap media_map