Skip to content

Commit

Permalink
add test case
Browse files Browse the repository at this point in the history
  • Loading branch information
@danieldouglas92
danieldouglas92 committed Nov 4, 2024
1 parent 30de055 commit 59be152
Show file tree
Hide file tree
Showing 9 changed files with 3,384 additions and 393 deletions.
390 changes: 0 additions & 390 deletions source/material_model/reaction_model/tian2019_parameterized_reactions.cc
View file

This file was deleted.

114 changes: 114 additions & 0 deletions source/material_model/reactive_fluid_transport.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,80 @@ namespace aspect
}
}



template <int dim>
std::vector<double>
ReactiveFluidTransport<dim>::
tian_equilibrium_bound_water_content (const MaterialModel::MaterialModelInputs<dim> &in,
unsigned int q) const
{
// Create arrays that will store the values of the polynomials at the current pressure
std::vector<double> LR_values(4);
std::vector<double> csat_values(4);
std::vector<double> Td_values(4);

// Loop over the four rock types (peridotite, gabbro, MORB, sediment) and the polynomial
// coefficients to fill the vectors defined above. The polynomials for LR are defined in
// equations 13, B2, B10, and B18. csat polynomials are defined in equations 14, B1, B9, and B17.
// Td polynomials are defined in equations 15, B3, B11, and B19.
for (unsigned int i = 0; i<devolatilization_enthalpy_changes.size(); ++i)
{
// Pressure, which must be in GPa for the parametrization, or GPa^-1. The polynomials for each lithology
// breaks down above certain pressures, make sure that we cap the pressure just before this break down.
// Introduce minimum pressure to avoid a division by 0.
const double minimum_pressure = 1e-12;
const double pressure = std::min(std::max(minimum_pressure, in.pressure[q]/1.e9), pressure_cutoffs[i]);
const double inverse_pressure = 1.0/pressure;
for (unsigned int j = 0; j<devolatilization_enthalpy_changes[i].size(); ++j)
{
#if DEAL_II_VERSION_GTE(9, 6, 0)
LR_values[i] += devolatilization_enthalpy_changes[i][j] * Utilities::pow(inverse_pressure, devolatilization_enthalpy_changes[i].size() - 1 - j);
#else
LR_values[i] += devolatilization_enthalpy_changes[i][j] * std::pow(inverse_pressure, devolatilization_enthalpy_changes[i].size() - 1 - j);
#endif
}

for (unsigned int j = 0; j<water_mass_fractions[i].size(); ++j)
{
#if DEAL_II_VERSION_GTE(9, 6, 0)
csat_values[i] += i==3 ? water_mass_fractions[i][j] * Utilities::pow(std::log10(pressure), water_mass_fractions[i].size() - 1 - j) :\
water_mass_fractions[i][j] * Utilities::pow(pressure, water_mass_fractions[i].size() - 1 - j);
#else
csat_values[i] += i==3 ? water_mass_fractions[i][j] * std::pow(std::log10(pressure), water_mass_fractions[i].size() - 1 - j) :\
water_mass_fractions[i][j] * std::pow(pressure, water_mass_fractions[i].size() - 1 - j);
#endif
}

for (unsigned int j = 0; j<devolatilization_onset_temperatures[i].size(); ++j)
{
#if DEAL_II_VERSION_GTE(9, 6, 0)
Td_values[i] += devolatilization_onset_temperatures[i][j] * Utilities::pow(pressure, devolatilization_onset_temperatures[i].size() - 1 - j);
#else
Td_values[i] += devolatilization_onset_temperatures[i][j] * std::pow(pressure, devolatilization_onset_temperatures[i].size() - 1 - j);
#endif
}
}

// Create an array for the equilibrium bound water content that is calculated from these polynomials
std::vector<double> eq_bound_water_content(4);

// Define the maximum bound water content allowed for the four different rock compositions
std::vector<double> max_bound_water_content = {tian_max_peridotite_water, tian_max_gabbro_water, tian_max_MORB_water, tian_max_sediment_water};

// Loop over all rock compositions and fill the equilibrium bound water content, divide by 100 to convert
// from percentage to fraction (equation 1)
for (unsigned int k = 0; k<LR_values.size(); ++k)
{
eq_bound_water_content[k] = (std::min(std::exp(csat_values[k]) * \
std::exp(std::exp(LR_values[k]) * (1/in.temperature[q] - 1/Td_values[k])), \
max_bound_water_content[k]) / 100.0);
}
return eq_bound_water_content;
}



template <int dim>
void
ReactiveFluidTransport<dim>::
Expand Down Expand Up @@ -86,6 +160,29 @@ namespace aspect
}
case tian_approximation:
{
// The bound fluid content is calculated using parametrized phase
// diagrams for four different rock types: sediment, MORB, gabbro, and
// peridotite.
const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
const unsigned int sediment_idx = this->introspection().compositional_index_for_name("sediment");
const unsigned int MORB_idx = this->introspection().compositional_index_for_name("MORB");
const unsigned int gabbro_idx = this->introspection().compositional_index_for_name("gabbro");
const unsigned int peridotite_idx = this->introspection().compositional_index_for_name("peridotite");

// Initialize a vector that stores the compositions (mass fractions) for
// the four different rock compositions,
std::vector<double> tracked_rock_mass_fractions(4);
tracked_rock_mass_fractions[0] = (in.composition[q][peridotite_idx]);
tracked_rock_mass_fractions[1] = (in.composition[q][gabbro_idx]);
tracked_rock_mass_fractions[2] = (in.composition[q][MORB_idx]);
tracked_rock_mass_fractions[3] = (in.composition[q][sediment_idx]);

// The bound water content (water within the solid phase) for the four different rock types
std::vector<double> tian_eq_bound_water_content = tian_equilibrium_bound_water_content(in, q);

// average the water content between the four different rock types
double average_eq_bound_water_content = MaterialUtilities::average_value (tracked_rock_mass_fractions, tian_eq_bound_water_content, MaterialUtilities::arithmetic);

// The fluid volume fraction in equilibrium with the solid (stored in the melt_fractions vector)
// is equal to the sum of the porosity and the change in bound fluid content
// (current bound fluid - updated average bound fluid).
Expand Down Expand Up @@ -294,6 +391,18 @@ namespace aspect
"computed. If the model does not use operator splitting, this parameter is not used. "
"Units: yr or s, depending on the ``Use years "
"in output instead of seconds'' parameter.");
prm.declare_entry ("Maximum weight percent water in sediment", "3",
Patterns::Double (0),
"The maximum allowed weight percent that the sediment composition can hold.");
prm.declare_entry ("Maximum weight percent water in MORB", "2",
Patterns::Double (0),
"The maximum allowed weight percent that the sediment composition can hold.");
prm.declare_entry ("Maximum weight percent water in gabbro", "1",
Patterns::Double (0),
"The maximum allowed weight percent that the sediment composition can hold.");
prm.declare_entry ("Maximum weight percent water in peridotite", "8",
Patterns::Double (0),
"The maximum allowed weight percent that the sediment composition can hold.");
prm.declare_entry ("Fluid-solid reaction scheme", "no reaction",
Patterns::Selection("no reaction|zero solubility|tian approximation|katz2003"),
"Select what type of scheme to use for reactions between fluid and solid phases. "
Expand Down Expand Up @@ -340,6 +449,11 @@ namespace aspect
fluid_reaction_time_scale = prm.get_double ("Fluid reaction time scale for operator splitting");
reference_T = prm.get_double ("Reference temperature");

tian_max_peridotite_water = prm.get_double ("Maximum weight percent water in peridotite");
tian_max_gabbro_water = prm.get_double ("Maximum weight percent water in gabbro");
tian_max_MORB_water = prm.get_double ("Maximum weight percent water in MORB");
tian_max_sediment_water = prm.get_double ("Maximum weight percent water in sediment");

// Create the base model and initialize its SimulatorAccess base
// class; it will get a chance to read its parameters below after we
// leave the current section.
Expand Down
7 changes: 4 additions & 3 deletions source/simulator/helper_functions.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2331,7 +2331,7 @@ namespace aspect
MaterialModel::MaterialModelInputs<dim> in(fe_face_values.n_quadrature_points, introspection.n_compositional_fields);
MaterialModel::MaterialModelOutputs<dim> out(fe_face_values.n_quadrature_points, introspection.n_compositional_fields);
MeltHandler<dim>::create_material_model_outputs(out);
MaterialModel::MeltOutputs<dim> *fluid_out = out.template get_additional_output<MaterialModel::MeltOutputs<dim>>();
MaterialModel::MeltOutputs<dim> *fluid_out = nullptr;

const auto &tangential_velocity_boundaries =
boundary_velocity_manager.get_tangential_boundary_velocity_indicators();
Expand Down Expand Up @@ -2366,12 +2366,13 @@ namespace aspect
{
in.reinit(fe_face_values, cell, introspection, solution);
material_model->evaluate(in, out);
fluid_out = out.template get_additional_output<MaterialModel::MeltOutputs<dim>>();
}

if (consider_fluid_velocity)
{
const FEValuesExtractors::Vector ex_u_f = introspection.variable("fluid velocity").extractor_vector();
fe_face_values[ex_u_f].get_function_values(current_linearization_point, face_current_fluid_velocity_values);
fe_face_values[introspection.variable("fluid velocity").extractor_vector()].get_function_values(current_linearization_point,
face_current_fluid_velocity_values);
}

// ... check if the face is an outflow boundary by integrating the normal velocities
Expand Down
188 changes: 188 additions & 0 deletions tests/darcy_outflow.prm
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,188 @@
# This test defines an initial blob of porosity implaced in a solid within a 2D cartesian box.
# Because the porosity blob is less dense, the blob will rise to the surface of the model and
# flow out of the top of the box. This tests the implementation of composition dependent outflow
# boundary conditions for when a composition is advected with the Darcy field advection method.
set Dimension = 2
set Use years in output instead of seconds = true
set End time = 10e3
set Pressure normalization = surface
set Maximum time step = 100
set CFL number = 0.5
set Surface pressure = 0
set Use operator splitting = true
set Output directory = output-darcy
set Nonlinear solver scheme = iterated Advection and Stokes
set Max nonlinear iterations = 3

# Define a function which
subsection Boundary velocity model
set Prescribed velocity boundary indicators = top:function, bottom:function, right:function, left:function

subsection Function
set Variable names = x,y,
set Function expression = 0;-0.1
end
end

subsection Gravity model
set Model name = vertical

subsection Vertical
set Magnitude = 10
end
end

subsection Formulation
set Formulation = Boussinesq approximation
end

subsection Solver parameters
set Temperature solver tolerance = 1e-10
subsection Stokes solver parameters
set Stokes solver type = block GMG
set GMRES solver restart length = 1000
set Number of cheap Stokes solver steps = 20000
set Use full A block as preconditioner = true
set Linear solver tolerance = 1e-7
set Maximum number of expensive Stokes solver steps = 0
end
end

# 10 km x 10 km box
subsection Geometry model
set Model name = box

subsection Box
set X extent = 7.5e3
set Y extent = 7.5e3

set X repetitions = 10
set Y repetitions = 10
end
end

# Uniform temperature of 293 K
subsection Initial temperature model
set Model name = function

subsection Function
set Function expression = 293
end
end

subsection Adiabatic conditions model
set Model name = function

subsection Function
set Variable names = z
set Function expression = 293; 3.3e4*z; 3300
end
end

subsection Boundary temperature model
set Allow fixed temperature on outflow boundaries = false
set List of model names = box
set Fixed temperature boundary indicators = top, bottom, right, left

subsection Box
set Bottom temperature = 293
set Top temperature = 293
set Left temperature = 293
set Right temperature = 293
end
end

subsection Compositional fields
set Number of fields = 3
set Names of fields = porosity, bound_fluid, solid_phase
set Compositional field methods = darcy field, field, field
end

subsection Melt settings
set Include melt transport = false
end

# Initialize a porosity of 1% (0.01) everywhere.
subsection Initial composition model
set Model name = function

subsection Function
set Variable names = x,y,t
set Function constants = pi=3.1415926, x0=3750, y0=3750, c=1000
set Function expression = 0.01 * exp(-((x-x0)*(x-x0)+(y-y0)*(y-y0))/(2*c*c)); 0.0; if(y>=7000, 1, 0)
end
end

subsection Boundary composition model
set Fixed composition boundary indicators = bottom, top
set Allow fixed composition on outflow boundaries = false
set List of model names = initial composition
end

subsection Material model
set Model name = reactive fluid transport
set Material averaging = geometric average only viscosity

subsection Reactive Fluid Transport Model
set Base model = visco plastic
set Reference fluid density = 1000 # density of water
set Shear to bulk viscosity ratio = 1
set Reference fluid viscosity = 10
set Reference permeability = 1e-6
set Exponential fluid weakening factor = 0
set Fluid compressibility = 0
set Fluid reaction time scale for operator splitting = 1e50
set Fluid-solid reaction scheme = zero solubility
end

subsection Visco Plastic
set Viscosity averaging scheme = geometric
set Viscous flow law = diffusion
set Prefactors for diffusion creep = 5e-21
set Stress exponents for diffusion creep = 1.0
set Activation energies for diffusion creep = 0
set Activation volumes for diffusion creep = 0
set Densities = 3000

set Angles of internal friction = 0
set Cohesions = 1e50

set Minimum viscosity = 1e21
set Maximum viscosity = 1e21
set Thermal expansivities = 0
end
end

# Set the global refinement to 0, 1 km x 1 km mesh.
subsection Mesh refinement
set Coarsening fraction = 0.1
set Refinement fraction = 0.9
set Initial adaptive refinement = 2
set Initial global refinement = 0
set Strategy = composition threshold, minimum refinement function
set Time steps between mesh refinement = 1

# Minimum of 4 global refinements
subsection Minimum refinement function
set Function expression = 0
end

# Refine where the porosity is bigger than 1e-6. Other compositions
# are set to 1e50 to ensure that we do not refine based on these
# compositions.
subsection Composition threshold
set Compositional field thresholds = 1e-3, 1e50, 1e50
end
end


# Output the darcy velocity
subsection Postprocess
set List of postprocessors = velocity statistics, visualization

subsection Visualization
set List of output variables = darcy velocity
set Time between graphical output = 250
set Output format = vtu
end
end
Loading

0 comments on commit 59be152

Please sign in to comment.