OGS
PhaseFieldProcess.cpp
Go to the documentation of this file.
1
11#include "PhaseFieldProcess.h"
12
13#include <cassert>
14
17#include "PhaseFieldFEM.h"
18#include "ProcessLib/Process.h"
20
21namespace ProcessLib
22{
23namespace PhaseField
24{
25template <int DisplacementDim>
27 std::string name,
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,
35 SecondaryVariableCollection&& secondary_variables,
36 bool const use_monolithic_scheme)
37 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
38 integration_order, std::move(process_variables),
39 std::move(secondary_variables), use_monolithic_scheme),
40 _process_data(std::move(process_data))
41{
42 if (use_monolithic_scheme)
43 {
45 "Monolithic scheme is not implemented for the PhaseField process.");
46 }
47
48 _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
49 mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
50}
51
52template <int DisplacementDim>
54{
55 return false;
56}
57
58template <int DisplacementDim>
61 const int process_id) const
62{
63 // For the M process (deformation) in the staggered scheme.
64 if (process_id == 0)
65 {
66 auto const& l = *_local_to_global_index_map;
67 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
68 &l.getGhostIndices(), &this->_sparsity_pattern};
69 }
70
71 // For staggered scheme and phase field process.
72 auto const& l = *_local_to_global_index_map_single_component;
73 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
74 &l.getGhostIndices(), &_sparsity_pattern_with_single_component};
75}
76
77template <int DisplacementDim>
80{
81 // For the M process (deformation) in the staggered scheme.
82 if (process_id == 0)
83 {
84 return *_local_to_global_index_map;
85 }
86
87 // For the equation of phasefield
88 return *_local_to_global_index_map_single_component;
89}
90
91template <int DisplacementDim>
93{
94 // For displacement equation.
95 const int mechanics_process_id = 0;
96 constructDofTableOfSpecifiedProcessStaggeredScheme(mechanics_process_id);
97
98 // TODO move the two data members somewhere else.
99 // for extrapolation of secondary variables of stress or strain
100 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
101 *_mesh_subset_all_nodes};
102 _local_to_global_index_map_single_component =
103 std::make_unique<NumLib::LocalToGlobalIndexMap>(
104 std::move(all_mesh_subsets_single_component),
105 // by location order is needed for output
107
108 assert(_local_to_global_index_map_single_component);
109
110 // For phase field equation.
111 _sparsity_pattern_with_single_component = NumLib::computeSparsityPattern(
112 *_local_to_global_index_map_single_component, _mesh);
113}
114
115template <int DisplacementDim>
117 NumLib::LocalToGlobalIndexMap const& dof_table,
118 MeshLib::Mesh const& mesh,
119 unsigned const integration_order)
120{
122 DisplacementDim, PhaseFieldLocalAssembler>(
123 mesh.getElements(), dof_table, _local_assemblers,
124 NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
125 _process_data);
126
127 _secondary_variables.addSecondaryVariable(
128 "sigma",
130 DisplacementDim>::RowsAtCompileTime,
131 getExtrapolator(), _local_assemblers,
132 &LocalAssemblerInterface::getIntPtSigma));
133
134 _secondary_variables.addSecondaryVariable(
135 "epsilon",
137 DisplacementDim>::RowsAtCompileTime,
138 getExtrapolator(), _local_assemblers,
139 &LocalAssemblerInterface::getIntPtEpsilon));
140
141 _secondary_variables.addSecondaryVariable(
142 "sigma_tensile",
144 DisplacementDim>::RowsAtCompileTime,
145 getExtrapolator(), _local_assemblers,
146 &LocalAssemblerInterface::getIntPtSigmaTensile));
147
148 _secondary_variables.addSecondaryVariable(
149 "sigma_compressive",
151 DisplacementDim>::RowsAtCompileTime,
152 getExtrapolator(), _local_assemblers,
153 &LocalAssemblerInterface::getIntPtSigmaCompressive));
154
155 _secondary_variables.addSecondaryVariable(
156 "eps_tensile",
158 DisplacementDim>::RowsAtCompileTime,
159 getExtrapolator(), _local_assemblers,
160 &LocalAssemblerInterface::getIntPtEpsilonTensile));
161
162 // Initialize local assemblers after all variables have been set.
164 &LocalAssemblerInterface::initialize, _local_assemblers,
165 *_local_to_global_index_map);
166}
167
168template <int DisplacementDim>
170 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
171{
172 // Staggered scheme:
173 // for the equations of deformation.
174 const int mechanical_process_id = 0;
175 initializeProcessBoundaryConditionsAndSourceTerms(
176 *_local_to_global_index_map, mechanical_process_id, media);
177 // for the phase field
178 const int phasefield_process_id = 1;
179 initializeProcessBoundaryConditionsAndSourceTerms(
180 *_local_to_global_index_map_single_component, phasefield_process_id,
181 media);
182}
183
184template <int DisplacementDim>
186 const double t, double const dt, std::vector<GlobalVector*> const& x,
187 std::vector<GlobalVector*> const& x_prev, int const process_id,
189{
190 DBUG("Assemble PhaseFieldProcess.");
191
192 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
193
194 // For the staggered scheme
195 if (process_id == 1)
196 {
197 DBUG(
198 "Assemble the equations of phase field in "
199 "PhaseFieldProcess for the staggered scheme.");
200 }
201 else
202 {
203 DBUG(
204 "Assemble the equations of deformation in "
205 "PhaseFieldProcess for the staggered scheme.");
206 }
207 dof_tables.emplace_back(_local_to_global_index_map_single_component.get());
208 dof_tables.emplace_back(_local_to_global_index_map.get());
209
210 // Call global assembler for each local assembly item.
212 _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
213 getActiveElementIDs(), dof_tables, t, dt, x, x_prev, process_id, &M, &K,
214 &b);
215}
216
217template <int DisplacementDim>
219 const double t, double const dt, std::vector<GlobalVector*> const& x,
220 std::vector<GlobalVector*> const& x_prev, int const process_id,
221 GlobalVector& b, GlobalMatrix& Jac)
222{
223 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
224
225 // For the staggered scheme
226 if (process_id == 1)
227 {
228 DBUG(
229 "Assemble the Jacobian equations of phase field in "
230 "PhaseFieldProcess for the staggered scheme.");
231 }
232 else
233 {
234 DBUG(
235 "Assemble the Jacobian equations of deformation in "
236 "PhaseFieldProcess for the staggered scheme.");
237 }
238 dof_tables.emplace_back(_local_to_global_index_map.get());
239 dof_tables.emplace_back(_local_to_global_index_map_single_component.get());
240
241 // Call global assembler for each local assembly item.
244 _local_assemblers, getActiveElementIDs(), dof_tables, t, dt, x, x_prev,
245 process_id, &b, &Jac);
246
247 if (process_id == 0)
248 {
249 b.copyValues(*_nodal_forces);
250 std::transform(_nodal_forces->begin(), _nodal_forces->end(),
251 _nodal_forces->begin(), [](double val) { return -val; });
252 }
253}
254
255template <int DisplacementDim>
257 std::vector<GlobalVector*> const& x, double const t, double const dt,
258 const int process_id)
259{
260 DBUG("PreTimestep PhaseFieldProcess {:d}.", process_id);
261
262 _process_data.injected_volume = t;
263
264 _x_previous_timestep =
266
268 &LocalAssemblerInterface::preTimestep, _local_assemblers,
269 getActiveElementIDs(), getDOFTable(process_id), *x[process_id], t, dt);
270}
271
272template <int DisplacementDim>
274 std::vector<GlobalVector*> const& x,
275 std::vector<GlobalVector*> const& /*x_prev*/, const double t,
276 const double /*delta_t*/, int const process_id)
277{
278 if (isPhaseFieldProcess(process_id))
279 {
280 DBUG("PostTimestep PhaseFieldProcess.");
281
282 _process_data.elastic_energy = 0.0;
283 _process_data.surface_energy = 0.0;
284 _process_data.pressure_work = 0.0;
285
286 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
287
288 dof_tables.emplace_back(_local_to_global_index_map.get());
289 dof_tables.emplace_back(
290 _local_to_global_index_map_single_component.get());
291
293 &LocalAssemblerInterface::computeEnergy, _local_assemblers,
294 getActiveElementIDs(), dof_tables, x, t,
295 _process_data.elastic_energy, _process_data.surface_energy,
296 _process_data.pressure_work);
297
298#ifdef USE_PETSC
299 double const elastic_energy = _process_data.elastic_energy;
300 MPI_Allreduce(&elastic_energy, &_process_data.elastic_energy, 1,
301 MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
302 double const surface_energy = _process_data.surface_energy;
303 MPI_Allreduce(&surface_energy, &_process_data.surface_energy, 1,
304 MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
305 double const pressure_work = _process_data.pressure_work;
306 MPI_Allreduce(&pressure_work, &_process_data.pressure_work, 1,
307 MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
308#endif
309
310 INFO(
311 "Elastic energy: {:g} Surface energy: {:g} Pressure work: {:g} at "
312 "time: {:g} ",
313 _process_data.elastic_energy, _process_data.surface_energy,
314 _process_data.pressure_work, t);
315 if (_process_data.propagating_pressurized_crack)
316 {
317 INFO("Pressure: {:g} at time: {:g} ", _process_data.pressure, t);
318 }
319 }
320}
321
322template <int DisplacementDim>
324 std::vector<GlobalVector*> const& x,
325 std::vector<GlobalVector*> const& /*x_prev*/, const double t,
326 double const /*dt*/, const int process_id)
327{
328 _process_data.crack_volume = 0.0;
329
330 if (isPhaseFieldProcess(process_id))
331 {
332 if (_process_data.propagating_pressurized_crack)
333 {
334 auto& u = *x[0];
335 MathLib::LinAlg::scale(const_cast<GlobalVector&>(u),
336 1 / _process_data.pressure);
337 }
338 return;
339 }
340
341 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
342
343 dof_tables.emplace_back(_local_to_global_index_map.get());
344 dof_tables.emplace_back(_local_to_global_index_map_single_component.get());
345
346 DBUG("PostNonLinearSolver crack volume computation.");
347
349 &LocalAssemblerInterface::computeCrackIntegral, _local_assemblers,
350 getActiveElementIDs(), dof_tables, x, t, _process_data.crack_volume);
351
352#ifdef USE_PETSC
353 double const crack_volume = _process_data.crack_volume;
354 MPI_Allreduce(&crack_volume, &_process_data.crack_volume, 1, MPI_DOUBLE,
355 MPI_SUM, PETSC_COMM_WORLD);
356#endif
357
358 INFO("Integral of crack: {:g}", _process_data.crack_volume);
359
360 if (_process_data.propagating_pressurized_crack)
361 {
362 _process_data.pressure_old = _process_data.pressure;
363 _process_data.pressure =
364 _process_data.injected_volume / _process_data.crack_volume;
365 _process_data.pressure_error =
366 std::abs(_process_data.pressure_old - _process_data.pressure) /
367 _process_data.pressure;
368 INFO("Internal pressure: {:g} and Pressure error: {:.4e}",
369 _process_data.pressure, _process_data.pressure_error);
370
371 auto& u = *x[0];
372 MathLib::LinAlg::scale(const_cast<GlobalVector&>(u),
373 _process_data.pressure);
374 }
375}
376
377template <int DisplacementDim>
379 GlobalVector& lower, GlobalVector& upper, int const /*process_id*/)
380{
381 lower.setZero();
382 MathLib::LinAlg::setLocalAccessibleVector(*_x_previous_timestep);
383 MathLib::LinAlg::copy(*_x_previous_timestep, upper);
384
385 GlobalIndexType const x_begin = _x_previous_timestep->getRangeBegin();
386 GlobalIndexType const x_end = _x_previous_timestep->getRangeEnd();
387
388 for (GlobalIndexType i = x_begin; i < x_end; i++)
389 {
390 if ((*_x_previous_timestep)[i] > _process_data.irreversible_threshold)
391 {
392 upper.set(i, 1.0);
393 }
394 }
395}
396
397template <int DisplacementDim>
399 int const process_id) const
400{
401 return process_id == 1;
402}
403
404template class PhaseFieldProcess<2>;
405template class PhaseFieldProcess<3>;
406
407} // namespace PhaseField
408} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
GlobalMatrix::IndexType GlobalIndexType
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Global vector based on Eigen vector.
Definition EigenVector.h:25
void copyValues(std::vector< double > &u) const
void set(IndexType rowId, double v)
set entry
Definition EigenVector.h:73
bool isAxiallySymmetric() const
Definition Mesh.h:137
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:109
virtual void preTimestep(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void initialize(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
bool isPhaseFieldProcess(int const process_id) const
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
void updateConstraints(GlobalVector &lower, GlobalVector &upper, int const process_id) override
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
MeshLib::PropertyVector< double > * _nodal_forces
void postNonLinearSolverConcreteProcess(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, const double t, double const dt, int const process_id) override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, const double t, const double delta_t, int const process_id) override
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
void initializeBoundaryConditions(std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media) override
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
PhaseFieldProcess(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, PhaseFieldProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector &b, GlobalMatrix &Jac) override
Handles configuration of several secondary variables from the project file.
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix *M, GlobalMatrix *K, GlobalVector *b)
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector *b, GlobalMatrix *Jac)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
void copy(PETScVector const &x, PETScVector &y)
Definition LinAlg.cpp:37
void setLocalAccessibleVector(PETScVector const &x)
Definition LinAlg.cpp:27
void scale(PETScVector &x, PetscScalar const a)
Definition LinAlg.cpp:44
@ BY_LOCATION
Ordering data by spatial location.
GlobalSparsityPattern computeSparsityPattern(LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
Computes a sparsity pattern for the given inputs.
void createLocalAssemblers(std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface > > &local_assemblers, ProviderOrOrder const &provider_or_order, ExtraCtorArgs &&... extra_ctor_args)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)