35 inline double operator() (
double alpha,
double density,
double gravity,
double pressure)
37 return -alpha*density*gravity*pressure;
63 inline double operator() (
double alpha,
double beta,
double pressure,
double old_pressure,
double div_u,
double old_div_u)
65 return (beta*(pressure-old_pressure) + alpha*(old_div_u - div_u)) /
time_->
dt();
75 std::string equation_name =
"Coupling_Iterative";
77 "Record with data for iterative coupling of flow and mechanics.\n")
84 "Flow equation, provides the velocity field as a result.")
87 "Mechanics, provides the displacement field.")
90 .make_field_descriptor_type(equation_name)),
92 "Input fields of the HM coupling.")
94 "Tuning parameter for iterative splitting. Its default value "
95 "corresponds to a theoretically optimal value with fastest convergence." )
97 "Maximal count of HM iterations." )
99 "Minimal count of HM iterations." )
101 "Absolute tolerance for difference in HM iteration." )
103 "Relative tolerance for difference in HM iteration." )
108 const int HM_Iterative::registrar = Input::register_class< HM_Iterative, Mesh &, const Input::Record >(
"Coupling_Iterative")
121 .
description(
"Volumetric mass density of the fluid.")
127 .
description(
"Gravitational acceleration constant.")
132 *
this +=
beta.
name(
"relaxation_beta")
133 .
description(
"Parameter of numerical method for iterative solution of hydro-mechanical coupling.")
138 .
description(
"Coupling term entering the mechanics equation.")
148 .
description(
"Pressure from last computed iteration.")
153 .
description(
"Displacement divergence from last computed time.")
158 .
description(
"Pressure potential on boundary (taking into account the flow boundary condition.")
163 .
description(
"Coupling term entering the flow equation.")
181 eq_data.
flow_->eq_fields().field_edge_pressure
185 ref_pressure_potential.set(ref_potential_ptr_, 0.0);
196 auto old_pressure_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(eq_data.
flow_->eq_data().dh_cr_, &eq_data.
flow_->eq_data().p_edge_solution_previous_time);
197 old_pressure.set(old_pressure_ptr_, 0.0);
200 old_iter_pressure_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(eq_data.
flow_->eq_fields().field_ele_pressure.get_field_fe()->get_dofhandler(),
nullptr, 1);
201 old_iter_pressure.set(old_iter_pressure_ptr_, 0.0);
203 old_div_u_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(eq_data.
mechanics_->eq_fields().output_div_ptr->get_dofhandler());
204 old_div_u.set(old_div_u_ptr_, 0.0);
210 eq_data.
flow_->eq_fields().field_edge_pressure,
212 eq_data.
mechanics_->eq_fields().output_divergence,
226 using namespace Input;
229 eq_data_ = std::make_shared<EqData>();
245 eq_data_->mechanics_->eq_fields()[
"cross_section"].copy_from(
eq_data_->flow_->eq_fields()[
"cross_section"]);
246 eq_data_->flow_->eq_fields() +=
eq_data_->mechanics_->eq_fields()[
"cross_section_updated"];
247 eq_data_->flow_->eq_fields() +=
eq_data_->mechanics_->eq_fields()[
"stress"];
248 eq_data_->flow_->eq_fields() +=
eq_data_->mechanics_->eq_fields()[
"von_mises_stress"];
249 eq_data_->flow_->eq_fields() +=
eq_data_->mechanics_->eq_fields()[
"mean_stress"];
251 std::stringstream ss;
282 template<
int dim,
class Value>
294 std::stringstream ss;
298 eq_data_->mechanics_->update_output_fields();
301 eq_data_->mechanics_->zero_time_step();
305 eq_fields_->old_iter_pressure.set_time_result_changed();
323 eq_data_->flow_->solve_time_step(
false);
327 eq_data_->mechanics_->solve_linear_system();
333 eq_data_->mechanics_->update_output_fields();
335 eq_fields_->old_iter_pressure.set_time_result_changed();
341 eq_data_->flow_->accept_time_step();
343 eq_data_->mechanics_->output_data();
351 auto ref_potential_vec_ =
eq_fields_->ref_potential_ptr_->vec();
352 auto dh =
eq_fields_->ref_potential_ptr_->get_dofhandler();
354 ref_potential_vec_.zero_entries();
357 ref_potential_vec_.local_to_ghost_begin();
358 ref_potential_vec_.local_to_ghost_end();
359 eq_fields_->pressure_potential.set_time_result_changed();
360 eq_fields_->ref_pressure_potential.set_time_result_changed();
368 eq_fields_->flow_source.set_time_result_changed();
376 auto dh =
eq_fields_->old_iter_pressure_ptr_->get_dofhandler();
385 abs_error = sqrt(recv_data[0]);
386 rel_error = abs_error / (sqrt(recv_data[1]) + std::numeric_limits<double>::min());
388 MessageOut().fmt(
"HM Iteration {} abs. difference: {} rel. difference: {}\n"
389 "--------------------------------------------------------",
393 THROW(ExcSolverDiverge() << EI_Reason(
"Reached max_it."));
static Input::Type::Abstract & get_input_type()
static const Input::Type::Record & get_input_type()
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportDG.
Input::Record input_record_
static Input::Type::Record & record_template()
Template Record with common keys for derived equations.
std::shared_ptr< FieldSet > eq_fieldset_
void init_user_fields(Input::Array user_fields, FieldSet &output_fields)
static Input::Type::Record & user_fields_template(std::string equation_name)
Template Record with common key user_fields for derived equations.
static bool print_message_table(ostream &stream, std::string equation_name)
FieldCommon & description(const string &description)
FieldCommon & flags_add(FieldFlag::Flags::Mask mask)
FieldCommon & flags(FieldFlag::Flags::Mask mask)
FieldCommon & name(const string &name)
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
FieldCommon & input_default(const string &input_default)
std::shared_ptr< DOFHandlerMultiDim > get_dofhandler() const
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
static constexpr Mask equation_external_output
Match an output field, that can be also copy of other field.
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
Container for various descendants of FieldCommonBase.
void set_default_fieldset()
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
std::shared_ptr< DarcyLMH > flow_
steady or unsteady water flow simulator based on MH scheme
std::shared_ptr< Elasticity > mechanics_
solute transport with chemistry through operator splitting
Field< 3, FieldValue< 3 >::Scalar > flow_source
Field< 3, FieldValue< 3 >::Scalar > density
Density of fluid.
Field< 3, FieldValue< 3 >::Scalar > old_pressure
Field< 3, FieldValue< 3 >::Scalar > old_iter_pressure
Field< 3, FieldValue< 3 >::Scalar > beta
Field< 3, FieldValue< 3 >::Scalar > ref_pressure_potential
Potential of reference (prescribed) pressure from flow b.c. TODO: Swith to BCField when possible.
Field< 3, FieldValue< 3 >::Scalar > alpha
Biot coefficient.
Field< 3, FieldValue< 3 >::Scalar > old_div_u
void initialize(Mesh &mesh, HM_Iterative::EqData &eq_data, const TimeGovernor *time_, double beta_)
Field< 3, FieldValue< 3 >::Scalar > pressure_potential
Potential -alpha*pressure whose gradient is passed to mechanics as additional load.
Field< 3, FieldValue< 3 >::Scalar > gravity
Standard gravity.
void solve_iteration() override
Solve equations and update data (fields).
GenericAssembly< FlowPotentialAssemblyHM > * flow_potential_assembly_
GenericAssembly< ResidualAssemblyHM > * residual_assembly_
void initialize() override
static const Input::Type::Record & get_input_type()
Define input record.
static const int registrar
void update_after_iteration() override
Save data (e.g. solution fields) for the next iteration.
std::shared_ptr< EqData > eq_data_
void update_solution() override
void compute_iteration_error(double &abs_error, double &rel_error) override
Compute absolute and relative error in the solution.
void zero_time_step() override
void update_flow_fields()
HM_Iterative(Mesh &mesh, Input::Record in_record)
void update_after_converged() override
Save data after iterations have finished.
std::shared_ptr< EqFields > eq_fields_
double r_tol_
Relative tolerance for difference between two succeeding iterations.
double a_tol_
Absolute tolerance for difference between two succeeding iterations.
unsigned int max_it_
Maximal number of iterations.
static const Input::Type::Record & record_template()
Basic time management functionality for unsteady (and steady) solvers (class Equation).
void view(const char *name="") const
const TimeStep & step(int index=-1) const
void next_time()
Proceed to the next time according to current estimated time step.
Class for representation SI units of Fields.
void copy_from(const VectorMPI &other)
Lumped mixed-hybrid model of linear Darcy flow, possibly unsteady.
double lame_mu(double young, double poisson)
double lame_lambda(double young, double poisson)
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
void copy_field(const FieldFE< dim, Value > &from_field, FieldFE< dim, Value > &to_field)
#define WarningOut()
Macro defining 'warning' record of log.
#define MessageOut()
Macro defining 'message' record of log.
#define MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm)
fn_fluid_source(const TimeGovernor *time)
double operator()(double alpha, double beta, double pressure, double old_pressure, double div_u, double old_div_u)
const TimeGovernor * time_
fn_hm_coupling_beta(double beta_f)
double operator()(double alpha, double lame_mu, double lame_lambda, double density, double gravity)
const double beta_factor
User-defined factor for iteration parameter.
double operator()(double alpha, double density, double gravity, double pressure)
#define START_TIMER(tag)
Starts a timer with specified tag.