OGS
LiquidFlowLocalAssembler-impl.h
Go to the documentation of this file.
1
13#pragma once
14
21
22namespace ProcessLib
23{
24namespace LiquidFlow
25{
26template <typename ShapeFunction, int GlobalDim>
28 double const t, double const dt, std::vector<double> const& local_x,
29 std::vector<double> const& /*local_x_prev*/,
30 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
31 std::vector<double>& local_b_data)
32{
34 pos.setElementID(_element.getID());
35
36 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
38 vars.temperature =
40 .template value<double>(vars, pos, t, dt);
41 vars.liquid_phase_pressure = std::numeric_limits<double>::quiet_NaN();
42 GlobalDimMatrixType const permeability =
45 vars, pos, t, dt));
46 // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
47 // the assert must be changed to permeability.rows() ==
48 // _element->getDimension()
49 assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
50
51 if (permeability.size() == 1)
52 { // isotropic or 1D problem.
53 assembleMatrixAndVector<IsotropicCalculator>(
54 t, dt, local_x, local_M_data, local_K_data, local_b_data);
55 }
56 else
57 {
58 assembleMatrixAndVector<AnisotropicCalculator>(
59 t, dt, local_x, local_M_data, local_K_data, local_b_data);
60 }
61}
62
63template <typename ShapeFunction, int GlobalDim>
65 MathLib::Point3d const& p_local_coords, double const t,
66 std::vector<double> const& local_x) const
67{
68 // TODO (tf) Temporary value not used by current material models. Need
69 // extension of getFlux interface
70 double const dt = std::numeric_limits<double>::quiet_NaN();
71
72 // Note: Axial symmetry is set to false here, because we only need dNdx
73 // here, which is not affected by axial symmetry.
74 auto const shape_matrices =
76 GlobalDim>(_element,
77 false /*is_axially_symmetric*/,
78 std::array{p_local_coords})[0];
79
80 // create pos object to access the correct media property
82 pos.setElementID(_element.getID());
83
84 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
85 auto const& fluid_phase = fluidPhase(medium);
86
88
89 double pressure = 0.0;
90 NumLib::shapeFunctionInterpolate(local_x, shape_matrices.N, pressure);
91 vars.liquid_phase_pressure = pressure;
92
93 GlobalDimMatrixType const intrinsic_permeability =
96 vars, pos, t, dt));
97 auto const viscosity =
99 .template value<double>(vars, pos, t, dt);
100
101 Eigen::Vector3d flux(0.0, 0.0, 0.0);
102 flux.head<GlobalDim>() =
103 -intrinsic_permeability / viscosity * shape_matrices.dNdx *
104 Eigen::Map<const NodalVectorType>(local_x.data(), local_x.size());
105
106 return flux;
107}
108
109template <typename ShapeFunction, int GlobalDim>
110template <typename LaplacianGravityVelocityCalculator>
112 assembleMatrixAndVector(double const t, double const dt,
113 std::vector<double> const& local_x,
114 std::vector<double>& local_M_data,
115 std::vector<double>& local_K_data,
116 std::vector<double>& local_b_data)
117{
118 auto const local_matrix_size = local_x.size();
119 assert(local_matrix_size == ShapeFunction::NPOINTS);
120
122 local_M_data, local_matrix_size, local_matrix_size);
124 local_K_data, local_matrix_size, local_matrix_size);
126 local_b_data, local_matrix_size);
127 const auto local_p_vec =
128 MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
129
130 unsigned const n_integration_points =
131 _integration_method.getNumberOfPoints();
132
134 pos.setElementID(_element.getID());
135
136 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
137 auto const& fluid_phase = fluidPhase(medium);
138
140 auto& phase_pressure = medium.hasPhase("Gas") ? vars.gas_phase_pressure
142
143 vars.temperature =
145 .template value<double>(vars, pos, t, dt);
146
147 GlobalDimVectorType const projected_body_force_vector =
148 _process_data.element_rotation_matrices[_element.getID()] *
149 _process_data.element_rotation_matrices[_element.getID()].transpose() *
150 _process_data.specific_body_force;
151
152 auto const& Ns = _process_data.shape_matrix_cache
153 .NsHigherOrder<typename ShapeFunction::MeshElement>();
154
155 for (unsigned ip = 0; ip < n_integration_points; ip++)
156 {
157 auto const& ip_data = _ip_data[ip];
158 auto const& N = Ns[ip];
159
160 phase_pressure = N.dot(local_p_vec);
161
162 auto const [fluid_density, viscosity] =
164 fluid_phase, vars);
165
166 auto const porosity =
168 .template value<double>(vars, pos, t, dt);
169 auto const specific_storage =
171 .template value<double>(vars, pos, t, dt);
172
173 auto const get_drho_dp = [&]()
174 {
176 .template dValue<double>(vars, _process_data.phase_variable,
177 pos, t, dt);
178 };
179 auto const storage =
180 _process_data.equation_balance_type == EquationBalanceType::volume
181 ? specific_storage
182 : specific_storage + porosity * get_drho_dp() / fluid_density;
183
184 double const scaling_factor =
185 _process_data.equation_balance_type == EquationBalanceType::volume
186 ? 1.0
187 : fluid_density;
188 // Assemble mass matrix, M
189 local_M.noalias() += scaling_factor * storage * N.transpose() * N *
190 ip_data.integration_weight;
191
192 GlobalDimMatrixType const permeability =
195 vars, pos, t, dt));
196
197 // Assemble Laplacian, K, and RHS by the gravitational term
198 LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
199 local_K, local_b, ip_data, scaling_factor * permeability, viscosity,
200 fluid_density, projected_body_force_vector,
201 _process_data.has_gravity);
202 }
203}
204
205template <typename ShapeFunction, int GlobalDim>
206template <typename VelocityCacheType>
208 bool const is_scalar_permeability, const double t, const double dt,
209 std::vector<double> const& local_x,
211 VelocityCacheType& darcy_velocity_at_ips) const
212{
213 if (is_scalar_permeability)
214 { // isotropic or 1D problem.
215 computeProjectedDarcyVelocity<IsotropicCalculator>(
216 t, dt, local_x, pos, darcy_velocity_at_ips);
217 }
218 else
219 {
220 computeProjectedDarcyVelocity<AnisotropicCalculator>(
221 t, dt, local_x, pos, darcy_velocity_at_ips);
222 }
223}
224
225template <typename ShapeFunction, int GlobalDim>
226std::vector<double> const&
228 const double t,
229 std::vector<GlobalVector*> const& x,
230 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
231 std::vector<double>& velocity_cache) const
232{
233 // TODO (tf) Temporary value not used by current material models. Need
234 // extension of secondary variable interface.
235 double const dt = std::numeric_limits<double>::quiet_NaN();
236
237 constexpr int process_id = 0;
238 auto const indices =
239 NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
240 auto const local_x = x[process_id]->get(indices);
241 auto const n_integration_points = _integration_method.getNumberOfPoints();
242 velocity_cache.clear();
243
245 pos.setElementID(_element.getID());
246
247 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
249 vars.temperature =
251 .template value<double>(vars, pos, t, dt);
252 vars.liquid_phase_pressure = std::numeric_limits<double>::quiet_NaN();
253
254 GlobalDimMatrixType const permeability =
257 vars, pos, t, dt));
258
259 assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
260
261 bool const is_scalar_permeability = (permeability.size() == 1);
262
263 auto velocity_cache_vectors = MathLib::createZeroedMatrix<
264 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
265 velocity_cache, GlobalDim, n_integration_points);
266
267 computeDarcyVelocity(is_scalar_permeability, t, dt, local_x, pos,
268 velocity_cache_vectors);
269
270 return velocity_cache;
271}
272
273template <typename ShapeFunction, int GlobalDim>
274template <typename LaplacianGravityVelocityCalculator,
275 typename VelocityCacheType>
278 const double t, const double dt, std::vector<double> const& local_x,
280 VelocityCacheType& darcy_velocity_at_ips) const
281{
282 auto const local_matrix_size = local_x.size();
283 assert(local_matrix_size == ShapeFunction::NPOINTS);
284
285 const auto local_p_vec =
286 MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
287
288 unsigned const n_integration_points =
289 _integration_method.getNumberOfPoints();
290
291 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
292 auto const& fluid_phase = fluidPhase(medium);
293
295 auto& phase_pressure = medium.hasPhase("Gas") ? vars.gas_phase_pressure
297
298 vars.temperature =
300 .template value<double>(vars, pos, t, dt);
301
302 GlobalDimVectorType const projected_body_force_vector =
303 _process_data.element_rotation_matrices[_element.getID()] *
304 _process_data.element_rotation_matrices[_element.getID()].transpose() *
305 _process_data.specific_body_force;
306
307 auto const& Ns = _process_data.shape_matrix_cache
308 .NsHigherOrder<typename ShapeFunction::MeshElement>();
309
310 for (unsigned ip = 0; ip < n_integration_points; ip++)
311 {
312 auto const& ip_data = _ip_data[ip];
313 auto const& N = Ns[ip];
314
315 phase_pressure = N.dot(local_p_vec);
316
317 auto const [fluid_density, viscosity] =
319 fluid_phase, vars);
320
321 GlobalDimMatrixType const permeability =
324 vars, pos, t, dt));
325
326 darcy_velocity_at_ips.col(ip) =
327 LaplacianGravityVelocityCalculator::calculateVelocity(
328 local_p_vec, ip_data, permeability, viscosity, fluid_density,
329 projected_body_force_vector, _process_data.has_gravity);
330 }
331}
332
333template <typename ShapeFunction, int GlobalDim>
336 Eigen::Map<NodalMatrixType>& local_K,
337 Eigen::Map<NodalVectorType>& local_b,
339 GlobalDimMatrixType const& permeability_with_density_factor,
340 double const mu, double const rho_L,
341 GlobalDimVectorType const& specific_body_force, bool const has_gravity)
342{
343 const double K = permeability_with_density_factor(0, 0) / mu;
344 const double fac = K * ip_data.integration_weight;
345 local_K.noalias() += fac * ip_data.dNdx.transpose() * ip_data.dNdx;
346
347 if (has_gravity)
348 {
349 local_b.noalias() +=
350 (fac * rho_L) * ip_data.dNdx.transpose() * specific_body_force;
351 }
352}
353
354template <typename ShapeFunction, int GlobalDim>
355Eigen::Matrix<double, GlobalDim, 1>
358 Eigen::Map<const NodalVectorType> const& local_p,
360 GlobalDimMatrixType const& permeability, double const mu,
361 double const rho_L, GlobalDimVectorType const& specific_body_force,
362 bool const has_gravity)
363{
364 const double K = permeability(0, 0) / mu;
365 // Compute the velocity
366 GlobalDimVectorType velocity = -K * ip_data.dNdx * local_p;
367 // gravity term
368 if (has_gravity)
369 {
370 velocity += (K * rho_L) * specific_body_force;
371 }
372 return velocity;
373}
374
375template <typename ShapeFunction, int GlobalDim>
378 Eigen::Map<NodalMatrixType>& local_K,
379 Eigen::Map<NodalVectorType>& local_b,
381 GlobalDimMatrixType const& permeability_with_density_factor,
382 double const mu, double const rho_L,
383 GlobalDimVectorType const& specific_body_force, bool const has_gravity)
384{
385 const double fac = ip_data.integration_weight / mu;
386 local_K.noalias() += fac * ip_data.dNdx.transpose() *
387 permeability_with_density_factor * ip_data.dNdx;
388
389 if (has_gravity)
390 {
391 local_b.noalias() += (fac * rho_L) * ip_data.dNdx.transpose() *
392 permeability_with_density_factor *
393 specific_body_force;
394 }
395}
396
397template <typename ShapeFunction, int GlobalDim>
398Eigen::Matrix<double, GlobalDim, 1>
401 Eigen::Map<const NodalVectorType> const& local_p,
403 GlobalDimMatrixType const& permeability, double const mu,
404 double const rho_L, GlobalDimVectorType const& specific_body_force,
405 bool const has_gravity)
406{
407 // Compute the velocity
408 GlobalDimVectorType velocity = -permeability * ip_data.dNdx * local_p / mu;
409
410 // gravity term
411 if (has_gravity)
412 {
413 velocity += (rho_L / mu) * permeability * specific_body_force;
414 }
415 return velocity;
416}
417
418} // namespace LiquidFlow
419} // namespace ProcessLib
void setElementID(std::size_t element_id)
Eigen::Vector3d getFlux(MathLib::Point3d const &p_local_coords, double const t, std::vector< double > const &local_x) const override
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
void computeProjectedDarcyVelocity(const double t, const double dt, std::vector< double > const &local_x, ParameterLib::SpatialPosition const &pos, VelocityCacheType &darcy_velocity_at_ips) const
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
void assembleMatrixAndVector(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
void computeDarcyVelocity(bool const is_scalar_permeability, const double t, const double dt, std::vector< double > const &local_x, ParameterLib::SpatialPosition const &pos, VelocityCacheType &darcy_velocity_at_ips) const
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &velocity_cache) const override
std::tuple< double, double > getFluidDensityAndViscosity(double const t, double const dt, ParameterLib::SpatialPosition const &pos, MaterialPropertyLib::Phase const &fluid_phase, MaterialPropertyLib::VariableArray &vars)
It computes fluid density and viscosity for single phase flow model.
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
static void calculateLaplacianAndGravityTerm(Eigen::Map< NodalMatrixType > &local_K, Eigen::Map< NodalVectorType > &local_b, IntegrationPointData< GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability_with_density_factor, double const mu, double const rho_L, GlobalDimVectorType const &specific_body_force, bool const has_gravity)
static Eigen::Matrix< double, GlobalDim, 1 > calculateVelocity(Eigen::Map< const NodalVectorType > const &local_p, IntegrationPointData< GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability, double const mu, double const rho_L, GlobalDimVectorType const &specific_body_force, bool const has_gravity)
static void calculateLaplacianAndGravityTerm(Eigen::Map< NodalMatrixType > &local_K, Eigen::Map< NodalVectorType > &local_b, IntegrationPointData< GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability_with_density_factor, double const mu, double const rho_L, GlobalDimVectorType const &specific_body_force, bool const has_gravity)
static Eigen::Matrix< double, GlobalDim, 1 > calculateVelocity(Eigen::Map< const NodalVectorType > const &local_p, IntegrationPointData< GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability, double const mu, double const rho_L, GlobalDimVectorType const &specific_body_force, bool const has_gravity)