From d1cb2c6e6ca4f7ea62f3bf0e348584c0f40bd3de Mon Sep 17 00:00:00 2001 From: danieldouglas92 Date: Sat, 8 Jun 2024 21:18:51 -0600 Subject: [PATCH] Add darcy field convection step --- doc/modules/changes/20240609_danieldouglas92 | 4 + source/time_stepping/convection_time_step.cc | 57 ++++++- tests/darcy_convection_step.prm | 155 +++++++++++++++++++ tests/darcy_convection_step/screen-output | 43 +++++ tests/darcy_convection_step/statistics | 17 ++ tests/darcy_velocity/screen-output | 52 ++++++- tests/darcy_velocity/statistics | 4 +- 7 files changed, 320 insertions(+), 12 deletions(-) create mode 100644 doc/modules/changes/20240609_danieldouglas92 create mode 100644 tests/darcy_convection_step.prm create mode 100644 tests/darcy_convection_step/screen-output create mode 100644 tests/darcy_convection_step/statistics diff --git a/doc/modules/changes/20240609_danieldouglas92 b/doc/modules/changes/20240609_danieldouglas92 new file mode 100644 index 00000000000..cbbd82231cb --- /dev/null +++ b/doc/modules/changes/20240609_danieldouglas92 @@ -0,0 +1,4 @@ +Added: Functionality to limit the timestep based on the +Darcy velocity. +
+(Daniel Douglas, 2024/06/09) diff --git a/source/time_stepping/convection_time_step.cc b/source/time_stepping/convection_time_step.cc index 9460a30d9b3..14c42d536c7 100644 --- a/source/time_stepping/convection_time_step.cc +++ b/source/time_stepping/convection_time_step.cc @@ -19,8 +19,9 @@ */ -#include +#include #include +#include namespace aspect { @@ -33,16 +34,40 @@ namespace aspect const QIterated 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::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 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> velocity_values(n_q_points); std::vector> fluid_velocity_values(n_q_points); + MaterialModel::MaterialModelInputs in(fe_values.n_quadrature_points, this->n_compositional_fields()); + MaterialModel::MaterialModelOutputs out(fe_values.n_quadrature_points, this->n_compositional_fields()); + MeltHandler::create_material_model_outputs(out); + MaterialModel::MeltOutputs *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>(); + } + double max_local_velocity = 0; for (unsigned int q=0; q 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) { diff --git a/tests/darcy_convection_step.prm b/tests/darcy_convection_step.prm new file mode 100644 index 00000000000..5a08c5e6b13 --- /dev/null +++ b/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 diff --git a/tests/darcy_convection_step/screen-output b/tests/darcy_convection_step/screen-output new file mode 100644 index 00000000000..6d6be06c379 --- /dev/null +++ b/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 + + + diff --git a/tests/darcy_convection_step/statistics b/tests/darcy_convection_step/statistics new file mode 100644 index 00000000000..8dfb07706f3 --- /dev/null +++ b/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 diff --git a/tests/darcy_velocity/screen-output b/tests/darcy_velocity/screen-output index 6a775b35d2f..3421897d6bc 100644 --- a/tests/darcy_velocity/screen-output +++ b/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 diff --git a/tests/darcy_velocity/statistics b/tests/darcy_velocity/statistics index 2927a55aaaf..0bae53133a0 100644 --- a/tests/darcy_velocity/statistics +++ b/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