OGS 6.2.0-97-g4a610c866
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<ParameterLib::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  const int heat_transport_process_id,
40  const int hydraulic_process_id)
41  : Process(mesh, std::move(jacobian_assembler), parameters,
42  integration_order, std::move(process_variables),
43  std::move(secondary_variables), std::move(named_function_caller),
44  use_monolithic_scheme),
45  _material_properties(std::move(material_properties)),
46  _surfaceflux(std::move(surfaceflux)),
47  _heat_transport_process_id(heat_transport_process_id),
48  _hydraulic_process_id(hydraulic_process_id)
49 {
50 }
51 
53  NumLib::LocalToGlobalIndexMap const& dof_table,
54  MeshLib::Mesh const& mesh,
55  unsigned const integration_order)
56 {
57  // For the staggered scheme, both processes are assumed to use the same
58  // element order. Therefore the order of shape function can be fetched from
59  // any set of the sets of process variables of the coupled processes. Here,
60  // we take the one from the first process by setting process_id = 0.
61  const int process_id = 0;
62  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
63 
65  {
66  ProcessLib::createLocalAssemblers<MonolithicHTFEM>(
67  mesh.getDimension(), mesh.getElements(), dof_table,
69  mesh.isAxiallySymmetric(), integration_order,
71  }
72  else
73  {
74  ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
75  mesh.getDimension(), mesh.getElements(), dof_table,
77  mesh.isAxiallySymmetric(), integration_order, *_material_properties,
79  }
80 
82  "darcy_velocity",
86 }
87 
89  GlobalVector const& x,
90  GlobalMatrix& M,
91  GlobalMatrix& K,
92  GlobalVector& b)
93 {
94  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
95  dof_tables;
97  {
98  DBUG("Assemble HTProcess.");
99  dof_tables.emplace_back(*_local_to_global_index_map);
100  }
101  else
102  {
104  {
105  DBUG(
106  "Assemble the equations of heat transport process within "
107  "HTProcess.");
108  }
109  else
110  {
111  DBUG(
112  "Assemble the equations of single phase fully saturated "
113  "fluid flow process within HTProcess.");
114  }
116  dof_tables.emplace_back(*_local_to_global_index_map);
117  dof_tables.emplace_back(*_local_to_global_index_map);
118  }
119 
120  const int process_id =
122  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
123  // Call global assembler for each local assembly item.
126  pv.getActiveElementIDs(), dof_tables, t, x, M, K, b,
128 }
129 
131  const double t, GlobalVector const& x, GlobalVector const& xdot,
132  const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
133  GlobalVector& b, GlobalMatrix& Jac)
134 {
135  DBUG("AssembleWithJacobian HTProcess.");
136 
137  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
138  dof_tables;
140  {
142  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
143  }
144  else
145  {
146  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
147  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
148  }
149 
150  // Call global assembler for each local assembly item.
151  const int process_id =
153  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
156  _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, x, xdot,
157  dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
158 }
159 
160 void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
161  const double /*t*/,
162  const double /*delta_t*/,
163  const int process_id)
164 {
165  assert(process_id < 2);
166 
168  {
169  return;
170  }
171 
172  if (!_xs_previous_timestep[process_id])
173  {
174  _xs_previous_timestep[process_id] =
176  }
177  else
178  {
179  auto& x0 = *_xs_previous_timestep[process_id];
180  MathLib::LinAlg::copy(x, x0);
181  }
182 
183  auto& x0 = *_xs_previous_timestep[process_id];
184  MathLib::LinAlg::setLocalAccessibleVector(x0);
185 }
186 
188 {
189  DBUG("Set the coupled term for the staggered scheme to local assembers.");
190 
191  const int process_id =
193  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  const int process_id)
225 {
226  const auto& x_t0 = _xs_previous_timestep[process_id];
227  if (x_t0 == nullptr)
228  {
229  OGS_FATAL(
230  "Memory is not allocated for the global vector of the solution of "
231  "the previous time step for the staggered scheme.\n It can be done "
232  "by overriding Process::preTimestepConcreteProcess (ref. "
233  "HTProcess::preTimestepConcreteProcess) ");
234  }
235 
236  _coupled_solutions->coupled_xs_t0[process_id] = x_t0.get();
237 }
238 
240 {
244 }
245 
246 Eigen::Vector3d HTProcess::getFlux(std::size_t element_id,
247  MathLib::Point3d const& p,
248  double const t,
249  GlobalVector const& x) const
250 {
251  // fetch local_x from primary variable
252  std::vector<GlobalIndexType> indices_cache;
253  auto const r_c_indices = NumLib::getRowColumnIndices(
254  element_id, *_local_to_global_index_map, indices_cache);
255  std::vector<double> local_x(x.get(r_c_indices.rows));
256 
257  return _local_assemblers[element_id]->getFlux(p, t, local_x);
258 }
259 
260 // this is almost a copy of the implementation in the GroundwaterFlow
261 void HTProcess::postTimestepConcreteProcess(GlobalVector const& x,
262  const double t,
263  const double /*delta_t*/,
264  int const process_id)
265 {
266  // For the monolithic scheme, process_id is always zero.
267  if (_use_monolithic_scheme && process_id != 0)
268  {
269  OGS_FATAL(
270  "The condition of process_id = 0 must be satisfied for "
271  "monolithic HTProcess, which is a single process.");
272  }
273  if (!_use_monolithic_scheme && process_id != _hydraulic_process_id)
274  {
275  DBUG("This is the thermal part of the staggered HTProcess.");
276  return;
277  }
278  if (!_surfaceflux) // computing the surfaceflux is optional
279  {
280  return;
281  }
282 
283  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
284 
285  _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
286  pv.getActiveElementIDs());
287  _surfaceflux->save(t);
288 }
289 
290 } // namespace HT
291 } // namespace ProcessLib
std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const override
Definition: HTProcess.cpp:200
MeshLib::Mesh & _mesh
Definition: Process.h:262
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:160
void assembleConcreteProcess(const double t, GlobalVector const &x, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
Definition: HTProcess.cpp:88
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition: Process.h:263
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:125
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:124
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:285
void setCoupledSolutionsOfPreviousTimeStep()
Definition: HTProcess.cpp:239
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)
const int _hydraulic_process_id
Definition: HTProcess.h:128
NumLib::Extrapolator & getExtrapolator() const
Definition: Process.h:151
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:265
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:267
Eigen::Vector3d getFlux(std::size_t element_id, MathLib::Point3d const &p, double const t, GlobalVector const &x) const override
Definition: HTProcess.cpp:246
HTProcess(MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::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, const int heat_transport_process_id, const int hydraulic_process_id)
Definition: HTProcess.cpp:27
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:130
void postTimestepConcreteProcess(GlobalVector const &x, const double t, const double delta_t, int const process_id) override
Definition: HTProcess.cpp:261
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:52
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:108
CoupledSolutionsForStaggeredScheme * _coupled_solutions
Definition: Process.h:280
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
const int _heat_transport_process_id
Definition: HTProcess.h:127
void setCoupledSolutionsOfPreviousTimeStepPerProcess(const int process_id)
Definition: HTProcess.cpp:223
Ordering data by spatial location.
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
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:118
std::array< std::unique_ptr< GlobalVector >, 2 > _xs_previous_timestep
Solutions of the previous time step.
Definition: HTProcess.h:123
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
const bool _use_monolithic_scheme
Definition: Process.h:276
void setCoupledTermForTheStaggeredSchemeToLocalAssemblers() override
Definition: HTProcess.cpp:187
VectorMatrixAssembler _global_assembler
Definition: Process.h:274
std::vector< std::unique_ptr< HTLocalAssemblerInterface > > _local_assemblers
Definition: HTProcess.h:120
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