Merge pull request #5856 from danieldouglas92/limit_timestep_by_darcy…

…_velocity

Add Darcy field convection timestep

doc/modules/changes/20240609_danieldouglas92

@@ -0,0 +1,4 @@
+Added: Functionality to limit the timestep based on the
+Darcy velocity.
+<br>
+(Daniel Douglas, 2024/06/09)

source/time_stepping/convection_time_step.cc

@@ -19,8 +19,9 @@
 */
 
 
-#include <aspect/global.h>
+#include <aspect/simulator.h>
 #include <aspect/time_stepping/convection_time_step.h>
+#include <aspect/melt.h>
 
 namespace aspect
 {
@@ -33,16 +34,40 @@ namespace aspect
       const QIterated<dim> quadrature_formula (QTrapezoid<1>(),
                                                this->get_parameters().stokes_velocity_degree);
 
+      bool consider_darcy_timestep = false;
+      unsigned int porosity_idx = numbers::invalid_unsigned_int;
+      if (this->introspection().composition_type_exists(CompositionalFieldDescription::porosity))
+        {
+          porosity_idx = this->introspection().find_composition_type(CompositionalFieldDescription::porosity);
+          if (this->get_parameters().compositional_field_methods[porosity_idx] == Parameters<dim>::AdvectionFieldMethod::fem_darcy_field)
+            consider_darcy_timestep = true;
+        }
+
+      const UpdateFlags update_flags
+        = UpdateFlags(
+            consider_darcy_timestep
+            ?
+            update_values |
+            update_gradients |
+            update_quadrature_points |
+            update_JxW_values
+            :
+            update_values);
+
       FEValues<dim> fe_values (this->get_mapping(),
                                this->get_fe(),
                                quadrature_formula,
-                               update_values);
+                               update_flags);
 
       const unsigned int n_q_points = quadrature_formula.size();
-
       std::vector<Tensor<1,dim>> velocity_values(n_q_points);
       std::vector<Tensor<1,dim>> fluid_velocity_values(n_q_points);
 
+      MaterialModel::MaterialModelInputs<dim> in(fe_values.n_quadrature_points, this->n_compositional_fields());
+      MaterialModel::MaterialModelOutputs<dim> out(fe_values.n_quadrature_points, this->n_compositional_fields());
+      MeltHandler<dim>::create_material_model_outputs(out);
+      MaterialModel::MeltOutputs<dim> *fluid_out = nullptr;
+
       double max_local_speed_over_meshsize = 0;
 
       for (const auto &cell : this->get_dof_handler().active_cell_iterators())
@@ -52,10 +77,32 @@ namespace aspect
             fe_values[this->introspection().extractors.velocities].get_function_values (this->get_solution(),
                                                                                         velocity_values);
 
+            if (consider_darcy_timestep == true)
+              {
+                in.reinit(fe_values, cell, this->introspection(), this->get_solution());
+                this->get_material_model().evaluate(in, out);
+                fluid_out = out.template get_additional_output<MaterialModel::MeltOutputs<dim>>();
+              }
+
             double max_local_velocity = 0;
             for (unsigned int q=0; q<n_q_points; ++q)
-              max_local_velocity = std::max (max_local_velocity,
-                                             velocity_values[q].norm());
+              {
+                if (consider_darcy_timestep)
+                  {
+                    const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(fe_values.quadrature_point(q));
+                    const double porosity = std::max(in.composition[q][porosity_idx], 1e-10);
+                    const double solid_density = out.densities[q];
+                    const double fluid_density = fluid_out->fluid_densities[q];
+                    const double fluid_viscosity = fluid_out->fluid_viscosities[q];
+                    const double permeability = fluid_out->permeabilities[q];
+                    const Tensor<1, dim> fluid_velocity = velocity_values[q] -
+                                                          (permeability / fluid_viscosity / porosity) *
+                                                          gravity * (solid_density - fluid_density);
+                    max_local_velocity = std::max(max_local_velocity, fluid_velocity.norm());
+                  }
+                max_local_velocity = std::max (max_local_velocity,
+                                               velocity_values[q].norm());
+              }
 
             if (this->get_parameters().include_melt_transport)
               {

tests/darcy_convection_step.prm

@@ -0,0 +1,155 @@
+#########################################################
+# This is a test for limiting the time step when advecting
+# a field with the 'darcy field' method without setting
+# Include melt transport = true. With the 'darcy field'
+# advection method, the field is advected using the Darcy
+# velocity, which is expressed as:
+#
+# u_f = u_s - K_D / phi * (rho_s * g - rho_f * g)
+# u_f = fluid velocity # m/yr
+# u_s = solid velocity # m/yr
+# K_D = Darcy Coefficient = k / eta_f
+# k = permeability # m2 / s
+# eta_f = fluid viscosity # Pa s
+# phi = porosity
+# rho_f = fluid density # kg/m3
+# rhos_s = solid density # kg/m3
+# g = gravity # m/s2
+#
+# Additionally, K_D and k can be expanded into:
+# K_D = k / eta_f
+# k = k_0 * phi**3 * (1 - phi)**2
+#
+# Where k_0 is the reference permeability. The test is a 10 km x 10 km 2D
+# box with a uniform porosity of 0.02 (2%) everywhere. The fluid has a
+# density of 1000 kg/m3, and a viscosity of 10 Pa s, and the solid has a
+# density of 3000 kg/m3 and a uniform velocity of 0 m/yr. The reference
+# pemeability is 1e-6, and gravity is set to -10. This results in a fluid
+# velocity of 24.25 m/yr. The mesh is 2.5 km x 2.5 km, and with CFL = 1,
+# this means that the time step should be 2500 m / 24.25 m/yr = 103.11
+
+############### Global parameters
+
+set Dimension                              = 2
+set Start time                             = 0
+set End time                               = 200
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme                = iterated Advection and Stokes
+set Max nonlinear iterations               = 1
+set CFL number                             = 1.0
+set Output directory                       = darcy_convection_step
+set Pressure normalization                 = surface
+set Surface pressure                       = 0
+set Maximum time step                      = 200
+set Use operator splitting                 = true
+
+# 10 km x 10 km box
+subsection Geometry model
+  set Model name = box
+  subsection Box
+    set X extent = 10e3
+    set Y extent = 10e3
+  end
+end
+
+# Uniform temperature of 293 K
+subsection Initial temperature model
+  set Model name = function
+  subsection Function
+    set Function expression = 293
+  end
+end
+
+subsection Boundary temperature model
+  set List of model names = box
+  set Fixed temperature boundary indicators   = top, bottom, left, right
+  subsection Box
+    set Bottom temperature = 293
+    set Top temperature    = 293
+    set Left temperature   = 293
+    set Right temperature  = 293
+  end
+end
+
+# Free slip boundary on all sides
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = left, right, top, bottom
+end
+
+# porosity and bound_fluid are required compositional fields when
+# using the reactive fluid transport. Set the porosity field method
+# to 'darcy field' so that the fluid is adveted with the Darcy
+# velocity.
+subsection Compositional fields
+  set Number of fields = 2
+  set Names of fields = porosity, bound_fluid
+  set Compositional field methods = darcy field, field
+end
+
+# Initialize a porosity of 2% (0.02) everywhere.
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Variable names      = x,y
+    set Function constants  =
+    set Function expression = 0.02; 0.0
+  end
+end
+
+# 10 m/s2 vertical gravity
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 10.0
+  end
+end
+
+# The reactive fluid transport model allows us to set the parameters which
+# influence fluid velocity, such as the fluid viscosity, fluid density, and
+# the reference permeability. We set the fluid weakening factors (shear to
+# bulk viscosity ratio, exponential fluid weakening factor) to 0 for
+# simplicity. We use the zero solubility fluid-solid reaction scheme to prevent
+# the fluid from partitioning into the solid.
+subsection Material model
+  set Model name = reactive fluid transport
+  subsection Reactive Fluid Transport Model
+    set Base model                                       = visco plastic
+    set Reference fluid density                          = 1000
+    set Shear to bulk viscosity ratio                    = 0.
+    set Reference fluid viscosity                        = 10
+    set Reference permeability                           = 1e-6
+    set Exponential fluid weakening factor               = 0
+    set Fluid-solid reaction scheme                      = zero solubility
+  end
+
+# Set the solid density to 3000 kg/m3, and set the minimum/maximum viscosity
+# to 1e21 Pa s for an isoviscous model.
+  subsection Visco Plastic
+    set Reference temperature                   = 1600
+    set Prefactors for diffusion creep          = 5e-21
+    set Viscous flow law                        = diffusion
+    set Densities                               = 3000
+    set Viscosity averaging scheme              = harmonic
+    set Minimum viscosity                       = 1e21
+    set Maximum viscosity                       = 1e21
+    set Thermal expansivities                   = 0
+  end
+end
+
+# Set the global refinement to 1, bringing the global mesh to 2.5 km x 2.5 km.
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 1
+  set Time steps between mesh refinement = 0
+end
+
+# Melt transport = false
+subsection Melt settings
+  set Include melt transport             = false
+end
+
+subsection Postprocess
+  set List of postprocessors =
+end

tests/darcy_convection_step/screen-output

@@ -0,0 +1,43 @@
+
+Number of active cells: 4 (on 2 levels)
+Number of degrees of freedom: 134 (50+9+25+25+25)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Skipping bound_fluid composition solve because RHS is zero.
+   Solving Stokes system... 7+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 6.2054e-17, 1.34414e-16, 0, 5.09696e-16
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 5.09696e-16
+
+
+   Postprocessing:
+
+*** Timestep 1:  t=103.11 years, dt=103.11 years
+   Solving composition reactions... in 1 substep(s).
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Skipping bound_fluid composition solve because RHS is zero.
+   Solving Stokes system... 6+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 3.51395e-16, 1.5726e-16, 0, 2.2185e-16
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 3.51395e-16
+
+
+   Postprocessing:
+
+*** Timestep 2:  t=200 years, dt=96.8895 years
+   Solving composition reactions... in 1 substep(s).
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Skipping bound_fluid composition solve because RHS is zero.
+   Solving Stokes system... 7+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.14854e-16, 1.32969e-16, 0, 3.48944e-16
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 3.48944e-16
+
+
+   Postprocessing:
+
+Termination requested by criterion: end time
+
+
+

tests/darcy_convection_step/statistics

@@ -0,0 +1,17 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Number of nonlinear iterations
+# 9: Iterations for temperature solver
+# 10: Iterations for composition solver 1
+# 11: Iterations for composition solver 2
+# 12: Iterations for Stokes solver
+# 13: Velocity iterations in Stokes preconditioner
+# 14: Schur complement iterations in Stokes preconditioner
+0 0.000000000000e+00 0.000000000000e+00 4 59 25 50 1 0 0 0 6 8 8 
+1 1.031104829590e+02 1.031104829590e+02 4 59 25 50 1 0 0 0 5 7 7 
+2 2.000000000000e+02 9.688951704104e+01 4 59 25 50 1 0 0 0 6 8 8 

tests/darcy_velocity/screen-output

@@ -24,26 +24,66 @@ Number of degrees of freedom: 44,070 (16,770+2,145+8,385+8,385+8,385)
      Compositions min/max/mass:  0.2/0.2/4e+09 // 0/1/4.872e+08
      Compositions max depth [m]: -1.798e+308 // 6.215e+04
 
-*** Timestep 1:  t=500 years, dt=500 years
+*** Timestep 1:  t=187.745 years, dt=187.745 years
    Solving temperature system... 0 iterations.
    Solving porosity system ... 0 iterations.
-   Solving fluid system ... 35 iterations.
+   Solving fluid system ... 17 iterations.
    Solving Stokes system... 0+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.81143e-16, 2.24127e-16, 0.407873, 3.71636e-09
-      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.407873
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.76026e-16, 1.62502e-16, 0.153152, 3.71636e-09
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.153152
 
    Solving temperature system... 0 iterations.
    Solving porosity system ... 0 iterations.
    Solving fluid system ... 0 iterations.
    Solving Stokes system... 0+0 iterations.
-      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.81143e-16, 2.24127e-16, 6.91442e-13, 3.71636e-09
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.76026e-16, 1.62502e-16, 8.478e-13, 3.71636e-09
       Relative nonlinear residual (total system) after nonlinear iteration 2: 3.71636e-09
 
 
    Postprocessing:
-     Compositions min/max/mass:  0.2/0.2/4e+09 // -0.1127/1.098/4.872e+08
+     Compositions min/max/mass:  0.2/0.2/4e+09 // -0.1357/1.131/4.872e+08
      Compositions max depth [m]: -1.798e+308 // 6.094e+04
 
+*** Timestep 2:  t=375.49 years, dt=187.745 years
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Solving fluid system ... 16 iterations.
+   Solving Stokes system... 0+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.01835e-16, 1.32239e-16, 0.0482778, 3.71636e-09
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.0482778
+
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Solving fluid system ... 0 iterations.
+   Solving Stokes system... 0+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 2.01835e-16, 1.32239e-16, 4.35337e-13, 3.71636e-09
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 3.71636e-09
+
+
+   Postprocessing:
+     Compositions min/max/mass:  0.2/0.2/4e+09 // -0.1747/1.142/4.872e+08
+     Compositions max depth [m]: -1.798e+308 // 6.094e+04
+
+*** Timestep 3:  t=500 years, dt=124.51 years
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Solving fluid system ... 14 iterations.
+   Solving Stokes system... 0+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.40836e-16, 2.12626e-16, 0.0286959, 3.71636e-09
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 0.0286959
+
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Solving fluid system ... 0 iterations.
+   Solving Stokes system... 0+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.40836e-16, 2.12487e-16, 4.06506e-13, 3.71636e-09
+      Relative nonlinear residual (total system) after nonlinear iteration 2: 3.71636e-09
+
+
+   Postprocessing:
+     Compositions min/max/mass:  0.2/0.2/4e+09 // -0.1528/1.137/4.872e+08
+     Compositions max depth [m]: -1.798e+308 // 5.973e+04
+
 Termination requested by criterion: end time
 
 

tests/darcy_velocity/statistics

@@ -21,4 +21,6 @@
 # 21: Max depth [m] for composition porosity
 # 22: Max depth [m] for composition fluid
 0 0.000000000000e+00 0.000000000000e+00 2048 18915 8385 16770 2 0 0  0          6 9 9 2.00000000e-01 2.00000000e-01 4.00000000e+09  0.00000000e+00 1.00000000e+00 4.87196181e+08 -1.79769313e+308 6.21478073e+04 
-1 5.000000000000e+02 5.000000000000e+02 2048 18915 8385 16770 2 0 0 35 4294967294 0 0 2.00000000e-01 2.00000000e-01 4.00000000e+09 -1.12670620e-01 1.09808535e+00 4.87198685e+08 -1.79769313e+308 6.09375000e+04 
+1 1.877448208391e+02 1.877448208391e+02 2048 18915 8385 16770 2 0 0 17 4294967294 0 0 2.00000000e-01 2.00000000e-01 4.00000000e+09 -1.35663795e-01 1.13147565e+00 4.87196208e+08 -1.79769313e+308 6.09375000e+04 
+2 3.754896416781e+02 1.877448208391e+02 2048 18915 8385 16770 2 0 0 16 4294967294 0 0 2.00000000e-01 2.00000000e-01 4.00000000e+09 -1.74738673e-01 1.14198327e+00 4.87196334e+08 -1.79769313e+308 6.09375000e+04 
+3 5.000000000000e+02 1.245103583219e+02 2048 18915 8385 16770 2 0 0 14 4294967294 0 0 2.00000000e-01 2.00000000e-01 4.00000000e+09 -1.52800550e-01 1.13661511e+00 4.87196548e+08 -1.79769313e+308 5.97271927e+04