OGS
TH2MFEM-impl.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/LU>
14
25#include "TH2MProcessData.h"
26
27namespace ProcessLib
28{
29namespace TH2M
30{
31namespace MPL = MaterialPropertyLib;
32
33template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
34 int DisplacementDim>
35TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
36 DisplacementDim>::
37 TH2MLocalAssembler(
38 MeshLib::Element const& e,
39 std::size_t const /*local_matrix_size*/,
40 NumLib::GenericIntegrationMethod const& integration_method,
41 bool const is_axially_symmetric,
43 : LocalAssemblerInterface<DisplacementDim>(
44 e, integration_method, is_axially_symmetric, process_data)
45{
46 unsigned const n_integration_points =
48
49 _ip_data.resize(n_integration_points);
50 _secondary_data.N_u.resize(n_integration_points);
51
52 auto const shape_matrices_u =
53 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
55 DisplacementDim>(e, is_axially_symmetric,
57
58 auto const shape_matrices_p =
59 NumLib::initShapeMatrices<ShapeFunctionPressure,
60 ShapeMatricesTypePressure, DisplacementDim>(
61 e, is_axially_symmetric, this->integration_method_);
62
63 for (unsigned ip = 0; ip < n_integration_points; ip++)
64 {
65 auto& ip_data = _ip_data[ip];
66 auto const& sm_u = shape_matrices_u[ip];
67 ip_data.integration_weight =
69 sm_u.integralMeasure * sm_u.detJ;
70
71 ip_data.N_u = sm_u.N;
72 ip_data.dNdx_u = sm_u.dNdx;
73
74 ip_data.N_p = shape_matrices_p[ip].N;
75 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
76
77 _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
78 }
79}
80
81template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
82 int DisplacementDim>
83std::tuple<
84 std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>,
85 std::vector<ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>>>
86TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
87 DisplacementDim>::
88 updateConstitutiveVariables(
89 Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
90 double const t, double const dt,
92 models)
93{
94 [[maybe_unused]] auto const matrix_size =
95 gas_pressure_size + capillary_pressure_size + temperature_size +
96 displacement_size;
97
98 assert(local_x.size() == matrix_size);
99
100 auto const gas_pressure =
101 local_x.template segment<gas_pressure_size>(gas_pressure_index);
102 auto const capillary_pressure =
103 local_x.template segment<capillary_pressure_size>(
104 capillary_pressure_index);
105
106 auto const temperature =
107 local_x.template segment<temperature_size>(temperature_index);
108 auto const temperature_prev =
109 local_x_prev.template segment<temperature_size>(temperature_index);
110
111 auto const displacement =
112 local_x.template segment<displacement_size>(displacement_index);
113 auto const displacement_prev =
114 local_x_prev.template segment<displacement_size>(displacement_index);
115
116 auto const& medium =
117 *this->process_data_.media_map.getMedium(this->element_.getID());
118 ConstitutiveRelations::MediaData media_data{medium};
119
120 unsigned const n_integration_points =
121 this->integration_method_.getNumberOfPoints();
122
123 std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>
124 ip_constitutive_data(n_integration_points);
125 std::vector<ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>>
126 ip_constitutive_variables(n_integration_points);
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 auto& ip_cd = ip_constitutive_data[ip];
133 auto& ip_out = this->output_data_[ip];
134 auto& current_state = this->current_states_[ip];
135 auto& prev_state = this->prev_states_[ip];
136
137 auto const& Np = ip_data.N_p;
138 auto const& NT = Np;
139 auto const& Nu = ip_data.N_u;
140 auto const& gradNu = ip_data.dNdx_u;
141 auto const& gradNp = ip_data.dNdx_p;
143 std::nullopt, this->element_.getID(), ip,
145 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
147 this->element_, Nu))};
148 auto const x_coord =
149 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
151 this->element_, Nu);
152
153 double const T = NT.dot(temperature);
154 double const T_prev = NT.dot(temperature_prev);
155 ConstitutiveRelations::TemperatureData const T_data{T, T_prev};
157 Np.dot(gas_pressure)};
159 Np.dot(capillary_pressure)};
161 this->process_data_.reference_temperature(t, pos)[0]};
163 grad_p_GR{gradNp * gas_pressure};
165 DisplacementDim> const grad_p_cap{gradNp * capillary_pressure};
167 grad_T{gradNp * temperature};
168
169 // medium properties
170 models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
171 ip_cv.C_el_data);
172
173 models.biot_model.eval({pos, t, dt}, media_data, ip_cv.biot_data);
174
175 auto const Bu =
176 LinearBMatrix::computeBMatrix<DisplacementDim,
177 ShapeFunctionDisplacement::NPOINTS,
179 gradNu, Nu, x_coord, this->is_axially_symmetric_);
180
181 ip_out.eps_data.eps.noalias() = Bu * displacement;
182 models.S_L_model.eval({pos, t, dt}, media_data, pCap_data,
183 current_state.S_L_data);
184
185 models.chi_S_L_model.eval({pos, t, dt}, media_data,
186 current_state.S_L_data, ip_cv.chi_S_L);
187
188 // solid phase compressibility
189 models.beta_p_SR_model.eval({pos, t, dt}, ip_cv.biot_data,
190 ip_cv.C_el_data, ip_cv.beta_p_SR);
191
192 // If there is swelling stress rate, compute swelling stress.
193 models.swelling_model.eval(
194 {pos, t, dt}, media_data, ip_cv.C_el_data, current_state.S_L_data,
195 prev_state.S_L_data, prev_state.swelling_data,
196 current_state.swelling_data, ip_cv.swelling_data);
197
198 // solid phase linear thermal expansion coefficient
199 models.s_therm_exp_model.eval({pos, t, dt}, media_data, T_data, T0,
200 ip_cv.s_therm_exp_data);
201
202 models.mechanical_strain_model.eval(
203 T_data, ip_cv.s_therm_exp_data, ip_out.eps_data,
204 Bu * displacement_prev, prev_state.mechanical_strain_data,
205 ip_cv.swelling_data, current_state.mechanical_strain_data);
206
207 models.s_mech_model.eval(
208 {pos, t, dt}, T_data, current_state.mechanical_strain_data,
209 prev_state.mechanical_strain_data, prev_state.eff_stress_data,
210 current_state.eff_stress_data, this->material_states_[ip],
211 ip_cd.s_mech_data, ip_cv.equivalent_plastic_strain_data);
212
213 models.total_stress_model.eval(current_state.eff_stress_data,
214 ip_cv.biot_data, ip_cv.chi_S_L, pGR_data,
215 pCap_data, ip_cv.total_stress_data);
216
217 models.permeability_model.eval(
218 {pos, t, dt}, media_data, current_state.S_L_data, pCap_data, T_data,
219 ip_cv.total_stress_data, ip_out.eps_data,
220 ip_cv.equivalent_plastic_strain_data, ip_out.permeability_data);
221
222 models.pure_liquid_density_model.eval({pos, t, dt}, media_data,
223 pGR_data, pCap_data, T_data,
224 current_state.rho_W_LR);
225
227 {pos, t, dt}, media_data, pGR_data, pCap_data, T_data,
228 current_state.rho_W_LR, ip_out.fluid_enthalpy_data,
229 ip_out.mass_mole_fractions_data, ip_out.fluid_density_data,
230 ip_out.vapour_pressure_data, current_state.constituent_density_data,
231 ip_cv.phase_transition_data);
232
233 models.viscosity_model.eval({pos, t, dt}, media_data, T_data,
234 ip_out.mass_mole_fractions_data,
235 ip_cv.viscosity_data);
236
237 models.porosity_model.eval({pos, t, dt}, media_data,
238#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
239 ip_cv.biot_data, ip_out.eps_data,
240 ip_cv.s_therm_exp_data,
241#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
242 ip_out.porosity_data);
243
244 models.solid_density_model.eval({pos, t, dt}, media_data, T_data,
245#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
246 ip_cv.biot_data, ip_out.eps_data,
247 ip_cv.s_therm_exp_data,
248#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
249 ip_out.solid_density_data);
250
251 models.solid_heat_capacity_model.eval({pos, t, dt}, media_data, T_data,
252 ip_cv.solid_heat_capacity_data);
253
254 models.thermal_conductivity_model.eval(
255 {pos, t, dt}, media_data, T_data, ip_out.porosity_data,
256 current_state.S_L_data, ip_cv.thermal_conductivity_data);
257
258 models.advection_model.eval(current_state.constituent_density_data,
259 ip_out.permeability_data,
260 current_state.rho_W_LR,
261 ip_cv.viscosity_data,
262 ip_cv.advection_data);
263
264 models.gravity_model.eval(
265 ip_out.fluid_density_data,
266 ip_out.porosity_data,
267 current_state.S_L_data,
268 ip_out.solid_density_data,
270 this->process_data_.specific_body_force},
271 ip_cv.volumetric_body_force);
272
273 models.diffusion_velocity_model.eval(grad_p_cap,
274 grad_p_GR,
275 ip_out.mass_mole_fractions_data,
276 ip_cv.phase_transition_data,
277 ip_out.porosity_data,
278 current_state.S_L_data,
279 grad_T,
280 ip_out.diffusion_velocity_data);
281
282 models.solid_enthalpy_model.eval(ip_cv.solid_heat_capacity_data, T_data,
283 ip_out.solid_enthalpy_data);
284
285 models.internal_energy_model.eval(ip_out.fluid_density_data,
286 ip_cv.phase_transition_data,
287 ip_out.porosity_data,
288 current_state.S_L_data,
289 ip_out.solid_density_data,
290 ip_out.solid_enthalpy_data,
291 current_state.internal_energy_data);
292
294 ip_out.fluid_density_data,
295 ip_out.fluid_enthalpy_data,
296 ip_out.porosity_data,
297 current_state.S_L_data,
298 ip_out.solid_density_data,
299 ip_out.solid_enthalpy_data,
300 ip_cv.effective_volumetric_enthalpy_data);
301
302 models.fC_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data,
303 ip_cv.fC_1);
304
305 if (!this->process_data_.apply_mass_lumping)
306 {
307 models.fC_2a_model.eval(ip_cv.biot_data,
308 pCap_data,
309 current_state.constituent_density_data,
310 ip_out.porosity_data,
311 current_state.S_L_data,
312 ip_cv.beta_p_SR,
313 ip_cv.fC_2a);
314 }
315 models.fC_3a_model.eval(dt,
316 current_state.constituent_density_data,
317 prev_state.constituent_density_data,
318 current_state.S_L_data,
319 ip_cv.fC_3a);
320
321 models.fC_4_LCpG_model.eval(ip_cv.advection_data,
322 ip_out.fluid_density_data,
323 ip_cv.phase_transition_data,
324 ip_out.porosity_data,
325 current_state.S_L_data,
326 ip_cv.fC_4_LCpG);
327
328 models.fC_4_LCpC_model.eval(ip_cv.advection_data,
329 ip_out.fluid_density_data,
330 ip_cv.phase_transition_data,
331 ip_out.porosity_data,
332 current_state.S_L_data,
333 ip_cv.fC_4_LCpC);
334
335 models.fC_4_LCT_model.eval(ip_out.fluid_density_data,
336 ip_cv.phase_transition_data,
337 ip_out.porosity_data,
338 current_state.S_L_data,
339 ip_cv.fC_4_LCT);
340
341 models.fC_4_MCpG_model.eval(ip_cv.biot_data,
342 current_state.constituent_density_data,
343 ip_out.porosity_data,
344 current_state.S_L_data,
345 ip_cv.beta_p_SR,
346 ip_cv.fC_4_MCpG);
347
348 models.fC_4_MCpC_model.eval(ip_cv.biot_data,
349 pCap_data,
350 current_state.constituent_density_data,
351 ip_out.porosity_data,
352 prev_state.S_L_data,
353 current_state.S_L_data,
354 ip_cv.beta_p_SR,
355 ip_cv.fC_4_MCpC);
356
357 models.fC_4_MCT_model.eval(ip_cv.biot_data,
358 current_state.constituent_density_data,
359 ip_out.porosity_data,
360 current_state.S_L_data,
361 ip_cv.s_therm_exp_data,
362 ip_cv.fC_4_MCT);
363
364 models.fC_4_MCu_model.eval(ip_cv.biot_data,
365 current_state.constituent_density_data,
366 current_state.S_L_data,
367 ip_cv.fC_4_MCu);
368
369 models.fW_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data,
370 ip_cv.fW_1);
371
372 if (!this->process_data_.apply_mass_lumping)
373 {
374 models.fW_2_model.eval(ip_cv.biot_data,
375 pCap_data,
376 current_state.constituent_density_data,
377 ip_out.porosity_data,
378 current_state.rho_W_LR,
379 current_state.S_L_data,
380 ip_cv.beta_p_SR,
381 ip_cv.fW_2);
382 }
383 models.fW_3a_model.eval(dt,
384 current_state.constituent_density_data,
385 prev_state.constituent_density_data,
386 prev_state.rho_W_LR,
387 current_state.rho_W_LR,
388 current_state.S_L_data,
389 ip_cv.fW_3a);
390
391 models.fW_4_LWpG_model.eval(ip_cv.advection_data,
392 ip_out.fluid_density_data,
393 ip_cv.phase_transition_data,
394 ip_out.porosity_data,
395 current_state.S_L_data,
396 ip_cv.fW_4_LWpG);
397
398 models.fW_4_LWpC_model.eval(ip_cv.advection_data,
399 ip_out.fluid_density_data,
400 ip_cv.phase_transition_data,
401 ip_out.porosity_data,
402 current_state.S_L_data,
403 ip_cv.fW_4_LWpC);
404
405 models.fW_4_LWT_model.eval(ip_out.fluid_density_data,
406 ip_cv.phase_transition_data,
407 ip_out.porosity_data,
408 current_state.S_L_data,
409 ip_cv.fW_4_LWT);
410
411 models.fW_4_MWpG_model.eval(ip_cv.biot_data,
412 current_state.constituent_density_data,
413 ip_out.porosity_data,
414 current_state.rho_W_LR,
415 current_state.S_L_data,
416 ip_cv.beta_p_SR,
417 ip_cv.fW_4_MWpG);
418
419 models.fW_4_MWpC_model.eval(ip_cv.biot_data,
420 pCap_data,
421 current_state.constituent_density_data,
422 ip_out.porosity_data,
423 prev_state.S_L_data,
424 current_state.rho_W_LR,
425 current_state.S_L_data,
426 ip_cv.beta_p_SR,
427 ip_cv.fW_4_MWpC);
428
429 models.fW_4_MWT_model.eval(ip_cv.biot_data,
430 current_state.constituent_density_data,
431 ip_out.porosity_data,
432 current_state.rho_W_LR,
433 current_state.S_L_data,
434 ip_cv.s_therm_exp_data,
435 ip_cv.fW_4_MWT);
436
437 models.fW_4_MWu_model.eval(ip_cv.biot_data,
438 current_state.constituent_density_data,
439 current_state.rho_W_LR,
440 current_state.S_L_data,
441 ip_cv.fW_4_MWu);
442
443 models.fT_1_model.eval(dt,
444 current_state.internal_energy_data,
445 prev_state.internal_energy_data,
446 ip_cv.fT_1);
447
448 // ---------------------------------------------------------------------
449 // Derivatives for Jacobian
450 // ---------------------------------------------------------------------
451
452 models.darcy_velocity_model.eval(
453 grad_p_cap,
454 ip_out.fluid_density_data,
455 grad_p_GR,
456 ip_out.permeability_data,
458 this->process_data_.specific_body_force},
459 ip_cv.viscosity_data,
460 ip_out.darcy_velocity_data);
461
462 models.fT_2_model.eval(ip_out.darcy_velocity_data,
463 ip_out.fluid_density_data,
464 ip_out.fluid_enthalpy_data,
465 ip_cv.fT_2);
466
467 models.fT_3_model.eval(
468 current_state.constituent_density_data,
469 ip_out.darcy_velocity_data,
470 ip_out.diffusion_velocity_data,
471 ip_out.fluid_density_data,
472 ip_cv.phase_transition_data,
474 this->process_data_.specific_body_force},
475 ip_cv.fT_3);
476
477 models.fu_2_KupC_model.eval(ip_cv.biot_data, ip_cv.chi_S_L,
478 ip_cv.fu_2_KupC);
479 }
480
481 return {ip_constitutive_data, ip_constitutive_variables};
482}
483
484template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
485 int DisplacementDim>
486std::vector<ConstitutiveRelations::DerivativesData<DisplacementDim>>
487TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
488 DisplacementDim>::
489 updateConstitutiveVariablesDerivatives(
490 Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
491 double const t, double const dt,
492 std::vector<
494 ip_constitutive_data,
495 std::vector<
497 ip_constitutive_variables,
499 models)
500{
501 [[maybe_unused]] auto const matrix_size =
502 gas_pressure_size + capillary_pressure_size + temperature_size +
503 displacement_size;
504
505 assert(local_x.size() == matrix_size);
506
507 auto const temperature =
508 local_x.template segment<temperature_size>(temperature_index);
509 auto const temperature_prev =
510 local_x_prev.template segment<temperature_size>(temperature_index);
511
512 auto const capillary_pressure =
513 local_x.template segment<capillary_pressure_size>(
514 capillary_pressure_index);
515
516 auto const& medium =
517 *this->process_data_.media_map.getMedium(this->element_.getID());
518 ConstitutiveRelations::MediaData media_data{medium};
519
520 unsigned const n_integration_points =
521 this->integration_method_.getNumberOfPoints();
522
523 std::vector<ConstitutiveRelations::DerivativesData<DisplacementDim>>
524 ip_d_data(n_integration_points);
525
526 for (unsigned ip = 0; ip < n_integration_points; ip++)
527 {
528 auto const& ip_data = _ip_data[ip];
529 auto& ip_dd = ip_d_data[ip];
530 auto const& ip_cd = ip_constitutive_data[ip];
531 auto const& ip_cv = ip_constitutive_variables[ip];
532 auto const& ip_out = this->output_data_[ip];
533 auto const& current_state = this->current_states_[ip];
534 auto const& prev_state = this->prev_states_[ip];
535
536 auto const& Nu = ip_data.N_u;
537 auto const& Np = ip_data.N_p;
538 auto const& NT = Np;
539
541 std::nullopt, this->element_.getID(), ip,
543 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
545 this->element_, Nu))};
546
547 double const T = NT.dot(temperature);
548 double const T_prev = NT.dot(temperature_prev);
549 ConstitutiveRelations::TemperatureData const T_data{T, T_prev};
551 Np.dot(capillary_pressure)};
552
553 models.S_L_model.dEval({pos, t, dt}, media_data, pCap_data,
554 ip_dd.dS_L_dp_cap);
555
556 models.advection_model.dEval(current_state.constituent_density_data,
557 ip_out.permeability_data,
558 ip_cv.viscosity_data,
559 ip_dd.dS_L_dp_cap,
560 ip_cv.phase_transition_data,
561 ip_dd.advection_d_data);
562
563 models.porosity_model.dEval(
564 {pos, t, dt}, media_data, ip_out.porosity_data, ip_dd.dS_L_dp_cap,
565#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
566 ip_cv.biot_data, ip_out.eps_data, ip_cv.s_therm_exp_data,
567#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
568 ip_dd.porosity_d_data);
569
570 models.thermal_conductivity_model.dEval(
571 {pos, t, dt}, media_data, T_data, ip_out.porosity_data,
572 ip_dd.porosity_d_data, current_state.S_L_data,
573 ip_dd.thermal_conductivity_d_data);
574
575 models.solid_density_model.dEval({pos, t, dt}, media_data, T_data,
576#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
577 ip_cv.biot_data, ip_out.eps_data,
578 ip_cv.s_therm_exp_data,
579#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
580 ip_dd.solid_density_d_data);
581
583 ip_out.fluid_density_data,
584 ip_cv.phase_transition_data,
585 ip_out.porosity_data,
586 ip_dd.porosity_d_data,
587 current_state.S_L_data,
588 ip_out.solid_density_data,
589 ip_dd.solid_density_d_data,
590 ip_out.solid_enthalpy_data,
591 ip_cv.solid_heat_capacity_data,
592 ip_dd.effective_volumetric_internal_energy_d_data);
593
595 ip_out.fluid_density_data,
596 ip_out.fluid_enthalpy_data,
597 ip_cv.phase_transition_data,
598 ip_out.porosity_data,
599 ip_dd.porosity_d_data,
600 current_state.S_L_data,
601 ip_out.solid_density_data,
602 ip_dd.solid_density_d_data,
603 ip_out.solid_enthalpy_data,
604 ip_cv.solid_heat_capacity_data,
605 ip_dd.effective_volumetric_enthalpy_d_data);
606 if (!this->process_data_.apply_mass_lumping)
607 {
608 models.fC_2a_model.dEval(ip_cv.biot_data,
609 pCap_data,
610 current_state.constituent_density_data,
611 ip_cv.phase_transition_data,
612 ip_out.porosity_data,
613 ip_dd.porosity_d_data,
614 current_state.S_L_data,
615 ip_dd.dS_L_dp_cap,
616 ip_cv.beta_p_SR,
617 ip_dd.dfC_2a);
618 }
619 models.fC_3a_model.dEval(dt,
620 current_state.constituent_density_data,
621 prev_state.constituent_density_data,
622 ip_cv.phase_transition_data,
623 current_state.S_L_data,
624 ip_dd.dS_L_dp_cap,
625 ip_dd.dfC_3a);
626
627 models.fC_4_LCpG_model.dEval(ip_out.permeability_data,
628 ip_cv.viscosity_data,
629 ip_cv.phase_transition_data,
630 ip_dd.advection_d_data,
631 ip_dd.dfC_4_LCpG);
632
633 models.fC_4_LCpC_model.dEval(current_state.constituent_density_data,
634 ip_out.permeability_data,
635 ip_cv.phase_transition_data,
636 ip_dd.dS_L_dp_cap,
637 ip_cv.viscosity_data,
638 ip_dd.dfC_4_LCpC);
639
640 models.fC_4_MCpG_model.dEval(ip_cv.biot_data,
641 current_state.constituent_density_data,
642 ip_cv.phase_transition_data,
643 ip_out.porosity_data,
644 ip_dd.porosity_d_data,
645 current_state.S_L_data,
646 ip_cv.beta_p_SR,
647 ip_dd.dfC_4_MCpG);
648
649 models.fC_4_MCT_model.dEval(ip_cv.biot_data,
650 current_state.constituent_density_data,
651 ip_cv.phase_transition_data,
652 ip_out.porosity_data,
653 ip_dd.porosity_d_data,
654 current_state.S_L_data,
655 ip_cv.s_therm_exp_data,
656 ip_dd.dfC_4_MCT);
657
658 models.fC_4_MCu_model.dEval(ip_cv.biot_data,
659 ip_cv.phase_transition_data,
660 current_state.S_L_data,
661 ip_dd.dfC_4_MCu);
662
663 if (!this->process_data_.apply_mass_lumping)
664 {
665 models.fW_2_model.dEval(ip_cv.biot_data,
666 pCap_data,
667 current_state.constituent_density_data,
668 ip_cv.phase_transition_data,
669 ip_out.porosity_data,
670 ip_dd.porosity_d_data,
671 current_state.rho_W_LR,
672 current_state.S_L_data,
673 ip_dd.dS_L_dp_cap,
674 ip_cv.beta_p_SR,
675 ip_dd.dfW_2);
676 }
677
678 models.fW_3a_model.dEval(dt,
679 current_state.constituent_density_data,
680 ip_cv.phase_transition_data,
681 prev_state.constituent_density_data,
682 prev_state.rho_W_LR,
683 current_state.rho_W_LR,
684 current_state.S_L_data,
685 ip_dd.dS_L_dp_cap,
686 ip_dd.dfW_3a);
687
688 models.fW_4_LWpG_model.dEval(current_state.constituent_density_data,
689 ip_out.permeability_data,
690 ip_cv.phase_transition_data,
691 current_state.rho_W_LR,
692 ip_dd.dS_L_dp_cap,
693 ip_cv.viscosity_data,
694 ip_dd.dfW_4_LWpG);
695
696 models.fW_4_LWpC_model.dEval(ip_cv.advection_data,
697 ip_out.fluid_density_data,
698 ip_out.permeability_data,
699 ip_cv.phase_transition_data,
700 ip_out.porosity_data,
701 current_state.rho_W_LR,
702 current_state.S_L_data,
703 ip_dd.dS_L_dp_cap,
704 ip_cv.viscosity_data,
705 ip_dd.dfW_4_LWpC);
706
707 models.fT_1_model.dEval(
708 dt, ip_dd.effective_volumetric_internal_energy_d_data, ip_dd.dfT_1);
709
710 models.fT_2_model.dEval(
711 ip_out.darcy_velocity_data,
712 ip_out.fluid_density_data,
713 ip_out.fluid_enthalpy_data,
714 ip_out.permeability_data,
715 ip_cv.phase_transition_data,
717 this->process_data_.specific_body_force},
718 ip_cv.viscosity_data,
719 ip_dd.dfT_2);
720
721 models.fu_1_KuT_model.dEval(ip_cd.s_mech_data, ip_cv.s_therm_exp_data,
722 ip_dd.dfu_1_KuT);
723
724 models.fu_2_KupC_model.dEval(ip_cv.biot_data,
725 ip_cv.chi_S_L,
726 pCap_data,
727 ip_dd.dS_L_dp_cap,
728 ip_dd.dfu_2_KupC);
729 }
730
731 return ip_d_data;
732}
733
734template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
735 int DisplacementDim>
736std::size_t TH2MLocalAssembler<
737 ShapeFunctionDisplacement, ShapeFunctionPressure,
738 DisplacementDim>::setIPDataInitialConditions(std::string_view name,
739 double const* values,
740 int const integration_order)
741{
742 if (integration_order !=
743 static_cast<int>(this->integration_method_.getIntegrationOrder()))
744 {
745 OGS_FATAL(
746 "Setting integration point initial conditions; The integration "
747 "order of the local assembler for element {:d} is different "
748 "from the integration order in the initial condition.",
749 this->element_.getID());
750 }
751
752 if (name == "sigma" && this->process_data_.initial_stress.value)
753 {
754 OGS_FATAL(
755 "Setting initial conditions for stress from integration "
756 "point data and from a parameter '{:s}' is not possible "
757 "simultaneously.",
758 this->process_data_.initial_stress.value->name);
759 }
760
761 if (name.starts_with("material_state_variable_"))
762 {
763 name.remove_prefix(24);
764 DBUG("Setting material state variable '{:s}'", name);
765
766 auto const& internal_variables =
767 this->solid_material_.getInternalVariables();
768 if (auto const iv = std::find_if(
769 begin(internal_variables), end(internal_variables),
770 [&name](auto const& iv) { return iv.name == name; });
771 iv != end(internal_variables))
772 {
773 DBUG("Setting material state variable '{:s}'", name);
775 values, this->material_states_,
777 DisplacementDim>::material_state_variables,
778 iv->reference);
779 }
780
781 WARN(
782 "Could not find variable {:s} in solid material model's "
783 "internal variables.",
784 name);
785 return 0;
786 }
787
788 // TODO this logic could be pulled out of the local assembler into the
789 // process. That might lead to a slightly better performance due to less
790 // string comparisons.
791 return ProcessLib::Reflection::reflectSetIPData<DisplacementDim>(
792 name, values, this->current_states_);
793}
794
795template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
796 int DisplacementDim>
797void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
798 DisplacementDim>::
799 setInitialConditionsConcrete(Eigen::VectorXd const local_x,
800 double const t,
801 int const /*process_id*/)
802{
803 [[maybe_unused]] auto const matrix_size =
804 gas_pressure_size + capillary_pressure_size + temperature_size +
805 displacement_size;
806
807 assert(local_x.size() == matrix_size);
808
809 auto const capillary_pressure =
810 local_x.template segment<capillary_pressure_size>(
811 capillary_pressure_index);
812
813 auto const p_GR =
814 local_x.template segment<gas_pressure_size>(gas_pressure_index);
815
816 auto const temperature =
817 local_x.template segment<temperature_size>(temperature_index);
818
819 auto const displacement =
820 local_x.template segment<displacement_size>(displacement_index);
821
822 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
823 auto const& medium =
824 *this->process_data_.media_map.getMedium(this->element_.getID());
825 auto const& solid_phase = medium.phase("Solid");
826
828 this->solid_material_, *this->process_data_.phase_transition_model_};
829
830 unsigned const n_integration_points =
831 this->integration_method_.getNumberOfPoints();
832
833 for (unsigned ip = 0; ip < n_integration_points; ip++)
834 {
836
837 auto& ip_data = _ip_data[ip];
838 auto& ip_out = this->output_data_[ip];
839 auto& prev_state = this->prev_states_[ip];
840 auto const& Np = ip_data.N_p;
841 auto const& NT = Np;
842 auto const& Nu = ip_data.N_u;
843 auto const& gradNu = ip_data.dNdx_u;
844 auto const x_coord =
845 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
847 this->element_, Nu);
849 std::nullopt, this->element_.getID(), ip,
851 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
853 this->element_, ip_data.N_u))};
854
855 double const pCap = Np.dot(capillary_pressure);
856 vars.capillary_pressure = pCap;
857
858 double const T = NT.dot(temperature);
860 T, T}; // T_prev = T in initialization.
861 vars.temperature = T;
862
863 auto const Bu =
864 LinearBMatrix::computeBMatrix<DisplacementDim,
865 ShapeFunctionDisplacement::NPOINTS,
867 gradNu, Nu, x_coord, this->is_axially_symmetric_);
868
869 auto& eps = ip_out.eps_data.eps;
870 eps.noalias() = Bu * displacement;
871
872 // Set volumetric strain rate for the general case without swelling.
873 vars.volumetric_strain = Invariants::trace(eps);
874
875 double const S_L =
876 medium.property(MPL::PropertyType::saturation)
877 .template value<double>(
878 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
879 this->prev_states_[ip].S_L_data->S_L = S_L;
880
881 // TODO (naumov) Double computation of C_el might be avoided if
882 // updateConstitutiveVariables is called before. But it might interfere
883 // with eps_m initialization.
885 C_el_data;
886 models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
887 C_el_data);
888 auto const& C_el = C_el_data.stiffness_tensor;
889
890 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
891 // restart.
892 auto const& sigma_sw = this->current_states_[ip].swelling_data.sigma_sw;
893 prev_state.mechanical_strain_data->eps_m.noalias() =
894 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
895 ? eps + C_el.inverse() * sigma_sw
896 : eps;
897
898 if (this->process_data_.initial_stress.isTotalStress())
899 {
900 auto const alpha_b =
901 medium.property(MPL::PropertyType::biot_coefficient)
902 .template value<double>(vars, pos, t, 0.0 /*dt*/);
903
904 vars.liquid_saturation = S_L;
905 double const bishop =
906 medium.property(MPL::PropertyType::bishops_effective_stress)
907 .template value<double>(vars, pos, t, 0.0 /*dt*/);
908
909 this->current_states_[ip].eff_stress_data.sigma.noalias() +=
910 alpha_b * Np.dot(p_GR - bishop * capillary_pressure) *
911 Invariants::identity2;
912 this->prev_states_[ip].eff_stress_data =
913 this->current_states_[ip].eff_stress_data;
914 }
915 }
916
917 // local_x_prev equal to local_x s.t. the local_x_dot is zero.
918 updateConstitutiveVariables(local_x, local_x, t, 0, models);
919
920 for (unsigned ip = 0; ip < n_integration_points; ip++)
921 {
922 this->material_states_[ip].pushBackState();
923 this->prev_states_[ip] = this->current_states_[ip];
924 }
925}
926
927template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
928 int DisplacementDim>
930 ShapeFunctionDisplacement, ShapeFunctionPressure,
931 DisplacementDim>::assemble(double const t, double const dt,
932 std::vector<double> const& local_x,
933 std::vector<double> const& local_x_prev,
934 std::vector<double>& local_M_data,
935 std::vector<double>& local_K_data,
936 std::vector<double>& local_rhs_data)
937{
938 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
939 temperature_size + displacement_size;
940 assert(local_x.size() == matrix_size);
941
942 auto const capillary_pressure =
943 Eigen::Map<VectorType<capillary_pressure_size> const>(
944 local_x.data() + capillary_pressure_index, capillary_pressure_size);
945
946 auto const capillary_pressure_prev =
947 Eigen::Map<VectorType<capillary_pressure_size> const>(
948 local_x_prev.data() + capillary_pressure_index,
949 capillary_pressure_size);
950
951 // pointer to local_M_data vector
952 auto local_M =
953 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
954 local_M_data, matrix_size, matrix_size);
955
956 // pointer to local_K_data vector
957 auto local_K =
958 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
959 local_K_data, matrix_size, matrix_size);
960
961 // pointer to local_rhs_data vector
962 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
963 local_rhs_data, matrix_size);
964
965 // component-formulation
966 // W - liquid phase main component
967 // C - gas phase main component
968 // pointer-matrices to the mass matrix - C component equation
969 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
970 C_index, gas_pressure_index);
971 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
972 C_index, capillary_pressure_index);
973 auto MCT = local_M.template block<C_size, temperature_size>(
974 C_index, temperature_index);
975 auto MCu = local_M.template block<C_size, displacement_size>(
976 C_index, displacement_index);
977
978 // pointer-matrices to the stiffness matrix - C component equation
979 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
980 C_index, gas_pressure_index);
981 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
982 C_index, capillary_pressure_index);
983 auto LCT = local_K.template block<C_size, temperature_size>(
984 C_index, temperature_index);
985
986 // pointer-matrices to the mass matrix - W component equation
987 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
988 W_index, gas_pressure_index);
989 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
990 W_index, capillary_pressure_index);
991 auto MWT = local_M.template block<W_size, temperature_size>(
992 W_index, temperature_index);
993 auto MWu = local_M.template block<W_size, displacement_size>(
994 W_index, displacement_index);
995
996 // pointer-matrices to the stiffness matrix - W component equation
997 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
998 W_index, gas_pressure_index);
999 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
1000 W_index, capillary_pressure_index);
1001 auto LWT = local_K.template block<W_size, temperature_size>(
1002 W_index, temperature_index);
1003
1004 // pointer-matrices to the mass matrix - temperature equation
1005 auto MTu = local_M.template block<temperature_size, displacement_size>(
1006 temperature_index, displacement_index);
1007
1008 // pointer-matrices to the stiffness matrix - temperature equation
1009 auto KTT = local_K.template block<temperature_size, temperature_size>(
1010 temperature_index, temperature_index);
1011
1012 // pointer-matrices to the stiffness matrix - displacement equation
1013 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1014 displacement_index, gas_pressure_index);
1015 auto KUpC =
1016 local_K.template block<displacement_size, capillary_pressure_size>(
1017 displacement_index, capillary_pressure_index);
1018
1019 // pointer-vectors to the right hand side terms - C-component equation
1020 auto fC = local_f.template segment<C_size>(C_index);
1021 // pointer-vectors to the right hand side terms - W-component equation
1022 auto fW = local_f.template segment<W_size>(W_index);
1023 // pointer-vectors to the right hand side terms - temperature equation
1024 auto fT = local_f.template segment<temperature_size>(temperature_index);
1025 // pointer-vectors to the right hand side terms - displacement equation
1026 auto fU = local_f.template segment<displacement_size>(displacement_index);
1027
1028 unsigned const n_integration_points =
1029 this->integration_method_.getNumberOfPoints();
1030
1032 this->solid_material_, *this->process_data_.phase_transition_model_};
1033
1034 auto const [ip_constitutive_data, ip_constitutive_variables] =
1035 updateConstitutiveVariables(
1036 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1037 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1038 local_x_prev.size()),
1039 t, dt, models);
1040
1041 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1042 {
1043 auto& ip = _ip_data[int_point];
1044 auto& ip_cv = ip_constitutive_variables[int_point];
1045 auto& ip_out = this->output_data_[int_point];
1046 auto& current_state = this->current_states_[int_point];
1047 auto const& prev_state = this->prev_states_[int_point];
1048
1049 auto const& Np = ip.N_p;
1050 auto const& NT = Np;
1051 auto const& Nu = ip.N_u;
1053 std::nullopt, this->element_.getID(), int_point,
1055 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1057 this->element_, Nu))};
1058
1059 auto const& NpT = Np.transpose().eval();
1060 auto const& NTT = NT.transpose().eval();
1061
1062 auto const& gradNp = ip.dNdx_p;
1063 auto const& gradNT = gradNp;
1064 auto const& gradNu = ip.dNdx_u;
1065
1066 auto const& gradNpT = gradNp.transpose().eval();
1067 auto const& gradNTT = gradNT.transpose().eval();
1068
1069 auto const& w = ip.integration_weight;
1070
1071 auto const x_coord =
1072 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1074 this->element_, Nu);
1075
1076 auto const Bu =
1077 LinearBMatrix::computeBMatrix<DisplacementDim,
1078 ShapeFunctionDisplacement::NPOINTS,
1080 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1081
1082 auto const NTN = (Np.transpose() * Np).eval();
1083 auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval();
1084
1085 double const pCap = Np.dot(capillary_pressure);
1086 double const pCap_prev = Np.dot(capillary_pressure_prev);
1087
1088 auto const s_L = current_state.S_L_data.S_L;
1089 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1090
1091 auto const& b = this->process_data_.specific_body_force;
1092
1093 // ---------------------------------------------------------------------
1094 // C-component equation
1095 // ---------------------------------------------------------------------
1096
1097 MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w);
1098 MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w);
1099
1100 if (this->process_data_.apply_mass_lumping)
1101 {
1102 if (pCap - pCap_prev != 0.) // avoid division by Zero
1103 {
1104 MCpC.noalias() +=
1105 NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w);
1106 }
1107 }
1108
1109 MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w);
1110 MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w);
1111
1112 LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w;
1113
1114 LCpC.noalias() += gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w;
1115
1116 LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w;
1117
1118 fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w;
1119
1120 if (!this->process_data_.apply_mass_lumping)
1121 {
1122 fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w);
1123 }
1124 // fC_III
1125 fC.noalias() -= NpT * (ip_out.porosity_data.phi * ip_cv.fC_3a.a * w);
1126
1127 // ---------------------------------------------------------------------
1128 // W-component equation
1129 // ---------------------------------------------------------------------
1130
1131 MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w);
1132 MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w);
1133
1134 if (this->process_data_.apply_mass_lumping)
1135 {
1136 if (pCap - pCap_prev != 0.) // avoid division by Zero
1137 {
1138 MWpC.noalias() +=
1139 NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w);
1140 }
1141 }
1142
1143 MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w);
1144
1145 MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w);
1146
1147 LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w;
1148
1149 LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
1150
1151 LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w;
1152
1153 fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
1154
1155 if (!this->process_data_.apply_mass_lumping)
1156 {
1157 fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w);
1158 }
1159
1160 fW.noalias() -= NpT * (ip_out.porosity_data.phi * ip_cv.fW_3a.a * w);
1161
1162 // ---------------------------------------------------------------------
1163 // - temperature equation
1164 // ---------------------------------------------------------------------
1165
1166 MTu.noalias() +=
1167 BTI2N.transpose() *
1168 (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w);
1169
1170 KTT.noalias() +=
1171 gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
1172
1173 fT.noalias() -= NTT * (ip_cv.fT_1.m * w);
1174
1175 fT.noalias() += gradNTT * ip_cv.fT_2.A * w;
1176
1177 fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w;
1178
1179 fT.noalias() += NTT * (ip_cv.fT_3.N * w);
1180
1181 // ---------------------------------------------------------------------
1182 // - displacement equation
1183 // ---------------------------------------------------------------------
1184
1185 KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w);
1186
1187 KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w);
1188
1189 fU.noalias() -=
1190 (Bu.transpose() * current_state.eff_stress_data.sigma -
1191 N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
1192 w;
1193
1194 if (this->process_data_.apply_mass_lumping)
1195 {
1196 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1197 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1198 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1199 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1200 }
1201 } // int_point-loop
1202}
1203
1204// Assembles the local Jacobian matrix. So far, the linearisation of HT part is
1205// not considered as that in HT process.
1206template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1207 int DisplacementDim>
1208void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1209 DisplacementDim>::
1210 assembleWithJacobian(double const t, double const dt,
1211 std::vector<double> const& local_x,
1212 std::vector<double> const& local_x_prev,
1213 std::vector<double>& /*local_M_data*/,
1214 std::vector<double>& /*local_K_data*/,
1215 std::vector<double>& local_rhs_data,
1216 std::vector<double>& local_Jac_data)
1217{
1218 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1219 temperature_size + displacement_size;
1220 assert(local_x.size() == matrix_size);
1221
1222 auto const temperature = Eigen::Map<VectorType<temperature_size> const>(
1223 local_x.data() + temperature_index, temperature_size);
1224
1225 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size> const>(
1226 local_x.data() + gas_pressure_index, gas_pressure_size);
1227
1228 auto const capillary_pressure =
1229 Eigen::Map<VectorType<capillary_pressure_size> const>(
1230 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1231
1232 auto const displacement = Eigen::Map<VectorType<displacement_size> const>(
1233 local_x.data() + displacement_index, displacement_size);
1234
1235 auto const gas_pressure_prev =
1236 Eigen::Map<VectorType<gas_pressure_size> const>(
1237 local_x_prev.data() + gas_pressure_index, gas_pressure_size);
1238
1239 auto const capillary_pressure_prev =
1240 Eigen::Map<VectorType<capillary_pressure_size> const>(
1241 local_x_prev.data() + capillary_pressure_index,
1242 capillary_pressure_size);
1243
1244 auto const temperature_prev =
1245 Eigen::Map<VectorType<temperature_size> const>(
1246 local_x_prev.data() + temperature_index, temperature_size);
1247
1248 auto const displacement_prev =
1249 Eigen::Map<VectorType<displacement_size> const>(
1250 local_x_prev.data() + displacement_index, displacement_size);
1251
1252 auto local_Jac =
1253 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1254 local_Jac_data, matrix_size, matrix_size);
1255
1256 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1257 local_rhs_data, matrix_size);
1258
1259 // component-formulation
1260 // W - liquid phase main component
1261 // C - gas phase main component
1262
1263 // C component equation matrices
1265 MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1268 C_size, capillary_pressure_size);
1270 MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1272 MatrixType<C_size, displacement_size>::Zero(C_size, displacement_size);
1273
1275 MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1278 C_size, capillary_pressure_size);
1280 MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1281
1282 // mass matrix - W component equation
1284 MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1287 W_size, capillary_pressure_size);
1289 MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1291 MatrixType<W_size, displacement_size>::Zero(W_size, displacement_size);
1292
1293 // stiffness matrix - W component equation
1295 MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1298 W_size, capillary_pressure_size);
1300 MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1301
1302 // mass matrix - temperature equation
1305 temperature_size, displacement_size);
1306
1307 // stiffness matrix - temperature equation
1310 temperature_size);
1311
1312 // stiffness matrices - displacement equation coupling into pressures
1315 displacement_size, gas_pressure_size);
1318 displacement_size, capillary_pressure_size);
1319
1320 // pointer-vectors to the right hand side terms - C-component equation
1321 auto fC = local_f.template segment<C_size>(C_index);
1322 // pointer-vectors to the right hand side terms - W-component equation
1323 auto fW = local_f.template segment<W_size>(W_index);
1324 // pointer-vectors to the right hand side terms - temperature equation
1325 auto fT = local_f.template segment<temperature_size>(temperature_index);
1326 // pointer-vectors to the right hand side terms - displacement equation
1327 auto fU = local_f.template segment<displacement_size>(displacement_index);
1328
1329 unsigned const n_integration_points =
1330 this->integration_method_.getNumberOfPoints();
1331
1333 this->solid_material_, *this->process_data_.phase_transition_model_};
1334
1335 auto const [ip_constitutive_data, ip_constitutive_variables] =
1336 updateConstitutiveVariables(
1337 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1338 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1339 local_x_prev.size()),
1340 t, dt, models);
1341
1342 auto const ip_d_data = updateConstitutiveVariablesDerivatives(
1343 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1344 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1345 local_x_prev.size()),
1346 t, dt, ip_constitutive_data, ip_constitutive_variables, models);
1347
1348 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1349 {
1350 auto& ip = _ip_data[int_point];
1351 auto& ip_cd = ip_constitutive_data[int_point];
1352 auto& ip_dd = ip_d_data[int_point];
1353 auto& ip_cv = ip_constitutive_variables[int_point];
1354 auto& ip_out = this->output_data_[int_point];
1355 auto& current_state = this->current_states_[int_point];
1356 auto& prev_state = this->prev_states_[int_point];
1357
1358 auto const& Np = ip.N_p;
1359 auto const& NT = Np;
1360 auto const& Nu = ip.N_u;
1362 std::nullopt, this->element_.getID(), int_point,
1364 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1366 this->element_, Nu))};
1367
1368 auto const& NpT = Np.transpose().eval();
1369 auto const& NTT = NT.transpose().eval();
1370
1371 auto const& gradNp = ip.dNdx_p;
1372 auto const& gradNT = gradNp;
1373 auto const& gradNu = ip.dNdx_u;
1374
1375 auto const& gradNpT = gradNp.transpose().eval();
1376 auto const& gradNTT = gradNT.transpose().eval();
1377
1378 auto const& w = ip.integration_weight;
1379
1380 auto const x_coord =
1381 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
1383 this->element_, Nu);
1384
1385 auto const Bu =
1386 LinearBMatrix::computeBMatrix<DisplacementDim,
1387 ShapeFunctionDisplacement::NPOINTS,
1389 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1390
1391 auto const NTN = (Np.transpose() * Np).eval();
1392 auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval();
1393
1394 double const div_u_dot =
1395 Invariants::trace(Bu * (displacement - displacement_prev) / dt);
1396
1397 double const pGR = Np.dot(gas_pressure);
1398 double const pCap = Np.dot(capillary_pressure);
1399 double const T = NT.dot(temperature);
1400
1401 GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
1402 GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
1403 GlobalDimVectorType const gradT = gradNT * temperature;
1404
1405 double const pGR_prev = Np.dot(gas_pressure_prev);
1406 double const pCap_prev = Np.dot(capillary_pressure_prev);
1407 double const T_prev = NT.dot(temperature_prev);
1408
1409 auto const& s_L = current_state.S_L_data.S_L;
1410 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1411
1412 auto const& b = this->process_data_.specific_body_force;
1413
1414 // ---------------------------------------------------------------------
1415 // C-component equation
1416 // ---------------------------------------------------------------------
1417
1418 MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w);
1419 MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w);
1420
1421 if (this->process_data_.apply_mass_lumping)
1422 {
1423 if (pCap - pCap_prev != 0.) // avoid division by Zero
1424 {
1425 MCpC.noalias() +=
1426 NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w);
1427 }
1428 }
1429
1430 MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w);
1431 // d (fC_4_MCT * T_dot)/d T
1432 local_Jac
1433 .template block<C_size, temperature_size>(C_index,
1434 temperature_index)
1435 .noalias() += NTN * (ip_dd.dfC_4_MCT.dT * (T - T_prev) / dt * w);
1436
1437 MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w);
1438 // d (fC_4_MCu * u_dot)/d T
1439 local_Jac
1440 .template block<C_size, temperature_size>(C_index,
1441 temperature_index)
1442 .noalias() += NTN * (ip_dd.dfC_4_MCu.dT * div_u_dot * w);
1443
1444 LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w;
1445
1446 // d (fC_4_LCpG * grad p_GR)/d p_GR
1447 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1448 gradNpT * ip_dd.dfC_4_LCpG.dp_GR * gradpGR * Np * w;
1449
1450 // d (fC_4_LCpG * grad p_GR)/d p_cap
1451 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1452 gradNpT * ip_dd.dfC_4_LCpG.dp_cap * gradpGR * Np * w;
1453
1454 // d (fC_4_LCpG * grad p_GR)/d T
1455 local_Jac
1456 .template block<C_size, temperature_size>(C_index,
1457 temperature_index)
1458 .noalias() += gradNpT * ip_dd.dfC_4_LCpG.dT * gradpGR * NT * w;
1459
1460 // d (fC_4_MCpG * p_GR_dot)/d p_GR
1461 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1462 NTN * (ip_dd.dfC_4_MCpG.dp_GR * (pGR - pGR_prev) / dt * w);
1463
1464 // d (fC_4_MCpG * p_GR_dot)/d T
1465 local_Jac
1466 .template block<C_size, temperature_size>(C_index,
1467 temperature_index)
1468 .noalias() +=
1469 NTN * (ip_dd.dfC_4_MCpG.dT * (pGR - pGR_prev) / dt * w);
1470
1471 LCpC.noalias() -= gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w;
1472
1473 /* TODO (naumov) This part is not tested by any of the current ctests.
1474 // d (fC_4_LCpC * grad p_cap)/d p_GR
1475 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1476 gradNpT * ip_dd.dfC_4_LCpC.dp_GR * gradpCap * Np * w;
1477 // d (fC_4_LCpC * grad p_cap)/d p_cap
1478 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1479 gradNpT * ip_dd.dfC_4_LCpC.dp_cap * gradpCap * Np * w;
1480
1481 local_Jac
1482 .template block<C_size, temperature_size>(C_index,
1483 temperature_index)
1484 .noalias() += gradNpT * ip_dd.dfC_4_LCpC.dT * gradpCap * Np * w;
1485 */
1486
1487 LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w;
1488
1489 // fC_1
1490 fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w;
1491
1492 if (!this->process_data_.apply_mass_lumping)
1493 {
1494 // fC_2 = \int a * s_L_dot
1495 fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w);
1496
1497 local_Jac.template block<C_size, C_size>(C_index, C_index)
1498 .noalias() +=
1499 NTN * ((ip_dd.dfC_2a.dp_GR * s_L_dot
1500 /*- ip_cv.fC_2a.a * (ds_L_dp_GR = 0) / dt*/) *
1501 w);
1502
1503 local_Jac.template block<C_size, W_size>(C_index, W_index)
1504 .noalias() +=
1505 NTN * ((ip_dd.dfC_2a.dp_cap * s_L_dot +
1506 ip_cv.fC_2a.a * ip_dd.dS_L_dp_cap() / dt) *
1507 w);
1508
1509 local_Jac
1510 .template block<C_size, temperature_size>(C_index,
1511 temperature_index)
1512 .noalias() += NTN * (ip_dd.dfC_2a.dT * s_L_dot * w);
1513 }
1514 {
1515 // fC_3 = \int phi * a
1516 fC.noalias() -=
1517 NpT * (ip_out.porosity_data.phi * ip_cv.fC_3a.a * w);
1518
1519 local_Jac.template block<C_size, C_size>(C_index, C_index)
1520 .noalias() +=
1521 NTN * (ip_out.porosity_data.phi * ip_dd.dfC_3a.dp_GR * w);
1522
1523 local_Jac.template block<C_size, W_size>(C_index, W_index)
1524 .noalias() +=
1525 NTN * (ip_out.porosity_data.phi * ip_dd.dfC_3a.dp_cap * w);
1526
1527 local_Jac
1528 .template block<C_size, temperature_size>(C_index,
1529 temperature_index)
1530 .noalias() +=
1531 NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fC_3a.a +
1532 ip_out.porosity_data.phi * ip_dd.dfC_3a.dT) *
1533 w);
1534 }
1535 // ---------------------------------------------------------------------
1536 // W-component equation
1537 // ---------------------------------------------------------------------
1538
1539 MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w);
1540 MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w);
1541
1542 if (this->process_data_.apply_mass_lumping)
1543 {
1544 if (pCap - pCap_prev != 0.) // avoid division by Zero
1545 {
1546 MWpC.noalias() +=
1547 NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w);
1548 }
1549 }
1550
1551 MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w);
1552
1553 MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w);
1554
1555 LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w;
1556
1557 // fW_4 LWpG' parts; LWpG = \int grad (a + d) grad
1558 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1559 gradNpT * ip_dd.dfW_4_LWpG.dp_GR * gradpGR * Np * w;
1560
1561 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1562 gradNpT * ip_dd.dfW_4_LWpG.dp_cap * gradpGR * Np * w;
1563
1564 local_Jac
1565 .template block<W_size, temperature_size>(W_index,
1566 temperature_index)
1567 .noalias() += gradNpT * ip_dd.dfW_4_LWpG.dT * gradpGR * NT * w;
1568
1569 LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
1570
1571 // fW_4 LWp_cap' parts; LWpC = \int grad (a + d) grad
1572 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1573 gradNpT * ip_dd.dfW_4_LWpC.dp_GR * gradpCap * Np * w;
1574
1575 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1576 gradNpT * ip_dd.dfW_4_LWpC.dp_cap * gradpCap * Np * w;
1577
1578 local_Jac
1579 .template block<W_size, temperature_size>(W_index,
1580 temperature_index)
1581 .noalias() -= gradNpT * ip_dd.dfW_4_LWpC.dT * gradpCap * NT * w;
1582
1583 LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w;
1584
1585 // fW_1
1586 fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
1587
1588 // fW_2 = \int a * s_L_dot
1589 if (!this->process_data_.apply_mass_lumping)
1590 {
1591 fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w);
1592
1593 local_Jac.template block<W_size, C_size>(W_index, C_index)
1594 .noalias() += NTN * (ip_dd.dfW_2.dp_GR * s_L_dot * w);
1595
1596 // sign negated because of dp_cap = -dp_LR
1597 // TODO (naumov) Had to change the sign to get equal Jacobian WW
1598 // blocks in A2 Test. Where is the error?
1599 local_Jac.template block<W_size, W_size>(W_index, W_index)
1600 .noalias() += NTN * ((ip_dd.dfW_2.dp_cap * s_L_dot +
1601 ip_cv.fW_2.a * ip_dd.dS_L_dp_cap() / dt) *
1602 w);
1603
1604 local_Jac
1605 .template block<W_size, temperature_size>(W_index,
1606 temperature_index)
1607 .noalias() += NTN * (ip_dd.dfW_2.dT * s_L_dot * w);
1608 }
1609
1610 // fW_3 = \int phi * a
1611 fW.noalias() -= NpT * (ip_out.porosity_data.phi * ip_cv.fW_3a.a * w);
1612
1613 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1614 NTN * (ip_out.porosity_data.phi * ip_dd.dfW_3a.dp_GR * w);
1615
1616 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1617 NTN * (ip_out.porosity_data.phi * ip_dd.dfW_3a.dp_cap * w);
1618
1619 local_Jac
1620 .template block<W_size, temperature_size>(W_index,
1621 temperature_index)
1622 .noalias() +=
1623 NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fW_3a.a +
1624 ip_out.porosity_data.phi * ip_dd.dfW_3a.dT) *
1625 w);
1626
1627 // ---------------------------------------------------------------------
1628 // - temperature equation
1629 // ---------------------------------------------------------------------
1630
1631 MTu.noalias() +=
1632 BTI2N.transpose() *
1633 (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w);
1634
1635 // dfT_4/dp_GR
1636 // d (MTu * u_dot)/dp_GR
1637 local_Jac
1638 .template block<temperature_size, C_size>(temperature_index,
1639 C_index)
1640 .noalias() +=
1641 NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR *
1642 div_u_dot * w);
1643
1644 // dfT_4/dp_cap
1645 // d (MTu * u_dot)/dp_cap
1646 local_Jac
1647 .template block<temperature_size, W_size>(temperature_index,
1648 W_index)
1649 .noalias() -=
1650 NTN *
1651 (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap *
1652 div_u_dot * w);
1653
1654 // dfT_4/dT
1655 // d (MTu * u_dot)/dT
1656 local_Jac
1657 .template block<temperature_size, temperature_size>(
1658 temperature_index, temperature_index)
1659 .noalias() +=
1660 NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dT *
1661 div_u_dot * w);
1662
1663 KTT.noalias() +=
1664 gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
1665
1666 // d KTT/dp_GR * T
1667 // TODO (naumov) always zero if lambda_xR have no derivatives wrt. p_GR.
1668 // dlambda_dp_GR =
1669 // (dphi_G_dp_GR = 0) * lambdaGR + phi_G * dlambda_GR_dp_GR +
1670 // (dphi_L_dp_GR = 0) * lambdaLR + phi_L * dlambda_LR_dp_GR +
1671 // (dphi_S_dp_GR = 0) * lambdaSR + phi_S * dlambda_SR_dp_GR +
1672 // = 0
1673 //
1674 // Since dlambda/dp_GR is 0 the derivative is omitted:
1675 // local_Jac
1676 // .template block<temperature_size, C_size>(temperature_index,
1677 // C_index)
1678 // .noalias() += gradNTT * dlambda_dp_GR * gradT * Np * w;
1679
1680 // d KTT/dp_cap * T
1681 local_Jac
1682 .template block<temperature_size, W_size>(temperature_index,
1683 W_index)
1684 .noalias() += gradNTT *
1685 ip_dd.thermal_conductivity_d_data.dlambda_dp_cap *
1686 gradT * Np * w;
1687
1688 // d KTT/dT * T
1689 local_Jac
1690 .template block<temperature_size, temperature_size>(
1691 temperature_index, temperature_index)
1692 .noalias() += gradNTT *
1693 ip_dd.thermal_conductivity_d_data.dlambda_dT * gradT *
1694 NT * w;
1695
1696 // fT_1
1697 fT.noalias() -= NTT * (ip_cv.fT_1.m * w);
1698
1699 // dfT_1/dp_GR
1700 local_Jac
1701 .template block<temperature_size, C_size>(temperature_index,
1702 C_index)
1703 .noalias() += NTN * (ip_dd.dfT_1.dp_GR * w);
1704
1705 // dfT_1/dp_cap
1706 local_Jac
1707 .template block<temperature_size, W_size>(temperature_index,
1708 W_index)
1709 .noalias() += NTN * (ip_dd.dfT_1.dp_cap * w);
1710
1711 // dfT_1/dT
1712 // MTT
1713 local_Jac
1714 .template block<temperature_size, temperature_size>(
1715 temperature_index, temperature_index)
1716 .noalias() += NTN * (ip_dd.dfT_1.dT * w);
1717
1718 // fT_2
1719 fT.noalias() += gradNTT * ip_cv.fT_2.A * w;
1720
1721 // dfT_2/dp_GR
1722 local_Jac
1723 .template block<temperature_size, C_size>(temperature_index,
1724 C_index)
1725 .noalias() -=
1726 // dfT_2/dp_GR first part
1727 gradNTT * ip_dd.dfT_2.dp_GR_Npart * Np * w +
1728 // dfT_2/dp_GR second part
1729 gradNTT * ip_dd.dfT_2.dp_GR_gradNpart * gradNp * w;
1730
1731 // dfT_2/dp_cap
1732 local_Jac
1733 .template block<temperature_size, W_size>(temperature_index,
1734 W_index)
1735 .noalias() -=
1736 // first part of dfT_2/dp_cap
1737 gradNTT * (-ip_dd.dfT_2.dp_cap_Npart) * Np * w +
1738 // second part of dfT_2/dp_cap
1739 gradNTT * (-ip_dd.dfT_2.dp_cap_gradNpart) * gradNp * w;
1740
1741 // dfT_2/dT
1742 local_Jac
1743 .template block<temperature_size, temperature_size>(
1744 temperature_index, temperature_index)
1745 .noalias() -= gradNTT * ip_dd.dfT_2.dT * NT * w;
1746
1747 // fT_3
1748 fT.noalias() += NTT * (ip_cv.fT_3.N * w);
1749
1750 fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w;
1751
1752 // ---------------------------------------------------------------------
1753 // - displacement equation
1754 // ---------------------------------------------------------------------
1755
1756 KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w);
1757
1758 // dfU_2/dp_GR = dKUpG/dp_GR * p_GR + KUpG. The former is zero, the
1759 // latter is handled below.
1760
1761 KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w);
1762
1763 // dfU_2/dp_cap = dKUpC/dp_cap * p_cap + KUpC. The former is handled
1764 // here, the latter below.
1765 local_Jac
1766 .template block<displacement_size, W_size>(displacement_index,
1767 W_index)
1768 .noalias() += BTI2N * (ip_dd.dfu_2_KupC.dp_cap * w);
1769
1770 local_Jac
1771 .template block<displacement_size, displacement_size>(
1772 displacement_index, displacement_index)
1773 .noalias() +=
1774 Bu.transpose() * ip_cd.s_mech_data.stiffness_tensor * Bu * w;
1775
1776 // fU_1
1777 fU.noalias() -=
1778 (Bu.transpose() * current_state.eff_stress_data.sigma -
1779 N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
1780 w;
1781
1782 // KuT
1783 local_Jac
1784 .template block<displacement_size, temperature_size>(
1785 displacement_index, temperature_index)
1786 .noalias() -= Bu.transpose() * ip_dd.dfu_1_KuT.dT * NT * w;
1787
1788 /* TODO (naumov) Test with gravity needed to check this Jacobian part.
1789 local_Jac
1790 .template block<displacement_size, temperature_size>(
1791 displacement_index, temperature_index)
1792 .noalias() += N_u_op(Nu).transpose() * ip_cv.drho_dT * b *
1793 N_u_op(Nu).transpose() * w;
1794 */
1795
1796 if (this->process_data_.apply_mass_lumping)
1797 {
1798 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1799 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1800 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1801 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1802 }
1803 } // int_point-loop
1804
1805 // --- Gas ---
1806 // fC_4
1807 fC.noalias() -= LCpG * gas_pressure + LCpC * capillary_pressure +
1808 LCT * temperature +
1809 MCpG * (gas_pressure - gas_pressure_prev) / dt +
1810 MCpC * (capillary_pressure - capillary_pressure_prev) / dt +
1811 MCT * (temperature - temperature_prev) / dt +
1812 MCu * (displacement - displacement_prev) / dt;
1813
1814 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1815 LCpG + MCpG / dt;
1816 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1817 LCpC + MCpC / dt;
1818 local_Jac
1819 .template block<C_size, temperature_size>(C_index, temperature_index)
1820 .noalias() += LCT + MCT / dt;
1821 local_Jac
1822 .template block<C_size, displacement_size>(C_index, displacement_index)
1823 .noalias() += MCu / dt;
1824
1825 // --- Capillary pressure ---
1826 // fW_4
1827 fW.noalias() -= LWpG * gas_pressure + LWpC * capillary_pressure +
1828 LWT * temperature +
1829 MWpG * (gas_pressure - gas_pressure_prev) / dt +
1830 MWpC * (capillary_pressure - capillary_pressure_prev) / dt +
1831 MWT * (temperature - temperature_prev) / dt +
1832 MWu * (displacement - displacement_prev) / dt;
1833
1834 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1835 LWpC + MWpC / dt;
1836 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1837 LWpG + MWpG / dt;
1838 local_Jac
1839 .template block<W_size, temperature_size>(W_index, temperature_index)
1840 .noalias() += LWT + MWT / dt;
1841 local_Jac
1842 .template block<W_size, displacement_size>(W_index, displacement_index)
1843 .noalias() += MWu / dt;
1844
1845 // --- Temperature ---
1846 // fT_4
1847 fT.noalias() -=
1848 KTT * temperature + MTu * (displacement - displacement_prev) / dt;
1849
1850 local_Jac
1851 .template block<temperature_size, temperature_size>(temperature_index,
1852 temperature_index)
1853 .noalias() += KTT;
1854 local_Jac
1855 .template block<temperature_size, displacement_size>(temperature_index,
1856 displacement_index)
1857 .noalias() += MTu / dt;
1858
1859 // --- Displacement ---
1860 // fU_2
1861 fU.noalias() -= KUpG * gas_pressure + KUpC * capillary_pressure;
1862
1863 local_Jac
1864 .template block<displacement_size, C_size>(displacement_index, C_index)
1865 .noalias() += KUpG;
1866 local_Jac
1867 .template block<displacement_size, W_size>(displacement_index, W_index)
1868 .noalias() += KUpC;
1869}
1870
1871template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1872 int DisplacementDim>
1873void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1874 DisplacementDim>::
1875 computeSecondaryVariableConcrete(double const t, double const dt,
1876 Eigen::VectorXd const& local_x,
1877 Eigen::VectorXd const& local_x_prev)
1878{
1879 auto const gas_pressure =
1880 local_x.template segment<gas_pressure_size>(gas_pressure_index);
1881 auto const capillary_pressure =
1882 local_x.template segment<capillary_pressure_size>(
1883 capillary_pressure_index);
1884 auto const liquid_pressure = (gas_pressure - capillary_pressure).eval();
1885
1887 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1888 DisplacementDim>(this->element_, this->is_axially_symmetric_,
1889 gas_pressure,
1890 *this->process_data_.gas_pressure_interpolated);
1891
1893 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1894 DisplacementDim>(this->element_, this->is_axially_symmetric_,
1895 capillary_pressure,
1896 *this->process_data_.capillary_pressure_interpolated);
1897
1899 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1900 DisplacementDim>(this->element_, this->is_axially_symmetric_,
1901 liquid_pressure,
1902 *this->process_data_.liquid_pressure_interpolated);
1903
1904 auto const temperature =
1905 local_x.template segment<temperature_size>(temperature_index);
1906
1908 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1909 DisplacementDim>(this->element_, this->is_axially_symmetric_,
1910 temperature,
1911 *this->process_data_.temperature_interpolated);
1912
1913 unsigned const n_integration_points =
1914 this->integration_method_.getNumberOfPoints();
1915
1916 double saturation_avg = 0;
1917
1919 this->solid_material_, *this->process_data_.phase_transition_model_};
1920
1921 updateConstitutiveVariables(local_x, local_x_prev, t, dt, models);
1922
1923 for (unsigned ip = 0; ip < n_integration_points; ip++)
1924 {
1925 saturation_avg += this->current_states_[ip].S_L_data.S_L;
1926 }
1927 saturation_avg /= n_integration_points;
1928 (*this->process_data_.element_saturation)[this->element_.getID()] =
1929 saturation_avg;
1930}
1931
1932} // namespace TH2M
1933} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
double getWeight() const
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition TH2MFEM.h:51
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition TH2MFEM.h:223
typename ShapeMatricesTypePressure::template MatrixType< M, N > MatrixType
Definition TH2MFEM.h:62
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition TH2MFEM.h:54
typename ShapeMatricesTypePressure::GlobalDimVectorType GlobalDimVectorType
Definition TH2MFEM.h:68
std::vector< IpData > _ip_data
Definition TH2MFEM.h:219
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
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.
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u
void eval(SpaceTimeData const &x_t, MediaData const &media_data, BiotData &out) const
Definition Biot.cpp:16
void eval(SpaceTimeData const &x_t, MediaData const &media_data, SaturationData const &S_L_data, BishopsData &out) const
Definition Bishops.cpp:16
Data that is needed for the equation system assembly.
DiffusionVelocityModel< DisplacementDim > diffusion_velocity_model
ElasticTangentStiffnessModel< DisplacementDim > elastic_tangent_stiffness_model
ThermalConductivityModel< DisplacementDim > thermal_conductivity_model
SolidCompressibilityModel< DisplacementDim, SolidConstitutiveRelation< DisplacementDim > > beta_p_SR_model
SolidThermalExpansionModel< DisplacementDim > s_therm_exp_model
MechanicalStrainModel< DisplacementDim > mechanical_strain_model
void dEval(FluidDensityData const &fluid_density_data, FluidEnthalpyData const &fluid_enthalpy_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, SaturationData const &S_L_data, SolidDensityData const &solid_density_data, SolidDensityDerivativeData const &solid_density_d_data, SolidEnthalpyData const &solid_enthalpy_data, SolidHeatCapacityData const &solid_heat_capacity_data, EffectiveVolumetricEnthalpyDerivatives &effective_volumetric_enthalpy_d_data) const
Definition Enthalpy.cpp:35
void eval(FluidDensityData const &fluid_density_data, FluidEnthalpyData const &fluid_enthalpy_data, PorosityData const &porosity_data, SaturationData const &S_L_data, SolidDensityData const &solid_density_data, SolidEnthalpyData const &solid_enthalpy_data, EffectiveVolumetricEnthalpy &effective_volumetric_enthalpy_data) const
Definition Enthalpy.cpp:16
void dEval(BiotData const &biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, SaturationData const &S_L_data, SaturationDataDeriv const &dS_L_dp_cap, SolidCompressibilityData const &beta_p_SR, FC2aDerivativeData &dfC_2a) const
Definition CEquation.cpp:47
void eval(BiotData const biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, SaturationData const &S_L_data, SolidCompressibilityData const beta_p_SR, FC2aData &fC_2a) const
Definition CEquation.cpp:29
void eval(double const dt, ConstituentDensityData const &constituent_density_data, PrevState< ConstituentDensityData > const &constituent_density_data_prev, SaturationData const &S_L_data, FC3aData &fC_3a) const
void dEval(double const dt, ConstituentDensityData const &constituent_density_data, PrevState< ConstituentDensityData > const &constituent_density_data_prev, PhaseTransitionData const &phase_transition_data, SaturationData const &S_L_data, SaturationDataDeriv const &dS_L_dp_cap, FC3aDerivativeData &dfC_3a) const
void eval(BiotData const &biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, PrevState< SaturationData > const &S_L_data_prev, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FC4MCpCData &fC_4_MCpC) const
void eval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FC4MCpGData &fC_4_MCpG) const
void dEval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FC4MCpGDerivativeData &dfC_4_MCpG) const
void eval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, SaturationData const &S_L_data, FC4MCuData &fC_4_MCu) const
void dEval(BiotData const &biot_data, PhaseTransitionData const &phase_transition_data, SaturationData const &S_L_data, FC4MCuDerivativeData &dfC_4_MCu) const
void dEval(double const dt, EffectiveVolumetricInternalEnergyDerivatives const &effective_volumetric_internal_energy_d_data, FT1DerivativeData &dfT_1) const
Definition TEquation.cpp:33
void eval(double const dt, InternalEnergyData const &internal_energy_data, PrevState< InternalEnergyData > const &internal_energy_data_prev, FT1Data &fT_1) const
Definition TEquation.cpp:16
void dEval(BiotData const &biot_data, BishopsData const &chi_S_L, CapillaryPressureData const &p_cap, SaturationDataDeriv const &dS_L_dp_cap, FU2KUpCDerivativeData &dfu_2_KupC) const
Definition UEquation.cpp:36
void dEval(BiotData const &biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SaturationDataDeriv const &dS_L_dp_cap, SolidCompressibilityData const &beta_p_SR, FW2DerivativeData &dfW_2) const
Definition WEquation.cpp:48
void eval(BiotData const biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SolidCompressibilityData const beta_p_SR, FW2Data &fW_2) const
Definition WEquation.cpp:29
void dEval(double const dt, ConstituentDensityData const &constituent_density_data, PhaseTransitionData const &phase_transition_data, PrevState< ConstituentDensityData > const &constituent_density_data_prev, PrevState< PureLiquidDensityData > const &rho_W_LR_prev, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SaturationDataDeriv const &dS_L_dp_cap, FW3aDerivativeData &dfW_3a) const
void eval(double const dt, ConstituentDensityData const &constituent_density_data, PrevState< ConstituentDensityData > const &constituent_density_data_prev, PrevState< PureLiquidDensityData > const &rho_W_LR_prev, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, FW3aData &fW_3a) const
void eval(BiotData const &biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, PrevState< SaturationData > const &S_L_data_prev, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FW4MWpCData &fW_4_MWpC) const
void eval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FW4MWpGData &fW_4_MWpG) const
void eval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, FW4MWuData &fW_4_MWu) const
void eval(FluidDensityData const &fluid_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, SaturationData const &S_L_data, SolidDensityData const &solid_density_data, SolidEnthalpyData const &solid_enthalpy_data, InternalEnergyData &internal_energy_data) const
void dEval(FluidDensityData const &fluid_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, SaturationData const &S_L_data, SolidDensityData const &solid_density_data, SolidDensityDerivativeData const &solid_density_d_data, SolidEnthalpyData const &solid_enthalpy_data, SolidHeatCapacityData const &solid_heat_capacity_data, EffectiveVolumetricInternalEnergyDerivatives &effective_volumetric_internal_energy_d_data) const
virtual void eval(SpaceTimeData const &x_t, MediaData const &media_data, GasPressureData const &p_GR, CapillaryPressureData const &p_cap, TemperatureData const &T_data, PureLiquidDensityData const &rho_W_LR, FluidEnthalpyData &fluid_enthalpy_data, MassMoleFractionsData &mass_mole_fractions_data, FluidDensityData &fluid_density_data, VapourPartialPressureData &vapour_pressure_data, ConstituentDensityData &constituent_density_data, PhaseTransitionData &cv) const =0
void dEval(SpaceTimeData const &x_t, MediaData const &media_data, PorosityData const &porosity_data, SaturationDataDeriv const &dS_L_dp_cap, PorosityDerivativeData &porosity_d_data) const
Definition Porosity.cpp:30
void eval(SpaceTimeData const &x_t, MediaData const &media_data, PorosityData &porosity_data) const
Definition Porosity.cpp:17
void eval(SpaceTimeData const &x_t, MediaData const &media_data, GasPressureData const &p_GR, CapillaryPressureData const &p_cap, TemperatureData const &T_data, PureLiquidDensityData &out) const
void dEval(SpaceTimeData const &x_t, MediaData const &media_data, CapillaryPressureData const &p_cap, SaturationDataDeriv &dS_L_data) const
void eval(SpaceTimeData const &x_t, MediaData const &media_data, CapillaryPressureData const &p_cap, SaturationData &S_L_data) const
void dEval(SpaceTimeData const &x_t, MediaData const &media_data, TemperatureData const &T_data, SolidDensityDerivativeData &solid_density_d_data) const
void eval(SpaceTimeData const &x_t, MediaData const &media_data, TemperatureData const &T_data, SolidDensityData &solid_density_data) const
void eval(SolidHeatCapacityData const &solid_heat_capacity_data, TemperatureData const &T_data, SolidEnthalpyData &solid_enthalpy_data) const
Definition Enthalpy.cpp:92
void eval(SpaceTimeData const &x_t, MediaData const &media_data, TemperatureData const &T_data, SolidHeatCapacityData &solid_heat_capacity) const
void eval(SpaceTimeData const &x_t, MediaData const &media_data, TemperatureData const &T_data, MassMoleFractionsData const &mass_mole_fractions_data, ViscosityData &viscosity_data) const
Definition Viscosity.cpp:16
NumLib::GenericIntegrationMethod const & integration_method_