OGS
ComponentTransportFEM.h
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
2// SPDX-License-Identifier: BSD-3-Clause
3
4#pragma once
5
6#include <numeric>
7#include <vector>
8
9#include "ChemistryLib/ChemicalSolverInterface.h"
10#include "ComponentTransportProcessData.h"
11#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
12#include "MaterialLib/MPL/Medium.h"
13#include "MaterialLib/MPL/Property.h"
14#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
15#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
16#include "NumLib/DOF/DOFTableUtil.h"
17#include "NumLib/Extrapolation/ExtrapolatableElement.h"
18#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
19#include "NumLib/Fem/InitShapeMatrices.h"
20#include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
21#include "NumLib/Fem/Interpolation.h"
22#include "NumLib/Fem/ShapeMatrixCache.h"
23#include "NumLib/Fem/ShapeMatrixPolicy.h"
24#include "NumLib/NumericalStability/AdvectionMatrixAssembler.h"
25#include "NumLib/NumericalStability/HydrodynamicDispersion.h"
26#include "ParameterLib/Parameter.h"
27#include "ProcessLib/LocalAssemblerInterface.h"
28#include "ProcessLib/ProcessVariable.h"
29
30namespace ProcessLib
31{
32namespace ComponentTransport
33{
34template <typename GlobalDimNodalMatrixType>
35struct IntegrationPointData final
36{
37 IntegrationPointData(GlobalDimNodalMatrixType const& dNdx_,
38 double const& integration_weight_)
39 : dNdx(dNdx_), integration_weight(integration_weight_)
40 {
41 }
42
43 void pushBackState() { porosity_prev = porosity; }
44 GlobalDimNodalMatrixType const dNdx;
45 double const integration_weight;
46
47 // -1 indicates that no chemical reaction takes place in the element to
48 // which the integration point belongs.
49 GlobalIndexType chemical_system_id = -1;
50
51 double porosity = std::numeric_limits<double>::quiet_NaN();
52 double porosity_prev = std::numeric_limits<double>::quiet_NaN();
53 EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
54};
55
56class ComponentTransportLocalAssemblerInterface
57 : public ProcessLib::LocalAssemblerInterface,
58 public NumLib::ExtrapolatableElement
59{
60public:
61 ComponentTransportLocalAssemblerInterface() = default;
62
63 virtual void setChemicalSystemID(std::size_t const /*mesh_item_id*/) = 0;
64
65 void initializeChemicalSystem(
66 std::size_t const mesh_item_id,
67 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
68 std::vector<GlobalVector*> const& x, double const t)
69 {
70 std::vector<double> local_x_vec;
71
72 auto const n_processes = x.size();
73 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
74 {
75 auto const indices =
76 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
77 assert(!indices.empty());
78 auto const local_solution = x[process_id]->get(indices);
79 local_x_vec.insert(std::end(local_x_vec),
80 std::begin(local_solution),
81 std::end(local_solution));
82 }
83 auto const local_x = MathLib::toVector(local_x_vec);
84
85 initializeChemicalSystemConcrete(local_x, t);
86 }
87
88 void setChemicalSystem(
89 std::size_t const mesh_item_id,
90 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
91 std::vector<GlobalVector*> const& x, double const t, double const dt)
92 {
93 std::vector<double> local_x_vec;
94
95 auto const n_processes = x.size();
96 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
97 {
98 auto const indices =
99 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
100 assert(!indices.empty());
101 auto const local_solution = x[process_id]->get(indices);
102 local_x_vec.insert(std::end(local_x_vec),
103 std::begin(local_solution),
104 std::end(local_solution));
105 }
106 auto const local_x = MathLib::toVector(local_x_vec);
107
108 setChemicalSystemConcrete(local_x, t, dt);
109 }
110
111 void assembleReactionEquation(
112 std::size_t const mesh_item_id,
113 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
114 std::vector<GlobalVector*> const& x, double const t, double const dt,
115 GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, int const process_id)
116 {
117 std::vector<double> local_x_vec;
118
119 auto const n_processes = x.size();
120 for (std::size_t pcs_id = 0; pcs_id < n_processes; ++pcs_id)
121 {
122 auto const indices =
123 NumLib::getIndices(mesh_item_id, *dof_tables[pcs_id]);
124 assert(!indices.empty());
125 auto const local_solution = x[pcs_id]->get(indices);
126 local_x_vec.insert(std::end(local_x_vec),
127 std::begin(local_solution),
128 std::end(local_solution));
129 }
130 auto const local_x = MathLib::toVector(local_x_vec);
131
132 auto const indices =
133 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
134 auto const num_r_c = indices.size();
135
136 std::vector<double> local_M_data;
137 local_M_data.reserve(num_r_c * num_r_c);
138 std::vector<double> local_K_data;
139 local_K_data.reserve(num_r_c * num_r_c);
140 std::vector<double> local_b_data;
141 local_b_data.reserve(num_r_c);
142
143 assembleReactionEquationConcrete(t, dt, local_x, local_M_data,
144 local_K_data, local_b_data,
145 process_id);
146
147 auto const r_c_indices =
148 NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices, indices);
149 if (!local_M_data.empty())
150 {
151 auto const local_M =
152 MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
153 M.add(r_c_indices, local_M);
154 }
155 if (!local_K_data.empty())
156 {
157 auto const local_K =
158 MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
159 K.add(r_c_indices, local_K);
160 }
161 if (!local_b_data.empty())
162 {
163 b.add(indices, local_b_data);
164 }
165 }
166
167 virtual void postSpeciationCalculation(std::size_t const ele_id,
168 double const t, double const dt) = 0;
169
170 virtual void computeReactionRelatedSecondaryVariable(
171 std::size_t const ele_id) = 0;
172
173 virtual std::vector<double> const& getIntPtLiquidDensity(
174 const double t,
175 std::vector<GlobalVector*> const& x,
176 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
177 std::vector<double>& cache) const = 0;
178
179 virtual std::vector<double> const& getIntPtDarcyVelocity(
180 const double t,
181 std::vector<GlobalVector*> const& x,
182 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
183 std::vector<double>& cache) const = 0;
184
185 virtual std::vector<double> const& getIntPtMolarFlux(
186 const double t, std::vector<GlobalVector*> const& x,
187 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
188 std::vector<double>& cache, int const component_id) const = 0;
189
190private:
191 virtual void initializeChemicalSystemConcrete(
192 Eigen::VectorXd const& /*local_x*/, double const /*t*/) = 0;
193
194 virtual void setChemicalSystemConcrete(Eigen::VectorXd const& /*local_x*/,
195 double const /*t*/,
196 double const /*dt*/) = 0;
197
198 virtual void assembleReactionEquationConcrete(
199 double const t, double const dt, Eigen::VectorXd const& local_x,
200 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
201 std::vector<double>& local_b_data, int const transport_process_id) = 0;
202};
203
204template <typename ShapeFunction, int GlobalDim>
205class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
206{
207 // When monolithic scheme is adopted, nodal pressure and nodal concentration
208 // are accessed by vector index.
209 static const int pressure_index = 0;
210 const int temperature_index = -1;
211 const int first_concentration_index = -1;
212
213 static const int pressure_size = ShapeFunction::NPOINTS;
214 static const int temperature_size = ShapeFunction::NPOINTS;
215 static const int concentration_size =
216 ShapeFunction::NPOINTS; // per component
217
218 using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
219 using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
220
221 using LocalBlockMatrixType =
222 typename ShapeMatricesType::template MatrixType<pressure_size,
223 pressure_size>;
224 using LocalSegmentVectorType =
225 typename ShapeMatricesType::template VectorType<pressure_size>;
226
227 using LocalMatrixType =
228 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
229 using LocalVectorType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
230
231 using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
232 using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
233
234 using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
235 using GlobalDimNodalMatrixType =
236 typename ShapeMatricesType::GlobalDimNodalMatrixType;
237 using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
238
239public:
240 LocalAssemblerData(
241 MeshLib::Element const& element,
242 std::size_t const local_matrix_size,
243 NumLib::GenericIntegrationMethod const& integration_method,
244 bool is_axially_symmetric,
245 ComponentTransportProcessData const& process_data,
246 std::vector<std::reference_wrapper<ProcessVariable>> const&
247 transport_process_variables)
248 : temperature_index(process_data.isothermal ? -1
249 : ShapeFunction::NPOINTS),
250 first_concentration_index(process_data.isothermal
251 ? ShapeFunction::NPOINTS
252 : 2 * ShapeFunction::NPOINTS),
253 _element(element),
254 _process_data(process_data),
255 _integration_method(integration_method),
256 _transport_process_variables(transport_process_variables)
257 {
258 (void)local_matrix_size;
259
260 unsigned const n_integration_points =
261 _integration_method.getNumberOfPoints();
262 _ip_data.reserve(n_integration_points);
263
264 ParameterLib::SpatialPosition pos;
265 pos.setElementID(_element.getID());
266
267 double const aperture_size = _process_data.aperture_size(0.0, pos)[0];
268
269 auto const shape_matrices =
270 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
271 GlobalDim>(element, is_axially_symmetric,
272 _integration_method);
273 auto const& medium =
274 *_process_data.media_map.getMedium(_element.getID());
275 for (unsigned ip = 0; ip < n_integration_points; ip++)
276 {
277 _ip_data.emplace_back(
278 shape_matrices[ip].dNdx,
279 _integration_method.getWeightedPoint(ip).getWeight() *
280 shape_matrices[ip].integralMeasure *
281 shape_matrices[ip].detJ * aperture_size);
282
283 _ip_data[ip].porosity =
284 medium[MaterialPropertyLib::PropertyType::porosity]
285 .template initialValue<double>(
286 pos, std::numeric_limits<double>::quiet_NaN() /*t*/);
287
288 _ip_data[ip].pushBackState();
289 }
290 }
291
292 void setChemicalSystemID(std::size_t const /*mesh_item_id*/) override
293 {
294 assert(_process_data.chemical_solver_interface);
295 // chemical system index map
296 auto& chemical_system_index_map =
297 _process_data.chemical_solver_interface->chemical_system_index_map;
298
299 unsigned const n_integration_points =
300 _integration_method.getNumberOfPoints();
301 for (unsigned ip = 0; ip < n_integration_points; ip++)
302 {
303 _ip_data[ip].chemical_system_id =
304 chemical_system_index_map.empty()
305 ? 0
306 : chemical_system_index_map.back() + 1;
307 chemical_system_index_map.push_back(
308 _ip_data[ip].chemical_system_id);
309 }
310 }
311
312 void initializeChemicalSystemConcrete(Eigen::VectorXd const& local_x,
313 double const t) override
314 {
315 assert(_process_data.chemical_solver_interface);
316
317 auto const& medium =
318 *_process_data.media_map.getMedium(_element.getID());
319
320 ParameterLib::SpatialPosition pos;
321 pos.setElementID(_element.getID());
322
323 auto const& Ns =
324 _process_data.shape_matrix_cache
325 .NsHigherOrder<typename ShapeFunction::MeshElement>();
326
327 unsigned const n_integration_points =
328 _integration_method.getNumberOfPoints();
329
330 for (unsigned ip = 0; ip < n_integration_points; ip++)
331 {
332 auto& ip_data = _ip_data[ip];
333 auto const& N = Ns[ip];
334 auto const& chemical_system_id = ip_data.chemical_system_id;
335
336 // set position with N as the shape matrix at the current
337 // integration point
338 pos.setCoordinates(MathLib::Point3d(
339 NumLib::interpolateCoordinates<ShapeFunction,
340 ShapeMatricesType>(_element,
341 N)));
342
343 auto const n_component = _transport_process_variables.size();
344 std::vector<double> C_int_pt(n_component);
345 for (unsigned component_id = 0; component_id < n_component;
346 ++component_id)
347 {
348 auto const concentration_index =
349 first_concentration_index +
350 component_id * concentration_size;
351 auto const local_C =
352 local_x.template segment<concentration_size>(
353 concentration_index);
354
355 NumLib::shapeFunctionInterpolate(local_C, N,
356 C_int_pt[component_id]);
357 }
358
359 _process_data.chemical_solver_interface
360 ->initializeChemicalSystemConcrete(C_int_pt, chemical_system_id,
361 medium, pos, t);
362 }
363 }
364
365 void setChemicalSystemConcrete(Eigen::VectorXd const& local_x,
366 double const t, double dt) override
367 {
368 assert(_process_data.chemical_solver_interface);
369
370 auto const& medium =
371 _process_data.media_map.getMedium(_element.getID());
372
373 MaterialPropertyLib::VariableArray vars;
374 MaterialPropertyLib::VariableArray vars_prev;
375
376 ParameterLib::SpatialPosition pos;
377 pos.setElementID(_element.getID());
378
379 auto const& Ns =
380 _process_data.shape_matrix_cache
381 .NsHigherOrder<typename ShapeFunction::MeshElement>();
382
383 unsigned const n_integration_points =
384 _integration_method.getNumberOfPoints();
385
386 for (unsigned ip = 0; ip < n_integration_points; ip++)
387 {
388 auto& ip_data = _ip_data[ip];
389 auto const& N = Ns[ip];
390 auto& porosity = ip_data.porosity;
391 auto const& porosity_prev = ip_data.porosity_prev;
392 auto const& chemical_system_id = ip_data.chemical_system_id;
393
394 auto const n_component = _transport_process_variables.size();
395 std::vector<double> C_int_pt(n_component);
396 for (unsigned component_id = 0; component_id < n_component;
397 ++component_id)
398 {
399 auto const concentration_index =
400 first_concentration_index +
401 component_id * concentration_size;
402 auto const local_C =
403 local_x.template segment<concentration_size>(
404 concentration_index);
405
406 NumLib::shapeFunctionInterpolate(local_C, N,
407 C_int_pt[component_id]);
408 }
409
410 {
411 vars_prev.porosity = porosity_prev;
412
413 porosity =
414 _process_data.chemically_induced_porosity_change
415 ? porosity_prev
416 : medium
417 ->property(
418 MaterialPropertyLib::PropertyType::porosity)
419 .template value<double>(vars, vars_prev, pos, t,
420 dt);
421
422 vars.porosity = porosity;
423 }
424
425 _process_data.chemical_solver_interface->setChemicalSystemConcrete(
426 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
427 }
428 }
429
430 void postSpeciationCalculation(std::size_t const ele_id, double const t,
431 double const dt) override
432 {
433 if (!_process_data.chemically_induced_porosity_change)
434 {
435 return;
436 }
437
438 auto const& medium = *_process_data.media_map.getMedium(ele_id);
439
440 ParameterLib::SpatialPosition pos;
441 pos.setElementID(ele_id);
442
443 for (auto& ip_data : _ip_data)
444 {
445 ip_data.porosity = ip_data.porosity_prev;
446
447 _process_data.chemical_solver_interface
448 ->updateVolumeFractionPostReaction(ip_data.chemical_system_id,
449 medium, pos,
450 ip_data.porosity, t, dt);
451
452 _process_data.chemical_solver_interface->updatePorosityPostReaction(
453 ip_data.chemical_system_id, medium, ip_data.porosity);
454 }
455 }
456
457 void assemble(double const t, double const dt,
458 std::vector<double> const& local_x,
459 std::vector<double> const& /*local_x_prev*/,
460 std::vector<double>& local_M_data,
461 std::vector<double>& local_K_data,
462 std::vector<double>& local_b_data) override
463 {
464 auto const local_matrix_size = local_x.size();
465 // Nodal DOFs include pressure
466 int const num_nodal_dof = 1 + _transport_process_variables.size();
467 // This assertion is valid only if all nodal d.o.f. use the same shape
468 // matrices.
469 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
470
471 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
472 local_M_data, local_matrix_size, local_matrix_size);
473 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
474 local_K_data, local_matrix_size, local_matrix_size);
475 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
476 local_b_data, local_matrix_size);
477
478 // Get block matrices
479 auto Kpp = local_K.template block<pressure_size, pressure_size>(
480 pressure_index, pressure_index);
481 auto Mpp = local_M.template block<pressure_size, pressure_size>(
482 pressure_index, pressure_index);
483 auto Bp = local_b.template segment<pressure_size>(pressure_index);
484
485 auto local_p = Eigen::Map<const NodalVectorType>(
486 &local_x[pressure_index], pressure_size);
487
488 auto const& b =
489 _process_data
490 .projected_specific_body_force_vectors[_element.getID()];
491
492 auto const number_of_components = num_nodal_dof - 1;
493 for (int component_id = 0; component_id < number_of_components;
494 ++component_id)
495 {
496 /* Partitioned assembler matrix
497 * | pp | pc1 | pc2 | pc3 |
498 * |-----|-----|-----|-----|
499 * | c1p | c1c1| 0 | 0 |
500 * |-----|-----|-----|-----|
501 * | c2p | 0 | c2c2| 0 |
502 * |-----|-----|-----|-----|
503 * | c3p | 0 | 0 | c3c3|
504 */
505 auto concentration_index =
506 pressure_size + component_id * concentration_size;
507
508 auto KCC =
509 local_K.template block<concentration_size, concentration_size>(
510 concentration_index, concentration_index);
511 auto MCC =
512 local_M.template block<concentration_size, concentration_size>(
513 concentration_index, concentration_index);
514 auto MCp =
515 local_M.template block<concentration_size, pressure_size>(
516 concentration_index, pressure_index);
517 auto MpC =
518 local_M.template block<pressure_size, concentration_size>(
519 pressure_index, concentration_index);
520
521 auto local_C = Eigen::Map<const NodalVectorType>(
522 &local_x[concentration_index], concentration_size);
523
524 assembleBlockMatrices(b, component_id, t, dt, local_C, local_p, KCC,
525 MCC, MCp, MpC, Kpp, Mpp, Bp);
526
527 if (_process_data.chemical_solver_interface)
528 {
529 auto const stoichiometric_matrix =
530 _process_data.chemical_solver_interface
531 ->getStoichiometricMatrix();
532
533 assert(stoichiometric_matrix);
534
535 for (Eigen::SparseMatrix<double>::InnerIterator it(
536 *stoichiometric_matrix, component_id);
537 it;
538 ++it)
539 {
540 auto const stoichiometric_coefficient = it.value();
541 auto const coupled_component_id = it.row();
542 auto const kinetic_prefactor =
543 _process_data.chemical_solver_interface
544 ->getKineticPrefactor(coupled_component_id);
545
546 auto const concentration_index =
547 pressure_size + component_id * concentration_size;
548 auto const coupled_concentration_index =
549 pressure_size +
550 coupled_component_id * concentration_size;
551 auto KCmCn = local_K.template block<concentration_size,
552 concentration_size>(
553 concentration_index, coupled_concentration_index);
554
555 // account for the coupling between components
556 assembleKCmCn(component_id, t, dt, KCmCn,
557 stoichiometric_coefficient,
558 kinetic_prefactor);
559 }
560 }
561 }
562 }
563
564 void assembleBlockMatrices(
565 GlobalDimVectorType const& b, int const component_id, double const t,
566 double const dt,
567 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
568 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
569 Eigen::Ref<LocalBlockMatrixType> KCC,
570 Eigen::Ref<LocalBlockMatrixType> MCC,
571 Eigen::Ref<LocalBlockMatrixType> MCp,
572 Eigen::Ref<LocalBlockMatrixType> MpC,
573 Eigen::Ref<LocalBlockMatrixType> Kpp,
574 Eigen::Ref<LocalBlockMatrixType> Mpp,
575 Eigen::Ref<LocalSegmentVectorType> Bp)
576 {
577 unsigned const n_integration_points =
578 _integration_method.getNumberOfPoints();
579
580 ParameterLib::SpatialPosition pos;
581 pos.setElementID(_element.getID());
582
583 MaterialPropertyLib::VariableArray vars;
584
585 // Get material properties
586 auto const& medium =
587 *_process_data.media_map.getMedium(_element.getID());
588 // Select the only valid for component transport liquid phase.
589 auto const& phase = medium.phase("AqueousLiquid");
590
591 // Assume that the component name is the same as the process variable
592 // name. Components are shifted by one because the first one is always
593 // pressure.
594 auto const& component = phase.component(
595 _transport_process_variables[component_id].get().getName());
596
597 LocalBlockMatrixType KCC_Laplacian =
598 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
599
600 std::vector<GlobalDimVectorType> ip_flux_vector;
601 double average_velocity_norm = 0.0;
602 if (!_process_data.non_advective_form)
603 {
604 ip_flux_vector.reserve(n_integration_points);
605 }
606
607 auto const& Ns =
608 _process_data.shape_matrix_cache
609 .NsHigherOrder<typename ShapeFunction::MeshElement>();
610
611 for (unsigned ip(0); ip < n_integration_points; ++ip)
612 {
613 auto& ip_data = _ip_data[ip];
614 auto const& dNdx = ip_data.dNdx;
615 auto const& N = Ns[ip];
616 auto const& w = ip_data.integration_weight;
617 auto& porosity = ip_data.porosity;
618
619 double C_int_pt = 0.0;
620 double p_int_pt = 0.0;
621
622 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
623 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
624
625 // set position with N as the shape matrix at the current
626 // integration point
627 pos.setCoordinates(MathLib::Point3d(
628 NumLib::interpolateCoordinates<ShapeFunction,
629 ShapeMatricesType>(_element,
630 N)));
631
632 vars.concentration = C_int_pt;
633 vars.liquid_phase_pressure = p_int_pt;
634
635 // update according to a particular porosity model
636 porosity = medium[MaterialPropertyLib::PropertyType::porosity]
637 .template value<double>(vars, pos, t, dt);
638 vars.porosity = porosity;
639
640 auto const& retardation_factor =
641 component[MaterialPropertyLib::PropertyType::retardation_factor]
642 .template value<double>(vars, pos, t, dt);
643
644 auto const& solute_dispersivity_transverse = medium.template value<
645 double>(
646 MaterialPropertyLib::PropertyType::transversal_dispersivity);
647
648 auto const& solute_dispersivity_longitudinal =
649 medium.template value<double>(
650 MaterialPropertyLib::PropertyType::
651 longitudinal_dispersivity);
652
653 // Use the fluid density model to compute the density
654 // TODO (renchao): concentration of which component as the argument
655 // for calculation of fluid density
656 auto const density =
657 phase[MaterialPropertyLib::PropertyType::density]
658 .template value<double>(vars, pos, t, dt);
659
660 auto const decay_rate =
661 component[MaterialPropertyLib::PropertyType::decay_rate]
662 .template value<double>(vars, pos, t, dt);
663
664 auto const& pore_diffusion_coefficient =
665 MaterialPropertyLib::formEigenTensor<GlobalDim>(
666 component[MaterialPropertyLib::PropertyType::pore_diffusion]
667 .value(vars, pos, t, dt));
668
669 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
670 medium[MaterialPropertyLib::PropertyType::permeability].value(
671 vars, pos, t, dt));
672
673 // Use the viscosity model to compute the viscosity
674 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
675 .template value<double>(vars, pos, t, dt);
676
677 double storage = 0;
678 if (medium.hasProperty(MaterialPropertyLib::PropertyType::storage))
679 {
680 storage = medium[MaterialPropertyLib::PropertyType::storage]
681 .template value<double>(vars, pos, t, dt);
682 }
683
684 GlobalDimMatrixType const K_over_mu = K / mu;
685 GlobalDimVectorType const velocity =
686 _process_data.has_gravity
687 ? GlobalDimVectorType(-K_over_mu *
688 (dNdx * p_nodal_values - density * b))
689 : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
690
691 const double drho_dp =
692 phase[MaterialPropertyLib::PropertyType::density]
693 .template dValue<double>(
694 vars,
695 MaterialPropertyLib::Variable::liquid_phase_pressure,
696 pos, t, dt);
697
698 const double drho_dC =
699 phase[MaterialPropertyLib::PropertyType::density]
700 .template dValue<double>(
701 vars, MaterialPropertyLib::Variable::concentration, pos,
702 t, dt);
703
704 GlobalDimMatrixType const hydrodynamic_dispersion =
705 NumLib::computeHydrodynamicDispersion(
706 _process_data.stabilizer, _element.getID(),
707 pore_diffusion_coefficient, velocity, porosity,
708 solute_dispersivity_transverse,
709 solute_dispersivity_longitudinal);
710
711 const double R_times_phi(retardation_factor * porosity);
712 GlobalDimVectorType const mass_density_flow = velocity * density;
713 auto const N_t_N = (N.transpose() * N).eval();
714
715 if (_process_data.non_advective_form)
716 {
717 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
718 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
719 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
720 }
721 else
722 {
723 ip_flux_vector.emplace_back(mass_density_flow);
724 average_velocity_norm += velocity.norm();
725 }
726 MCC.noalias() += N_t_N * (R_times_phi * density * w);
727 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
728 KCC_Laplacian.noalias() +=
729 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
730
731 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
732
733 // Calculate Mpp, Kpp, and bp in the first loop over components
734 if (component_id == 0)
735 {
736 Mpp.noalias() +=
737 N_t_N * (porosity * drho_dp * w + density * storage * w);
738 Kpp.noalias() +=
739 dNdx.transpose() * K_over_mu * dNdx * (density * w);
740
741 if (_process_data.has_gravity)
742 {
743 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
744 (density * density * w);
745 }
746 }
747 }
748
749 if (!_process_data.non_advective_form)
750 {
751 NumLib::assembleAdvectionMatrix<
752 typename ShapeFunction::MeshElement>(
753 _process_data.stabilizer,
754 _ip_data,
755 _process_data.shape_matrix_cache,
756 ip_flux_vector,
757 average_velocity_norm /
758 static_cast<double>(n_integration_points),
759 KCC_Laplacian);
760 }
761
762 KCC.noalias() += KCC_Laplacian;
763 }
764
765 void assembleKCmCn(int const component_id, double const t, double const dt,
766 Eigen::Ref<LocalBlockMatrixType> KCmCn,
767 double const stoichiometric_coefficient,
768 double const kinetic_prefactor)
769 {
770 unsigned const n_integration_points =
771 _integration_method.getNumberOfPoints();
772
773 ParameterLib::SpatialPosition pos;
774 pos.setElementID(_element.getID());
775
776 MaterialPropertyLib::VariableArray vars;
777
778 auto const& medium =
779 *_process_data.media_map.getMedium(_element.getID());
780 auto const& phase = medium.phase("AqueousLiquid");
781 auto const& component = phase.component(
782 _transport_process_variables[component_id].get().getName());
783
784 auto const& Ns =
785 _process_data.shape_matrix_cache
786 .NsHigherOrder<typename ShapeFunction::MeshElement>();
787
788 for (unsigned ip(0); ip < n_integration_points; ++ip)
789 {
790 auto& ip_data = _ip_data[ip];
791 auto const& w = ip_data.integration_weight;
792 auto const& N = Ns[ip];
793 auto& porosity = ip_data.porosity;
794
795 // set position with N as the shape matrix at the current
796 // integration point
797 pos.setCoordinates(MathLib::Point3d(
798 NumLib::interpolateCoordinates<ShapeFunction,
799 ShapeMatricesType>(_element,
800 N)));
801
802 auto const retardation_factor =
803 component[MaterialPropertyLib::PropertyType::retardation_factor]
804 .template value<double>(vars, pos, t, dt);
805
806 porosity = medium[MaterialPropertyLib::PropertyType::porosity]
807 .template value<double>(vars, pos, t, dt);
808
809 auto const density =
810 phase[MaterialPropertyLib::PropertyType::density]
811 .template value<double>(vars, pos, t, dt);
812
813 KCmCn.noalias() -= N.transpose() * N *
814 (stoichiometric_coefficient * kinetic_prefactor *
815 retardation_factor * porosity * density * w);
816 }
817 }
818
819 void assembleForStaggeredScheme(double const t, double const dt,
820 Eigen::VectorXd const& local_x,
821 Eigen::VectorXd const& local_x_prev,
822 int const process_id,
823 std::vector<double>& local_M_data,
824 std::vector<double>& local_K_data,
825 std::vector<double>& local_b_data) override
826 {
827 if (process_id == _process_data.hydraulic_process_id)
828 {
829 assembleHydraulicEquation(t, dt, local_x, local_x_prev,
830 local_M_data, local_K_data, local_b_data);
831 }
832 else if (process_id == _process_data.thermal_process_id)
833 {
834 assembleHeatTransportEquation(t, dt, local_x, local_x_prev,
835 local_M_data, local_K_data,
836 local_b_data);
837 }
838 else
839 {
840 // Go for assembling in an order of transport process id.
841 assembleComponentTransportEquation(t, dt, local_x, local_x_prev,
842 local_M_data, local_K_data,
843 local_b_data, process_id);
844 }
845 }
846
847 void assembleHydraulicEquation(double const t,
848 double const dt,
849 Eigen::VectorXd const& local_x,
850 Eigen::VectorXd const& local_x_prev,
851 std::vector<double>& local_M_data,
852 std::vector<double>& local_K_data,
853 std::vector<double>& local_b_data)
854 {
855 auto const local_p =
856 local_x.template segment<pressure_size>(pressure_index);
857 auto const local_C = local_x.template segment<concentration_size>(
858 first_concentration_index);
859 auto const local_C_prev =
860 local_x_prev.segment<concentration_size>(first_concentration_index);
861
862 NodalVectorType local_T = getLocalTemperature(t, local_x);
863
864 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
865 local_M_data, pressure_size, pressure_size);
866 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
867 local_K_data, pressure_size, pressure_size);
868 auto local_b = MathLib::createZeroedVector<LocalSegmentVectorType>(
869 local_b_data, pressure_size);
870
871 unsigned const n_integration_points =
872 _integration_method.getNumberOfPoints();
873
874 ParameterLib::SpatialPosition pos;
875 pos.setElementID(_element.getID());
876
877 auto const& b =
878 _process_data
879 .projected_specific_body_force_vectors[_element.getID()];
880
881 auto const& medium =
882 *_process_data.media_map.getMedium(_element.getID());
883 auto const& phase = medium.phase("AqueousLiquid");
884
885 MaterialPropertyLib::VariableArray vars;
886 MaterialPropertyLib::VariableArray vars_prev;
887
888 auto const& Ns =
889 _process_data.shape_matrix_cache
890 .NsHigherOrder<typename ShapeFunction::MeshElement>();
891
892 for (unsigned ip(0); ip < n_integration_points; ++ip)
893 {
894 auto& ip_data = _ip_data[ip];
895 auto const& dNdx = ip_data.dNdx;
896 auto const& w = ip_data.integration_weight;
897 auto const& N = Ns[ip];
898 auto& porosity = ip_data.porosity;
899 auto const& porosity_prev = ip_data.porosity_prev;
900
901 double const C_int_pt = N.dot(local_C);
902 double const p_int_pt = N.dot(local_p);
903 double const T_int_pt = N.dot(local_T);
904
905 vars.concentration = C_int_pt;
906 vars.liquid_phase_pressure = p_int_pt;
907 vars.temperature = T_int_pt;
908
909 // porosity
910 {
911 vars_prev.porosity = porosity_prev;
912
913 porosity =
914 _process_data.chemically_induced_porosity_change
915 ? porosity_prev
916 : medium[MaterialPropertyLib::PropertyType::porosity]
917 .template value<double>(vars, vars_prev, pos, t,
918 dt);
919
920 vars.porosity = porosity;
921 }
922
923 // Use the fluid density model to compute the density
924 // TODO: Concentration of which component as one of arguments for
925 // calculation of fluid density
926 auto const density =
927 phase[MaterialPropertyLib::PropertyType::density]
928 .template value<double>(vars, pos, t, dt);
929
930 double storage = 0;
931 if (medium.hasProperty(MaterialPropertyLib::PropertyType::storage))
932 {
933 storage = medium[MaterialPropertyLib::PropertyType::storage]
934 .template value<double>(vars, pos, t, dt);
935 }
936
937 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
938 medium[MaterialPropertyLib::PropertyType::permeability].value(
939 vars, pos, t, dt));
940
941 // Use the viscosity model to compute the viscosity
942 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
943 .template value<double>(vars, pos, t, dt);
944
945 GlobalDimMatrixType const K_over_mu = K / mu;
946
947 const double drho_dp =
948 phase[MaterialPropertyLib::PropertyType::density]
949 .template dValue<double>(
950 vars,
951 MaterialPropertyLib::Variable::liquid_phase_pressure,
952 pos, t, dt);
953 const double drho_dC =
954 phase[MaterialPropertyLib::PropertyType::density]
955 .template dValue<double>(
956 vars, MaterialPropertyLib::Variable::concentration, pos,
957 t, dt);
958
959 // matrix assembly
960 local_M.noalias() +=
961 N.transpose() * N *
962 (porosity * drho_dp * w + density * storage * w);
963 local_K.noalias() +=
964 w * dNdx.transpose() * density * K_over_mu * dNdx;
965
966 if (_process_data.has_gravity)
967 {
968 local_b.noalias() +=
969 w * density * density * dNdx.transpose() * K_over_mu * b;
970 }
971
972 // coupling term
973 {
974 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
975
976 local_b.noalias() -=
977 N.transpose() * (porosity * drho_dC * C_dot * w);
978 }
979 }
980 }
981
982 void assembleHeatTransportEquation(double const t, double const dt,
983 Eigen::VectorXd const& local_x,
984 Eigen::VectorXd const& /*local_x_prev*/,
985 std::vector<double>& local_M_data,
986 std::vector<double>& local_K_data,
987 std::vector<double>& /*local_b_data*/)
988 {
989 // In the staggered HTC process, number of components might be non-zero.
990 assert(local_x.size() ==
991 pressure_size + temperature_size +
992 concentration_size *
993 static_cast<int>(_transport_process_variables.size()));
994
995 auto const local_p =
996 local_x.template segment<pressure_size>(pressure_index);
997 auto const local_T = getLocalTemperature(t, local_x);
998
999 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1000 local_M_data, temperature_size, temperature_size);
1001 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1002 local_K_data, temperature_size, temperature_size);
1003
1004 ParameterLib::SpatialPosition pos;
1005 pos.setElementID(this->_element.getID());
1006
1007 auto const& process_data = this->_process_data;
1008 auto const& medium =
1009 *process_data.media_map.getMedium(this->_element.getID());
1010 auto const& liquid_phase = medium.phase("AqueousLiquid");
1011
1012 auto const& b =
1013 _process_data
1014 .projected_specific_body_force_vectors[_element.getID()];
1015
1016 MaterialPropertyLib::VariableArray vars;
1017
1018 unsigned const n_integration_points =
1019 this->_integration_method.getNumberOfPoints();
1020
1021 std::vector<GlobalDimVectorType> ip_flux_vector;
1022 double average_velocity_norm = 0.0;
1023 ip_flux_vector.reserve(n_integration_points);
1024
1025 auto const& Ns =
1026 _process_data.shape_matrix_cache
1027 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1028
1029 for (unsigned ip(0); ip < n_integration_points; ip++)
1030 {
1031 auto const& ip_data = this->_ip_data[ip];
1032 auto const& dNdx = ip_data.dNdx;
1033 auto const& w = ip_data.integration_weight;
1034 auto const& N = Ns[ip];
1035
1036 ParameterLib::SpatialPosition const pos(
1037 {}, this->_element.getID(),
1038 MathLib::Point3d(
1039 NumLib::interpolateCoordinates<ShapeFunction,
1040 ShapeMatricesType>(_element,
1041 N)));
1042
1043 double p_at_xi = 0.;
1044 NumLib::shapeFunctionInterpolate(local_p, N, p_at_xi);
1045 double T_at_xi = 0.;
1046 NumLib::shapeFunctionInterpolate(local_T, N, T_at_xi);
1047
1048 vars.temperature = T_at_xi;
1049 vars.liquid_phase_pressure = p_at_xi;
1050
1051 vars.liquid_saturation = 1.0;
1052
1053 auto const porosity =
1054 medium.property(MaterialPropertyLib::PropertyType::porosity)
1055 .template value<double>(vars, pos, t, dt);
1056 vars.porosity = porosity;
1057
1058 // Use the fluid density model to compute the density
1059 auto const fluid_density =
1060 liquid_phase
1061 .property(MaterialPropertyLib::PropertyType::density)
1062 .template value<double>(vars, pos, t, dt);
1063 vars.density = fluid_density;
1064 auto const specific_heat_capacity_fluid =
1065 liquid_phase
1066 .property(MaterialPropertyLib::specific_heat_capacity)
1067 .template value<double>(vars, pos, t, dt);
1068
1069 // Assemble mass matrix
1070 local_M.noalias() +=
1071 N.transpose() * N *
1072 (this->getHeatEnergyCoefficient(vars, porosity, fluid_density,
1073 specific_heat_capacity_fluid,
1074 pos, t, dt) *
1075 w);
1076
1077 // Assemble Laplace matrix
1078 auto const viscosity =
1079 liquid_phase
1080 .property(MaterialPropertyLib::PropertyType::viscosity)
1081 .template value<double>(vars, pos, t, dt);
1082
1083 auto const intrinsic_permeability =
1084 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1085 medium
1086 .property(
1087 MaterialPropertyLib::PropertyType::permeability)
1088 .value(vars, pos, t, dt));
1089
1090 GlobalDimMatrixType const K_over_mu =
1091 intrinsic_permeability / viscosity;
1092 GlobalDimVectorType const velocity =
1093 process_data.has_gravity
1094 ? GlobalDimVectorType(-K_over_mu *
1095 (dNdx * local_p - fluid_density * b))
1096 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1097
1098 GlobalDimMatrixType const thermal_conductivity_dispersivity =
1099 this->getThermalConductivityDispersivity(
1100 vars, fluid_density, specific_heat_capacity_fluid, velocity,
1101 pos, t, dt);
1102
1103 local_K.noalias() +=
1104 w * dNdx.transpose() * thermal_conductivity_dispersivity * dNdx;
1105
1106 ip_flux_vector.emplace_back(velocity * fluid_density *
1107 specific_heat_capacity_fluid);
1108 average_velocity_norm += velocity.norm();
1109 }
1110
1111 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1112 process_data.stabilizer, this->_ip_data,
1113 _process_data.shape_matrix_cache, ip_flux_vector,
1114 average_velocity_norm / static_cast<double>(n_integration_points),
1115 local_K);
1116 }
1117
1118 void assembleComponentTransportEquation(
1119 double const t, double const dt, Eigen::VectorXd const& local_x,
1120 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_M_data,
1121 std::vector<double>& local_K_data,
1122 std::vector<double>& /*local_b_data*/, int const transport_process_id)
1123 {
1124 assert(static_cast<int>(local_x.size()) ==
1125 pressure_size +
1126 concentration_size *
1127 static_cast<int>(_transport_process_variables.size()) +
1128 (_process_data.isothermal ? 0 : temperature_size));
1129
1130 auto const local_p =
1131 local_x.template segment<pressure_size>(pressure_index);
1132
1133 NodalVectorType local_T = getLocalTemperature(t, local_x);
1134
1135 auto const local_C = local_x.template segment<concentration_size>(
1136 first_concentration_index +
1137 (transport_process_id - (_process_data.isothermal ? 1 : 2)) *
1138 concentration_size);
1139 auto const local_p_prev =
1140 local_x_prev.segment<pressure_size>(pressure_index);
1141
1142 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1143 local_M_data, concentration_size, concentration_size);
1144 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1145 local_K_data, concentration_size, concentration_size);
1146
1147 LocalBlockMatrixType KCC_Laplacian =
1148 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1149
1150 unsigned const n_integration_points =
1151 _integration_method.getNumberOfPoints();
1152
1153 std::vector<GlobalDimVectorType> ip_flux_vector;
1154 double average_velocity_norm = 0.0;
1155 if (!_process_data.non_advective_form)
1156 {
1157 ip_flux_vector.reserve(n_integration_points);
1158 }
1159
1160 ParameterLib::SpatialPosition pos;
1161 pos.setElementID(_element.getID());
1162
1163 auto const& b =
1164 _process_data
1165 .projected_specific_body_force_vectors[_element.getID()];
1166
1167 MaterialPropertyLib::VariableArray vars;
1168 MaterialPropertyLib::VariableArray vars_prev;
1169
1170 auto const& medium =
1171 *_process_data.media_map.getMedium(_element.getID());
1172 auto const& phase = medium.phase("AqueousLiquid");
1173 auto const component_id =
1174 transport_process_id - (_process_data.isothermal ? 1 : 2);
1175 auto const& component = phase.component(
1176 _transport_process_variables[component_id].get().getName());
1177
1178 auto const& Ns =
1179 _process_data.shape_matrix_cache
1180 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1181
1182 for (unsigned ip(0); ip < n_integration_points; ++ip)
1183 {
1184 auto& ip_data = _ip_data[ip];
1185 auto const& dNdx = ip_data.dNdx;
1186 auto const& w = ip_data.integration_weight;
1187 auto const& N = Ns[ip];
1188 auto& porosity = ip_data.porosity;
1189 auto const& porosity_prev = ip_data.porosity_prev;
1190
1191 double const C_int_pt = N.dot(local_C);
1192 double const p_int_pt = N.dot(local_p);
1193 double const T_int_pt = N.dot(local_T);
1194
1195 vars.concentration = C_int_pt;
1196 vars.liquid_phase_pressure = p_int_pt;
1197 vars.temperature = T_int_pt;
1198
1199 if (_process_data.temperature)
1200 {
1201 vars.temperature = N.dot(local_T);
1202 }
1203
1204 // porosity
1205 {
1206 vars_prev.porosity = porosity_prev;
1207
1208 porosity =
1209 _process_data.chemically_induced_porosity_change
1210 ? porosity_prev
1211 : medium[MaterialPropertyLib::PropertyType::porosity]
1212 .template value<double>(vars, vars_prev, pos, t,
1213 dt);
1214
1215 vars.porosity = porosity;
1216 }
1217
1218 auto const& retardation_factor =
1219 component[MaterialPropertyLib::PropertyType::retardation_factor]
1220 .template value<double>(vars, pos, t, dt);
1221
1222 auto const& solute_dispersivity_transverse = medium.template value<
1223 double>(
1224 MaterialPropertyLib::PropertyType::transversal_dispersivity);
1225 auto const& solute_dispersivity_longitudinal =
1226 medium.template value<double>(
1227 MaterialPropertyLib::PropertyType::
1228 longitudinal_dispersivity);
1229
1230 // Use the fluid density model to compute the density
1231 auto const density =
1232 phase[MaterialPropertyLib::PropertyType::density]
1233 .template value<double>(vars, pos, t, dt);
1234 auto const decay_rate =
1235 component[MaterialPropertyLib::PropertyType::decay_rate]
1236 .template value<double>(vars, pos, t, dt);
1237
1238 auto const& pore_diffusion_coefficient =
1239 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1240 component[MaterialPropertyLib::PropertyType::pore_diffusion]
1241 .value(vars, pos, t, dt));
1242
1243 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1244 medium[MaterialPropertyLib::PropertyType::permeability].value(
1245 vars, pos, t, dt));
1246 // Use the viscosity model to compute the viscosity
1247 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1248 .template value<double>(vars, pos, t, dt);
1249
1250 GlobalDimMatrixType const K_over_mu = K / mu;
1251 GlobalDimVectorType const velocity =
1252 _process_data.has_gravity
1253 ? GlobalDimVectorType(-K_over_mu *
1254 (dNdx * local_p - density * b))
1255 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1256
1257 GlobalDimMatrixType const hydrodynamic_dispersion =
1258 NumLib::computeHydrodynamicDispersion(
1259 _process_data.stabilizer, _element.getID(),
1260 pore_diffusion_coefficient, velocity, porosity,
1261 solute_dispersivity_transverse,
1262 solute_dispersivity_longitudinal);
1263
1264 double const R_times_phi = retardation_factor * porosity;
1265 auto const N_t_N = (N.transpose() * N).eval();
1266
1267 if (_process_data.non_advective_form)
1268 {
1269 const double drho_dC =
1270 phase[MaterialPropertyLib::PropertyType::density]
1271 .template dValue<double>(
1272 vars, MaterialPropertyLib::Variable::concentration,
1273 pos, t, dt);
1274 local_M.noalias() +=
1275 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1276 }
1277
1278 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1279
1280 // coupling term
1281 if (_process_data.non_advective_form)
1282 {
1283 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1284
1285 const double drho_dp =
1286 phase[MaterialPropertyLib::PropertyType::density]
1287 .template dValue<double>(vars,
1288 MaterialPropertyLib::Variable::
1289 liquid_phase_pressure,
1290 pos, t, dt);
1291
1292 local_K.noalias() +=
1293 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1294 dNdx.transpose() * velocity * N * (density * w);
1295 }
1296 else
1297 {
1298 ip_flux_vector.emplace_back(velocity * density);
1299 average_velocity_norm += velocity.norm();
1300 }
1301 local_K.noalias() +=
1302 N_t_N * (decay_rate * R_times_phi * density * w);
1303
1304 KCC_Laplacian.noalias() += dNdx.transpose() *
1305 hydrodynamic_dispersion * dNdx *
1306 (density * w);
1307 }
1308
1309 if (!_process_data.non_advective_form)
1310 {
1311 NumLib::assembleAdvectionMatrix<
1312 typename ShapeFunction::MeshElement>(
1313 _process_data.stabilizer, _ip_data,
1314 _process_data.shape_matrix_cache, ip_flux_vector,
1315 average_velocity_norm /
1316 static_cast<double>(n_integration_points),
1317 KCC_Laplacian);
1318 }
1319 local_K.noalias() += KCC_Laplacian;
1320 }
1321
1322 void assembleWithJacobianForStaggeredScheme(
1323 double const t, double const dt, Eigen::VectorXd const& local_x,
1324 Eigen::VectorXd const& local_x_prev, int const process_id,
1325 std::vector<double>& local_b_data,
1326 std::vector<double>& local_Jac_data) override
1327 {
1328 if (process_id == _process_data.hydraulic_process_id)
1329 {
1330 assembleWithJacobianHydraulicEquation(t, dt, local_x, local_x_prev,
1331 local_b_data, local_Jac_data);
1332 }
1333 else
1334 {
1335 int const component_id = process_id - 1;
1336 assembleWithJacobianComponentTransportEquation(
1337 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1338 component_id);
1339 }
1340 }
1341
1342 void assembleWithJacobianHydraulicEquation(
1343 double const t, double const dt, Eigen::VectorXd const& local_x,
1344 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1345 std::vector<double>& local_Jac_data)
1346 {
1347 auto const p = local_x.template segment<pressure_size>(pressure_index);
1348 auto const c = local_x.template segment<concentration_size>(
1349 first_concentration_index);
1350
1351 auto const p_prev = local_x_prev.segment<pressure_size>(pressure_index);
1352 auto const c_prev =
1353 local_x_prev.segment<concentration_size>(first_concentration_index);
1354
1355 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1356 local_Jac_data, pressure_size, pressure_size);
1357 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1358 local_b_data, pressure_size);
1359
1360 unsigned const n_integration_points =
1361 _integration_method.getNumberOfPoints();
1362
1363 ParameterLib::SpatialPosition pos;
1364 pos.setElementID(_element.getID());
1365 auto const& b =
1366 _process_data
1367 .projected_specific_body_force_vectors[_element.getID()];
1368
1369 auto const& medium =
1370 *_process_data.media_map.getMedium(_element.getID());
1371 auto const& phase = medium.phase("AqueousLiquid");
1372
1373 MaterialPropertyLib::VariableArray vars;
1374 MaterialPropertyLib::VariableArray vars_prev;
1375
1376 auto const& Ns =
1377 _process_data.shape_matrix_cache
1378 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1379
1380 for (unsigned ip(0); ip < n_integration_points; ++ip)
1381 {
1382 auto& ip_data = _ip_data[ip];
1383 auto const& dNdx = ip_data.dNdx;
1384 auto const& w = ip_data.integration_weight;
1385 auto const& N = Ns[ip];
1386 auto& phi = ip_data.porosity;
1387 auto const& phi_prev = ip_data.porosity_prev;
1388
1389 double const p_ip = N.dot(p);
1390 double const c_ip = N.dot(c);
1391
1392 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1393
1394 vars.liquid_phase_pressure = p_ip;
1395 vars.concentration = c_ip;
1396
1397 // porosity
1398 {
1399 vars_prev.porosity = phi_prev;
1400
1401 phi = _process_data.chemically_induced_porosity_change
1402 ? phi_prev
1403 : medium[MaterialPropertyLib::PropertyType::porosity]
1404 .template value<double>(vars, vars_prev, pos, t,
1405 dt);
1406
1407 vars.porosity = phi;
1408 }
1409
1410 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1411 .template value<double>(vars, pos, t, dt);
1412
1413 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1414 medium[MaterialPropertyLib::PropertyType::permeability].value(
1415 vars, pos, t, dt));
1416
1417 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1418 .template value<double>(vars, pos, t, dt);
1419
1420 auto const drho_dp =
1421 phase[MaterialPropertyLib::PropertyType::density]
1422 .template dValue<double>(
1423 vars,
1424 MaterialPropertyLib::Variable::liquid_phase_pressure,
1425 pos, t, dt);
1426 auto const drho_dc =
1427 phase[MaterialPropertyLib::PropertyType::density]
1428 .template dValue<double>(
1429 vars, MaterialPropertyLib::Variable::concentration, pos,
1430 t, dt);
1431
1432 // matrix assembly
1433 local_Jac.noalias() +=
1434 N.transpose() * N * (phi * drho_dp / dt * w) +
1435 w * dNdx.transpose() * rho * k / mu * dNdx;
1436
1437 local_rhs.noalias() -=
1438 N.transpose() * (drho_dp * N * p_prev + drho_dc * cdot_ip) *
1439 (phi * w) +
1440 dNdx.transpose() * k / mu * dNdx * p * (rho * w);
1441
1442 if (_process_data.has_gravity)
1443 {
1444 local_rhs.noalias() +=
1445 w * rho * dNdx.transpose() * k / mu * rho * b;
1446 }
1447 }
1448 }
1449
1450 void assembleWithJacobianComponentTransportEquation(
1451 double const t, double const dt, Eigen::VectorXd const& local_x,
1452 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1453 std::vector<double>& local_Jac_data, int const component_id)
1454 {
1455 auto const concentration_index =
1456 first_concentration_index + component_id * concentration_size;
1457
1458 auto const p = local_x.template segment<pressure_size>(pressure_index);
1459 auto const c =
1460 local_x.template segment<concentration_size>(concentration_index);
1461 auto const c_prev =
1462 local_x_prev.segment<concentration_size>(concentration_index);
1463
1464 NodalVectorType T;
1465 if (_process_data.temperature)
1466 {
1467 T = _process_data.temperature->getNodalValuesOnElement(_element, t);
1468 }
1469
1470 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1471 local_Jac_data, concentration_size, concentration_size);
1472 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1473 local_b_data, concentration_size);
1474
1475 LocalBlockMatrixType KCC_Laplacian =
1476 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1477
1478 unsigned const n_integration_points =
1479 _integration_method.getNumberOfPoints();
1480
1481 std::vector<GlobalDimVectorType> ip_flux_vector;
1482 double average_velocity_norm = 0.0;
1483 ip_flux_vector.reserve(n_integration_points);
1484
1485 ParameterLib::SpatialPosition pos;
1486 pos.setElementID(_element.getID());
1487
1488 auto const& b =
1489 _process_data
1490 .projected_specific_body_force_vectors[_element.getID()];
1491
1492 MaterialPropertyLib::VariableArray vars;
1493 MaterialPropertyLib::VariableArray vars_prev;
1494
1495 auto const& medium =
1496 *_process_data.media_map.getMedium(_element.getID());
1497 auto const& phase = medium.phase("AqueousLiquid");
1498 auto const& component = phase.component(
1499 _transport_process_variables[component_id].get().getName());
1500
1501 auto const& Ns =
1502 _process_data.shape_matrix_cache
1503 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1504
1505 for (unsigned ip(0); ip < n_integration_points; ++ip)
1506 {
1507 auto& ip_data = _ip_data[ip];
1508 auto const& dNdx = ip_data.dNdx;
1509 auto const& w = ip_data.integration_weight;
1510 auto const& N = Ns[ip];
1511 auto& phi = ip_data.porosity;
1512 auto const& phi_prev = ip_data.porosity_prev;
1513
1514 double const p_ip = N.dot(p);
1515 double const c_ip = N.dot(c);
1516
1517 vars.liquid_phase_pressure = p_ip;
1518 vars.concentration = c_ip;
1519
1520 if (_process_data.temperature)
1521 {
1522 vars.temperature = N.dot(T);
1523 }
1524
1525 // porosity
1526 {
1527 vars_prev.porosity = phi_prev;
1528
1529 phi = _process_data.chemically_induced_porosity_change
1530 ? phi_prev
1531 : medium[MaterialPropertyLib::PropertyType::porosity]
1532 .template value<double>(vars, vars_prev, pos, t,
1533 dt);
1534
1535 vars.porosity = phi;
1536 }
1537
1538 auto const R =
1539 component[MaterialPropertyLib::PropertyType::retardation_factor]
1540 .template value<double>(vars, pos, t, dt);
1541
1542 auto const alpha_T = medium.template value<double>(
1543 MaterialPropertyLib::PropertyType::transversal_dispersivity);
1544 auto const alpha_L = medium.template value<double>(
1545 MaterialPropertyLib::PropertyType::longitudinal_dispersivity);
1546
1547 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1548 .template value<double>(vars, pos, t, dt);
1549 // first-order decay constant
1550 auto const alpha =
1551 component[MaterialPropertyLib::PropertyType::decay_rate]
1552 .template value<double>(vars, pos, t, dt);
1553
1554 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1555 component[MaterialPropertyLib::PropertyType::pore_diffusion]
1556 .value(vars, pos, t, dt));
1557
1558 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1559 medium[MaterialPropertyLib::PropertyType::permeability].value(
1560 vars, pos, t, dt));
1561 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1562 .template value<double>(vars, pos, t, dt);
1563 // Darcy flux
1564 GlobalDimVectorType const q =
1565 _process_data.has_gravity
1566 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1567 : GlobalDimVectorType(-k / mu * dNdx * p);
1568
1569 GlobalDimMatrixType const D = NumLib::computeHydrodynamicDispersion(
1570 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1571 alpha_L);
1572
1573 // matrix assembly
1574 local_Jac.noalias() +=
1575 N.transpose() * N * (rho * phi * R * (alpha + 1 / dt) * w);
1576
1577 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1578
1579 auto const cdot = (c - c_prev) / dt;
1580 local_rhs.noalias() -=
1581 N.transpose() * N * (cdot + alpha * c) * (rho * phi * R * w);
1582
1583 ip_flux_vector.emplace_back(q * rho);
1584 average_velocity_norm += q.norm();
1585 }
1586
1587 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1588 _process_data.stabilizer, _ip_data,
1589 _process_data.shape_matrix_cache, ip_flux_vector,
1590 average_velocity_norm / static_cast<double>(n_integration_points),
1591 KCC_Laplacian);
1592
1593 local_rhs.noalias() -= KCC_Laplacian * c;
1594
1595 local_Jac.noalias() += KCC_Laplacian;
1596 }
1597
1598 void assembleReactionEquationConcrete(
1599 double const t, double const dt, Eigen::VectorXd const& local_x,
1600 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1601 std::vector<double>& local_b_data,
1602 int const transport_process_id) override
1603 {
1604 auto const local_C = local_x.template segment<concentration_size>(
1605 first_concentration_index +
1606 (transport_process_id - 1) * concentration_size);
1607
1608 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1609 local_M_data, concentration_size, concentration_size);
1610 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1611 local_K_data, concentration_size, concentration_size);
1612 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1613 local_b_data, concentration_size);
1614
1615 unsigned const n_integration_points =
1616 _integration_method.getNumberOfPoints();
1617
1618 ParameterLib::SpatialPosition pos;
1619 pos.setElementID(_element.getID());
1620
1621 MaterialPropertyLib::VariableArray vars;
1622 MaterialPropertyLib::VariableArray vars_prev;
1623
1624 auto const& medium =
1625 *_process_data.media_map.getMedium(_element.getID());
1626 auto const component_id = transport_process_id - 1;
1627
1628 auto const& Ns =
1629 _process_data.shape_matrix_cache
1630 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1631
1632 for (unsigned ip(0); ip < n_integration_points; ++ip)
1633 {
1634 auto& ip_data = _ip_data[ip];
1635 auto const w = ip_data.integration_weight;
1636 auto const& N = Ns[ip];
1637 auto& porosity = ip_data.porosity;
1638 auto const& porosity_prev = ip_data.porosity_prev;
1639 auto const chemical_system_id = ip_data.chemical_system_id;
1640
1641 double C_int_pt = 0.0;
1642 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
1643
1644 vars.concentration = C_int_pt;
1645
1646 auto const porosity_dot = (porosity - porosity_prev) / dt;
1647
1648 // porosity
1649 {
1650 vars_prev.porosity = porosity_prev;
1651
1652 porosity =
1653 _process_data.chemically_induced_porosity_change
1654 ? porosity_prev
1655 : medium[MaterialPropertyLib::PropertyType::porosity]
1656 .template value<double>(vars, vars_prev, pos, t,
1657 dt);
1658 }
1659
1660 local_M.noalias() += w * N.transpose() * porosity * N;
1661
1662 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1663
1664 if (chemical_system_id == -1)
1665 {
1666 continue;
1667 }
1668
1669 auto const C_post_int_pt =
1670 _process_data.chemical_solver_interface->getConcentration(
1671 component_id, chemical_system_id);
1672
1673 local_b.noalias() += N.transpose() * ((C_post_int_pt - C_int_pt) /
1674 dt * porosity * w);
1675 }
1676 }
1677
1678 std::vector<double> const& getIntPtLiquidDensity(
1679 const double t,
1680 std::vector<GlobalVector*> const& x,
1681 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
1682 std::vector<double>& cache) const override
1683 {
1684 assert(x.size() == dof_table.size());
1685
1686 auto const n_processes = x.size();
1687 std::vector<std::vector<double>> local_x;
1688 local_x.reserve(n_processes);
1689
1690 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1691 {
1692 auto const indices =
1693 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1694 assert(!indices.empty());
1695 local_x.push_back(x[process_id]->get(indices));
1696 }
1697
1698 // only one process, must be monolithic.
1699 if (n_processes == 1)
1700 {
1701 auto const local_p = Eigen::Map<const NodalVectorType>(
1702 &local_x[0][pressure_index], pressure_size);
1703 auto const local_C = Eigen::Map<const NodalVectorType>(
1704 &local_x[0][first_concentration_index], concentration_size);
1705 NodalVectorType local_T;
1706 local_T.setConstant(ShapeFunction::NPOINTS,
1707 std::numeric_limits<double>::quiet_NaN());
1708 int const temperature_index =
1709 _process_data.isothermal ? -1 : ShapeFunction::NPOINTS;
1710 if (temperature_index != -1)
1711 {
1712 local_T = Eigen::Map<const NodalVectorType>(
1713 &local_x[0][temperature_index], temperature_size);
1714 }
1715 return calculateIntPtLiquidDensity(t, local_p, local_C, local_T,
1716 cache);
1717 }
1718
1719 // multiple processes, must be staggered.
1720 {
1721 constexpr int pressure_process_id = 0;
1722 int concentration_process_id = 1;
1723 // Normally temperature process is not there,
1724 // hence set the default temperature index to -1
1725 int temperature_process_id = -1;
1726
1727 // check whether temperature process exists
1728 if (!_process_data.isothermal)
1729 {
1730 // if temperature process exists, its id is 1
1731 temperature_process_id = 1;
1732 // then the concentration index shifts to 2
1733 concentration_process_id = 2;
1734 }
1735
1736 auto const local_p = Eigen::Map<const NodalVectorType>(
1737 &local_x[pressure_process_id][0], pressure_size);
1738 auto const local_C = Eigen::Map<const NodalVectorType>(
1739 &local_x[concentration_process_id][0], concentration_size);
1740 NodalVectorType local_T;
1741 local_T.setConstant(ShapeFunction::NPOINTS,
1742 std::numeric_limits<double>::quiet_NaN());
1743 if (temperature_process_id != -1)
1744 {
1745 local_T = Eigen::Map<const NodalVectorType>(
1746 &local_x[temperature_process_id][0], temperature_size);
1747 }
1748 return calculateIntPtLiquidDensity(t, local_p, local_C, local_T,
1749 cache);
1750 }
1751 }
1752
1753 std::vector<double> const& calculateIntPtLiquidDensity(
1754 const double t,
1755 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
1756 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
1757 Eigen::Ref<const NodalVectorType> const& T_nodal_values,
1758 std::vector<double>& cache) const
1759 {
1760 auto const n_integration_points =
1761 _integration_method.getNumberOfPoints();
1762
1763 cache.clear();
1764 auto cache_vec = MathLib::createZeroedVector<Eigen::VectorXd>(
1765 cache, n_integration_points);
1766
1767 ParameterLib::SpatialPosition pos;
1768 pos.setElementID(_element.getID());
1769
1770 MaterialPropertyLib::VariableArray vars;
1771
1772 auto const& medium =
1773 *_process_data.media_map.getMedium(_element.getID());
1774 auto const& phase = medium.phase("AqueousLiquid");
1775
1776 auto const& Ns =
1777 _process_data.shape_matrix_cache
1778 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1779
1780 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1781 {
1782 auto const& N = Ns[ip];
1783
1784 double C_int_pt = 0.0;
1785 double p_int_pt = 0.0;
1786 double T_int_pt = 0.0;
1787
1788 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1789 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1790 NumLib::shapeFunctionInterpolate(T_nodal_values, N, T_int_pt);
1791
1792 vars.concentration = C_int_pt;
1793 vars.liquid_phase_pressure = p_int_pt;
1794 vars.temperature = T_int_pt;
1795
1796 // TODO (naumov) Temporary value not used by current material
1797 // models. Need extension of secondary variables interface.
1798 double const dt = std::numeric_limits<double>::quiet_NaN();
1799
1800 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
1801 .template value<double>(vars, pos, t, dt);
1802 cache_vec[ip] = rho_w;
1803 }
1804
1805 return cache;
1806 }
1807
1808 std::vector<double> const& getIntPtDarcyVelocity(
1809 const double t,
1810 std::vector<GlobalVector*> const& x,
1811 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
1812 std::vector<double>& cache) const override
1813 {
1814 assert(x.size() == dof_table.size());
1815
1816 auto const n_processes = x.size();
1817 std::vector<std::vector<double>> local_x;
1818 local_x.reserve(n_processes);
1819
1820 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1821 {
1822 auto const indices =
1823 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1824 assert(!indices.empty());
1825 local_x.push_back(x[process_id]->get(indices));
1826 }
1827
1828 // only one process, must be monolithic.
1829 if (n_processes == 1)
1830 {
1831 auto const local_p = Eigen::Map<const NodalVectorType>(
1832 &local_x[0][pressure_index], pressure_size);
1833 auto const local_C = Eigen::Map<const NodalVectorType>(
1834 &local_x[0][first_concentration_index], concentration_size);
1835 NodalVectorType local_T;
1836 local_T.setConstant(ShapeFunction::NPOINTS,
1837 std::numeric_limits<double>::quiet_NaN());
1838 int const temperature_index =
1839 _process_data.isothermal ? -1 : ShapeFunction::NPOINTS;
1840 if (temperature_index != -1)
1841 {
1842 local_T = Eigen::Map<const NodalVectorType>(
1843 &local_x[0][temperature_index], temperature_size);
1844 }
1845 return calculateIntPtDarcyVelocity(t, local_p, local_C, local_T,
1846 cache);
1847 }
1848
1849 // multiple processes, must be staggered.
1850 {
1851 constexpr int pressure_process_id = 0;
1852 int concentration_process_id = 1;
1853 // Normally temperature process is not there,
1854 // hence set the default temperature index to -1
1855 int temperature_process_id = -1;
1856
1857 // check whether temperature process exists
1858 if (!_process_data.isothermal)
1859 {
1860 // if temperature process exists, its id is 1
1861 temperature_process_id = 1;
1862 // then the concentration index shifts to 2
1863 concentration_process_id = 2;
1864 }
1865
1866 auto const local_p = Eigen::Map<const NodalVectorType>(
1867 &local_x[pressure_process_id][0], pressure_size);
1868 auto const local_C = Eigen::Map<const NodalVectorType>(
1869 &local_x[concentration_process_id][0], concentration_size);
1870 NodalVectorType local_T;
1871 local_T.setConstant(ShapeFunction::NPOINTS,
1872 std::numeric_limits<double>::quiet_NaN());
1873 if (temperature_process_id != -1)
1874 {
1875 local_T = Eigen::Map<const NodalVectorType>(
1876 &local_x[temperature_process_id][0], temperature_size);
1877 }
1878 return calculateIntPtDarcyVelocity(t, local_p, local_C, local_T,
1879 cache);
1880 }
1881 }
1882
1883 std::vector<double> const& calculateIntPtDarcyVelocity(
1884 const double t,
1885 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
1886 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
1887 Eigen::Ref<const NodalVectorType> const& T_nodal_values,
1888 std::vector<double>& cache) const
1889 {
1890 auto const n_integration_points =
1891 _integration_method.getNumberOfPoints();
1892
1893 cache.clear();
1894 auto cache_mat = MathLib::createZeroedMatrix<
1895 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1896 cache, GlobalDim, n_integration_points);
1897
1898 ParameterLib::SpatialPosition pos;
1899 pos.setElementID(_element.getID());
1900
1901 auto const& b =
1902 _process_data
1903 .projected_specific_body_force_vectors[_element.getID()];
1904
1905 MaterialPropertyLib::VariableArray vars;
1906
1907 auto const& medium =
1908 *_process_data.media_map.getMedium(_element.getID());
1909 auto const& phase = medium.phase("AqueousLiquid");
1910
1911 auto const& Ns =
1912 _process_data.shape_matrix_cache
1913 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1914
1915 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1916 {
1917 auto const& ip_data = _ip_data[ip];
1918 auto const& dNdx = ip_data.dNdx;
1919 auto const& N = Ns[ip];
1920 auto const& porosity = ip_data.porosity;
1921
1922 double C_int_pt = 0.0;
1923 double p_int_pt = 0.0;
1924 double T_int_pt = 0.0;
1925
1926 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1927 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1928 NumLib::shapeFunctionInterpolate(T_nodal_values, N, T_int_pt);
1929
1930 vars.concentration = C_int_pt;
1931 vars.liquid_phase_pressure = p_int_pt;
1932 vars.porosity = porosity;
1933 vars.temperature = T_int_pt;
1934
1935 // TODO (naumov) Temporary value not used by current material
1936 // models. Need extension of secondary variables interface.
1937 double const dt = std::numeric_limits<double>::quiet_NaN();
1938 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1939 medium[MaterialPropertyLib::PropertyType::permeability].value(
1940 vars, pos, t, dt));
1941 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1942 .template value<double>(vars, pos, t, dt);
1943 GlobalDimMatrixType const K_over_mu = K / mu;
1944
1945 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1946 if (_process_data.has_gravity)
1947 {
1948 auto const rho_w =
1949 phase[MaterialPropertyLib::PropertyType::density]
1950 .template value<double>(vars, pos, t, dt);
1951 // here it is assumed that the vector b is directed 'downwards'
1952 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1953 }
1954 }
1955
1956 return cache;
1957 }
1958
1959 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
1960 const unsigned integration_point) const override
1961 {
1962 auto const& N = _process_data.shape_matrix_cache.NsHigherOrder<
1963 typename ShapeFunction::MeshElement>()[integration_point];
1964
1965 // assumes N is stored contiguously in memory
1966 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1967 }
1968
1969 Eigen::Vector3d getFlux(MathLib::Point3d const& pnt_local_coords,
1970 double const t,
1971 std::vector<double> const& local_x) const override
1972 {
1973 auto const local_p = Eigen::Map<const NodalVectorType>(
1974 &local_x[pressure_index], pressure_size);
1975 auto const local_C = Eigen::Map<const NodalVectorType>(
1976 &local_x[first_concentration_index], concentration_size);
1977
1978 // Eval shape matrices at given point
1979 // Note: Axial symmetry is set to false here, because we only need dNdx
1980 // here, which is not affected by axial symmetry.
1981 auto const shape_matrices =
1982 NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
1983 GlobalDim>(
1984 _element, false /*is_axially_symmetric*/,
1985 std::array{pnt_local_coords})[0];
1986
1987 ParameterLib::SpatialPosition pos;
1988 pos.setElementID(_element.getID());
1989 auto const& b =
1990 _process_data
1991 .projected_specific_body_force_vectors[_element.getID()];
1992
1993 MaterialPropertyLib::VariableArray vars;
1994
1995 auto const& medium =
1996 *_process_data.media_map.getMedium(_element.getID());
1997 auto const& phase = medium.phase("AqueousLiquid");
1998
1999 // local_x contains the local concentration and pressure values
2000 double c_int_pt;
2001 NumLib::shapeFunctionInterpolate(local_C, shape_matrices.N, c_int_pt);
2002 vars.concentration = c_int_pt;
2003
2004 double p_int_pt;
2005 NumLib::shapeFunctionInterpolate(local_p, shape_matrices.N, p_int_pt);
2006 vars.liquid_phase_pressure = p_int_pt;
2007
2008 // TODO (naumov) Temporary value not used by current material models.
2009 // Need extension of secondary variables interface.
2010 double const dt = std::numeric_limits<double>::quiet_NaN();
2011 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
2012 medium[MaterialPropertyLib::PropertyType::permeability].value(
2013 vars, pos, t, dt));
2014
2015 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
2016 .template value<double>(vars, pos, t, dt);
2017 GlobalDimMatrixType const K_over_mu = K / mu;
2018
2019 GlobalDimVectorType q = -K_over_mu * shape_matrices.dNdx * local_p;
2020 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
2021 .template value<double>(vars, pos, t, dt);
2022 if (_process_data.has_gravity)
2023 {
2024 q += K_over_mu * rho_w * b;
2025 }
2026 Eigen::Vector3d flux(0.0, 0.0, 0.0);
2027 flux.head<GlobalDim>() = rho_w * q;
2028 return flux;
2029 }
2030
2031 void computeSecondaryVariableConcrete(
2032 double const t,
2033 double const /*dt*/,
2034 Eigen::VectorXd const& local_x,
2035 Eigen::VectorXd const& /*local_x_prev*/) override
2036 {
2037 auto const local_p =
2038 local_x.template segment<pressure_size>(pressure_index);
2039 auto const local_C = local_x.template segment<concentration_size>(
2040 first_concentration_index);
2041 auto const local_T = getLocalTemperature(t, local_x);
2042
2043 std::vector<double> ele_velocity;
2044 calculateIntPtDarcyVelocity(t, local_p, local_C, local_T, ele_velocity);
2045
2046 auto const n_integration_points =
2047 _integration_method.getNumberOfPoints();
2048 auto const ele_velocity_mat =
2049 MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
2050
2051 auto const ele_id = _element.getID();
2052 Eigen::Map<LocalVectorType>(
2053 &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
2054 GlobalDim) =
2055 ele_velocity_mat.rowwise().sum() / n_integration_points;
2056 }
2057
2058 void computeReactionRelatedSecondaryVariable(
2059 std::size_t const ele_id) override
2060 {
2061 auto const n_integration_points =
2062 _integration_method.getNumberOfPoints();
2063
2064 if (_process_data.chemically_induced_porosity_change)
2065 {
2066 auto const& medium = *_process_data.media_map.getMedium(ele_id);
2067
2068 for (auto& ip_data : _ip_data)
2069 {
2070 ip_data.porosity = ip_data.porosity_prev;
2071
2072 _process_data.chemical_solver_interface
2073 ->updatePorosityPostReaction(ip_data.chemical_system_id,
2074 medium, ip_data.porosity);
2075 }
2076
2077 (*_process_data.mesh_prop_porosity)[ele_id] =
2078 std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
2079 [](double const s, auto const& ip)
2080 { return s + ip.porosity; }) /
2081 n_integration_points;
2082 }
2083
2084 std::vector<GlobalIndexType> chemical_system_indices;
2085 chemical_system_indices.reserve(n_integration_points);
2086 std::transform(_ip_data.begin(), _ip_data.end(),
2087 std::back_inserter(chemical_system_indices),
2088 [](auto const& ip_data)
2089 { return ip_data.chemical_system_id; });
2090
2091 _process_data.chemical_solver_interface->computeSecondaryVariable(
2092 ele_id, chemical_system_indices);
2093 }
2094
2095 std::vector<double> const& getIntPtMolarFlux(
2096 const double t, std::vector<GlobalVector*> const& x,
2097 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
2098 std::vector<double>& cache, int const component_id) const override
2099 {
2100 std::vector<double> local_x_vec;
2101
2102 auto const n_processes = x.size();
2103 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
2104 {
2105 auto const indices =
2106 NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
2107 assert(!indices.empty());
2108 auto const local_solution = x[process_id]->get(indices);
2109 local_x_vec.insert(std::end(local_x_vec),
2110 std::begin(local_solution),
2111 std::end(local_solution));
2112 }
2113 auto const local_x = MathLib::toVector(local_x_vec);
2114
2115 auto const p = local_x.template segment<pressure_size>(pressure_index);
2116 auto const c = local_x.template segment<concentration_size>(
2117 first_concentration_index + component_id * concentration_size);
2118
2119 auto const n_integration_points =
2120 _integration_method.getNumberOfPoints();
2121
2122 cache.clear();
2123 auto cache_mat = MathLib::createZeroedMatrix<
2124 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
2125 cache, GlobalDim, n_integration_points);
2126
2127 ParameterLib::SpatialPosition pos;
2128 pos.setElementID(_element.getID());
2129
2130 auto const& b =
2131 _process_data
2132 .projected_specific_body_force_vectors[_element.getID()];
2133
2134 MaterialPropertyLib::VariableArray vars;
2135
2136 auto const& medium =
2137 *_process_data.media_map.getMedium(_element.getID());
2138 auto const& phase = medium.phase("AqueousLiquid");
2139
2140 auto const& component = phase.component(
2141 _transport_process_variables[component_id].get().getName());
2142
2143 auto const& Ns =
2144 _process_data.shape_matrix_cache
2145 .NsHigherOrder<typename ShapeFunction::MeshElement>();
2146
2147 for (unsigned ip = 0; ip < n_integration_points; ++ip)
2148 {
2149 auto const& ip_data = _ip_data[ip];
2150 auto const& dNdx = ip_data.dNdx;
2151 auto const& N = Ns[ip];
2152 auto const& phi = ip_data.porosity;
2153
2154 double const p_ip = N.dot(p);
2155 double const c_ip = N.dot(c);
2156
2157 vars.concentration = c_ip;
2158 vars.liquid_phase_pressure = p_ip;
2159 vars.porosity = phi;
2160
2161 double const dt = std::numeric_limits<double>::quiet_NaN();
2162
2163 auto const& k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
2164 medium[MaterialPropertyLib::PropertyType::permeability].value(
2165 vars, pos, t, dt));
2166 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
2167 .template value<double>(vars, pos, t, dt);
2168 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
2169 .template value<double>(vars, pos, t, dt);
2170
2171 // Darcy flux
2172 GlobalDimVectorType const q =
2173 _process_data.has_gravity
2174 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
2175 : GlobalDimVectorType(-k / mu * dNdx * p);
2176
2177 auto const alpha_T = medium.template value<double>(
2178 MaterialPropertyLib::PropertyType::transversal_dispersivity);
2179 auto const alpha_L = medium.template value<double>(
2180 MaterialPropertyLib::PropertyType::longitudinal_dispersivity);
2181 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
2182 component[MaterialPropertyLib::PropertyType::pore_diffusion]
2183 .value(vars, pos, t, dt));
2184
2185 // Hydrodynamic dispersion
2186 GlobalDimMatrixType const D = NumLib::computeHydrodynamicDispersion(
2187 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
2188 alpha_L);
2189
2190 cache_mat.col(ip).noalias() = q * c_ip - D * dNdx * c;
2191 }
2192
2193 return cache;
2194 }
2195
2196 void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
2197 Eigen::VectorXd const& /*local_x_prev*/,
2198 double const /*t*/, double const /*dt*/,
2199 int const /*process_id*/) override
2200 {
2201 unsigned const n_integration_points =
2202 _integration_method.getNumberOfPoints();
2203
2204 for (unsigned ip = 0; ip < n_integration_points; ip++)
2205 {
2206 _ip_data[ip].pushBackState();
2207 }
2208 }
2209
2210private:
2211 MeshLib::Element const& _element;
2212 ComponentTransportProcessData const& _process_data;
2213
2214 NumLib::GenericIntegrationMethod const& _integration_method;
2215 std::vector<std::reference_wrapper<ProcessVariable>> const
2216 _transport_process_variables;
2217
2218 std::vector<IntegrationPointData<GlobalDimNodalMatrixType>> _ip_data;
2219
2220 double getHeatEnergyCoefficient(
2221 MaterialPropertyLib::VariableArray const& vars, const double porosity,
2222 const double fluid_density, const double specific_heat_capacity_fluid,
2223 ParameterLib::SpatialPosition const& pos, double const t,
2224 double const dt)
2225 {
2226 auto const& medium =
2227 *_process_data.media_map.getMedium(this->_element.getID());
2228 auto const& solid_phase = medium.phase("Solid");
2229
2230 auto const specific_heat_capacity_solid =
2231 solid_phase
2232 .property(
2233 MaterialPropertyLib::PropertyType::specific_heat_capacity)
2234 .template value<double>(vars, pos, t, dt);
2235
2236 auto const solid_density =
2237 solid_phase.property(MaterialPropertyLib::PropertyType::density)
2238 .template value<double>(vars, pos, t, dt);
2239
2240 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2241 fluid_density * specific_heat_capacity_fluid * porosity;
2242 }
2243
2244 GlobalDimMatrixType getThermalConductivityDispersivity(
2245 MaterialPropertyLib::VariableArray const& vars,
2246 const double fluid_density, const double specific_heat_capacity_fluid,
2247 const GlobalDimVectorType& velocity,
2248 ParameterLib::SpatialPosition const& pos, double const t,
2249 double const dt)
2250 {
2251 auto const& medium =
2252 *_process_data.media_map.getMedium(_element.getID());
2253
2254 auto thermal_conductivity =
2255 MaterialPropertyLib::formEigenTensor<GlobalDim>(
2256 medium
2257 .property(
2258 MaterialPropertyLib::PropertyType::thermal_conductivity)
2259 .value(vars, pos, t, dt));
2260
2261 auto const thermal_dispersivity_transversal =
2262 medium
2263 .property(MaterialPropertyLib::PropertyType::
2264 thermal_transversal_dispersivity)
2265 .template value<double>();
2266
2267 auto const thermal_dispersivity_longitudinal =
2268 medium
2269 .property(MaterialPropertyLib::PropertyType::
2270 thermal_longitudinal_dispersivity)
2271 .template value<double>();
2272
2273 // Thermal conductivity is moved outside and zero matrix is passed
2274 // instead due to multiplication with fluid's density times specific
2275 // heat capacity.
2276 return thermal_conductivity +
2277 fluid_density * specific_heat_capacity_fluid *
2278 NumLib::computeHydrodynamicDispersion(
2279 _process_data.stabilizer, _element.getID(),
2280 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2281 velocity, 0 /* phi */, thermal_dispersivity_transversal,
2282 thermal_dispersivity_longitudinal);
2283 }
2284
2285 NodalVectorType getLocalTemperature(double const t,
2286 Eigen::VectorXd const& local_x) const
2287 {
2288 NodalVectorType local_T;
2289 if (_process_data.isothermal)
2290 {
2291 if (_process_data.temperature)
2292 {
2293 local_T = _process_data.temperature->getNodalValuesOnElement(
2294 _element, t);
2295 }
2296 else
2297 {
2298 local_T = NodalVectorType::Zero(temperature_size);
2299 }
2300 }
2301 else
2302 {
2303 local_T =
2304 local_x.template segment<temperature_size>(temperature_index);
2305 }
2306 return local_T;
2307 }
2308};
2309
2310} // namespace ComponentTransport
2311} // namespace ProcessLib
ChemicalSolverInterface.h
Interface for coupling OpenGeoSys with an external geochemical solver.
GlobalMatrix
MathLib::EigenMatrix GlobalMatrix
Definition GlobalMatrixVectorTypes.h:17
GlobalVector
MathLib::EigenVector GlobalVector
Definition GlobalMatrixVectorTypes.h:16
GlobalIndexType
GlobalMatrix::IndexType GlobalIndexType
Definition GlobalMatrixVectorTypes.h:44
ShapeMatrixPolicyType
EigenFixedShapeMatrixPolicy< ShapeFunction, GlobalDim > ShapeMatrixPolicyType
Definition ShapeMatrixPolicy.h:134
getName
std::string getName(std::string const &line)
Returns the name/title from the "Zone"-description.
Definition TecPlotTools.cpp:56
MathLib::EigenMatrix::add
int add(IndexType row, IndexType col, double val)
Definition EigenMatrix.h:80
MathLib::EigenVector::add
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:70
MeshLib::Element::getID
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:80
NumLib::LocalToGlobalIndexMap::RowColumnIndices
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
Definition LocalToGlobalIndexMap.h:44
ParameterLib::SpatialPosition::setCoordinates
void setCoordinates(MathLib::Point3d const &coordinates)
Definition SpatialPosition.h:76
ParameterLib::SpatialPosition::setElementID
void setElementID(std::size_t element_id)
Definition SpatialPosition.h:69
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::getIntPtLiquidDensity
virtual std::vector< double > const & getIntPtLiquidDensity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::getIntPtDarcyVelocity
virtual std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::setChemicalSystemConcrete
virtual void setChemicalSystemConcrete(Eigen::VectorXd const &, double const, double const)=0
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::postSpeciationCalculation
virtual void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt)=0
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::computeReactionRelatedSecondaryVariable
virtual void computeReactionRelatedSecondaryVariable(std::size_t const ele_id)=0
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::initializeChemicalSystem
void initializeChemicalSystem(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t)
Definition ComponentTransportFEM.h:65
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::setChemicalSystem
void setChemicalSystem(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
Definition ComponentTransportFEM.h:88
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::assembleReactionEquationConcrete
virtual void assembleReactionEquationConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, int const transport_process_id)=0
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::assembleReactionEquation
void assembleReactionEquation(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, int const process_id)
Definition ComponentTransportFEM.h:111
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::initializeChemicalSystemConcrete
virtual void initializeChemicalSystemConcrete(Eigen::VectorXd const &, double const)=0
ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::getIntPtMolarFlux
virtual std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache, int const component_id) const =0
ProcessLib::ComponentTransport::LocalAssemblerData::getLocalTemperature
NodalVectorType getLocalTemperature(double const t, Eigen::VectorXd const &local_x) const
Definition ComponentTransportFEM.h:2285
ProcessLib::ComponentTransport::LocalAssemblerData::GlobalDimVectorType
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
Definition ComponentTransportFEM.h:234
ProcessLib::ComponentTransport::LocalAssemblerData::assembleForStaggeredScheme
void assembleForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
Definition ComponentTransportFEM.h:819
ProcessLib::ComponentTransport::LocalAssemblerData::ShapeMatrices
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
Definition ComponentTransportFEM.h:219
ProcessLib::ComponentTransport::LocalAssemblerData::getIntPtLiquidDensity
std::vector< double > const & getIntPtLiquidDensity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
Definition ComponentTransportFEM.h:1678
ProcessLib::ComponentTransport::LocalAssemblerData::assembleWithJacobianHydraulicEquation
void assembleWithJacobianHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition ComponentTransportFEM.h:1342
ProcessLib::ComponentTransport::LocalAssemblerData::assembleHeatTransportEquation
void assembleHeatTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &)
Definition ComponentTransportFEM.h:982
ProcessLib::ComponentTransport::LocalAssemblerData::LocalMatrixType
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > LocalMatrixType
Definition ComponentTransportFEM.h:227
ProcessLib::ComponentTransport::LocalAssemblerData::postTimestepConcrete
void postTimestepConcrete(Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
Definition ComponentTransportFEM.h:2196
ProcessLib::ComponentTransport::LocalAssemblerData::getIntPtMolarFlux
std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< double > &cache, int const component_id) const override
Definition ComponentTransportFEM.h:2095
ProcessLib::ComponentTransport::LocalAssemblerData::postSpeciationCalculation
void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt) override
Definition ComponentTransportFEM.h:430
ProcessLib::ComponentTransport::LocalAssemblerData::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition ComponentTransportFEM.h:2214
ProcessLib::ComponentTransport::LocalAssemblerData::calculateIntPtDarcyVelocity
std::vector< double > const & calculateIntPtDarcyVelocity(const double t, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &T_nodal_values, std::vector< double > &cache) const
Definition ComponentTransportFEM.h:1883
ProcessLib::ComponentTransport::LocalAssemblerData::initializeChemicalSystemConcrete
void initializeChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t) override
Definition ComponentTransportFEM.h:312
ProcessLib::ComponentTransport::LocalAssemblerData::computeReactionRelatedSecondaryVariable
void computeReactionRelatedSecondaryVariable(std::size_t const ele_id) override
Definition ComponentTransportFEM.h:2058
ProcessLib::ComponentTransport::LocalAssemblerData::assembleBlockMatrices
void assembleBlockMatrices(GlobalDimVectorType const &b, int const component_id, double const t, double const dt, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< LocalBlockMatrixType > KCC, Eigen::Ref< LocalBlockMatrixType > MCC, Eigen::Ref< LocalBlockMatrixType > MCp, Eigen::Ref< LocalBlockMatrixType > MpC, Eigen::Ref< LocalBlockMatrixType > Kpp, Eigen::Ref< LocalBlockMatrixType > Mpp, Eigen::Ref< LocalSegmentVectorType > Bp)
Definition ComponentTransportFEM.h:564
ProcessLib::ComponentTransport::LocalAssemblerData::LocalBlockMatrixType
typename ShapeMatricesType::template MatrixType< pressure_size, pressure_size > LocalBlockMatrixType
Definition ComponentTransportFEM.h:221
ProcessLib::ComponentTransport::LocalAssemblerData::GlobalDimMatrixType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
Definition ComponentTransportFEM.h:237
ProcessLib::ComponentTransport::LocalAssemblerData::assembleWithJacobianForStaggeredScheme
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
Definition ComponentTransportFEM.h:1322
ProcessLib::ComponentTransport::LocalAssemblerData::_ip_data
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > > _ip_data
Definition ComponentTransportFEM.h:2218
ProcessLib::ComponentTransport::LocalAssemblerData::_process_data
ComponentTransportProcessData const & _process_data
Definition ComponentTransportFEM.h:2212
ProcessLib::ComponentTransport::LocalAssemblerData::getShapeMatrix
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
Definition ComponentTransportFEM.h:1959
ProcessLib::ComponentTransport::LocalAssemblerData::ShapeMatricesType
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
Definition ComponentTransportFEM.h:218
ProcessLib::ComponentTransport::LocalAssemblerData::computeSecondaryVariableConcrete
void computeSecondaryVariableConcrete(double const t, double const, Eigen::VectorXd const &local_x, Eigen::VectorXd const &) override
Definition ComponentTransportFEM.h:2031
ProcessLib::ComponentTransport::LocalAssemblerData::_transport_process_variables
std::vector< std::reference_wrapper< ProcessVariable > > const _transport_process_variables
Definition ComponentTransportFEM.h:2216
ProcessLib::ComponentTransport::LocalAssemblerData::assembleHydraulicEquation
void assembleHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
Definition ComponentTransportFEM.h:847
ProcessLib::ComponentTransport::LocalAssemblerData::GlobalDimNodalMatrixType
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
Definition ComponentTransportFEM.h:235
ProcessLib::ComponentTransport::LocalAssemblerData::NodalRowVectorType
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
Definition ComponentTransportFEM.h:232
ProcessLib::ComponentTransport::LocalAssemblerData::getIntPtDarcyVelocity
std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
Definition ComponentTransportFEM.h:1808
ProcessLib::ComponentTransport::LocalAssemblerData::assembleComponentTransportEquation
void assembleComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &, int const transport_process_id)
Definition ComponentTransportFEM.h:1118
ProcessLib::ComponentTransport::LocalAssemblerData::calculateIntPtLiquidDensity
std::vector< double > const & calculateIntPtLiquidDensity(const double t, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &T_nodal_values, std::vector< double > &cache) const
Definition ComponentTransportFEM.h:1753
ProcessLib::ComponentTransport::LocalAssemblerData::getFlux
Eigen::Vector3d getFlux(MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
Definition ComponentTransportFEM.h:1969
ProcessLib::ComponentTransport::LocalAssemblerData::getThermalConductivityDispersivity
GlobalDimMatrixType getThermalConductivityDispersivity(MaterialPropertyLib::VariableArray const &vars, const double fluid_density, const double specific_heat_capacity_fluid, const GlobalDimVectorType &velocity, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
Definition ComponentTransportFEM.h:2244
ProcessLib::ComponentTransport::LocalAssemblerData::assemble
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
Definition ComponentTransportFEM.h:457
ProcessLib::ComponentTransport::LocalAssemblerData::assembleWithJacobianComponentTransportEquation
void assembleWithJacobianComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, int const component_id)
Definition ComponentTransportFEM.h:1450
ProcessLib::ComponentTransport::LocalAssemblerData::LocalVectorType
Eigen::Matrix< double, Eigen::Dynamic, 1 > LocalVectorType
Definition ComponentTransportFEM.h:229
ProcessLib::ComponentTransport::LocalAssemblerData::getHeatEnergyCoefficient
double getHeatEnergyCoefficient(MaterialPropertyLib::VariableArray const &vars, const double porosity, const double fluid_density, const double specific_heat_capacity_fluid, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
Definition ComponentTransportFEM.h:2220
ProcessLib::ComponentTransport::LocalAssemblerData::assembleReactionEquationConcrete
void assembleReactionEquationConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, int const transport_process_id) override
Definition ComponentTransportFEM.h:1598
ProcessLib::ComponentTransport::LocalAssemblerData::LocalAssemblerData
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool is_axially_symmetric, ComponentTransportProcessData const &process_data, std::vector< std::reference_wrapper< ProcessVariable > > const &transport_process_variables)
Definition ComponentTransportFEM.h:240
ProcessLib::ComponentTransport::LocalAssemblerData::LocalSegmentVectorType
typename ShapeMatricesType::template VectorType< pressure_size > LocalSegmentVectorType
Definition ComponentTransportFEM.h:224
ProcessLib::ComponentTransport::LocalAssemblerData::assembleKCmCn
void assembleKCmCn(int const component_id, double const t, double const dt, Eigen::Ref< LocalBlockMatrixType > KCmCn, double const stoichiometric_coefficient, double const kinetic_prefactor)
Definition ComponentTransportFEM.h:765
ProcessLib::ComponentTransport::LocalAssemblerData::setChemicalSystemConcrete
void setChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t, double dt) override
Definition ComponentTransportFEM.h:365
ProcessLib::ComponentTransport::LocalAssemblerData::NodalVectorType
typename ShapeMatricesType::NodalVectorType NodalVectorType
Definition ComponentTransportFEM.h:231
MaterialPropertyLib::formEigenTensor
constexpr Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition FormEigenTensor.h:163
MaterialPropertyLib::longitudinal_dispersivity
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
Definition PropertyType.h:50
MaterialPropertyLib::transversal_dispersivity
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
Definition PropertyType.h:101
MaterialPropertyLib::retardation_factor
@ retardation_factor
specify retardation factor used in component transport process.
Definition PropertyType.h:74
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:146
MathLib::toVector
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Definition EigenMapTools.h:160
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:25
MathLib::toMatrix
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:65
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:20
NumLib::detail::assembleAdvectionMatrix
void assembleAdvectionMatrix(IPData const &ip_data_vector, NumLib::ShapeMatrixCache const &shape_matrix_cache, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
Definition AdvectionMatrixAssembler.h:22
NumLib::getIndices
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
Definition DOFTableUtil.cpp:105
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:48
NumLib::computeShapeMatrices
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
Definition InitShapeMatrices.h:20
NumLib::computeHydrodynamicDispersion
Eigen::MatrixXd computeHydrodynamicDispersion(NumericalStabilization const &stabilizer, std::size_t const element_id, Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
Definition HydrodynamicDispersion.h:65
NumLib::interpolateCoordinates
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:93
EigenFixedShapeMatrixPolicy::ShapeMatrices
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
Definition ShapeMatrixPolicy.h:80
EigenFixedShapeMatrixPolicy::GlobalDimNodalMatrixType
MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType
Definition ShapeMatrixPolicy.h:75
EigenFixedShapeMatrixPolicy::GlobalDimMatrixType
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
Definition ShapeMatrixPolicy.h:77
EigenFixedShapeMatrixPolicy::GlobalDimVectorType
VectorType< GlobalDim > GlobalDimVectorType
Definition ShapeMatrixPolicy.h:78
EigenFixedShapeMatrixPolicy::NodalVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
Definition ShapeMatrixPolicy.h:69
EigenFixedShapeMatrixPolicy::NodalRowVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
Definition ShapeMatrixPolicy.h:71
ProcessLib::ComponentTransport::IntegrationPointData::IntegrationPointData
IntegrationPointData(GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_)
Definition ComponentTransportFEM.h:37