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