OGS
TH2MFEM-impl.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/LU>
14
15#include "MaterialLib/MPL/Medium.h"
16#include "MaterialLib/MPL/Property.h"
17#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
18#include "MaterialLib/PhysicalConstant.h"
19#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
20#include "MathLib/EigenBlockMatrixView.h"
21#include "MathLib/KelvinVector.h"
22#include "NumLib/Function/Interpolation.h"
23#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
24
25namespace ProcessLib
26{
27namespace TH2M
28{
29namespace MPL = MaterialPropertyLib;
30
31template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
32 int DisplacementDim>
33TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
34 DisplacementDim>::
35 TH2MLocalAssembler(
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,
40 TH2MProcessData<DisplacementDim>& process_data)
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 =
47 _integration_method.getNumberOfPoints();
48
49 _ip_data.reserve(n_integration_points);
50 _secondary_data.N_u.resize(n_integration_points);
51
52 auto const shape_matrices_u =
53 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
54 ShapeMatricesTypeDisplacement,
55 DisplacementDim>(e, is_axially_symmetric,
56 _integration_method);
57
58 auto const shape_matrices_p =
59 NumLib::initShapeMatrices<ShapeFunctionPressure,
60 ShapeMatricesTypePressure, DisplacementDim>(
61 e, is_axially_symmetric, _integration_method);
62
63 auto const& solid_material =
64 MaterialLib::Solids::selectSolidConstitutiveRelation(
65 _process_data.solid_materials, _process_data.material_ids,
66 e.getID());
67
68 for (unsigned ip = 0; ip < n_integration_points; ip++)
69 {
70 _ip_data.emplace_back(solid_material);
71 auto& ip_data = _ip_data[ip];
72 auto const& sm_u = shape_matrices_u[ip];
73 ip_data.integration_weight =
74 _integration_method.getWeightedPoint(ip).getWeight() *
75 sm_u.integralMeasure * sm_u.detJ;
76
77 ip_data.N_u = sm_u.N;
78 ip_data.dNdx_u = sm_u.dNdx;
79
80 ip_data.N_p = shape_matrices_p[ip].N;
81 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
82
83 _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
84 }
85}
86
87template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
88 int DisplacementDim>
89std::vector<ConstitutiveVariables<DisplacementDim>> TH2MLocalAssembler<
90 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
91 updateConstitutiveVariables(Eigen::VectorXd const& local_x,
92 Eigen::VectorXd const& local_x_prev,
93 double const t, double const dt)
94{
95 [[maybe_unused]] auto const matrix_size =
96 gas_pressure_size + capillary_pressure_size + temperature_size +
97 displacement_size;
98
99 assert(local_x.size() == matrix_size);
100
101 auto const gas_pressure =
102 local_x.template segment<gas_pressure_size>(gas_pressure_index);
103 auto const capillary_pressure =
104 local_x.template segment<capillary_pressure_size>(
105 capillary_pressure_index);
106
107 auto const temperature =
108 local_x.template segment<temperature_size>(temperature_index);
109 auto const temperature_prev =
110 local_x_prev.template segment<temperature_size>(temperature_index);
111
112 auto const displacement =
113 local_x.template segment<displacement_size>(displacement_index);
114
115 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
116 auto const& gas_phase = medium.phase("Gas");
117 auto const& liquid_phase = medium.phase("AqueousLiquid");
118 auto const& solid_phase = medium.phase("Solid");
119
120 unsigned const n_integration_points =
121 _integration_method.getNumberOfPoints();
122
123 std::vector<ConstitutiveVariables<DisplacementDim>>
124 ip_constitutive_variables(n_integration_points);
125
126 PhaseTransitionModelVariables ptmv;
127
128 for (unsigned ip = 0; ip < n_integration_points; ip++)
129 {
130 auto& ip_data = _ip_data[ip];
131 auto& ip_cv = ip_constitutive_variables[ip];
132
133 auto const& Np = ip_data.N_p;
134 auto const& NT = Np;
135 auto const& Nu = ip_data.N_u;
136 auto const& gradNu = ip_data.dNdx_u;
137 auto const& gradNp = ip_data.dNdx_p;
138 ParameterLib::SpatialPosition const pos{
139 std::nullopt, _element.getID(), ip,
140 MathLib::Point3d(
141 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
142 ShapeMatricesTypeDisplacement>(
143 _element, Nu))};
144 auto const x_coord =
145 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
146 ShapeMatricesTypeDisplacement>(
147 _element, Nu);
148
149 double const T = NT.dot(temperature);
150 double const T_prev = NT.dot(temperature_prev);
151 double const pGR = Np.dot(gas_pressure);
152 double const pCap = Np.dot(capillary_pressure);
153 double const pLR = pGR - pCap;
154 GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
155 GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
156 GlobalDimVectorType const gradT = gradNp * temperature;
157
158 MPL::VariableArray vars;
159 MPL::VariableArray vars_prev;
160 vars.temperature = T;
161 vars.phase_pressure = pGR;
162 vars.capillary_pressure = pCap;
163 vars.liquid_phase_pressure = pLR;
164
165 // medium properties
166 auto const C_el =
167 ip_data.computeElasticTangentStiffness(t, pos, dt, T_prev, T);
168 auto const K_S = ip_data.solid_material.getBulkModulus(t, pos, &C_el);
169
170 ip_data.alpha_B = medium.property(MPL::PropertyType::biot_coefficient)
171 .template value<double>(vars, pos, t, dt);
172
173 auto const Bu =
174 LinearBMatrix::computeBMatrix<DisplacementDim,
175 ShapeFunctionDisplacement::NPOINTS,
176 typename BMatricesType::BMatrixType>(
177 gradNu, Nu, x_coord, _is_axially_symmetric);
178
179 auto& eps = ip_data.eps;
180 eps.noalias() = Bu * displacement;
181
182 // Set volumetric strain for the general case without swelling.
183 vars.volumetric_strain = Invariants::trace(eps);
184
185 ip_data.s_L =
186 medium.property(MPL::PropertyType::saturation)
187 .template value<double>(
188 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
189
190 vars.liquid_saturation = ip_data.s_L;
191 vars_prev.liquid_saturation = ip_data.s_L_prev;
192
193 auto const chi = [&medium, pos, t, dt](double const s_L)
194 {
195 MPL::VariableArray vs;
196 vs.liquid_saturation = s_L;
197 return medium.property(MPL::PropertyType::bishops_effective_stress)
198 .template value<double>(vs, pos, t, dt);
199 };
200 ip_cv.chi_s_L = chi(ip_data.s_L);
201 ip_cv.dchi_ds_L =
202 medium.property(MPL::PropertyType::bishops_effective_stress)
203 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
204 pos, t, dt);
205
206 // relative permeability
207 // Set mechanical variables for the intrinsic permeability model
208 // For stress dependent permeability.
209 {
210 auto const sigma_total =
211 (_ip_data[ip].sigma_eff - ip_data.alpha_B *
212 (pGR - ip_cv.chi_s_L * pCap) *
213 Invariants::identity2)
214 .eval();
215
216 vars.total_stress.emplace<SymmetricTensor>(
217 MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
218 sigma_total));
219 }
220 vars.equivalent_plastic_strain =
221 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
222 vars.mechanical_strain
223 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
224 eps);
225
226 // intrinsic permeability
227 ip_data.k_S = MPL::formEigenTensor<DisplacementDim>(
228 medium.property(MPL::PropertyType::permeability)
229 .value(vars, pos, t, dt));
230
231 ip_data.k_rel_G =
232 medium
233 .property(
234 MPL::PropertyType::relative_permeability_nonwetting_phase)
235 .template value<double>(vars, pos, t, dt);
236
237 auto const dk_rel_G_ds_L =
238 medium[MPL::PropertyType::relative_permeability_nonwetting_phase]
239 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
240 pos, t, dt);
241
242 ip_data.k_rel_L =
243 medium.property(MPL::PropertyType::relative_permeability)
244 .template value<double>(vars, pos, t, dt);
245
246 auto const dk_rel_L_ds_L =
247 medium[MPL::PropertyType::relative_permeability]
248 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
249 pos, t, dt);
250
251 // solid phase compressibility
252 ip_data.beta_p_SR = (1. - ip_data.alpha_B) / K_S;
253
254 // If there is swelling stress rate, compute swelling stress.
255 if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
256 {
257 auto& sigma_sw = ip_data.sigma_sw;
258 auto const& sigma_sw_prev = ip_data.sigma_sw_prev;
259
260 sigma_sw = sigma_sw_prev;
261
262 auto const sigma_sw_dot =
263 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
264 MPL::formEigenTensor<3>(
265 solid_phase[MPL::PropertyType::swelling_stress_rate]
266 .value(vars, vars_prev, pos, t, dt)));
267 sigma_sw += sigma_sw_dot * dt;
268 }
269
270 // solid phase linear thermal expansion coefficient
271 ip_data.alpha_T_SR = MathLib::KelvinVector::tensorToKelvin<
272 DisplacementDim>(MaterialPropertyLib::formEigenTensor<3>(
273 solid_phase
274 .property(
275 MaterialPropertyLib::PropertyType::thermal_expansivity)
276 .value(vars, pos, t, dt)));
277
278 // isotropic solid phase volumetric thermal expansion coefficient
279 ip_data.beta_T_SR = Invariants::trace(ip_data.alpha_T_SR);
280
281 MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
282 dthermal_strain = ip_data.alpha_T_SR * (T - T_prev);
283
284 auto& eps_prev = ip_data.eps_prev;
285 auto& eps_m = ip_data.eps_m;
286 auto& eps_m_prev = ip_data.eps_m_prev;
287
288 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
289
290 if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
291 {
292 eps_m.noalias() +=
293 C_el.inverse() * (ip_data.sigma_sw - ip_data.sigma_sw_prev);
294 }
295
296 vars.mechanical_strain
297 .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
298 eps_m);
299
300 auto const rho_ref_SR =
301 solid_phase.property(MPL::PropertyType::density)
302 .template value<double>(
303 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
304
305 double const T0 = _process_data.reference_temperature(t, pos)[0];
306 double const delta_T(T - T0);
307 ip_data.thermal_volume_strain = ip_data.beta_T_SR * delta_T;
308
309 // initial porosity
310 auto const phi_0 = medium.property(MPL::PropertyType::porosity)
311 .template value<double>(vars, pos, t, dt);
312
313 auto const phi_S_0 = 1. - phi_0;
314
315#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
316 auto const& m = Invariants::identity2;
317 double const div_u = m.transpose() * eps;
318
319 const double phi_S = phi_S_0 * (1. + ip_data.thermal_volume_strain -
320 ip_data.alpha_B * div_u);
321#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
322 const double phi_S = phi_S_0;
323#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
324
325 // porosity
326 ip_data.phi = 1. - phi_S;
327 vars.porosity = ip_data.phi;
328
329 // solid phase density
330#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
331 auto const rhoSR = rho_ref_SR * (1. - ip_data.thermal_volume_strain +
332 (ip_data.alpha_B - 1.) * div_u);
333#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
334 auto const rhoSR = rho_ref_SR;
335#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
336
337 ip_cv.C = ip_data.updateConstitutiveRelation(vars, t, pos, dt, T_prev);
338
339 // constitutive model object as specified in process creation
340 auto& ptm = *_process_data.phase_transition_model_;
341 ptmv = ptm.updateConstitutiveVariables(ptmv, &medium, vars, pos, t, dt);
342 auto const& c = ptmv;
343
344 auto const phi_L = ip_data.s_L * ip_data.phi;
345 auto const phi_G = (1. - ip_data.s_L) * ip_data.phi;
346
347 // thermal conductivity
348 ip_data.lambda = MaterialPropertyLib::formEigenTensor<DisplacementDim>(
349 medium
350 .property(
351 MaterialPropertyLib::PropertyType::thermal_conductivity)
352 .value(vars, pos, t, dt));
353
354 auto const cpS =
355 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
356 .template value<double>(vars, pos, t, dt);
357 ip_data.h_S = cpS * T;
358 auto const u_S = ip_data.h_S;
359
360 ip_data.rho_u_eff = phi_G * c.rhoGR * c.uG + phi_L * c.rhoLR * c.uL +
361 phi_S * rhoSR * u_S;
362
363 ip_data.rho_G_h_G = phi_G * c.rhoGR * c.hG;
364 ip_data.rho_L_h_L = phi_L * c.rhoLR * c.hL;
365 ip_data.rho_S_h_S = phi_S * rhoSR * ip_data.h_S;
366
367 ip_data.muGR = c.muGR;
368 ip_data.muLR = c.muLR;
369
370 ip_data.rhoGR = c.rhoGR;
371 ip_data.rhoLR = c.rhoLR;
372 ip_data.rhoSR = rhoSR;
373
374 ip_data.rhoCGR = c.rhoCGR;
375 ip_data.rhoCLR = c.rhoCLR;
376 ip_data.rhoWGR = c.rhoWGR;
377 ip_data.rhoWLR = c.rhoWLR;
378
379 ip_data.dxmWG_dpGR = c.dxmWG_dpGR;
380 ip_data.dxmWG_dpCap = c.dxmWG_dpCap;
381 ip_data.dxmWG_dT = c.dxmWG_dT;
382
383 ip_data.dxmWL_dpGR = c.dxmWL_dpGR;
384 ip_data.dxmWL_dpCap = c.dxmWL_dpCap;
385 ip_data.dxmWL_dT = c.dxmWL_dT;
386
387 ip_data.dxmWL_dpLR = c.dxmWL_dpLR;
388
389 // for variable output
390 ip_data.xnCG = 1. - c.xnWG;
391 ip_data.xmCG = 1. - c.xmWG;
392 ip_data.xmWG = c.xmWG;
393 ip_data.xmWL = c.xmWL;
394 auto const xmCL = 1. - c.xmWL;
395
396 ip_data.diffusion_coefficient_vapour = c.diffusion_coefficient_vapour;
397 ip_data.diffusion_coefficient_solute = c.diffusion_coefficient_solute;
398
399 ip_data.h_G = c.hG;
400 ip_data.h_CG = c.hCG;
401 ip_data.h_WG = c.hWG;
402 ip_data.h_L = c.hL;
403 ip_data.pWGR = c.pWGR;
404
405 const GlobalDimVectorType gradxmWG = ip_data.dxmWG_dpGR * gradpGR +
406 ip_data.dxmWG_dpCap * gradpCap +
407 ip_data.dxmWG_dT * gradT;
408 const GlobalDimVectorType gradxmCG = -gradxmWG;
409
410 const GlobalDimVectorType gradxmWL = ip_data.dxmWL_dpGR * gradpGR +
411 ip_data.dxmWL_dpCap * gradpCap +
412 ip_data.dxmWL_dT * gradT;
413 const GlobalDimVectorType gradxmCL = -gradxmWL;
414
415 // Todo: factor -phiAlpha / xmZetaAlpha * DZetaAlpha can be evaluated in
416 // the respective phase transition model, here only the multiplication
417 // with the gradient of the mass fractions should take place.
418
419 ip_data.d_CG = ip_data.xmCG == 0.
420 ? 0. * gradxmCG // Keep d_CG's dimension and prevent
421 // division by zero
422 : -phi_G / ip_data.xmCG *
423 ip_data.diffusion_coefficient_vapour *
424 gradxmCG;
425
426 ip_data.d_WG = ip_data.xmWG == 0.
427 ? 0. * gradxmWG // Keep d_WG's dimension and prevent
428 // division by zero
429 : -phi_G / ip_data.xmWG *
430 ip_data.diffusion_coefficient_vapour *
431 gradxmWG;
432
433 ip_data.d_CL = xmCL == 0. ? 0. * gradxmCL // Keep d_CL's dimension and
434 // prevent division by zero
435 : -phi_L / xmCL *
436 ip_data.diffusion_coefficient_solute *
437 gradxmCL;
438
439 ip_data.d_WL = ip_data.xmWL == 0.
440 ? 0. * gradxmWL // Keep d_WG's dimension and prevent
441 // division by zero
442 : -phi_L / ip_data.xmWL *
443 ip_data.diffusion_coefficient_solute *
444 gradxmWL;
445
446 // ---------------------------------------------------------------------
447 // Derivatives for Jacobian
448 // ---------------------------------------------------------------------
449 auto const drho_LR_dT =
450 liquid_phase.property(MPL::PropertyType::density)
451 .template dValue<double>(vars, MPL::Variable::temperature, pos,
452 t, dt);
453 auto const drho_SR_dT =
454 solid_phase.property(MPL::PropertyType::density)
455 .template dValue<double>(vars, MPL::Variable::temperature,
456 pos, t, dt)
457#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
458 * (1. - ip_data.thermal_volume_strain +
459 (ip_data.alpha_B - 1.) * div_u) -
460 rho_ref_SR * ip_data.beta_T_SR
461#endif
462 ;
463
464 // porosity
465 auto const dphi_0_dT =
466 medium[MPL::PropertyType::porosity].template dValue<double>(
467 vars, MPL::Variable::temperature, pos, t, dt);
468
469 auto const dphi_S_0_dT = -dphi_0_dT;
470 const double dphi_S_dT = dphi_S_0_dT
471#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
472 * (1. + ip_data.thermal_volume_strain -
473 ip_data.alpha_B * div_u) +
474 phi_S_0 * ip_data.beta_T_SR
475#endif
476 ;
477
478 ip_cv.drho_u_eff_dT =
479 phi_G * c.drho_GR_dT * c.uG + phi_G * c.rhoGR * c.du_G_dT +
480 phi_L * drho_LR_dT * c.uL + phi_L * c.rhoLR * c.du_L_dT +
481 phi_S * drho_SR_dT * u_S + phi_S * rhoSR * cpS +
482 dphi_S_dT * rhoSR * u_S;
483
484 ip_cv.ds_L_dp_cap =
485 medium[MPL::PropertyType::saturation].template dValue<double>(
486 vars, MPL::Variable::capillary_pressure, pos, t, dt);
487
488 // ds_L_dp_GR = 0;
489 // ds_G_dp_GR = -ds_L_dp_GR;
490 double const ds_G_dp_cap = -ip_cv.ds_L_dp_cap;
491
492 // dphi_G_dp_GR = -ds_L_dp_GR * ip_data.phi = 0;
493 double const dphi_G_dp_cap = -ip_cv.ds_L_dp_cap * ip_data.phi;
494 // dphi_L_dp_GR = ds_L_dp_GR * ip_data.phi = 0;
495 double const dphi_L_dp_cap = ip_cv.ds_L_dp_cap * ip_data.phi;
496
497 auto const lambdaGR =
498 gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
499 ? MPL::formEigenTensor<DisplacementDim>(
500 gas_phase
501 .property(MPL::PropertyType::thermal_conductivity)
502 .value(vars, pos, t, dt))
503 : MPL::formEigenTensor<DisplacementDim>(0.);
504
505 auto const dlambda_GR_dT =
506 gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
507 ? MPL::formEigenTensor<DisplacementDim>(
508 gas_phase[MPL::PropertyType::thermal_conductivity].dValue(
509 vars, MPL::Variable::temperature, pos, t, dt))
510 : MPL::formEigenTensor<DisplacementDim>(0.);
511
512 auto const lambdaLR =
513 liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
514 ? MPL::formEigenTensor<DisplacementDim>(
515 liquid_phase
516 .property(MPL::PropertyType::thermal_conductivity)
517 .value(vars, pos, t, dt))
518 : MPL::formEigenTensor<DisplacementDim>(0.);
519
520 auto const dlambda_LR_dT =
521 liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
522 ? MPL::formEigenTensor<DisplacementDim>(
523 liquid_phase[MPL::PropertyType::thermal_conductivity]
524 .dValue(vars, MPL::Variable::temperature, pos, t, dt))
525 : MPL::formEigenTensor<DisplacementDim>(0.);
526
527 auto const lambdaSR =
528 solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
529 ? MPL::formEigenTensor<DisplacementDim>(
530 solid_phase
531 .property(MPL::PropertyType::thermal_conductivity)
532 .value(vars, pos, t, dt))
533 : MPL::formEigenTensor<DisplacementDim>(0.);
534
535 auto const dlambda_SR_dT =
536 solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
537 ? MPL::formEigenTensor<DisplacementDim>(
538 solid_phase[MPL::PropertyType::thermal_conductivity]
539 .dValue(vars, MPL::Variable::temperature, pos, t, dt))
540 : MPL::formEigenTensor<DisplacementDim>(0.);
541
542 ip_cv.dlambda_dp_cap =
543 dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
544
545 ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
546 phi_S * dlambda_SR_dT + dphi_S_dT * lambdaSR;
547
548 // From p_LR = p_GR - p_cap it follows for
549 // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR
550 // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR)
551 // = drho_LR/dp_LR * (1 - 0)
552 double const drho_LR_dp_GR = c.drho_LR_dp_LR;
553 double const drho_LR_dp_cap = -c.drho_LR_dp_LR;
554 // drho_GR_dp_cap = 0;
555
556 ip_cv.drho_h_eff_dp_GR =
557 /*(dphi_G_dp_GR = 0) * c.rhoGR * c.hG +*/ phi_G * c.drho_GR_dp_GR *
558 c.hG +
559 /*(dphi_L_dp_GR = 0) * c.rhoLR * c.hL +*/ phi_L * drho_LR_dp_GR *
560 c.hL;
561 ip_cv.drho_h_eff_dp_cap = dphi_G_dp_cap * c.rhoGR * c.hG +
562 /*phi_G * (drho_GR_dp_cap = 0) * c.hG +*/
563 dphi_L_dp_cap * c.rhoLR * c.hL +
564 phi_L * drho_LR_dp_cap * c.hL;
565
566 // TODO (naumov) Extend for temperature dependent porosities.
567 constexpr double dphi_G_dT = 0;
568 constexpr double dphi_L_dT = 0;
569 ip_cv.drho_h_eff_dT =
570 dphi_G_dT * c.rhoGR * c.hG + phi_G * c.drho_GR_dT * c.hG +
571 phi_G * c.rhoGR * c.dh_G_dT + dphi_L_dT * c.rhoLR * c.hL +
572 phi_L * drho_LR_dT * c.hL + phi_L * c.rhoLR * c.dh_L_dT +
573 dphi_S_dT * rhoSR * ip_data.h_S + phi_S * drho_SR_dT * ip_data.h_S +
574 phi_S * rhoSR * cpS;
575
576 ip_cv.drho_u_eff_dp_GR =
577 /*(dphi_G_dp_GR = 0) * c.rhoGR * c.uG +*/
578 phi_G * c.drho_GR_dp_GR * c.uG + phi_G * c.rhoGR * c.du_G_dp_GR +
579 /*(dphi_L_dp_GR = 0) * c.rhoLR * c.uL +*/
580 phi_L * drho_LR_dp_GR * c.uL + phi_L * c.rhoLR * c.du_L_dp_GR;
581
582 ip_cv.drho_u_eff_dp_cap = dphi_G_dp_cap * c.rhoGR * c.uG +
583 /*phi_G * (drho_GR_dp_cap = 0) * c.uG +*/
584 dphi_L_dp_cap * c.rhoLR * c.uL +
585 phi_L * drho_LR_dp_cap * c.uL +
586 phi_L * c.rhoLR * c.du_L_dp_cap;
587
588 auto const& b = _process_data.specific_body_force;
589 GlobalDimMatrixType const k_over_mu_G =
590 ip_data.k_S * ip_data.k_rel_G / ip_data.muGR;
591 GlobalDimMatrixType const k_over_mu_L =
592 ip_data.k_S * ip_data.k_rel_L / ip_data.muLR;
593
594 // dk_over_mu_G_dp_GR =
595 // ip_data.k_S * dk_rel_G_ds_L * (ds_L_dp_GR = 0) / ip_data.muGR =
596 // 0;
597 // dk_over_mu_L_dp_GR =
598 // ip_data.k_S * dk_rel_L_ds_L * (ds_L_dp_GR = 0) / ip_data.muLR =
599 // 0;
600 ip_cv.dk_over_mu_G_dp_cap =
601 ip_data.k_S * dk_rel_G_ds_L * ip_cv.ds_L_dp_cap / ip_data.muGR;
602 ip_cv.dk_over_mu_L_dp_cap =
603 ip_data.k_S * dk_rel_L_ds_L * ip_cv.ds_L_dp_cap / ip_data.muLR;
604
605 ip_data.w_GS = k_over_mu_G * c.rhoGR * b - k_over_mu_G * gradpGR;
606 ip_data.w_LS = k_over_mu_L * gradpCap + k_over_mu_L * c.rhoLR * b -
607 k_over_mu_L * gradpGR;
608
609 ip_cv.drho_GR_h_w_eff_dp_GR_Npart =
610 c.drho_GR_dp_GR * c.hG * ip_data.w_GS +
611 c.rhoGR * c.hG * k_over_mu_G * c.drho_GR_dp_GR * b;
612 ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart =
613 -c.rhoGR * c.hG * k_over_mu_G - c.rhoLR * c.hL * k_over_mu_L;
614
615 ip_cv.drho_LR_h_w_eff_dp_cap_Npart =
616 -drho_LR_dp_cap * c.hL * ip_data.w_LS -
617 c.rhoLR * c.hL * k_over_mu_L * drho_LR_dp_cap * b;
618 ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart =
619 // TODO (naumov) why the minus sign??????
620 -c.rhoLR * c.hL * k_over_mu_L;
621
622 ip_cv.drho_GR_h_w_eff_dT = c.drho_GR_dT * c.hG * ip_data.w_GS +
623 c.rhoGR * c.dh_G_dT * ip_data.w_GS +
624 drho_LR_dT * c.hL * ip_data.w_LS +
625 c.rhoLR * c.dh_L_dT * ip_data.w_LS;
626 // TODO (naumov) + k_over_mu_G * drho_GR_dT * b + k_over_mu_L *
627 // drho_LR_dT * b
628
629 // Derivatives of s_G * rho_C_GR_dot + s_L * rho_C_LR_dot abbreviated
630 // here with S_rho_C_eff.
631 double const s_L = ip_data.s_L;
632 double const s_G = 1. - ip_data.s_L;
633 double const rho_C_FR = s_G * ip_data.rhoCGR + s_L * ip_data.rhoCLR;
634 double const rho_W_FR = s_G * ip_data.rhoWGR + s_L * ip_data.rhoWLR;
635 // TODO (naumov) Extend for partially saturated media.
636 constexpr double drho_C_GR_dp_cap = 0;
637 if (dt == 0.)
638 {
639 ip_cv.dfC_3a_dp_GR = 0.;
640 ip_cv.dfC_3a_dp_cap = 0.;
641 ip_cv.dfC_3a_dT = 0.;
642 }
643 else
644 {
645 double const rho_C_GR_dot =
646 (ip_data.rhoCGR - ip_data.rhoCGR_prev) / dt;
647 double const rho_C_LR_dot =
648 (ip_data.rhoCLR - ip_data.rhoCLR_prev) / dt;
649 ip_cv.dfC_3a_dp_GR =
650 /*(ds_G_dp_GR = 0) * rho_C_GR_dot +*/ s_G * c.drho_C_GR_dp_GR /
651 dt +
652 /*(ds_L_dp_GR = 0) * rho_C_LR_dot +*/ s_L * c.drho_C_LR_dp_GR /
653 dt;
654 ip_cv.dfC_3a_dp_cap =
655 ds_G_dp_cap * rho_C_GR_dot + s_G * drho_C_GR_dp_cap / dt +
656 ip_cv.ds_L_dp_cap * rho_C_LR_dot - s_L * c.drho_C_LR_dp_LR / dt;
657 ip_cv.dfC_3a_dT =
658 s_G * c.drho_C_GR_dT / dt + s_L * c.drho_C_LR_dT / dt;
659 }
660
661 double const drho_C_FR_dp_GR =
662 /*(ds_G_dp_GR = 0) * ip_data.rhoCGR +*/ s_G * c.drho_C_GR_dp_GR +
663 /*(ds_L_dp_GR = 0) * ip_data.rhoCLR +*/ s_L * c.drho_C_LR_dp_GR;
664 ip_cv.dfC_4_MCpG_dp_GR = drho_C_FR_dp_GR *
665 (ip_data.alpha_B - ip_data.phi) *
666 ip_data.beta_p_SR;
667
668 double const drho_C_FR_dT = s_G * c.drho_C_GR_dT + s_L * c.drho_C_LR_dT;
669 ip_cv.dfC_4_MCpG_dT =
670 drho_C_FR_dT * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR
671#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
672 - rho_C_FR * ip_data.dphi_dT * ip_data.beta_p_SR
673#endif
674 ;
675
676 ip_cv.dfC_4_MCT_dT =
677 drho_C_FR_dT * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_T_SR
678#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
679 + rho_C_FR * (ip_data.alpha_B - ip_data.dphi_dT) * ip_data.beta_T_SR
680#endif
681 ;
682
683 ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_data.alpha_B;
684
685 ip_cv.dfC_2a_dp_GR = -ip_data.phi * c.drho_C_GR_dp_GR -
686 drho_C_FR_dp_GR * pCap *
687 (ip_data.alpha_B - ip_data.phi) *
688 ip_data.beta_p_SR;
689
690 double const drho_C_FR_dp_cap =
691 ds_G_dp_cap * ip_data.rhoCGR + s_G * drho_C_GR_dp_cap +
692 ip_cv.ds_L_dp_cap * ip_data.rhoCLR - s_L * c.drho_C_LR_dp_LR;
693
694 ip_cv.dfC_2a_dp_cap =
695 ip_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
696 drho_C_FR_dp_cap * pCap * (ip_data.alpha_B - ip_data.phi) *
697 ip_data.beta_p_SR +
698 rho_C_FR * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR;
699
700 ip_cv.dfC_2a_dT =
701#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
702 ip_data.dphi_dT * (ip_data.rhoCLR - ip_data.rhoCGR) +
703#endif
704 ip_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) -
705 drho_C_FR_dT * pCap * (ip_data.alpha_B - ip_data.phi) *
706 ip_data.beta_p_SR
707#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
708 + rho_C_FR * pCap * ip_data.dphi_dT * ip_data.beta_p_SR
709#endif
710 ;
711
712 ip_cv.dadvection_C_dp_GR = c.drho_C_GR_dp_GR * k_over_mu_G
713 // + rhoCGR * (dk_over_mu_G_dp_GR = 0)
714 // + rhoCLR * (dk_over_mu_L_dp_GR = 0)
715 + c.drho_C_LR_dp_GR * k_over_mu_L;
716
717 ip_cv.dadvection_C_dp_cap =
718 //(drho_C_GR_dp_cap = 0) * k_over_mu_G
719 ip_data.rhoCGR * ip_cv.dk_over_mu_G_dp_cap +
720 (-c.drho_C_LR_dp_LR) * k_over_mu_L +
721 ip_data.rhoCLR * ip_cv.dk_over_mu_L_dp_cap;
722
723 ip_cv.dfC_4_LCpG_dT =
724 c.drho_C_GR_dT * k_over_mu_G + c.drho_C_LR_dT * k_over_mu_L
725 // + ip_cv.ddiffusion_C_p_dT TODO (naumov)
726 ;
727
728 double const drho_W_FR_dp_GR =
729 /*(ds_G_dp_GR = 0) * ip_data.rhoWGR +*/ s_G * c.drho_W_GR_dp_GR +
730 /*(ds_L_dp_GR = 0) * ip_data.rhoWLR +*/ s_L * c.drho_W_LR_dp_GR;
731 double const drho_W_FR_dp_cap =
732 ds_G_dp_cap * ip_data.rhoWGR + s_G * c.drho_W_GR_dp_cap +
733 ip_cv.ds_L_dp_cap * ip_data.rhoWLR - s_L * c.drho_W_LR_dp_LR;
734 double const drho_W_FR_dT = s_G * c.drho_W_GR_dT + s_L * c.drho_W_LR_dT;
735
736 ip_cv.dfW_2a_dp_GR =
737 ip_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR);
738 ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap *
739 (ip_data.alpha_B - ip_data.phi) *
740 ip_data.beta_p_SR;
741 ip_cv.dfW_2a_dp_cap =
742 ip_data.phi * (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap);
743 ip_cv.dfW_2b_dp_cap =
744 drho_W_FR_dp_cap * pCap * (ip_data.alpha_B - ip_data.phi) *
745 ip_data.beta_p_SR +
746 rho_W_FR * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR;
747
748 ip_cv.dfW_2a_dT =
749#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
750 ip_data.dphi_dT * (ip_data.rhoWLR - ip_data.rhoWGR) +
751#endif
752 ip_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT);
753 ip_cv.dfW_2b_dT =
754 drho_W_FR_dT * pCap * (ip_data.alpha_B - ip_data.phi) *
755 ip_data.beta_p_SR
756#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
757 - rho_W_FR * pCap * ip_data.dphi_dT * ip_data.beta_p_SR
758#endif
759 ;
760
761 if (dt == 0.)
762 {
763 ip_cv.dfW_3a_dp_GR = 0.;
764 ip_cv.dfW_3a_dp_cap = 0.;
765 ip_cv.dfW_3a_dT = 0.;
766 }
767 else
768 {
769 double const rho_W_GR_dot =
770 (ip_data.rhoWGR - ip_data.rhoWGR_prev) / dt;
771 double const rho_W_LR_dot =
772 (ip_data.rhoWLR - ip_data.rhoWLR_prev) / dt;
773
774 ip_cv.dfW_3a_dp_GR =
775 /*(ds_G_dp_GR = 0) * rho_W_GR_dot +*/ s_G * c.drho_W_GR_dp_GR /
776 dt +
777 /*(ds_L_dp_GR = 0) * rho_W_LR_dot +*/ s_L * c.drho_W_LR_dp_GR /
778 dt;
779 ip_cv.dfW_3a_dp_cap =
780 ds_G_dp_cap * rho_W_GR_dot + s_G * c.drho_W_GR_dp_cap / dt +
781 ip_cv.ds_L_dp_cap * rho_W_LR_dot - s_L * c.drho_W_LR_dp_LR / dt;
782 ip_cv.dfW_3a_dT =
783 s_G * c.drho_W_GR_dT / dt + s_L * c.drho_W_LR_dT / dt;
784 }
785
786 ip_cv.dfW_4_LWpG_a_dp_GR = c.drho_W_GR_dp_GR * k_over_mu_G
787 // + rhoWGR * (dk_over_mu_G_dp_GR = 0)
788 + c.drho_W_LR_dp_GR * k_over_mu_L
789 // + rhoWLR * (dk_over_mu_L_dp_GR = 0)
790 ;
791 ip_cv.dfW_4_LWpG_a_dp_cap = c.drho_W_GR_dp_cap * k_over_mu_G +
792 ip_data.rhoWGR * ip_cv.dk_over_mu_G_dp_cap +
793 -c.drho_W_LR_dp_LR * k_over_mu_L +
794 ip_data.rhoWLR * ip_cv.dk_over_mu_L_dp_cap;
795
796 ip_cv.dfW_4_LWpG_a_dT =
797 c.drho_W_GR_dT * k_over_mu_G
798 //+ rhoWGR * (dk_over_mu_G_dT != 0 TODO for mu_G(T))
799 + c.drho_W_LR_dT * k_over_mu_L
800 //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_G(T))
801 ;
802
803 // TODO (naumov) for dxmW*/d* != 0
804 ip_cv.dfW_4_LWpG_d_dp_GR =
805 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
806 ip_cv.dfW_4_LWpG_d_dp_cap =
807 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
808 ip_cv.dfW_4_LWpG_d_dT =
809 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
810
811 ip_cv.dfW_4_LWpC_a_dp_GR = c.drho_W_LR_dp_GR * k_over_mu_L
812 //+ rhoWLR * (dk_over_mu_L_dp_GR = 0)
813 ;
814 ip_cv.dfW_4_LWpC_a_dp_cap = -c.drho_W_LR_dp_LR * k_over_mu_L +
815 ip_data.rhoWLR * ip_cv.dk_over_mu_L_dp_cap;
816 ip_cv.dfW_4_LWpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
817 //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
818 ;
819
820 // TODO (naumov) for dxmW*/d* != 0
821 ip_cv.dfW_4_LWpC_d_dp_GR =
822 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
823 ip_cv.dfW_4_LWpC_d_dp_cap =
824 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
825 ip_cv.dfW_4_LWpC_d_dT =
826 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
827
828 ip_cv.dfC_4_LCpC_a_dp_GR = c.drho_C_LR_dp_GR * k_over_mu_L
829 //+ rhoCLR * (dk_over_mu_L_dp_GR = 0)
830 ;
831 ip_cv.dfC_4_LCpC_a_dp_cap = -c.drho_C_LR_dp_LR * k_over_mu_L +
832 ip_data.rhoCLR * ip_cv.dk_over_mu_L_dp_cap;
833 ip_cv.dfC_4_LCpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
834 //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T))
835 ;
836
837 // TODO (naumov) for dxmW*/d* != 0
838 ip_cv.dfC_4_LCpC_d_dp_GR =
839 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
840 ip_cv.dfC_4_LCpC_d_dp_cap =
841 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
842 ip_cv.dfC_4_LCpC_d_dT =
843 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
844 }
845
846 return ip_constitutive_variables;
847}
848
849template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
850 int DisplacementDim>
851std::size_t TH2MLocalAssembler<
852 ShapeFunctionDisplacement, ShapeFunctionPressure,
853 DisplacementDim>::setIPDataInitialConditions(std::string const& name,
854 double const* values,
855 int const integration_order)
856{
857 if (integration_order !=
858 static_cast<int>(_integration_method.getIntegrationOrder()))
859 {
860 OGS_FATAL(
861 "Setting integration point initial conditions; The integration "
862 "order of the local assembler for element {:d} is different "
863 "from the integration order in the initial condition.",
864 _element.getID());
865 }
866
867 if (name == "sigma_ip")
868 {
869 if (_process_data.initial_stress != nullptr)
870 {
871 OGS_FATAL(
872 "Setting initial conditions for stress from integration "
873 "point data and from a parameter '{:s}' is not possible "
874 "simultaneously.",
875 _process_data.initial_stress->name);
876 }
877 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
878 values, _ip_data, &IpData::sigma_eff);
879 }
880
881 if (name == "saturation_ip")
882 {
883 return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
884 &IpData::s_L);
885 }
886 if (name == "swelling_stress_ip")
887 {
888 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
889 values, _ip_data, &IpData::sigma_sw);
890 }
891 if (name == "epsilon_ip")
892 {
893 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
894 values, _ip_data, &IpData::eps);
895 }
896 if (name.starts_with("material_state_variable_") && name.ends_with("_ip"))
897 {
898 std::string const variable_name = name.substr(24, name.size() - 24 - 3);
899 DBUG("Setting material state variable '{:s}'", variable_name);
900
901 // Using first ip data for solid material. TODO (naumov) move solid
902 // material into element, store only material state in IPs.
903 auto const& internal_variables =
904 _ip_data[0].solid_material.getInternalVariables();
905 if (auto const iv =
906 std::find_if(begin(internal_variables), end(internal_variables),
907 [&variable_name](auto const& iv)
908 { return iv.name == variable_name; });
909 iv != end(internal_variables))
910 {
911 return ProcessLib::setIntegrationPointDataMaterialStateVariables(
912 values, _ip_data, &IpData::material_state_variables,
913 iv->reference);
914 }
915 }
916 return 0;
917}
918
919template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
920 int DisplacementDim>
921void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
922 DisplacementDim>::
923 setInitialConditionsConcrete(std::vector<double> const& local_x_data,
924 double const t,
925 bool const /*use_monolithic_scheme*/,
926 int const /*process_id*/)
927{
928 [[maybe_unused]] auto const matrix_size =
929 gas_pressure_size + capillary_pressure_size + temperature_size +
930 displacement_size;
931
932 assert(local_x_data.size() == matrix_size);
933
934 auto const local_x = Eigen::Map<Eigen::VectorXd const>(local_x_data.data(),
935 local_x_data.size());
936 auto const capillary_pressure =
937 local_x.template segment<capillary_pressure_size>(
938 capillary_pressure_index);
939
940 auto const temperature =
941 local_x.template segment<temperature_size>(temperature_index);
942
943 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
944 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
945 auto const& solid_phase = medium.phase("Solid");
946
947 unsigned const n_integration_points =
948 _integration_method.getNumberOfPoints();
949
950 for (unsigned ip = 0; ip < n_integration_points; ip++)
951 {
952 MPL::VariableArray vars;
953
954 auto& ip_data = _ip_data[ip];
955 auto const& Np = ip_data.N_p;
956 auto const& NT = Np;
957 ParameterLib::SpatialPosition const pos{
958 std::nullopt, _element.getID(), ip,
959 MathLib::Point3d(
960 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
961 ShapeMatricesTypeDisplacement>(
962 _element, ip_data.N_u))};
963
964 double const pCap = Np.dot(capillary_pressure);
965 vars.capillary_pressure = pCap;
966
967 double const T = NT.dot(temperature);
968 vars.temperature = T;
969
970 auto& eps = ip_data.eps;
971
972 // Set volumetric strain rate for the general case without swelling.
973 vars.volumetric_strain = Invariants::trace(eps);
974
975 ip_data.s_L_prev =
976 medium.property(MPL::PropertyType::saturation)
977 .template value<double>(
978 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
979
980 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
981 // restart.
982 auto const C_el =
983 ip_data.computeElasticTangentStiffness(t, pos, dt, T, T);
984 auto& sigma_sw = ip_data.sigma_sw;
985 ip_data.eps_m_prev.noalias() =
986 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
987 ? eps + C_el.inverse() * sigma_sw
988 : eps;
989 }
990 // local_x_prev equal to local_x s.t. the local_x_dot is zero.
991 updateConstitutiveVariables(local_x, local_x, t, 0);
992
993 for (unsigned ip = 0; ip < n_integration_points; ip++)
994 {
995 auto& ip_data = _ip_data[ip];
996 ip_data.pushBackState();
997 }
998}
999
1000template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1001 int DisplacementDim>
1002void TH2MLocalAssembler<
1003 ShapeFunctionDisplacement, ShapeFunctionPressure,
1004 DisplacementDim>::assemble(double const t, double const dt,
1005 std::vector<double> const& local_x,
1006 std::vector<double> const& local_x_prev,
1007 std::vector<double>& local_M_data,
1008 std::vector<double>& local_K_data,
1009 std::vector<double>& local_rhs_data)
1010{
1011 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1012 temperature_size + displacement_size;
1013 assert(local_x.size() == matrix_size);
1014
1015 auto const capillary_pressure =
1016 Eigen::Map<VectorType<capillary_pressure_size> const>(
1017 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1018
1019 auto const capillary_pressure_prev =
1020 Eigen::Map<VectorType<capillary_pressure_size> const>(
1021 local_x_prev.data() + capillary_pressure_index,
1022 capillary_pressure_size);
1023
1024 // pointer to local_M_data vector
1025 auto local_M =
1026 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1027 local_M_data, matrix_size, matrix_size);
1028
1029 // pointer to local_K_data vector
1030 auto local_K =
1031 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1032 local_K_data, matrix_size, matrix_size);
1033
1034 // pointer to local_rhs_data vector
1035 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1036 local_rhs_data, matrix_size);
1037
1038 // component-formulation
1039 // W - liquid phase main component
1040 // C - gas phase main component
1041 // pointer-matrices to the mass matrix - C component equation
1042 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
1043 C_index, gas_pressure_index);
1044 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
1045 C_index, capillary_pressure_index);
1046 auto MCT = local_M.template block<C_size, temperature_size>(
1047 C_index, temperature_index);
1048 auto MCu = local_M.template block<C_size, displacement_size>(
1049 C_index, displacement_index);
1050
1051 // pointer-matrices to the stiffness matrix - C component equation
1052 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
1053 C_index, gas_pressure_index);
1054 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
1055 C_index, capillary_pressure_index);
1056 auto LCT = local_K.template block<C_size, temperature_size>(
1057 C_index, temperature_index);
1058
1059 // pointer-matrices to the mass matrix - W component equation
1060 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
1061 W_index, gas_pressure_index);
1062 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
1063 W_index, capillary_pressure_index);
1064 auto MWT = local_M.template block<W_size, temperature_size>(
1065 W_index, temperature_index);
1066 auto MWu = local_M.template block<W_size, displacement_size>(
1067 W_index, displacement_index);
1068
1069 // pointer-matrices to the stiffness matrix - W component equation
1070 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
1071 W_index, gas_pressure_index);
1072 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
1073 W_index, capillary_pressure_index);
1074 auto LWT = local_K.template block<W_size, temperature_size>(
1075 W_index, temperature_index);
1076
1077 // pointer-matrices to the mass matrix - temperature equation
1078 auto MTu = local_M.template block<temperature_size, displacement_size>(
1079 temperature_index, displacement_index);
1080
1081 // pointer-matrices to the stiffness matrix - temperature equation
1082 auto KTT = local_K.template block<temperature_size, temperature_size>(
1083 temperature_index, temperature_index);
1084
1085 // pointer-matrices to the stiffness matrix - displacement equation
1086 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1087 displacement_index, gas_pressure_index);
1088 auto KUpC =
1089 local_K.template block<displacement_size, capillary_pressure_size>(
1090 displacement_index, capillary_pressure_index);
1091
1092 // pointer-vectors to the right hand side terms - C-component equation
1093 auto fC = local_f.template segment<C_size>(C_index);
1094 // pointer-vectors to the right hand side terms - W-component equation
1095 auto fW = local_f.template segment<W_size>(W_index);
1096 // pointer-vectors to the right hand side terms - temperature equation
1097 auto fT = local_f.template segment<temperature_size>(temperature_index);
1098 // pointer-vectors to the right hand side terms - displacement equation
1099 auto fU = local_f.template segment<displacement_size>(displacement_index);
1100
1101 unsigned const n_integration_points =
1102 _integration_method.getNumberOfPoints();
1103
1104 auto const ip_constitutive_variables = updateConstitutiveVariables(
1105 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1106 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1107 local_x_prev.size()),
1108 t, dt);
1109
1110 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1111 {
1112 auto& ip = _ip_data[int_point];
1113 auto& ip_cv = ip_constitutive_variables[int_point];
1114
1115 auto const& Np = ip.N_p;
1116 auto const& NT = Np;
1117 auto const& Nu = ip.N_u;
1118 ParameterLib::SpatialPosition const pos{
1119 std::nullopt, _element.getID(), int_point,
1120 MathLib::Point3d(
1121 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1122 ShapeMatricesTypeDisplacement>(
1123 _element, Nu))};
1124
1125 auto const& NpT = Np.transpose().eval();
1126 auto const& NTT = NT.transpose().eval();
1127
1128 auto const& gradNp = ip.dNdx_p;
1129 auto const& gradNT = gradNp;
1130 auto const& gradNu = ip.dNdx_u;
1131
1132 auto const& gradNpT = gradNp.transpose().eval();
1133 auto const& gradNTT = gradNT.transpose().eval();
1134
1135 auto const& w = ip.integration_weight;
1136
1137 auto const& m = Invariants::identity2;
1138
1139 auto const mT = m.transpose().eval();
1140
1141 auto const x_coord =
1142 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1143 ShapeMatricesTypeDisplacement>(
1144 _element, Nu);
1145
1146 auto const Bu =
1147 LinearBMatrix::computeBMatrix<DisplacementDim,
1148 ShapeFunctionDisplacement::NPOINTS,
1149 typename BMatricesType::BMatrixType>(
1150 gradNu, Nu, x_coord, _is_axially_symmetric);
1151
1152 auto const BuT = Bu.transpose().eval();
1153
1154 double const pCap = Np.dot(capillary_pressure);
1155 double const pCap_prev = Np.dot(capillary_pressure_prev);
1156
1157 auto& beta_T_SR = ip.beta_T_SR;
1158
1159 auto const I =
1160 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1161
1162 const double sD_G = ip.diffusion_coefficient_vapour;
1163 const double sD_L = ip.diffusion_coefficient_solute;
1164
1165 GlobalDimMatrixType const D_C_G = sD_G * I;
1166 GlobalDimMatrixType const D_W_G = sD_G * I;
1167 GlobalDimMatrixType const D_C_L = sD_L * I;
1168 GlobalDimMatrixType const D_W_L = sD_L * I;
1169
1170 auto& k_S = ip.k_S;
1171
1172 auto& s_L = ip.s_L;
1173 auto const s_G = 1. - s_L;
1174 auto const s_L_dot = (s_L - ip.s_L_prev) / dt;
1175
1176 auto& alpha_B = ip.alpha_B;
1177 auto& beta_p_SR = ip.beta_p_SR;
1178
1179 auto const& b = _process_data.specific_body_force;
1180
1181 // porosity
1182 auto& phi = ip.phi;
1183
1184 // volume fraction
1185 auto const phi_G = s_G * phi;
1186 auto const phi_L = s_L * phi;
1187 auto const phi_S = 1. - phi;
1188
1189 // solid phase density
1190 auto& rho_SR = ip.rhoSR;
1191 // effective density
1192 auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
1193
1194 // abbreviations
1195 const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
1196 const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
1197
1198 // phase specific enthalpies
1199 auto& h_G = ip.h_G;
1200 auto& h_L = ip.h_L;
1201
1202 auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
1203 auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
1204 auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
1205 auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
1206
1207 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1208
1209 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1210
1211 GlobalDimMatrixType const k_over_mu_G = k_S * ip.k_rel_G / ip.muGR;
1212 GlobalDimMatrixType const k_over_mu_L = k_S * ip.k_rel_L / ip.muLR;
1213
1214 // ---------------------------------------------------------------------
1215 // C-component equation
1216 // ---------------------------------------------------------------------
1217
1218 MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1219 MCpC.noalias() -=
1220 NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1221
1222 if (_process_data.apply_mass_lumping)
1223 {
1224 if (pCap - pCap_prev != 0.) // avoid division by Zero
1225 {
1226 MCpC.noalias() +=
1227 NpT *
1228 (phi * (ip.rhoCLR - ip.rhoCGR) -
1229 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1230 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1231 }
1232 }
1233
1234 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1235 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1236
1237 using DisplacementDimMatrix =
1238 Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
1239
1240 DisplacementDimMatrix const advection_C_G = ip.rhoCGR * k_over_mu_G;
1241 DisplacementDimMatrix const advection_C_L = ip.rhoCLR * k_over_mu_L;
1242
1243 DisplacementDimMatrix const diffusion_CGpGR =
1244 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dpGR;
1245 DisplacementDimMatrix const diffusion_CLpGR =
1246 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dpGR;
1247
1248 DisplacementDimMatrix const diffusion_CGpCap =
1249 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dpCap;
1250 DisplacementDimMatrix const diffusion_CLpCap =
1251 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dpCap;
1252
1253 DisplacementDimMatrix const diffusion_CGT =
1254 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dT;
1255 DisplacementDimMatrix const diffusion_CLT =
1256 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dT;
1257
1258 DisplacementDimMatrix const advection_C = advection_C_G + advection_C_L;
1259 DisplacementDimMatrix const diffusion_C_pGR =
1260 diffusion_CGpGR + diffusion_CLpGR;
1261 DisplacementDimMatrix const diffusion_C_pCap =
1262 diffusion_CGpCap + diffusion_CLpCap;
1263
1264 DisplacementDimMatrix const diffusion_C_T =
1265 diffusion_CGT + diffusion_CLT;
1266
1267 LCpG.noalias() +=
1268 gradNpT * (advection_C + diffusion_C_pGR) * gradNp * w;
1269
1270 LCpC.noalias() +=
1271 gradNpT * (diffusion_C_pCap - advection_C_L) * gradNp * w;
1272
1273 LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
1274
1275 fC.noalias() += gradNpT *
1276 (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1277 b * w;
1278
1279 if (!_process_data.apply_mass_lumping)
1280 {
1281 fC.noalias() -= NpT *
1282 (phi * (ip.rhoCLR - ip.rhoCGR) -
1283 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1284 s_L_dot * w;
1285 }
1286 // fC_III
1287 fC.noalias() -=
1288 NpT * phi * (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
1289
1290 // ---------------------------------------------------------------------
1291 // W-component equation
1292 // ---------------------------------------------------------------------
1293
1294 MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1295 MWpC.noalias() -=
1296 NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1297
1298 if (_process_data.apply_mass_lumping)
1299 {
1300 if (pCap - pCap_prev != 0.) // avoid division by Zero
1301 {
1302 MWpC.noalias() +=
1303 NpT *
1304 (phi * (ip.rhoWLR - ip.rhoWGR) -
1305 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1306 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1307 }
1308 }
1309
1310 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1311
1312 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1313
1314 DisplacementDimMatrix const advection_W_G = ip.rhoWGR * k_over_mu_G;
1315 DisplacementDimMatrix const advection_W_L = ip.rhoWLR * k_over_mu_L;
1316
1317 DisplacementDimMatrix const diffusion_WGpGR =
1318 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR;
1319 DisplacementDimMatrix const diffusion_WLpGR =
1320 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpGR;
1321
1322 DisplacementDimMatrix const diffusion_WGpCap =
1323 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpCap;
1324 DisplacementDimMatrix const diffusion_WLpCap =
1325 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpCap;
1326
1327 DisplacementDimMatrix const diffusion_WGT =
1328 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT;
1329 DisplacementDimMatrix const diffusion_WLT =
1330 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT;
1331
1332 DisplacementDimMatrix const advection_W = advection_W_G + advection_W_L;
1333 DisplacementDimMatrix const diffusion_W_pGR =
1334 diffusion_WGpGR + diffusion_WLpGR;
1335 DisplacementDimMatrix const diffusion_W_pCap =
1336 diffusion_WGpCap + diffusion_WLpCap;
1337
1338 DisplacementDimMatrix const diffusion_W_T =
1339 diffusion_WGT + diffusion_WLT;
1340
1341 LWpG.noalias() +=
1342 gradNpT * (advection_W + diffusion_W_pGR) * gradNp * w;
1343
1344 LWpC.noalias() +=
1345 gradNpT * (diffusion_W_pCap - advection_W_L) * gradNp * w;
1346
1347 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1348
1349 fW.noalias() += gradNpT *
1350 (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1351 b * w;
1352
1353 if (!_process_data.apply_mass_lumping)
1354 {
1355 fW.noalias() -= NpT *
1356 (phi * (ip.rhoWLR - ip.rhoWGR) -
1357 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1358 s_L_dot * w;
1359 }
1360
1361 fW.noalias() -=
1362 NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1363
1364 // ---------------------------------------------------------------------
1365 // - temperature equation
1366 // ---------------------------------------------------------------------
1367
1368 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1369
1370 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1371
1372 fT.noalias() -= NTT * rho_u_eff_dot * w;
1373
1374 fT.noalias() +=
1375 gradNTT * (ip.rhoGR * h_G * ip.w_GS + ip.rhoLR * h_L * ip.w_LS) * w;
1376
1377 fT.noalias() +=
1378 gradNTT *
1379 (ip.rhoCGR * ip.h_CG * ip.d_CG + ip.rhoWGR * ip.h_WG * ip.d_WG) * w;
1380
1381 fT.noalias() +=
1382 NTT *
1383 (ip.rhoGR * ip.w_GS.transpose() + ip.rhoLR * ip.w_LS.transpose()) *
1384 b * w;
1385
1386 // ---------------------------------------------------------------------
1387 // - displacement equation
1388 // ---------------------------------------------------------------------
1389
1390 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1391
1392 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_s_L * m * Np) * w;
1393
1394 fU.noalias() -=
1395 (BuT * ip.sigma_eff - N_u_op(Nu).transpose() * rho * b) * w;
1396
1397 if (_process_data.apply_mass_lumping)
1398 {
1399 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1400 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1401 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1402 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1403 }
1404 } // int_point-loop
1405}
1406
1407// Assembles the local Jacobian matrix. So far, the linearisation of HT part is
1408// not considered as that in HT process.
1409template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1410 int DisplacementDim>
1411void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1412 DisplacementDim>::
1413 assembleWithJacobian(double const t, double const dt,
1414 std::vector<double> const& local_x,
1415 std::vector<double> const& local_x_prev,
1416 std::vector<double>& /*local_M_data*/,
1417 std::vector<double>& /*local_K_data*/,
1418 std::vector<double>& local_rhs_data,
1419 std::vector<double>& local_Jac_data)
1420{
1421 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1422 temperature_size + displacement_size;
1423 assert(local_x.size() == matrix_size);
1424
1425 auto const temperature = Eigen::Map<VectorType<temperature_size> const>(
1426 local_x.data() + temperature_index, temperature_size);
1427
1428 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size> const>(
1429 local_x.data() + gas_pressure_index, gas_pressure_size);
1430
1431 auto const capillary_pressure =
1432 Eigen::Map<VectorType<capillary_pressure_size> const>(
1433 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1434
1435 auto const displacement = Eigen::Map<VectorType<displacement_size> const>(
1436 local_x.data() + displacement_index, displacement_size);
1437
1438 auto const gas_pressure_prev =
1439 Eigen::Map<VectorType<gas_pressure_size> const>(
1440 local_x_prev.data() + gas_pressure_index, gas_pressure_size);
1441
1442 auto const capillary_pressure_prev =
1443 Eigen::Map<VectorType<capillary_pressure_size> const>(
1444 local_x_prev.data() + capillary_pressure_index,
1445 capillary_pressure_size);
1446
1447 auto const temperature_prev =
1448 Eigen::Map<VectorType<temperature_size> const>(
1449 local_x_prev.data() + temperature_index, temperature_size);
1450
1451 auto const displacement_prev =
1452 Eigen::Map<VectorType<displacement_size> const>(
1453 local_x_prev.data() + displacement_index, displacement_size);
1454
1455 auto local_Jac =
1456 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1457 local_Jac_data, matrix_size, matrix_size);
1458
1459 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1460 local_rhs_data, matrix_size);
1461
1462 // component-formulation
1463 // W - liquid phase main component
1464 // C - gas phase main component
1465
1466 // C component equation matrices
1467 MatrixType<C_size, gas_pressure_size> MCpG =
1468 MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1469 MatrixType<C_size, capillary_pressure_size> MCpC =
1470 MatrixType<C_size, capillary_pressure_size>::Zero(
1471 C_size, capillary_pressure_size);
1472 MatrixType<C_size, temperature_size> MCT =
1473 MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1474 MatrixType<C_size, displacement_size> MCu =
1475 MatrixType<C_size, displacement_size>::Zero(C_size, displacement_size);
1476
1477 MatrixType<C_size, gas_pressure_size> LCpG =
1478 MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1479 MatrixType<C_size, capillary_pressure_size> LCpC =
1480 MatrixType<C_size, capillary_pressure_size>::Zero(
1481 C_size, capillary_pressure_size);
1482 MatrixType<C_size, temperature_size> LCT =
1483 MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1484
1485 // mass matrix - W component equation
1486 MatrixType<W_size, gas_pressure_size> MWpG =
1487 MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1488 MatrixType<W_size, capillary_pressure_size> MWpC =
1489 MatrixType<W_size, capillary_pressure_size>::Zero(
1490 W_size, capillary_pressure_size);
1491 MatrixType<W_size, temperature_size> MWT =
1492 MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1493 MatrixType<W_size, displacement_size> MWu =
1494 MatrixType<W_size, displacement_size>::Zero(W_size, displacement_size);
1495
1496 // stiffness matrix - W component equation
1497 MatrixType<W_size, gas_pressure_size> LWpG =
1498 MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1499 MatrixType<W_size, capillary_pressure_size> LWpC =
1500 MatrixType<W_size, capillary_pressure_size>::Zero(
1501 W_size, capillary_pressure_size);
1502 MatrixType<W_size, temperature_size> LWT =
1503 MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1504
1505 // mass matrix - temperature equation
1506 MatrixType<temperature_size, displacement_size> MTu =
1507 MatrixType<temperature_size, displacement_size>::Zero(
1508 temperature_size, displacement_size);
1509
1510 // stiffness matrix - temperature equation
1511 MatrixType<temperature_size, temperature_size> KTT =
1512 MatrixType<temperature_size, temperature_size>::Zero(temperature_size,
1513 temperature_size);
1514
1515 // stiffness matrices - displacement equation coupling into pressures
1516 MatrixType<displacement_size, gas_pressure_size> KUpG =
1517 MatrixType<displacement_size, gas_pressure_size>::Zero(
1518 displacement_size, gas_pressure_size);
1519 MatrixType<displacement_size, capillary_pressure_size> KUpC =
1520 MatrixType<displacement_size, capillary_pressure_size>::Zero(
1521 displacement_size, capillary_pressure_size);
1522
1523 // pointer-vectors to the right hand side terms - C-component equation
1524 auto fC = local_f.template segment<C_size>(C_index);
1525 // pointer-vectors to the right hand side terms - W-component equation
1526 auto fW = local_f.template segment<W_size>(W_index);
1527 // pointer-vectors to the right hand side terms - temperature equation
1528 auto fT = local_f.template segment<temperature_size>(temperature_index);
1529 // pointer-vectors to the right hand side terms - displacement equation
1530 auto fU = local_f.template segment<displacement_size>(displacement_index);
1531
1532 unsigned const n_integration_points =
1533 _integration_method.getNumberOfPoints();
1534
1535 auto const ip_constitutive_variables = updateConstitutiveVariables(
1536 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1537 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1538 local_x_prev.size()),
1539 t, dt);
1540
1541 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1542 {
1543 auto& ip = _ip_data[int_point];
1544 auto& ip_cv = ip_constitutive_variables[int_point];
1545
1546 auto const& Np = ip.N_p;
1547 auto const& NT = Np;
1548 auto const& Nu = ip.N_u;
1549 ParameterLib::SpatialPosition const pos{
1550 std::nullopt, _element.getID(), int_point,
1551 MathLib::Point3d(
1552 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1553 ShapeMatricesTypeDisplacement>(
1554 _element, Nu))};
1555
1556 auto const& NpT = Np.transpose().eval();
1557 auto const& NTT = NT.transpose().eval();
1558
1559 auto const& gradNp = ip.dNdx_p;
1560 auto const& gradNT = gradNp;
1561 auto const& gradNu = ip.dNdx_u;
1562
1563 auto const& gradNpT = gradNp.transpose().eval();
1564 auto const& gradNTT = gradNT.transpose().eval();
1565
1566 auto const& w = ip.integration_weight;
1567
1568 auto const& m = Invariants::identity2;
1569 auto const mT = m.transpose().eval();
1570
1571 auto const x_coord =
1572 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1573 ShapeMatricesTypeDisplacement>(
1574 _element, Nu);
1575
1576 auto const Bu =
1577 LinearBMatrix::computeBMatrix<DisplacementDim,
1578 ShapeFunctionDisplacement::NPOINTS,
1579 typename BMatricesType::BMatrixType>(
1580 gradNu, Nu, x_coord, _is_axially_symmetric);
1581
1582 auto const BuT = Bu.transpose().eval();
1583
1584 double const div_u_dot =
1585 Invariants::trace(Bu * (displacement - displacement_prev) / dt);
1586
1587 double const pGR = Np.dot(gas_pressure);
1588 double const pCap = Np.dot(capillary_pressure);
1589 double const T = NT.dot(temperature);
1590
1591 GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
1592 GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
1593 GlobalDimVectorType const gradT = gradNT * temperature;
1594
1595 double const pGR_prev = Np.dot(gas_pressure_prev);
1596 double const pCap_prev = Np.dot(capillary_pressure_prev);
1597 double const T_prev = NT.dot(temperature_prev);
1598 auto& beta_T_SR = ip.beta_T_SR;
1599
1600 auto const I =
1601 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1602
1603 const double sD_G = ip.diffusion_coefficient_vapour;
1604 const double sD_L = ip.diffusion_coefficient_solute;
1605
1606 GlobalDimMatrixType const D_C_G = sD_G * I;
1607 GlobalDimMatrixType const D_W_G = sD_G * I;
1608 GlobalDimMatrixType const D_C_L = sD_L * I;
1609 GlobalDimMatrixType const D_W_L = sD_L * I;
1610
1611 auto& k_S = ip.k_S;
1612
1613 auto& s_L = ip.s_L;
1614 auto const s_G = 1. - s_L;
1615 auto const s_L_dot = (s_L - ip.s_L_prev) / dt;
1616
1617 auto& alpha_B = ip.alpha_B;
1618 auto& beta_p_SR = ip.beta_p_SR;
1619
1620 auto const& b = _process_data.specific_body_force;
1621
1622 // porosity
1623 auto& phi = ip.phi;
1624
1625 // volume fraction
1626 auto const phi_G = s_G * phi;
1627 auto const phi_L = s_L * phi;
1628 auto const phi_S = 1. - phi;
1629
1630 // solid phase density
1631 auto& rho_SR = ip.rhoSR;
1632 // effective density
1633 auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
1634
1635 // abbreviations
1636 const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
1637 const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
1638
1639 // phase specific enthalpies
1640 auto& h_G = ip.h_G;
1641 auto& h_L = ip.h_L;
1642
1643 auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
1644 auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
1645 auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
1646 auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
1647
1648 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1649
1650 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1651
1652 GlobalDimMatrixType const k_over_mu_G = k_S * ip.k_rel_G / ip.muGR;
1653 GlobalDimMatrixType const k_over_mu_L = k_S * ip.k_rel_L / ip.muLR;
1654
1655 // ---------------------------------------------------------------------
1656 // C-component equation
1657 // ---------------------------------------------------------------------
1658
1659 MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1660 MCpC.noalias() -=
1661 NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1662
1663 if (_process_data.apply_mass_lumping)
1664 {
1665 if (pCap - pCap_prev != 0.) // avoid division by Zero
1666 {
1667 MCpC.noalias() +=
1668 NpT *
1669 (phi * (ip.rhoCLR - ip.rhoCGR) -
1670 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1671 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1672 }
1673 }
1674
1675 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1676 // d (fC_4_MCT * T_dot)/d T
1677 local_Jac
1678 .template block<C_size, temperature_size>(C_index,
1679 temperature_index)
1680 .noalias() += NpT * ip_cv.dfC_4_MCT_dT * (T - T_prev) / dt * NT * w;
1681
1682 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1683 // d (fC_4_MCu * u_dot)/d T
1684 local_Jac
1685 .template block<C_size, temperature_size>(C_index,
1686 temperature_index)
1687 .noalias() += NpT * ip_cv.dfC_4_MCu_dT * div_u_dot * NT * w;
1688
1689 GlobalDimMatrixType const advection_C_G = ip.rhoCGR * k_over_mu_G;
1690 GlobalDimMatrixType const advection_C_L = ip.rhoCLR * k_over_mu_L;
1691 GlobalDimMatrixType const diffusion_C_G_p =
1692 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dpGR;
1693 GlobalDimMatrixType const diffusion_C_L_p =
1694 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dpLR;
1695 GlobalDimMatrixType const diffusion_C_G_T =
1696 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dT;
1697 GlobalDimMatrixType const diffusion_C_L_T =
1698 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dT;
1699
1700 GlobalDimMatrixType const advection_C = advection_C_G + advection_C_L;
1701 GlobalDimMatrixType const diffusion_C_p =
1702 diffusion_C_G_p + diffusion_C_L_p;
1703 GlobalDimMatrixType const diffusion_C_T =
1704 diffusion_C_G_T + diffusion_C_L_T;
1705
1706 LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
1707
1708 // d (fC_4_LCpG * grad p_GR)/d p_GR
1709 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1710 gradNpT *
1711 (ip_cv.dadvection_C_dp_GR
1712 // + ip_cv.ddiffusion_C_p_dp_GR TODO (naumov)
1713 ) *
1714 gradpGR * Np * w;
1715
1716 // d (fC_4_LCpG * grad p_GR)/d p_cap
1717 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1718 gradNpT *
1719 (ip_cv.dadvection_C_dp_cap
1720 // + ip_cv.ddiffusion_C_p_dp_GR TODO (naumov)
1721 ) *
1722 gradpGR * Np * w;
1723
1724 // d (fC_4_LCpG * grad p_GR)/d T
1725 local_Jac
1726 .template block<C_size, temperature_size>(C_index,
1727 temperature_index)
1728 .noalias() += gradNpT * ip_cv.dfC_4_LCpG_dT * gradpGR * NT * w;
1729
1730 // d (fC_4_MCpG * p_GR_dot)/d p_GR
1731 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1732 NpT * ip_cv.dfC_4_MCpG_dp_GR * (pGR - pGR_prev) / dt * Np * w;
1733
1734 // d (fC_4_MCpG * p_GR_dot)/d T
1735 local_Jac
1736 .template block<C_size, temperature_size>(C_index,
1737 temperature_index)
1738 .noalias() +=
1739 NpT * ip_cv.dfC_4_MCpG_dT * (pGR - pGR_prev) / dt * NT * w;
1740
1741 LCpC.noalias() -=
1742 gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
1743
1744 /* TODO (naumov) This part is not tested by any of the current ctests.
1745 // d (fC_4_LCpC * grad p_cap)/d p_GR
1746 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1747 gradNpT *
1748 (ip_cv.dfC_4_LCpC_a_dp_GR
1749 // + ip_cv.dfC_4_LCpC_d_dp_GR TODO (naumov)
1750 ) *
1751 gradpCap * Np * w;
1752 // d (fC_4_LCpC * grad p_cap)/d p_cap
1753 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1754 gradNpT *
1755 (ip_cv.dfC_4_LCpC_a_dp_cap
1756 // + ip_cv.dfC_4_LCpC_d_dp_cap TODO (naumov)
1757 ) *
1758 gradpCap * Np * w;
1759
1760 local_Jac
1761 .template block<C_size, temperature_size>(C_index,
1762 temperature_index)
1763 .noalias() += gradNpT *
1764 (ip_cv.dfC_4_LCpC_a_dT
1765 // + ip_cv.dfC_4_LCpC_d_dT TODO (naumov)
1766 ) *
1767 gradpCap * Np * w;
1768 */
1769
1770 LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w;
1771
1772 // fC_1
1773 fC.noalias() += gradNpT *
1774 (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1775 b * w;
1776
1777 if (!_process_data.apply_mass_lumping)
1778 {
1779 // fC_2 = \int a * s_L_dot
1780 auto const a = phi * (ip.rhoCLR - ip.rhoCGR) -
1781 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR;
1782 fC.noalias() -= NpT * a * s_L_dot * w;
1783
1784 local_Jac.template block<C_size, C_size>(C_index, C_index)
1785 .noalias() +=
1786 NpT *
1787 (ip_cv.dfC_2a_dp_GR * s_L_dot /*- a * (ds_L_dp_GR = 0) / dt*/) *
1788 Np * w;
1789
1790 local_Jac.template block<C_size, W_size>(C_index, W_index)
1791 .noalias() +=
1792 NpT *
1793 (ip_cv.dfC_2a_dp_cap * s_L_dot + a * ip_cv.ds_L_dp_cap / dt) *
1794 Np * w;
1795
1796 local_Jac
1797 .template block<C_size, temperature_size>(C_index,
1798 temperature_index)
1799 .noalias() += NpT * ip_cv.dfC_2a_dT * s_L_dot * NT * w;
1800 }
1801 {
1802 // fC_3 = \int phi * a
1803 double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot;
1804 fC.noalias() -= NpT * phi * a * w;
1805
1806 local_Jac.template block<C_size, C_size>(C_index, C_index)
1807 .noalias() += NpT * phi * ip_cv.dfC_3a_dp_GR * Np * w;
1808
1809 local_Jac.template block<C_size, W_size>(C_index, W_index)
1810 .noalias() += NpT * phi * ip_cv.dfC_3a_dp_cap * Np * w;
1811
1812 local_Jac
1813 .template block<C_size, temperature_size>(C_index,
1814 temperature_index)
1815 .noalias() += NpT *
1816 (
1817#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1818 ip.dphi_dT * a +
1819#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1820 phi * ip_cv.dfC_3a_dT) *
1821 NT * w;
1822 }
1823 // ---------------------------------------------------------------------
1824 // W-component equation
1825 // ---------------------------------------------------------------------
1826
1827 MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1828 MWpC.noalias() -=
1829 NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1830
1831 if (_process_data.apply_mass_lumping)
1832 {
1833 if (pCap - pCap_prev != 0.) // avoid division by Zero
1834 {
1835 MWpC.noalias() +=
1836 NpT *
1837 (phi * (ip.rhoWLR - ip.rhoWGR) -
1838 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1839 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1840 }
1841 }
1842
1843 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1844
1845 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1846
1847 GlobalDimMatrixType const advection_W_G = ip.rhoWGR * k_over_mu_G;
1848 GlobalDimMatrixType const advection_W_L = ip.rhoWLR * k_over_mu_L;
1849 GlobalDimMatrixType const diffusion_W_G_p =
1850 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR;
1851 GlobalDimMatrixType const diffusion_W_L_p =
1852 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpLR;
1853 GlobalDimMatrixType const diffusion_W_G_T =
1854 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT;
1855 GlobalDimMatrixType const diffusion_W_L_T =
1856 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT;
1857
1858 GlobalDimMatrixType const advection_W = advection_W_G + advection_W_L;
1859 GlobalDimMatrixType const diffusion_W_p =
1860 diffusion_W_G_p + diffusion_W_L_p;
1861 GlobalDimMatrixType const diffusion_W_T =
1862 diffusion_W_G_T + diffusion_W_L_T;
1863
1864 LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
1865
1866 // fW_4 LWpG' parts; LWpG = \int grad (a + d) grad
1867 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1868 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_GR + ip_cv.dfW_4_LWpG_d_dp_GR) *
1869 gradpGR * Np * w;
1870
1871 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1872 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_cap + ip_cv.dfW_4_LWpG_d_dp_cap) *
1873 gradpGR * Np * w;
1874
1875 local_Jac
1876 .template block<W_size, temperature_size>(W_index,
1877 temperature_index)
1878 .noalias() += gradNpT *
1879 (ip_cv.dfW_4_LWpG_a_dT + ip_cv.dfW_4_LWpG_d_dT) *
1880 gradpGR * NT * w;
1881
1882 LWpC.noalias() -=
1883 gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
1884
1885 // fW_4 LWp_cap' parts; LWpC = \int grad (a + d) grad
1886 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1887 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) *
1888 gradpCap * Np * w;
1889
1890 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1891 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) *
1892 gradpCap * Np * w;
1893
1894 local_Jac
1895 .template block<W_size, temperature_size>(W_index,
1896 temperature_index)
1897 .noalias() -= gradNpT *
1898 (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) *
1899 gradpCap * NT * w;
1900
1901 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1902
1903 // fW_1
1904 fW.noalias() += gradNpT *
1905 (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1906 b * w;
1907
1908 // fW_2 = \int (f - g) * s_L_dot
1909 if (!_process_data.apply_mass_lumping)
1910 {
1911 double const f = phi * (ip.rhoWLR - ip.rhoWGR);
1912 double const g = rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR;
1913
1914 fW.noalias() -= NpT * (f - g) * s_L_dot * w;
1915
1916 local_Jac.template block<W_size, C_size>(W_index, C_index)
1917 .noalias() += NpT * (ip_cv.dfW_2a_dp_GR - ip_cv.dfW_2b_dp_GR) *
1918 s_L_dot * Np * w;
1919
1920 // sign negated because of dp_cap = -dp_LR
1921 // TODO (naumov) Had to change the sign to get equal Jacobian WW
1922 // blocks in A2 Test. Where is the error?
1923 local_Jac.template block<W_size, W_size>(W_index, W_index)
1924 .noalias() +=
1925 NpT *
1926 ((ip_cv.dfW_2a_dp_cap - ip_cv.dfW_2b_dp_cap) * s_L_dot +
1927 (f - g) * ip_cv.ds_L_dp_cap / dt) *
1928 Np * w;
1929
1930 local_Jac
1931 .template block<W_size, temperature_size>(W_index,
1932 temperature_index)
1933 .noalias() +=
1934 NpT * (ip_cv.dfW_2a_dT - ip_cv.dfW_2b_dT) * s_L_dot * Np * w;
1935 }
1936
1937 // fW_3 = \int phi * a
1938 fW.noalias() -=
1939 NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1940
1941 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1942 NpT * phi * ip_cv.dfW_3a_dp_GR * Np * w;
1943
1944 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1945 NpT * phi * ip_cv.dfW_3a_dp_cap * Np * w;
1946
1947 local_Jac
1948 .template block<W_size, temperature_size>(W_index,
1949 temperature_index)
1950 .noalias() +=
1951 NpT *
1952 (
1953#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1954 ip.dphi_dT * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
1955#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1956 phi * ip_cv.dfW_3a_dT) *
1957 NT * w;
1958
1959 // ---------------------------------------------------------------------
1960 // - temperature equation
1961 // ---------------------------------------------------------------------
1962
1963 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1964
1965 // dfT_4/dp_GR
1966 // d (MTu * u_dot)/dp_GR
1967 local_Jac
1968 .template block<temperature_size, C_size>(temperature_index,
1969 C_index)
1970 .noalias() += NTT * ip_cv.drho_h_eff_dp_GR * div_u_dot * NT * w;
1971
1972 // dfT_4/dp_cap
1973 // d (MTu * u_dot)/dp_cap
1974 local_Jac
1975 .template block<temperature_size, W_size>(temperature_index,
1976 W_index)
1977 .noalias() -= NTT * ip_cv.drho_h_eff_dp_cap * div_u_dot * NT * w;
1978
1979 // dfT_4/dT
1980 // d (MTu * u_dot)/dT
1981 local_Jac
1982 .template block<temperature_size, temperature_size>(
1983 temperature_index, temperature_index)
1984 .noalias() += NTT * ip_cv.drho_h_eff_dT * div_u_dot * NT * w;
1985
1986 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1987
1988 // d KTT/dp_GR * T
1989 // TODO (naumov) always zero if lambda_xR have no derivatives wrt. p_GR.
1990 // dlambda_dp_GR =
1991 // (dphi_G_dp_GR = 0) * lambdaGR + phi_G * dlambda_GR_dp_GR +
1992 // (dphi_L_dp_GR = 0) * lambdaLR + phi_L * dlambda_LR_dp_GR +
1993 // (dphi_S_dp_GR = 0) * lambdaSR + phi_S * dlambda_SR_dp_GR +
1994 // = 0
1995 //
1996 // Since dlambda/dp_GR is 0 the derivative is omitted:
1997 // local_Jac
1998 // .template block<temperature_size, C_size>(temperature_index,
1999 // C_index)
2000 // .noalias() += gradNTT * dlambda_dp_GR * gradT * Np * w;
2001
2002 // d KTT/dp_cap * T
2003 local_Jac
2004 .template block<temperature_size, W_size>(temperature_index,
2005 W_index)
2006 .noalias() += gradNTT * ip_cv.dlambda_dp_cap * gradT * Np * w;
2007
2008 // d KTT/dT * T
2009 local_Jac
2010 .template block<temperature_size, temperature_size>(
2011 temperature_index, temperature_index)
2012 .noalias() += gradNTT * ip_cv.dlambda_dT * gradT * NT * w;
2013
2014 // fT_1
2015 fT.noalias() -= NTT * rho_u_eff_dot * w;
2016
2017 // dfT_1/dp_GR
2018 local_Jac
2019 .template block<temperature_size, C_size>(temperature_index,
2020 C_index)
2021 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_GR * Np * w;
2022
2023 // dfT_1/dp_cap
2024 local_Jac
2025 .template block<temperature_size, W_size>(temperature_index,
2026 W_index)
2027 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_cap * Np * w;
2028
2029 // dfT_1/dT
2030 // MTT
2031 local_Jac
2032 .template block<temperature_size, temperature_size>(
2033 temperature_index, temperature_index)
2034 .noalias() += NTT * ip_cv.drho_u_eff_dT / dt * NT * w;
2035
2036 // fT_2
2037 fT.noalias() +=
2038 gradNTT * (ip.rhoGR * h_G * ip.w_GS + ip.rhoLR * h_L * ip.w_LS) * w;
2039
2040 // dfT_2/dp_GR
2041 local_Jac
2042 .template block<temperature_size, C_size>(temperature_index,
2043 C_index)
2044 .noalias() -=
2045 // dfT_2/dp_GR first part
2046 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_Npart * Np * w +
2047 // dfT_2/dp_GR second part
2048 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart * gradNp * w;
2049
2050 // dfT_2/dp_cap
2051 local_Jac
2052 .template block<temperature_size, W_size>(temperature_index,
2053 W_index)
2054 .noalias() -=
2055 // first part of dfT_2/dp_cap
2056 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_Npart) * Np * w +
2057 // second part of dfT_2/dp_cap
2058 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart) * gradNp * w;
2059
2060 // dfT_2/dT
2061 local_Jac
2062 .template block<temperature_size, temperature_size>(
2063 temperature_index, temperature_index)
2064 .noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
2065
2066 // fT_3
2067 fT.noalias() +=
2068 NTT *
2069 (ip.rhoGR * ip.w_GS.transpose() + ip.rhoLR * ip.w_LS.transpose()) *
2070 b * w;
2071
2072 fT.noalias() +=
2073 gradNTT *
2074 (ip.rhoCGR * ip.h_CG * ip.d_CG + ip.rhoWGR * ip.h_WG * ip.d_WG) * w;
2075
2076 // ---------------------------------------------------------------------
2077 // - displacement equation
2078 // ---------------------------------------------------------------------
2079
2080 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
2081
2082 // dfU_2/dp_GR = dKUpG/dp_GR * p_GR + KUpG. The former is zero, the
2083 // latter is handled below.
2084
2085 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_s_L * m * Np) * w;
2086
2087 // dfU_2/dp_cap = dKUpC/dp_cap * p_cap + KUpC. The former is handled
2088 // here, the latter below.
2089 local_Jac
2090 .template block<displacement_size, W_size>(displacement_index,
2091 W_index)
2092 .noalias() += BuT * alpha_B * ip_cv.dchi_ds_L * ip_cv.ds_L_dp_cap *
2093 pCap * m * Np * w;
2094
2095 local_Jac
2096 .template block<displacement_size, displacement_size>(
2097 displacement_index, displacement_index)
2098 .noalias() += BuT * ip_cv.C * Bu * w;
2099
2100 // fU_1
2101 fU.noalias() -=
2102 (BuT * ip.sigma_eff - N_u_op(Nu).transpose() * rho * b) * w;
2103
2104 // KuT
2105 local_Jac
2106 .template block<displacement_size, temperature_size>(
2107 displacement_index, temperature_index)
2108 .noalias() -= BuT * (ip_cv.C * ip.alpha_T_SR) * NT * w;
2109
2110 /* TODO (naumov) Test with gravity needed to check this Jacobian part.
2111 local_Jac
2112 .template block<displacement_size, temperature_size>(
2113 displacement_index, temperature_index)
2114 .noalias() += N_u_op(Nu).transpose() * ip_cv.drho_dT * b *
2115 N_u_op(Nu).transpose() * w;
2116 */
2117
2118 if (_process_data.apply_mass_lumping)
2119 {
2120 MCpG = MCpG.colwise().sum().eval().asDiagonal();
2121 MCpC = MCpC.colwise().sum().eval().asDiagonal();
2122 MWpG = MWpG.colwise().sum().eval().asDiagonal();
2123 MWpC = MWpC.colwise().sum().eval().asDiagonal();
2124 }
2125 } // int_point-loop
2126
2127 // --- Gas ---
2128 // fC_4
2129 fC.noalias() -= LCpG * gas_pressure + LCpC * capillary_pressure +
2130 LCT * temperature +
2131 MCpG * (gas_pressure - gas_pressure_prev) / dt +
2132 MCpC * (capillary_pressure - capillary_pressure_prev) / dt +
2133 MCT * (temperature - temperature_prev) / dt +
2134 MCu * (displacement - displacement_prev) / dt;
2135
2136 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
2137 LCpG + MCpG / dt;
2138 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
2139 LCpC + MCpC / dt;
2140 local_Jac
2141 .template block<C_size, temperature_size>(C_index, temperature_index)
2142 .noalias() += LCT + MCT / dt;
2143 local_Jac
2144 .template block<C_size, displacement_size>(C_index, displacement_index)
2145 .noalias() += MCu / dt;
2146
2147 // --- Capillary pressure ---
2148 // fW_4
2149 fW.noalias() -= LWpG * gas_pressure + LWpC * capillary_pressure +
2150 LWT * temperature +
2151 MWpG * (gas_pressure - gas_pressure_prev) / dt +
2152 MWpC * (capillary_pressure - capillary_pressure_prev) / dt +
2153 MWT * (temperature - temperature_prev) / dt +
2154 MWu * (displacement - displacement_prev) / dt;
2155
2156 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
2157 LWpC + MWpC / dt;
2158 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
2159 LWpG + MWpG / dt;
2160 local_Jac
2161 .template block<W_size, temperature_size>(W_index, temperature_index)
2162 .noalias() += LWT + MWT / dt;
2163 local_Jac
2164 .template block<W_size, displacement_size>(W_index, displacement_index)
2165 .noalias() += MWu / dt;
2166
2167 // --- Temperature ---
2168 // fT_4
2169 fT.noalias() -=
2170 KTT * temperature + MTu * (displacement - displacement_prev) / dt;
2171
2172 local_Jac
2173 .template block<temperature_size, temperature_size>(temperature_index,
2174 temperature_index)
2175 .noalias() += KTT;
2176 local_Jac
2177 .template block<temperature_size, displacement_size>(temperature_index,
2178 displacement_index)
2179 .noalias() += MTu / dt;
2180
2181 // --- Displacement ---
2182 // fU_2
2183 fU.noalias() -= KUpG * gas_pressure + KUpC * capillary_pressure;
2184
2185 local_Jac
2186 .template block<displacement_size, C_size>(displacement_index, C_index)
2187 .noalias() += KUpG;
2188 local_Jac
2189 .template block<displacement_size, W_size>(displacement_index, W_index)
2190 .noalias() += KUpC;
2191}
2192
2193template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
2194 int DisplacementDim>
2195std::vector<double> const& TH2MLocalAssembler<
2196 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2197 getIntPtDarcyVelocityGas(
2198 const double /*t*/,
2199 std::vector<GlobalVector*> const& /*x*/,
2200 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
2201 std::vector<double>& cache) const
2202{
2203 unsigned const n_integration_points =
2204 _integration_method.getNumberOfPoints();
2205
2206 cache.clear();
2207 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2208 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2209 cache, DisplacementDim, n_integration_points);
2210
2211 for (unsigned ip = 0; ip < n_integration_points; ip++)
2212 {
2213 cache_matrix.col(ip).noalias() = _ip_data[ip].w_GS;
2214 }
2215
2216 return cache;
2217}
2218
2219template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
2220 int DisplacementDim>
2221std::vector<double> const& TH2MLocalAssembler<
2222 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2223 getIntPtDarcyVelocityLiquid(
2224 const double /*t*/,
2225 std::vector<GlobalVector*> const& /*x*/,
2226 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
2227 std::vector<double>& cache) const
2228{
2229 unsigned const n_integration_points =
2230 _integration_method.getNumberOfPoints();
2231
2232 cache.clear();
2233 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2234 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2235 cache, DisplacementDim, n_integration_points);
2236
2237 for (unsigned ip = 0; ip < n_integration_points; ip++)
2238 {
2239 cache_matrix.col(ip).noalias() = _ip_data[ip].w_LS;
2240 }
2241
2242 return cache;
2243}
2244
2245template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
2246 int DisplacementDim>
2247std::vector<double> const& TH2MLocalAssembler<
2248 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2249 getIntPtDiffusionVelocityVapourGas(
2250 const double /*t*/,
2251 std::vector<GlobalVector*> const& /*x*/,
2252 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
2253 std::vector<double>& cache) const
2254{
2255 unsigned const n_integration_points =
2256 _integration_method.getNumberOfPoints();
2257
2258 cache.clear();
2259 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2260 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2261 cache, DisplacementDim, n_integration_points);
2262
2263 for (unsigned ip = 0; ip < n_integration_points; ip++)
2264 {
2265 cache_matrix.col(ip).noalias() = _ip_data[ip].d_WG;
2266 }
2267
2268 return cache;
2269}
2270
2271template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
2272 int DisplacementDim>
2273std::vector<double> const& TH2MLocalAssembler<
2274 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2275 getIntPtDiffusionVelocityGasGas(
2276 const double /*t*/,
2277 std::vector<GlobalVector*> const& /*x*/,
2278 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
2279 std::vector<double>& cache) const
2280{
2281 unsigned const n_integration_points =
2282 _integration_method.getNumberOfPoints();
2283
2284 cache.clear();
2285 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2286 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2287 cache, DisplacementDim, n_integration_points);
2288
2289 for (unsigned ip = 0; ip < n_integration_points; ip++)
2290 {
2291 cache_matrix.col(ip).noalias() = _ip_data[ip].d_CG;
2292 }
2293
2294 return cache;
2295}
2296
2297template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
2298 int DisplacementDim>
2299std::vector<double> const& TH2MLocalAssembler<
2300 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2301 getIntPtDiffusionVelocitySoluteLiquid(
2302 const double /*t*/,
2303 std::vector<GlobalVector*> const& /*x*/,
2304 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
2305 std::vector<double>& cache) const
2306{
2307 unsigned const n_integration_points =
2308 _integration_method.getNumberOfPoints();
2309
2310 cache.clear();
2311 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2312 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2313 cache, DisplacementDim, n_integration_points);
2314
2315 for (unsigned ip = 0; ip < n_integration_points; ip++)
2316 {
2317 cache_matrix.col(ip).noalias() = _ip_data[ip].d_CL;
2318 }
2319
2320 return cache;
2321}
2322
2323template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
2324 int DisplacementDim>
2325std::vector<double> const& TH2MLocalAssembler<
2326 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2327 getIntPtDiffusionVelocityLiquidLiquid(
2328 const double /*t*/,
2329 std::vector<GlobalVector*> const& /*x*/,
2330 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
2331 std::vector<double>& cache) const
2332{
2333 unsigned const n_integration_points =
2334 _integration_method.getNumberOfPoints();
2335
2336 cache.clear();
2337 auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
2338 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2339 cache, DisplacementDim, n_integration_points);
2340
2341 for (unsigned ip = 0; ip < n_integration_points; ip++)
2342 {
2343 cache_matrix.col(ip).noalias() = _ip_data[ip].d_WL;
2344 }
2345
2346 return cache;
2347}
2348
2349template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
2350 int DisplacementDim>
2351void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
2352 DisplacementDim>::
2353 computeSecondaryVariableConcrete(double const t, double const dt,
2354 Eigen::VectorXd const& local_x,
2355 Eigen::VectorXd const& local_x_prev)
2356{
2357 auto const gas_pressure =
2358 local_x.template segment<gas_pressure_size>(gas_pressure_index);
2359 auto const capillary_pressure =
2360 local_x.template segment<capillary_pressure_size>(
2361 capillary_pressure_index);
2362 auto const liquid_pressure = (gas_pressure - capillary_pressure).eval();
2363
2364 NumLib::interpolateToHigherOrderNodes<
2365 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2366 DisplacementDim>(_element, _is_axially_symmetric, gas_pressure,
2367 *_process_data.gas_pressure_interpolated);
2368
2369 NumLib::interpolateToHigherOrderNodes<
2370 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2371 DisplacementDim>(_element, _is_axially_symmetric, capillary_pressure,
2372 *_process_data.capillary_pressure_interpolated);
2373
2374 NumLib::interpolateToHigherOrderNodes<
2375 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2376 DisplacementDim>(_element, _is_axially_symmetric, liquid_pressure,
2377 *_process_data.liquid_pressure_interpolated);
2378
2379 auto const temperature =
2380 local_x.template segment<temperature_size>(temperature_index);
2381
2382 NumLib::interpolateToHigherOrderNodes<
2383 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
2384 DisplacementDim>(_element, _is_axially_symmetric, temperature,
2385 *_process_data.temperature_interpolated);
2386
2387 unsigned const n_integration_points =
2388 _integration_method.getNumberOfPoints();
2389
2390 double saturation_avg = 0;
2391
2392 updateConstitutiveVariables(local_x, local_x_prev, t, dt);
2393
2394 for (unsigned ip = 0; ip < n_integration_points; ip++)
2395 {
2396 saturation_avg += _ip_data[ip].s_L;
2397 }
2398 saturation_avg /= n_integration_points;
2399 (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
2400}
2401
2402} // namespace TH2M
2403} // namespace ProcessLib
OGS_FATAL
#define OGS_FATAL(...)
Definition: Error.h:26
DBUG
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:30
MaterialPropertyLib::VariableArray::mechanical_strain
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
Definition: VariableType.h:182
MaterialPropertyLib::VariableArray::total_stress
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > total_stress
Definition: VariableType.h:201
MathLib::WeightedPoint::getWeight
double getWeight() const
Definition: WeightedPoint.h:80
MeshLib::Element::getID
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
NumLib::GenericIntegrationMethod::getWeightedPoint
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
Definition: GenericIntegrationMethod.h:36
ProcessLib::BMatrixPolicyType::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
ProcessLib::TH2M::TH2MLocalAssembler::MatrixType
typename ShapeMatricesTypePressure::template MatrixType< M, N > MatrixType
Definition: TH2MFEM.h:65
ProcessLib::TH2M::TH2MLocalAssembler::SymmetricTensor
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
Definition: TH2MFEM.h:75
ProcessLib::TH2M::TH2MLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition: TH2MFEM.h:57
ProcessLib::TH2M::TH2MLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition: TH2MFEM.h:496
ProcessLib::TH2M::TH2MLocalAssembler::_process_data
TH2MProcessData< DisplacementDim > & _process_data
Definition: TH2MFEM.h:481
ProcessLib::TH2M::TH2MLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition: TH2MFEM.h:54
ProcessLib::TH2M::TH2MLocalAssembler::GlobalDimMatrixType
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
Definition: TH2MFEM.h:68
ProcessLib::TH2M::TH2MLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition: TH2MFEM.h:489
ProcessLib::TH2M::TH2MLocalAssembler::GlobalDimVectorType
typename ShapeMatricesTypePressure::GlobalDimVectorType GlobalDimVectorType
Definition: TH2MFEM.h:71
ProcessLib::TH2M::TH2MLocalAssembler::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition: TH2MFEM.h:491
MaterialLib::Solids::selectSolidConstitutiveRelation
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Definition: SelectSolidConstitutiveRelation.h:36
MaterialPropertyLib::formEigenTensor< 3 >
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
MathLib::KelvinVector::KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
MathLib::KelvinVector::kelvinVectorToSymmetricTensor
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition: KelvinVector.cpp:142
MathLib::KelvinVector::tensorToKelvin
KelvinVectorType< DisplacementDim > tensorToKelvin(Eigen::Matrix< double, 3, 3 > const &m)
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
NumLib::interpolateXCoordinate
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition: InitShapeMatrices.h:71
NumLib::interpolateToHigherOrderNodes
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
Definition: Interpolation.h:107
NumLib::initShapeMatrices
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)
Definition: InitShapeMatrices.h:53
NumLib::interpolateCoordinates
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition: InitShapeMatrices.h:82
ProcessLib::LinearBMatrix::computeBMatrix
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:42
ProcessLib::setIntegrationPointDataMaterialStateVariables
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)
Definition: SetOrGetIntegrationPointData.h:270
ProcessLib::setIPDataInitialConditions
void setIPDataInitialConditions(std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > const &_integration_point_writer, MeshLib::Properties const &mesh_properties, LocalAssemblersVector &local_assemblers, bool const remove_name_suffix=false)
Definition: SetIPDataInitialConditions.h:34
ProcessLib::setIntegrationPointScalarData
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
Definition: SetOrGetIntegrationPointData.h:215
ProcessLib::HydroMechanics::SecondaryData::N_u
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
Definition: HydroMechanicsFEM.h:161