OGS 6.2.2-330-gf48c72f61.dirty.20200225212913
HTProcess.cpp
Go to the documentation of this file.
1 
11 #include "HTProcess.h"
12 
13 #include <cassert>
14 
19 
20 #include "MonolithicHTFEM.h"
21 #include "StaggeredHTFEM.h"
22 
23 namespace ProcessLib
24 {
25 namespace HT
26 {
28  HTProcessData const& process_data)
29 {
30  DBUG("Check the media properties of HT process ...");
31 
32  std::array const requiredPropertyMedium = {
35 
36  std::array const requiredPropertyLiquidPhase = {
40 
41  std::array const requiredPropertySolidPhase = {
45 
46  for (auto const& element : mesh.getElements())
47  {
48  auto const element_id = element->getID();
49 
50  auto const& medium = *process_data.media_map->getMedium(element_id);
52  medium, requiredPropertyMedium);
53 
55  medium.phase("AqueousLiquid"), requiredPropertyLiquidPhase);
56 
58  medium.phase("Solid"), requiredPropertySolidPhase);
59  }
60  DBUG("Media properties verified.");
61 }
62 
64  std::string name,
65  MeshLib::Mesh& mesh,
66  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
67  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
68  unsigned const integration_order,
69  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
70  process_variables,
71  HTProcessData&& process_data,
72  SecondaryVariableCollection&& secondary_variables,
73  bool const use_monolithic_scheme,
74  std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux,
75  const int heat_transport_process_id,
76  const int hydraulic_process_id)
77  : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
78  integration_order, std::move(process_variables),
79  std::move(secondary_variables), use_monolithic_scheme),
80  _process_data(std::move(process_data)),
81  _surfaceflux(std::move(surfaceflux)),
82  _heat_transport_process_id(heat_transport_process_id),
83  _hydraulic_process_id(hydraulic_process_id)
84 {
85 }
86 
88  NumLib::LocalToGlobalIndexMap const& dof_table,
89  MeshLib::Mesh const& mesh,
90  unsigned const integration_order)
91 {
93 
94  // For the staggered scheme, both processes are assumed to use the same
95  // element order. Therefore the order of shape function can be fetched from
96  // any set of the sets of process variables of the coupled processes. Here,
97  // we take the one from the first process by setting process_id = 0.
98  const int process_id = 0;
99  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
100 
102  {
103  ProcessLib::createLocalAssemblers<MonolithicHTFEM>(
104  mesh.getDimension(), mesh.getElements(), dof_table,
106  mesh.isAxiallySymmetric(), integration_order, _process_data);
107  }
108  else
109  {
110  ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
111  mesh.getDimension(), mesh.getElements(), dof_table,
113  mesh.isAxiallySymmetric(), integration_order, _process_data,
115  }
116 
118  "darcy_velocity",
122 }
123 
124 void HTProcess::assembleConcreteProcess(const double t, double const dt,
125  std::vector<GlobalVector*> const& x,
126  std::vector<GlobalVector*> const& xdot,
127  int const process_id, GlobalMatrix& M,
128  GlobalMatrix& K, GlobalVector& b)
129 {
130  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
131  dof_tables;
133  {
134  DBUG("Assemble HTProcess.");
135  dof_tables.emplace_back(*_local_to_global_index_map);
136  }
137  else
138  {
139  if (process_id == _heat_transport_process_id)
140  {
141  DBUG(
142  "Assemble the equations of heat transport process within "
143  "HTProcess.");
144  }
145  else
146  {
147  DBUG(
148  "Assemble the equations of single phase fully saturated "
149  "fluid flow process within HTProcess.");
150  }
152  dof_tables.emplace_back(*_local_to_global_index_map);
153  dof_tables.emplace_back(*_local_to_global_index_map);
154  }
155 
156  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
157  // Call global assembler for each local assembly item.
160  pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
161  b, _coupled_solutions);
162 }
163 
165  const double t, double const dt, std::vector<GlobalVector*> const& x,
166  GlobalVector const& xdot, const double dxdot_dx, const double dx_dx,
167  int const process_id, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
168  GlobalMatrix& Jac)
169 {
170  DBUG("AssembleWithJacobian HTProcess.");
171 
172  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
173  dof_tables;
175  {
177  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
178  }
179  else
180  {
181  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
182  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
183  }
184 
185  // Call global assembler for each local assembly item.
186  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
189  _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
190  dxdot_dx, dx_dx, process_id, M, K, b, Jac, _coupled_solutions);
191 }
192 
193 void HTProcess::preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
194  const double /*t*/,
195  const double /*delta_t*/,
196  const int process_id)
197 {
198  assert(process_id < 2);
199 
201  {
202  return;
203  }
204 
205  if (!_xs_previous_timestep[process_id])
206  {
207  _xs_previous_timestep[process_id] =
209  *x[process_id]);
210  }
211  else
212  {
213  auto& x0 = *_xs_previous_timestep[process_id];
214  MathLib::LinAlg::copy(*x[process_id], x0);
215  }
216 
217  auto& x0 = *_xs_previous_timestep[process_id];
218  MathLib::LinAlg::setLocalAccessibleVector(x0);
219 }
220 
222  int const process_id)
223 {
224  DBUG("Set the coupled term for the staggered scheme to local assembers.");
225 
226  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
230 }
231 
232 std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
234 {
236  {
237  // For single-variable-single-component processes reuse the existing DOF
238  // table.
239  const bool manage_storage = false;
240  return std::make_tuple(_local_to_global_index_map.get(),
241  manage_storage);
242  }
243 
244  // Otherwise construct a new DOF table.
245  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
247 
248  const bool manage_storage = true;
249  return std::make_tuple(new NumLib::LocalToGlobalIndexMap(
250  std::move(all_mesh_subsets_single_component),
251  // by location order is needed for output
253  manage_storage);
254 }
255 
257  const int process_id)
258 {
259  const auto& x_t0 = _xs_previous_timestep[process_id];
260  if (x_t0 == nullptr)
261  {
262  OGS_FATAL(
263  "Memory is not allocated for the global vector of the solution of "
264  "the previous time step for the staggered scheme.\n It can be done "
265  "by overriding Process::preTimestepConcreteProcess (ref. "
266  "HTProcess::preTimestepConcreteProcess) ");
267  }
268 
269  _coupled_solutions->coupled_xs_t0[process_id] = x_t0.get();
270 }
271 
273 {
277 }
278 
279 Eigen::Vector3d HTProcess::getFlux(std::size_t element_id,
280  MathLib::Point3d const& p,
281  double const t,
282  std::vector<GlobalVector*> const& x) const
283 {
284  // fetch local_x from primary variable
285  std::vector<GlobalIndexType> indices_cache;
286  auto const r_c_indices = NumLib::getRowColumnIndices(
287  element_id, *_local_to_global_index_map, indices_cache);
288  std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes{
289  x.size(), r_c_indices.rows};
290  auto const local_x =
291  getCoupledLocalSolutions(x, indices_of_all_coupled_processes);
292 
293  return _local_assemblers[element_id]->getFlux(p, t, local_x);
294 }
295 
296 // this is almost a copy of the implementation in the GroundwaterFlow
297 void HTProcess::postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
298  const double t,
299  const double /*delta_t*/,
300  int const process_id)
301 {
302  // For the monolithic scheme, process_id is always zero.
303  if (_use_monolithic_scheme && process_id != 0)
304  {
305  OGS_FATAL(
306  "The condition of process_id = 0 must be satisfied for "
307  "monolithic HTProcess, which is a single process.");
308  }
309  if (!_use_monolithic_scheme && process_id != _hydraulic_process_id)
310  {
311  DBUG("This is the thermal part of the staggered HTProcess.");
312  return;
313  }
314  if (!_surfaceflux) // computing the surfaceflux is optional
315  {
316  return;
317  }
318 
319  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
320 
321  _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
322  pv.getActiveElementIDs());
323  _surfaceflux->save(t);
324 }
325 } // namespace HT
326 } // namespace ProcessLib
std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const override
Definition: HTProcess.cpp:233
MeshLib::Mesh & _mesh
Definition: Process.h:295
void setStaggeredCoupledSolutions(std::size_t const, CoupledSolutionsForStaggeredScheme *const coupling_term)
Eigen::Vector3d getFlux(std::size_t element_id, MathLib::Point3d const &p, double const t, std::vector< GlobalVector *> const &x) const override
Definition: HTProcess.cpp:279
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector *> const &x, std::vector< GlobalVector *> const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
Definition: HTProcess.cpp:124
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, double const dt, std::vector< GlobalVector *> const &x, std::vector< GlobalVector *> const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, CoupledSolutionsForStaggeredScheme const *const cpl_xs)
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition: Process.h:296
virtual std::vector< double > const & getIntPtDarcyVelocity(const double, std::vector< GlobalVector *> const &x, std::vector< NumLib::LocalToGlobalIndexMap const *> const &, std::vector< double > &) const =0
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:133
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:129
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:313
void setCoupledSolutionsOfPreviousTimeStep()
Definition: HTProcess.cpp:272
const int _hydraulic_process_id
Definition: HTProcess.h:136
HTProcessData _process_data
Definition: HTProcess.h:126
void preTimestepConcreteProcess(std::vector< GlobalVector *> const &x, double const t, double const dt, const int process_id) override
Definition: HTProcess.cpp:193
NumLib::Extrapolator & getExtrapolator() const
Definition: Process.h:157
bool isAxiallySymmetric() const
Definition: Mesh.h:134
void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(int const process_id) override
Definition: HTProcess.cpp:221
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:298
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:300
void postTimestepConcreteProcess(std::vector< GlobalVector *> const &x, const double t, const double delta_t, int const process_id) override
Definition: HTProcess.cpp:297
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector *> const &x, GlobalVector const &xdot, const double dxdot_dx, const double dx_dx, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
Definition: HTProcess.cpp:164
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector *> const &global_solutions, std::vector< std::vector< GlobalIndexType >> const &indices)
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:87
std::unique_ptr< MaterialPropertyLib::MaterialSpatialDistributionMap > media_map
Definition: HTProcessData.h:28
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:105
CoupledSolutionsForStaggeredScheme * _coupled_solutions
Definition: Process.h:308
Handles configuration of several secondary variables from the project file.
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, double const dt, std::vector< GlobalVector *> const &x, GlobalVector const &xdot, const double dxdot_dx, const double dx_dx, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac, CoupledSolutionsForStaggeredScheme const *const cpl_xs)
unsigned getShapeFunctionOrder() const
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:78
const int _heat_transport_process_id
Definition: HTProcess.h:135
void checkRequiredProperties(Component const &c, Container const &required_properties)
Definition: Component.h:89
void setCoupledSolutionsOfPreviousTimeStepPerProcess(const int process_id)
Definition: HTProcess.cpp:256
Ordering data by spatial location.
#define OGS_FATAL(fmt,...)
Definition: Error.h:64
void checkMPLProperties(MeshLib::Mesh const &mesh, HTProcessData const &process_data)
Definition: HTProcess.cpp:27
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)
std::array< std::unique_ptr< GlobalVector >, 2 > _xs_previous_timestep
Solutions of the previous time step.
Definition: HTProcess.h:131
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:37
HTProcess(std::string name, 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, HTProcessData &&process_data, SecondaryVariableCollection &&secondary_variables, 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:63
const bool _use_monolithic_scheme
Definition: Process.h:304
VectorMatrixAssembler _global_assembler
Definition: Process.h:302
std::vector< std::unique_ptr< HTLocalAssemblerInterface > > _local_assemblers
Definition: HTProcess.h:128
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