hm_iterative.cc

Go to the documentation of this file.

/*!
 *
 * Copyright (C) 2015 Technical University of Liberec.  All rights reserved.
 * 
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License version 3 as published by the
 * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
 * 
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
 *
 * 
 * @file    hm_iterative.cc
 * @brief   
 * @author  Jan Stebel
 */
 
#include "hm_iterative.hh"
#include "system/sys_profiler.hh"
#include "input/input_type.hh"
#include "flow/richards_lmh.hh"
#include "fields/field_fe.hh"         // for create_field_fe()
#include "fields/field_model.hh"      // for Model
 
 
FLOW123D_FORCE_LINK_IN_CHILD(coupling_iterative)
 
 
namespace it = Input::Type;
 
 
struct fn_pressure_potential {
    inline double operator() (double alpha, double density, double gravity, double pressure)
    {
        return -alpha*density*gravity*pressure;
    }
};
 
 
struct fn_hm_coupling_beta {
 
    fn_hm_coupling_beta(double beta_f) : beta_factor(beta_f) {}
 
 
    inline double operator() (double alpha, double lame_mu, double lame_lambda)
    {
        return beta_factor*0.5*alpha*alpha/(2*lame_mu/3 + lame_lambda);
    }
 
private:
 
    const double beta_factor; ///< User-defined factor for iteration parameter.
 
};
 
 
struct fn_fluid_source {
 
    fn_fluid_source(TimeGovernor *time) : time_(time) {}
 
    inline double operator() (double alpha, double beta, double pressure, double old_pressure, double div_u, double old_div_u)
    {
        return (beta*(pressure-old_pressure) + alpha*(old_div_u - div_u)) / time_->dt();
    }
 
private:
 
    const TimeGovernor *time_;
};
 
 
const it::Record & HM_Iterative::get_input_type() {
    return it::Record("Coupling_Iterative",
            "Record with data for iterative coupling of flow and mechanics.\n")
        .derive_from( DarcyFlowInterface::get_input_type() )
        .copy_keys(EquationBase::record_template())
        .copy_keys(IterativeCoupling::record_template())
        .declare_key("flow_equation", RichardsLMH::get_input_type(),
                it::Default::obligatory(),
                "Flow equation, provides the velocity field as a result.")
        .declare_key("mechanics_equation", Elasticity::get_input_type(),
                "Mechanics, provides the displacement field.")
        .declare_key("input_fields", it::Array(
                HM_Iterative::EqFields()
                    .make_field_descriptor_type("Coupling_Iterative")),
                IT::Default::obligatory(),
                "Input fields of the HM coupling.")
        .declare_key( "iteration_parameter", it::Double(), it::Default("1"),
                "Tuning parameter for iterative splitting. Its default value "
                "corresponds to a theoretically optimal value with fastest convergence." )
        .declare_key( "max_it", it::Integer(0), it::Default("100"),
                "Maximal count of HM iterations." )
        .declare_key( "min_it", it::Integer(0), it::Default("1"),
                "Minimal count of HM iterations." )
        .declare_key( "a_tol", it::Double(0), it::Default("0"),
                "Absolute tolerance for difference in HM iteration." )
        .declare_key( "r_tol", it::Double(0), it::Default("1e-7"),
                "Relative tolerance for difference in HM iteration." )
        .close();
}
 
 
const int HM_Iterative::registrar = Input::register_class< HM_Iterative, Mesh &, const Input::Record >("Coupling_Iterative")
                                    + HM_Iterative::get_input_type().size();
 
 
HM_Iterative::EqFields::EqFields()
{
    *this += alpha.name("biot_alpha")
                     .description("Biot poroelastic coefficient.")
                     .units(UnitSI().dimensionless())
                     .input_default("0.0")
                     .flags_add(FieldFlag::in_rhs);
    
    *this += density.name("fluid_density")
                     .description("Volumetric mass density of the fluid.")
                     .units(UnitSI().kg().m(-3))
                     .input_default("0.0")
                     .flags_add(FieldFlag::in_rhs);
    
    *this += gravity.name("gravity")
                     .description("Gravitational acceleration constant.")
                     .units(UnitSI().m().s(-2))
                     .input_default("9.81")
                     .flags_add(FieldFlag::in_rhs);
    
    *this += beta.name("relaxation_beta")
                     .description("Parameter of numerical method for iterative solution of hydro-mechanical coupling.")
                     .units(UnitSI().dimensionless())
                     .flags(FieldFlag::equation_external_output);
    
    *this += pressure_potential.name("pressure_potential")
                     .description("Coupling term entering the mechanics equation.")
                     .units(UnitSI().m())
                     .flags(FieldFlag::equation_result);
    
    *this += old_pressure.name("old_pressure")
                     .description("Pressure from last computed time.")
                     .units(UnitSI().m())
                     .flags(FieldFlag::equation_external_output);
 
    *this += old_div_u.name("old_displacement_divergence")
                     .description("Displacement divergence from last computed time.")
                     .units(UnitSI().dimensionless())
                     .flags(FieldFlag::equation_external_output);
 
    *this += ref_pressure_potential.name("ref_pressure_potential")
                     .description("Pressure potential on boundary (taking into account the flow boundary condition.")
                     .units(UnitSI().m())
                     .flags(FieldFlag::equation_result);
    
    *this += flow_source.name("extra_flow_source")
                     .description("Coupling term entering the flow equation.")
                     .units(UnitSI().s(-1))
                     .flags(FieldFlag::equation_result);
}
 
 
void HM_Iterative::EqFields::initialize(Mesh &mesh, HM_Iterative &parent)
{
    // initialize coupling fields with FieldFE
    set_mesh(mesh);
    
    ref_potential_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(mesh, MixedPtr<FE_CR>());
    ref_pressure_potential.set(ref_potential_ptr_, 0.0);
    
    old_pressure_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(parent.flow_->eq_data().dh_cr_, &parent.flow_->eq_data().p_edge_solution_previous_time);
    old_iter_pressure_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(mesh, MixedPtr<FE_P_disc>(0));
    old_div_u_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(parent.mechanics_->eq_fields().output_div_ptr->get_dofhandler());
}
 
                                    
 
HM_Iterative::HM_Iterative(Mesh &mesh, Input::Record in_record)
: DarcyFlowInterface(mesh, in_record),
  IterativeCoupling(in_record)
{
    START_TIMER("HM constructor");
    using namespace Input;
 
    time_ = new TimeGovernor(in_record.val<Record>("time"));
    ASSERT( time_->is_default() == false ).error("Missing key 'time' in Coupling_Iterative.");
    
    // setup flow equation
    Record flow_rec = in_record.val<Record>("flow_equation");
    // Need explicit template types here, since reference is used (automatically passing by value)
    flow_ = std::make_shared<RichardsLMH>(*mesh_, flow_rec, time_);
    
    // setup mechanics
    Record mech_rec = in_record.val<Record>("mechanics_equation");
    mechanics_ = std::make_shared<Elasticity>(*mesh_, mech_rec, this->time_);
    mechanics_->initialize();
    
    // setup coupling fields and finish initialization of flow
    mechanics_->eq_fields()["cross_section"].copy_from(flow_->eq_fields()["cross_section"]);
    // flow_->eq_fields()["cross_section_updated"].copy_from(mechanics_->eq_fields()["cross_section_updated"]);
    // flow_->eq_fields()["stress"].copy_from(mechanics_->eq_fields()["stress"]);
    // flow_->eq_fields()["von_mises_stress"].copy_from(mechanics_->eq_fields()["von_mises_stress"]);
 
    flow_->eq_fields() += mechanics_->eq_fields()["cross_section_updated"];
    flow_->eq_fields() += mechanics_->eq_fields()["stress"];
    flow_->eq_fields() += mechanics_->eq_fields()["von_mises_stress"];
    flow_->eq_fields() += mechanics_->eq_fields()["mean_stress"];
    flow_->initialize();
    std::stringstream ss; // print warning message with table of uninitialized fields
    if ( FieldCommon::print_message_table(ss, "flow") )
        WarningOut() << ss.str();
 
    eq_fields_ += *flow_->eq_fields().field("pressure_edge");
 
    this->eq_fieldset_ = &eq_fields_;
    
    // setup input fields
    eq_fields_.set_input_list( in_record.val<Input::Array>("input_fields"), time() );
 
    eq_fields_.initialize(*mesh_, *this);
    mechanics_->set_potential_load(eq_fields_.pressure_potential, eq_fields_.ref_pressure_potential);
}
 
 
void HM_Iterative::initialize()
{
    eq_fields_.old_pressure.set(eq_fields_.old_pressure_ptr_, 0.0);
    eq_fields_.old_div_u.set(eq_fields_.old_div_u_ptr_, 0.0);
 
    eq_fields_.pressure_potential.set(Model<3, FieldValue<3>::Scalar>::create(
        fn_pressure_potential(),
        eq_fields_.alpha,
        eq_fields_.density,
        eq_fields_.gravity,
        flow_->eq_fields().field_edge_pressure
        ), 0.0);
    
    eq_fields_.beta.set(Model<3, FieldValue<3>::Scalar>::create(
        fn_hm_coupling_beta(input_record_.val<double>("iteration_parameter")),
        eq_fields_.alpha,
        mechanics_->eq_fields().lame_mu,
        mechanics_->eq_fields().lame_lambda
        ), 0.0);
 
    eq_fields_.flow_source.set(Model<3, FieldValue<3>::Scalar>::create(
        fn_fluid_source(time_),
        eq_fields_.alpha,
        eq_fields_.beta,
        flow_->eq_fields().field_edge_pressure,
        eq_fields_.old_pressure,
        mechanics_->eq_fields().output_divergence,
        eq_fields_.old_div_u
        ), 0.0);
}
 
 
template<int dim, class Value>
void copy_field(const FieldCommon &from_field_common, FieldFE<dim, Value> &to_field)
{
    auto dh = to_field.get_dofhandler();
    auto vec = to_field.vec();
    Field<dim,Value> from_field;
    from_field.copy_from(from_field_common);
    
    for ( auto cell : dh->own_range() )
        vec.set( cell.local_idx(), from_field.value(cell.elm().centre(), cell.elm()) );
}
 
 
 
void HM_Iterative::zero_time_step()
{
    eq_fields_.set_time(time_->step(), LimitSide::right);
    std::stringstream ss;
    if ( FieldCommon::print_message_table(ss, "coupling_iterative") )
        WarningOut() << ss.str();
    
    mechanics_->update_output_fields(); // init field values for use in flow
    flow_->zero_time_step();
    update_potential();
    mechanics_->zero_time_step();
    
    copy_field(*flow_->eq_fields().field("pressure_p0"), *eq_fields_.old_iter_pressure_ptr_);
    copy_field(mechanics_->eq_fields().output_divergence, *eq_fields_.old_div_u_ptr_);
}
 
 
void HM_Iterative::update_solution()
{
    time_->next_time();
    time_->view("HM");
    eq_fields_.set_time(time_->step(), LimitSide::right);
 
    solve_step();
}
 
void HM_Iterative::solve_iteration()
{
    // pass displacement (divergence) to flow
    // and solve flow problem
    update_flow_fields();
    flow_->solve_time_step(false);
    
    // pass pressure to mechanics and solve mechanics
    update_potential();
    mechanics_->solve_linear_system();
}
 
 
void HM_Iterative::update_after_iteration()
{
    mechanics_->update_output_fields();
    copy_field(*flow_->eq_fields().field("pressure_p0"), *eq_fields_.old_iter_pressure_ptr_);
}
 
 
void HM_Iterative::update_after_converged()
{
    flow_->accept_time_step();
    flow_->output_data();
    mechanics_->output_data();
    
    copy_field(mechanics_->eq_fields().output_divergence, *eq_fields_.old_div_u_ptr_);
}
 
 
void HM_Iterative::update_potential()
{
    auto ref_potential_vec_ = eq_fields_.ref_potential_ptr_->vec();
    auto dh = eq_fields_.ref_potential_ptr_->get_dofhandler();
    Field<3, FieldValue<3>::Scalar> field_edge_pressure;
    field_edge_pressure.copy_from(*flow_->eq_fields().field("pressure_edge"));
 
    for ( auto ele : dh->local_range() )
    {
        auto elm = ele.elm();
        LocDofVec dof_indices = ele.get_loc_dof_indices();
        for ( auto side : ele.side_range() )
        {
            double alpha = eq_fields_.alpha.value(side.centre(), elm);
            double density = eq_fields_.density.value(side.centre(), elm);
            double gravity = eq_fields_.gravity.value(side.centre(), elm);
 
            // The reference potential is applied only on dirichlet and total_flux b.c.,
            // i.e. where only mechanical traction is prescribed.
            if (side.side().is_boundary() &&
                    (flow_->eq_fields().bc_type.value(side.centre(), side.cond().element_accessor()) == DarcyMH::EqFields::dirichlet ||
                    flow_->eq_fields().bc_type.value(side.centre(), side.cond().element_accessor()) == DarcyMH::EqFields::total_flux)
                )
            {
                double bc_pressure = flow_->eq_fields().bc_pressure.value(side.centre(), side.cond().element_accessor());
                ref_potential_vec_.set(dof_indices[side.side_idx()], -alpha*density*gravity*bc_pressure);
            }
            else
                ref_potential_vec_.set(dof_indices[side.side_idx()], 0);
        }
    }
    
    ref_potential_vec_.local_to_ghost_begin();
    ref_potential_vec_.local_to_ghost_end();
    eq_fields_.pressure_potential.set_time_result_changed();
    eq_fields_.ref_pressure_potential.set_time_result_changed();
    mechanics_->set_potential_load(eq_fields_.pressure_potential, eq_fields_.ref_pressure_potential);
}
 
 
void HM_Iterative::update_flow_fields()
{
    eq_fields_.beta.set_time_result_changed();
    eq_fields_.flow_source.set_time_result_changed();
    flow_->set_extra_storativity(eq_fields_.beta);
    flow_->set_extra_source(eq_fields_.flow_source);
}
 
 
void HM_Iterative::compute_iteration_error(double& abs_error, double& rel_error)
{
    auto dh = eq_fields_.old_iter_pressure_ptr_->get_dofhandler();
    double p_dif2 = 0, p_norm2 = 0;
    Field<3,FieldValue<3>::Scalar> field_ele_pressure;
    field_ele_pressure.copy_from(*flow_->eq_fields().field("pressure_p0"));
    for (auto cell : dh->own_range())
    {
        auto elm = cell.elm();
        double new_p = field_ele_pressure.value(elm.centre(), elm);
        double old_p = eq_fields_.old_iter_pressure_ptr_->value(elm.centre(), elm);
        p_dif2 += pow(new_p - old_p, 2)*elm.measure();
        p_norm2 += pow(old_p, 2)*elm.measure();
    }
    
    double send_data[] = { p_dif2, p_norm2 };
    double recv_data[2];
    MPI_Allreduce(&send_data, &recv_data, 2, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
    abs_error = sqrt(recv_data[0]);
    rel_error = abs_error / (sqrt(recv_data[1]) + std::numeric_limits<double>::min());
    
    MessageOut().fmt("HM Iteration {} abs. difference: {}  rel. difference: {}\n"
                         "--------------------------------------------------------",
                         iteration(), abs_error, rel_error);
 
    if(iteration() >= max_it_ && (abs_error > a_tol_ || rel_error > r_tol_))
        MessageOut().fmt("HM solver did not converge in {} iterations.\n", iteration());
}
 
 
 
HM_Iterative::~HM_Iterative() {
    flow_.reset();
    mechanics_.reset();
}