OGS
LiquidFlowLocalAssembler-impl.h
Go to the documentation of this file.
1
13#pragma once
14
20
21namespace ProcessLib
22{
23namespace LiquidFlow
24{
25template <typename ShapeFunction, int GlobalDim>
27 double const t, double const dt, std::vector<double> const& local_x,
28 std::vector<double> const& /*local_x_prev*/,
29 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
30 std::vector<double>& local_b_data)
31{
33 pos.setElementID(_element.getID());
34
35 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
37 vars.temperature =
39 .template value<double>(vars, pos, t, dt);
40 vars.liquid_phase_pressure = std::numeric_limits<double>::quiet_NaN();
41 GlobalDimMatrixType const permeability =
42 MaterialPropertyLib::formEigenTensor<GlobalDim>(
44 vars, pos, t, dt));
45 // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
46 // the assert must be changed to permeability.rows() ==
47 // _element->getDimension()
48 assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
49
50 if (permeability.size() == 1)
51 { // isotropic or 1D problem.
52 assembleMatrixAndVector<IsotropicCalculator>(
53 t, dt, local_x, local_M_data, local_K_data, local_b_data);
54 }
55 else
56 {
57 assembleMatrixAndVector<AnisotropicCalculator>(
58 t, dt, local_x, local_M_data, local_K_data, local_b_data);
59 }
60}
61
62template <typename ShapeFunction, int GlobalDim>
64 MathLib::Point3d const& p_local_coords, double const t,
65 std::vector<double> const& local_x) const
66{
67 // TODO (tf) Temporary value not used by current material models. Need
68 // extension of getFlux interface
69 double const dt = std::numeric_limits<double>::quiet_NaN();
70
71 // Note: Axial symmetry is set to false here, because we only need dNdx
72 // here, which is not affected by axial symmetry.
73 auto const shape_matrices =
75 GlobalDim>(_element,
76 false /*is_axially_symmetric*/,
77 std::array{p_local_coords})[0];
78
79 // create pos object to access the correct media property
81 pos.setElementID(_element.getID());
82
83 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
84 auto const& liquid_phase = medium.phase("AqueousLiquid");
85
87
88 double pressure = 0.0;
89 NumLib::shapeFunctionInterpolate(local_x, shape_matrices.N, pressure);
90 vars.liquid_phase_pressure = pressure;
91
92 GlobalDimMatrixType const intrinsic_permeability =
93 MaterialPropertyLib::formEigenTensor<GlobalDim>(
95 vars, pos, t, dt));
96 auto const viscosity =
98 .template value<double>(vars, pos, t, dt);
99
100 Eigen::Vector3d flux(0.0, 0.0, 0.0);
101 flux.head<GlobalDim>() =
102 -intrinsic_permeability / viscosity * shape_matrices.dNdx *
103 Eigen::Map<const NodalVectorType>(local_x.data(), local_x.size());
104
105 return flux;
106}
107
108template <typename ShapeFunction, int GlobalDim>
109template <typename LaplacianGravityVelocityCalculator>
111 assembleMatrixAndVector(double const t, double const dt,
112 std::vector<double> const& local_x,
113 std::vector<double>& local_M_data,
114 std::vector<double>& local_K_data,
115 std::vector<double>& local_b_data)
116{
117 auto const local_matrix_size = local_x.size();
118 assert(local_matrix_size == ShapeFunction::NPOINTS);
119
120 auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
121 local_M_data, local_matrix_size, local_matrix_size);
122 auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
123 local_K_data, local_matrix_size, local_matrix_size);
124 auto local_b = MathLib::createZeroedVector<NodalVectorType>(
125 local_b_data, local_matrix_size);
126
127 unsigned const n_integration_points =
128 _integration_method.getNumberOfPoints();
129
131 pos.setElementID(_element.getID());
132
133 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
134 auto const& liquid_phase = medium.phase("AqueousLiquid");
135
137 vars.temperature =
139 .template value<double>(vars, pos, t, dt);
140
141 GlobalDimVectorType const projected_body_force_vector =
142 _process_data.element_rotation_matrices[_element.getID()] *
143 _process_data.element_rotation_matrices[_element.getID()].transpose() *
144 _process_data.specific_body_force;
145
146 auto const& N = _shape_matrix_cache
147 .NsHigherOrder<typename ShapeFunction::MeshElement>();
148
149 for (unsigned ip = 0; ip < n_integration_points; ip++)
150 {
151 auto const& ip_data = _ip_data[ip];
152
153 double p = 0.;
154 NumLib::shapeFunctionInterpolate(local_x, N[ip], p);
155 vars.liquid_phase_pressure = p;
156
157 // Compute density:
158 auto const fluid_density =
160 .template value<double>(vars, pos, t, dt);
161 assert(fluid_density > 0.);
162 vars.density = fluid_density;
163
164 auto const ddensity_dpressure =
166 .template dValue<double>(
168 pos, t, dt);
169
170 auto const porosity =
172 .template value<double>(vars, pos, t, dt);
173 auto const storage = medium[MaterialPropertyLib::PropertyType::storage]
174 .template value<double>(vars, pos, t, dt);
175
176 // Assemble mass matrix, M
177 local_M.noalias() +=
178 (porosity * ddensity_dpressure / fluid_density + storage) *
179 N[ip].transpose() * N[ip] * ip_data.integration_weight;
180
181 // Compute viscosity:
182 auto const viscosity =
184 .template value<double>(vars, pos, t, dt);
185
186 pos.setIntegrationPoint(ip);
187 GlobalDimMatrixType const permeability =
188 MaterialPropertyLib::formEigenTensor<GlobalDim>(
190 vars, pos, t, dt));
191
192 // Assemble Laplacian, K, and RHS by the gravitational term
193 LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
194 local_K, local_b, ip_data, permeability, viscosity, fluid_density,
195 projected_body_force_vector, _process_data.has_gravity);
196 }
197}
198
199template <typename ShapeFunction, int GlobalDim>
200template <typename VelocityCacheType>
202 bool const is_scalar_permeability, const double t, const double dt,
203 std::vector<double> const& local_x,
205 VelocityCacheType& darcy_velocity_at_ips) const
206{
207 if (is_scalar_permeability)
208 { // isotropic or 1D problem.
209 computeProjectedDarcyVelocity<IsotropicCalculator>(
210 t, dt, local_x, pos, darcy_velocity_at_ips);
211 }
212 else
213 {
214 computeProjectedDarcyVelocity<AnisotropicCalculator>(
215 t, dt, local_x, pos, darcy_velocity_at_ips);
216 }
217}
218
219template <typename ShapeFunction, int GlobalDim>
220std::vector<double> const&
222 const double t,
223 std::vector<GlobalVector*> const& x,
224 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
225 std::vector<double>& velocity_cache) const
226{
227 // TODO (tf) Temporary value not used by current material models. Need
228 // extension of secondary variable interface.
229 double const dt = std::numeric_limits<double>::quiet_NaN();
230
231 constexpr int process_id = 0;
232 auto const indices =
233 NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
234 auto const local_x = x[process_id]->get(indices);
235 auto const n_integration_points = _integration_method.getNumberOfPoints();
236 velocity_cache.clear();
237
239 pos.setElementID(_element.getID());
240
241 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
243 vars.temperature =
245 .template value<double>(vars, pos, t, dt);
246 vars.liquid_phase_pressure = std::numeric_limits<double>::quiet_NaN();
247
248 GlobalDimMatrixType const permeability =
249 MaterialPropertyLib::formEigenTensor<GlobalDim>(
251 vars, pos, t, dt));
252
253 assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
254
255 bool const is_scalar_permeability = (permeability.size() == 1);
256
257 auto velocity_cache_vectors = MathLib::createZeroedMatrix<
258 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
259 velocity_cache, GlobalDim, n_integration_points);
260
261 computeDarcyVelocity(is_scalar_permeability, t, dt, local_x, pos,
262 velocity_cache_vectors);
263
264 return velocity_cache;
265}
266
267template <typename ShapeFunction, int GlobalDim>
268template <typename LaplacianGravityVelocityCalculator,
269 typename VelocityCacheType>
272 const double t, const double dt, std::vector<double> const& local_x,
274 VelocityCacheType& darcy_velocity_at_ips) const
275{
276 auto const local_matrix_size = local_x.size();
277 assert(local_matrix_size == ShapeFunction::NPOINTS);
278
279 const auto local_p_vec =
280 MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
281
282 unsigned const n_integration_points =
283 _integration_method.getNumberOfPoints();
284
285 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
286 auto const& liquid_phase = medium.phase("AqueousLiquid");
287
289 vars.temperature =
291 .template value<double>(vars, pos, t, dt);
292
293 GlobalDimVectorType const projected_body_force_vector =
294 _process_data.element_rotation_matrices[_element.getID()] *
295 _process_data.element_rotation_matrices[_element.getID()].transpose() *
296 _process_data.specific_body_force;
297
298 auto const& N = _shape_matrix_cache
299 .NsHigherOrder<typename ShapeFunction::MeshElement>();
300
301 for (unsigned ip = 0; ip < n_integration_points; ip++)
302 {
303 auto const& ip_data = _ip_data[ip];
304 double p = 0.;
305 NumLib::shapeFunctionInterpolate(local_x, N[ip], p);
306 vars.liquid_phase_pressure = p;
307
308 // Compute density:
309 auto const fluid_density =
311 .template value<double>(vars, pos, t, dt);
312 vars.density = fluid_density;
313
314 // Compute viscosity:
315 auto const viscosity =
317 .template value<double>(vars, pos, t, dt);
318
319 GlobalDimMatrixType const permeability =
320 MaterialPropertyLib::formEigenTensor<GlobalDim>(
322 vars, pos, t, dt));
323
324 darcy_velocity_at_ips.col(ip) =
325 LaplacianGravityVelocityCalculator::calculateVelocity(
326 local_p_vec, ip_data, permeability, viscosity, fluid_density,
327 projected_body_force_vector, _process_data.has_gravity);
328 }
329}
330
331template <typename ShapeFunction, int GlobalDim>
334 Eigen::Map<NodalMatrixType>& local_K,
335 Eigen::Map<NodalVectorType>& local_b,
337 GlobalDimNodalMatrixType> const& ip_data,
338 GlobalDimMatrixType const& permeability, double const mu,
339 double const rho_L, GlobalDimVectorType const& specific_body_force,
340 bool const has_gravity)
341{
342 const double K = permeability(0, 0) / mu;
343 const double fac = K * ip_data.integration_weight;
344 local_K.noalias() += fac * ip_data.dNdx.transpose() * ip_data.dNdx;
345
346 if (has_gravity)
347 {
348 local_b.noalias() +=
349 (fac * rho_L) * ip_data.dNdx.transpose() * specific_body_force;
350 }
351}
352
353template <typename ShapeFunction, int GlobalDim>
354Eigen::Matrix<double, GlobalDim, 1>
357 Eigen::Map<const NodalVectorType> const& local_p,
359 GlobalDimNodalMatrixType> const& ip_data,
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 GlobalDimNodalMatrixType> const& ip_data,
382 GlobalDimMatrixType const& permeability, double const mu,
383 double const rho_L, GlobalDimVectorType const& specific_body_force,
384 bool const has_gravity)
385{
386 const double fac = ip_data.integration_weight / mu;
387 local_K.noalias() +=
388 fac * ip_data.dNdx.transpose() * permeability * ip_data.dNdx;
389
390 if (has_gravity)
391 {
392 local_b.noalias() += (fac * rho_L) * ip_data.dNdx.transpose() *
393 permeability * 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 GlobalDimNodalMatrixType> const& ip_data,
404 GlobalDimMatrixType const& permeability, double const mu,
405 double const rho_L, GlobalDimVectorType const& specific_body_force,
406 bool const has_gravity)
407{
408 // Compute the velocity
409 GlobalDimVectorType velocity = -permeability * ip_data.dNdx * local_p / mu;
410
411 // gravity term
412 if (has_gravity)
413 {
414 velocity += (rho_L / mu) * permeability * specific_body_force;
415 }
416 return velocity;
417}
418
419} // namespace LiquidFlow
420} // namespace ProcessLib
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
Eigen::Vector3d getFlux(MathLib::Point3d const &p_local_coords, double const t, std::vector< double > const &local_x) const override
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
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
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
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
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< NodalRowVectorType, GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability, 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< NodalRowVectorType, 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< NodalRowVectorType, GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability, 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< NodalRowVectorType, GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability, double const mu, double const rho_L, GlobalDimVectorType const &specific_body_force, bool const has_gravity)