OGS 6.1.0-1721-g6382411ad
HTProcess.cpp
Go to the documentation of this file.
1 
10 #include "HTProcess.h"
11 
12 #include <cassert>
13 
18 
19 #include "HTMaterialProperties.h"
20 #include "MonolithicHTFEM.h"
21 #include "StaggeredHTFEM.h"
22 
23 namespace ProcessLib
24 {
25 namespace HT
26 {
28  MeshLib::Mesh& mesh,
29  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
30  std::vector<std::unique_ptr<ParameterBase>> const& parameters,
31  unsigned const integration_order,
32  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
33  process_variables,
34  std::unique_ptr<HTMaterialProperties>&& material_properties,
35  SecondaryVariableCollection&& secondary_variables,
36  NumLib::NamedFunctionCaller&& named_function_caller,
37  bool const use_monolithic_scheme,
38  std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux)
39  : Process(mesh, std::move(jacobian_assembler), parameters,
40  integration_order, std::move(process_variables),
41  std::move(secondary_variables), std::move(named_function_caller),
42  use_monolithic_scheme),
43  _material_properties(std::move(material_properties)),
44  _surfaceflux(std::move(surfaceflux))
45 {
46 }
47 
49  NumLib::LocalToGlobalIndexMap const& dof_table,
50  MeshLib::Mesh const& mesh,
51  unsigned const integration_order)
52 {
53  // For the staggered scheme, both processes are assumed to use the same
54  // element order. Therefore the order of shape function can be fetched from
55  // any set of the sets of process variables of the coupled processes. Here,
56  // we take the one from the first process by setting process_id = 0.
57  const int process_id = 0;
58  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
59 
61  {
62  ProcessLib::createLocalAssemblers<MonolithicHTFEM>(
63  mesh.getDimension(), mesh.getElements(), dof_table,
65  mesh.isAxiallySymmetric(), integration_order,
67  }
68  else
69  {
70  const int heat_transport_process_id = 0;
71  const int hydraulic_process_id = 1;
72 
73  ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
74  mesh.getDimension(), mesh.getElements(), dof_table,
76  mesh.isAxiallySymmetric(), integration_order, *_material_properties,
77  heat_transport_process_id, hydraulic_process_id);
78  }
79 
81  "darcy_velocity",
85 }
86 
88  GlobalVector const& x,
89  GlobalMatrix& M,
90  GlobalMatrix& K,
91  GlobalVector& b)
92 {
93  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
94  dof_tables;
96  {
97  DBUG("Assemble HTProcess.");
98  dof_tables.emplace_back(*_local_to_global_index_map);
99  }
100  else
101  {
102  if (_coupled_solutions->process_id == 0)
103  {
104  DBUG(
105  "Assemble the equations of heat transport process within "
106  "HTProcess.");
107  }
108  else
109  {
110  DBUG(
111  "Assemble the equations of single phase fully saturated "
112  "fluid flow process within HTProcess.");
113  }
115  dof_tables.emplace_back(*_local_to_global_index_map);
116  dof_tables.emplace_back(*_local_to_global_index_map);
117  }
118 
119  const int process_id =
121  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
122  // Call global assembler for each local assembly item.
125  pv.getActiveElementIDs(), dof_tables, t, x, M, K, b,
127 }
128 
130  const double t, GlobalVector const& x, GlobalVector const& xdot,
131  const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
132  GlobalVector& b, GlobalMatrix& Jac)
133 {
134  DBUG("AssembleWithJacobian HTProcess.");
135 
136  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
137  dof_tables;
139  {
141  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
142  }
143  else
144  {
145  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
146  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
147  }
148 
149  // Call global assembler for each local assembly item.
150  const int process_id =
152  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
155  _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, x,
156  xdot, dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
157 }
158 
159 void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
160  const double /*t*/,
161  const double /*delta_t*/,
162  const int process_id)
163 {
164  assert(process_id < 2);
165 
167  {
168  return;
169  }
170 
171  if (!_xs_previous_timestep[process_id])
172  {
173  _xs_previous_timestep[process_id] =
175  }
176  else
177  {
178  auto& x0 = *_xs_previous_timestep[process_id];
179  MathLib::LinAlg::copy(x, x0);
180  }
181 
182  auto& x0 = *_xs_previous_timestep[process_id];
183  MathLib::LinAlg::setLocalAccessibleVector(x0);
184 }
185 
187 {
188  DBUG("Set the coupled term for the staggered scheme to local assembers.");
189 
190  const int process_id =
192  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
197 }
198 
199 std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
201 {
203  {
204  // For single-variable-single-component processes reuse the existing DOF
205  // table.
206  const bool manage_storage = false;
207  return std::make_tuple(_local_to_global_index_map.get(),
208  manage_storage);
209  }
210 
211  // Otherwise construct a new DOF table.
212  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
214 
215  const bool manage_storage = true;
216  return std::make_tuple(new NumLib::LocalToGlobalIndexMap(
217  std::move(all_mesh_subsets_single_component),
218  // by location order is needed for output
220  manage_storage);
221 }
222 
224 {
225  const auto number_of_coupled_solutions =
228  _coupled_solutions->coupled_xs_t0.reserve(number_of_coupled_solutions);
229  const int process_id = _coupled_solutions->process_id;
230  for (std::size_t i = 0; i < number_of_coupled_solutions; i++)
231  {
232  const auto& x_t0 = _xs_previous_timestep[process_id];
233  if (x_t0 == nullptr)
234  {
235  OGS_FATAL(
236  "Memory is not allocated for the global vector "
237  "of the solution of the previous time step for the ."
238  "staggered scheme.\n It can be done by overriding "
239  "Process::preTimestepConcreteProcess"
240  "(ref. HTProcess::preTimestepConcreteProcess) ");
241  }
242 
243  MathLib::LinAlg::setLocalAccessibleVector(*x_t0);
244  _coupled_solutions->coupled_xs_t0.emplace_back(x_t0.get());
245  }
246 }
247 
248 Eigen::Vector3d HTProcess::getFlux(std::size_t element_id,
249  MathLib::Point3d const& p,
250  double const t,
251  GlobalVector const& x) const
252 {
253  // fetch local_x from primary variable
254  std::vector<GlobalIndexType> indices_cache;
255  auto const r_c_indices = NumLib::getRowColumnIndices(
256  element_id, *_local_to_global_index_map, indices_cache);
257  std::vector<double> local_x(x.get(r_c_indices.rows));
258 
259  return _local_assemblers[element_id]->getFlux(p, t, local_x);
260 }
261 
262 // this is almost a copy of the implemention in the GroundwaterFlow
263 void HTProcess::postTimestepConcreteProcess(GlobalVector const& x,
264  const double t,
265  const double /*delta_t*/,
266  int const process_id)
267 {
268  // For the monolithic scheme, process_id is always zero.
269  if (_use_monolithic_scheme && process_id != 0)
270  {
271  OGS_FATAL(
272  "The condition of process_id = 0 must be satisfied for "
273  "monolithic HTProcess, which is a single process.");
274  }
275  if (!_use_monolithic_scheme && process_id != 1)
276  {
277  DBUG("This is the thermal part of the staggered HTProcess.");
278  return;
279  }
280  if (!_surfaceflux) // computing the surfaceflux is optional
281  {
282  return;
283  }
284 
285  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
286 
287  _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
288  pv.getActiveElementIDs());
289  _surfaceflux->save(t);
290 }
291 
292 } // namespace HT
293 } // namespace ProcessLib
std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const override
Definition: HTProcess.cpp:200
MeshLib::Mesh & _mesh
Definition: Process.h:261
void setStaggeredCoupledSolutions(std::size_t const, CoupledSolutionsForStaggeredScheme *const coupling_term)
void preTimestepConcreteProcess(GlobalVector const &x, double const t, double const dt, const int process_id) override
Definition: HTProcess.cpp:159
void assembleConcreteProcess(const double t, GlobalVector const &x, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
Definition: HTProcess.cpp:87
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition: Process.h:262
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
std::unique_ptr< ProcessLib::SurfaceFluxData > _surfaceflux
Definition: HTProcess.h:120
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:123
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
unsigned const _integration_order
Definition: Process.h:284
HTProcess(MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterBase >> const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable >>> &&process_variables, std::unique_ptr< HTMaterialProperties > &&material_properties, SecondaryVariableCollection &&secondary_variables, NumLib::NamedFunctionCaller &&named_function_caller, bool const use_monolithic_scheme, std::unique_ptr< ProcessLib::SurfaceFluxData > &&surfaceflux)
Definition: HTProcess.cpp:27
void setCoupledSolutionsOfPreviousTimeStep()
Definition: HTProcess.cpp:223
virtual std::vector< double > const & getIntPtDarcyVelocity(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &) const =0
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, double const t, GlobalVector const &x, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, CoupledSolutionsForStaggeredScheme const *const cpl_xs)
NumLib::Extrapolator & getExtrapolator() const
Definition: Process.h:150
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, const double t, GlobalVector const &x, GlobalVector const &xdot, const double dxdot_dx, const double dx_dx, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac, CoupledSolutionsForStaggeredScheme const *const cpl_xs)
bool isAxiallySymmetric() const
Definition: Mesh.h:137
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:264
std::vector< std::reference_wrapper< GlobalVector const > > const & coupled_xs
References to the current solutions of the coupled processes.
Builds expression trees of named functions dynamically at runtime.
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
SecondaryVariableCollection _secondary_variables
Definition: Process.h:266
Eigen::Vector3d getFlux(std::size_t element_id, MathLib::Point3d const &p, double const t, GlobalVector const &x) const override
Definition: HTProcess.cpp:248
void assembleWithJacobianConcreteProcess(const double t, GlobalVector const &x, GlobalVector const &xdot, const double dxdot_dx, const double dx_dx, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
Definition: HTProcess.cpp:129
void postTimestepConcreteProcess(GlobalVector const &x, const double t, const double delta_t, int const process_id) override
Definition: HTProcess.cpp:263
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
Definition: HTProcess.cpp:48
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:108
CoupledSolutionsForStaggeredScheme * _coupled_solutions
Definition: Process.h:279
Handles configuration of several secondary variables from the project file.
unsigned getShapeFunctionOrder() const
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:81
Ordering data by spatial location.
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
void addSecondaryVariable(std::string const &internal_name, SecondaryVariableFunctions &&fcts)
NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalIndexType > &indices)
const std::unique_ptr< HTMaterialProperties > _material_properties
Definition: HTProcess.h:113
std::array< std::unique_ptr< GlobalVector >, 2 > _xs_previous_timestep
Solutions of the previous time step.
Definition: HTProcess.h:118
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
const bool _use_monolithic_scheme
Definition: Process.h:275
void setCoupledTermForTheStaggeredSchemeToLocalAssemblers() override
Definition: HTProcess.cpp:186
VectorMatrixAssembler _global_assembler
Definition: Process.h:273
std::vector< std::unique_ptr< HTLocalAssemblerInterface > > _local_assemblers
Definition: HTProcess.h:115
std::vector< GlobalVector * > coupled_xs_t0
Pointers to the vector of the solutions of the previous time step.
std::vector< std::size_t > & getActiveElementIDs() const