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