A self-adaptive boundary velocity plugin ?

[XieJChris]

Hey guys,

In my current work, we assume a domain driven by left/right boundary, but the (horizontal) boundary velocity is required to change with this process to ensure a constant background strain rate, things like the opposite case of cookbooks/continental_extension. However, I am still quite new to ASPECT, so I refer to a few benchmarks (e.g., benchmarks/burstedde) and try to write such a plugin and test it on cookbooks/continental_extension, but the result (velocity field) seems to be weird. So I post the code in here and ask for some help.

Best, Xie

set Use years in output instead of seconds = true

subsection Boundary velocity model
  set Prescribed velocity boundary indicators = 0: ConstantStrainRateBoundaryVelocity, 1: ConstantStrainRateBoundaryVelocity
  set Zero velocity boundary indicators = 2 #bottom

  subsection Constant strain
    set The initial extents of the x-axis = 400e3
  end
end

subsection Material model
  set Model name = visco plastic

  subsection Visco plastic
    set Background strain rate = 5e-18
  end
end

template <int dim>
class ConstantStrainRateBoundaryVelocity : public BoundaryVelocity::Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:

     ConstantStrainRateBoundaryVelocity();

     void update() override;

     static void declare_parameters(ParameterHandler &prm);

     void parse_parameters(ParameterHandler &prm) override;

     Tensor<1, dim> boundary_velocity(const types::boundary_id boundary_indicator, const Point<dim> &position) const override;

     double get_boundary_velocity(const double strain_rate, double extents);

private:

     bool is_first_call; 
     double boundary_velocities_x; 
     double x_extents; // extent of the x-axis, initial value equals to get_extents()[0]
};

template <int dim>
ConstantStrainRateBoundaryVelocity<dim>::ConstantStrainRateBoundaryVelocity(): is_first_call(true) {}

template <int dim>
void ConstantStrainRateBoundaryVelocity<dim>::declare_parameters(ParameterHandler &prm)
{
   prm.enter_subsection("Boundary velocity model");
   {
      prm.enter_subsection("Constant strain");
      {
         prm.declare_entry("The initial extents of the x-axis","1100e3",Patterns::Double(),"Only horizontal velocity is considered here. Units:\\si{\\m}.");
      }
      prm.leave_subsection();    
   }
   prm.leave_subsection();
}

template <int dim>
void ConstantStrainRateBoundaryVelocity<dim>::parse_parameters(ParameterHandler &prm)
{
   prm.enter_subsection("Boundary velocity model");
   {
      prm.enter_subsection("Constant strain");
      {
         x_extents = prm.get_double("The initial extents of the x-axis");
      }
      prm.leave_subsection();    
   }
   prm.leave_subsection();
}

template <int dim>
double ConstantStrainRateBoundaryVelocity<dim>::get_boundary_velocity(const double strain_rate, double extents)
{
   double v_x = 0.0;

   if(this->convert_output_to_years())
   {
     v_x = strain_rate * extents / years_in_seconds;
   }
   else
   {
     v_x = strain_rate * extents;
   }
   return v_x;
}

template <>
Tensor<1,3>
ConstantStrainRateBoundaryVelocity<3>::boundary_velocity(const types::boundary_id, const Point<3> &) const
{
   Assert(false, ExcNotImplemented());

   return Tensor<1,3>();
}

template <>
Tensor<1,2>
ConstantStrainRateBoundaryVelocity<2>::boundary_velocity(const types::boundary_id, const Point<2> &) const
{
   Tensor<1,2> boundary_velocity_values;

    double v = 0.0;

    const ViscoPlastic<2> &vp_material_model = Plugins::get_plugin_as_type<const Viscoelastic<2>>(this->get_material_model());

    const double bg_sr = vp_material_model.get_background_strain_rate();// return background strain rate, read from prm file.

    if(boundary_indicator == 0)
   {
     if(this->convert_output_to_years())
     {
       v = bg_sr * x_extents / year_in_seconds;
     }
     else
     {
       v =  bg_sr * x_extents;
     }
   }
   else if(boundary_indicator == 1)
   {
     if(this->convert_output_to_years())
     {
       v = -1 * bg_sr * x_extents / year_in_seconds;
     }
     else
     {
       v =  -1 * bg_sr * x_extents;
     }
   }
   boundary_velocity_values[0] = v;

   boundary_velocity_values[1]  = 0.0;

   return boundary_velocity_values;
}

template <int dim>
void ConstantStrainRateBoundaryVelocity<dim>::update()
{
   const double dt = this->get_timestep();

   if(is_fist_call == true)
   {
     AssertThrow(Plugins::plugin_type_matches<const GeometryModel::Box<dim>>(this->get_geometry_model()),ExcMessage("This plugin is only implemented if the geometry is a box."));

     const ViscoPlastic<dim> &vp_material_model = Plugins::get_plugin_as_type<const Viscoelastic<dim>>(this->get_material_model());

     const double bg_strain_rate = vp_material_model.get_background_strain_rate(); 

     boundary_velocities_x = get_boundary_velocity(bg_strain_rate, x_extents);

     is_first_call = false;
   }
   x_extents -= 2 * boundary_velocities_x * dt;
}

[gassmoeller]

Hi Xie, your current description is missing some details that I would need to understand your problem better. In particular: What would you expect the velocity field to look like, and what does it do instead (what do you mean by weird)? Also you refer to the upper boundary in your message, but you do not prescribe any velocities at the upper boundary in the parameter file. What velocities do you expect at the upper boundary and what velocities do you see? Have you checked that the velocities at the left and right boundary (the ones you prescribe) are what you expect them to be? Also please post your complete parameter file, at the moment I do not quite understand what the upper boundary condition is (is it a free surface?).

[XieJChris]

Hi Rene,

Sorry for the missing details. My expected velocity filed is that the velocity at the left and right boundary is horizontal (the y-component is zero), the purpose of this is to create a compressed domain, and with the horizontal shortening, the velocity must be variable in order to make the background strain rate constant (see below). I think the similar case would be 'Brittle thrust wedges benchmark' or 'continental extension' provided in the manual, but the difference is that the former benchmarks set the boundary velocity constant. Therefore, the velocity filed I posted earlier is strange, because the direction of the boundary velocity becomes vertical.

As for the complete prm file, I haven't been in the lab recently, so I don't have the complete prm at present, I modify the continental_extension.prm, hope this still helpful.

set Dimension                              = 2
set Start time                             = 0
set End time                               = 5e6
set Use years in output instead of seconds = true
set Nonlinear solver scheme                = single Advection, iterated Stokes
set Nonlinear solver tolerance             = 1e-7
set Max nonlinear iterations               = 30
set CFL number                             = 0.5
set Output directory                       = output-continental_compression
set Timing output frequency                = 1
set Pressure normalization                 = no

subsection Solver parameters
  subsection Stokes solver parameters
    set Number of cheap Stokes solver steps = 200
    set Linear solver tolerance = 1e-7
  end
end

subsection Geometry model
  set Model name       = box
  subsection Box
    set X repetitions  = 200
    set Y repetitions  = 50
    set X extent       = 400e3
    set Y extent       = 100e3
  end
end

subsection Mesh refinement
  set Initial adaptive refinement        = 0
  set Initial global refinement          = 0
  set Time steps between mesh refinement = 0
end

subsection Mesh deformation
  set Mesh deformation boundary indicators  = top: free surface
  subsection Free surface
    set Surface velocity projection = vertical
  end
end

subsection Boundary velocity model
  set Prescribed velocity boundary indicators = 0: ConstantStrainRateBoundaryVelocity, 1: ConstantStrainRateBoundaryVelocity
  set Zero velocity boundary indicators = bottom

  subsection Constant strain
    set The initial extents of the x-axis = 400e3
  end
end

subsection Compositional fields
  set Number of fields = 7
  set Names of fields  = stress_xx, stress_yy, stress_xy, upper, lower, mantle, seed
end

subsection Initial composition model
  set Model name = function
  subsection Function
    set Variable names      = x,y
    set Function expression = 0;0;0;\
                              if(y>=80.e3, 1, 0); \
                              if(y<80.e3 && y>=70.e3, 1, 0); \
                              if(y<70.e3 && y>-100.e3, 1, 0); \ 
                              if(y<68.e3 && y>60.e3 && x>=198.e3 && x<=202.e3, 1, 0);
  end
end

subsection Boundary composition model
  set Fixed composition boundary indicators = bottom
  set List of model names = initial composition
end

subsection Boundary temperature model
  set Fixed temperature boundary indicators = bottom, top
  set List of model names  = box

  subsection Box
    set Bottom temperature = 1573
    set Top temperature    =  273
  end
end

subsection Initial temperature model
  set Model name = function
  subsection Function
    set Variable names = x,y
    set Function constants = h=100e3,ts1=273,ts2=681.5714,ts3=823., \
                                     k1=2.5,k2=2.5,k3=3.3,A=1.5e-6, \
                             qs1=0.0653571,qs2=0.035357,qs3=0.035357,qb3=0.035357
    set Function expression = if((h-y)<=20.e3, \
                                  ts1 + (qs1/k1)*(h-y) - (A*(h-y)*(h-y))/(2.0*k1), \
                                  if((h-y)>20.e3 && (h-y)<=30.e3, ts2 + (qs2/k2)*(h-y-20.e3),\
                                  ts3 + (qs3/k3)*(h-y-30.e3)));
  end
end

subsection Heating model
  set List of model names = compositional heating
  subsection Compositional heating
    set Compositional heating values = 0., 0., 0., 0., 1.5e-6, 0., 0., 0.
  end
end

subsection Formulation
  set Enable elasticity = true
end

subsection Material model
  set Model name = visco plastic

  subsection Visco Plastic

    set Include viscoelasticity  = true
    set Reference temperature    = 293
    set Reference viscosity      = 1e22

    set Minimum strain rate      = 1.e-20
    set Reference strain rate    = 1.e-16

    set Minimum viscosity        = 1e18
    set Maximum viscosity        = 1e26

    set Thermal diffusivities    = 1.333333e-6, 1.333333e-6, 1.333333e-6, 1.333333e-6, 1.190476e-6, 1.149425e-6, 1.333333e-6, 1.333333e-6
    set Heat capacities          = 750.
    set Densities                = 3300, 3300, 3300, 3300, 2800, 2900, 3300, 3300
    set Thermal expansivities    = 2e-5

    set Viscosity averaging scheme = harmonic

    set Viscous flow law = dislocation

    set Prefactors for dislocation creep          = 6.52e-16, 6.52e-16, 6.52e-16, 6.52e-16, 8.57e-28, 7.13e-18, 6.52e-16, 7.13e-18
    set Stress exponents for dislocation creep    = 3.5, 3.5, 3.5, 3.5, 4.0, 3.0, 3.5, 3.0
    set Activation energies for dislocation creep = 530.e3, 530.e3, 530.e3, 530.e3, 223.e3, 345.e3, 530.e3, 345.e3
    set Activation volumes for dislocation creep  = 18.e-6, 18.e-6, 18.e-6, 18.e-6, 0., 0., 18.e-6, 0.

    set Elastic shear moduli                      = 1.e10
    set Use fixed elastic time step               = false
    set Fixed elastic time step                   = 1e3
    set Use stress averaging                      = false

    set Yield mechanism                           = drucker
    set Angles of internal friction               = 30
    set Cohesions                                 = 1.0e7
    set Maximum yield stress                      = 1e12
    set Use plastic damper                        = false

  end
end

subsection Gravity model
  set Model name = vertical

  subsection Vertical
    set Magnitude = 9.8
  end
end

subsection Postprocess
  set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization
  subsection Visualization
    set List of output variables = density, viscosity, strain rate, error indicator
    set Output format = vtu
    set Time between graphical output = 1e4
  end
end

Best,

Xie