Skip to content

Commit

Permalink
clean up, improve test
Browse files Browse the repository at this point in the history
  • Loading branch information
@bobmyhill
bobmyhill committed Aug 13, 2023
1 parent 0d0bb8b commit 52d7a57
Show file tree
Hide file tree
Showing 6 changed files with 3,663 additions and 103 deletions.
10 changes: 5 additions & 5 deletions include/aspect/material_model/steinberger.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,7 @@ namespace aspect
double thermal_conductivity (const double temperature,
const double pressure,
const Point<dim> &position) const;
MaterialUtilities::CompositionalAveragingOperation viscosity_averaging;

/**
* Whether the compositional fields representing mass fractions
* should be normalized to one when computing their fractions
Expand Down Expand Up @@ -276,13 +276,13 @@ namespace aspect
std::vector<double> saturation_scaling;
double maximum_conductivity;


/**
* viscosity_prefactors are dimensionless quantities which
* are multiplied by the reference viscosity profile based on
* volume fraction
* Compositional prefactors with which to multiply the reference viscosity.
* Volume fractions are used to weight the prefactors according to the
* assigned viscosity averaging scheme.
*/
std::vector<double> viscosity_prefactors;
MaterialUtilities::CompositionalAveragingOperation viscosity_averaging_scheme;

/**
* Information about lateral temperature averages.
Expand Down
78 changes: 37 additions & 41 deletions source/material_model/steinberger.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,8 +177,6 @@ namespace aspect
const double depth = this->get_geometry_model().depth(position);
const double adiabatic_temperature = this->get_adiabatic_conditions().temperature(position);



double delta_temperature;
if (use_lateral_average_temperature)
{
Expand All @@ -192,21 +190,18 @@ namespace aspect
// We here compute the lateral variation of viscosity due to temperature (temperature_prefactor) as
// V_lT = exp [-(H/nR)*dT/(T_adiabatic*(T_adiabatic + dT))] as in Eq. 6 of the paper.
// We get H/nR from the lateral_viscosity_lookup->lateral_viscosity function.
const double log_temperature_prefactor = -1.0*lateral_viscosity_lookup->lateral_viscosity(depth)*delta_temperature/(temperature*adiabatic_temperature);
const double log_thermal_prefactor = -1.0 * lateral_viscosity_lookup->lateral_viscosity(depth) * delta_temperature / (temperature * adiabatic_temperature);

// Limit the lateral viscosity variation to a reasonable interval
const double temperature_prefactor = std::max(std::min(std::exp(log_temperature_prefactor),max_lateral_eta_variation),1/max_lateral_eta_variation);
const double thermal_prefactor = std::max(std::min(std::exp(log_thermal_prefactor), max_lateral_eta_variation), 1/max_lateral_eta_variation);

const double compositional_prefactor = MaterialUtilities::average_value (volume_fractions, viscosity_prefactors, viscosity_averaging_scheme);

//Visc_rT = exp[(H/nR)/T_adiabatic], Eq. 7 of the paper
const double eta_ref = radial_viscosity_lookup->radial_viscosity(depth);

// std::cout<<"size vf"<<volume_fractions.size()<<std::endl;
// std::cout<<"size vis_pre"<<viscosity_prefactors.size()<<std::endl;
const double compositional_prefactor = MaterialUtilities::average_value (volume_fractions, viscosity_prefactors, viscosity_averaging);
// std::cout<<"size vf"<<volume_fractions.size()<<std::endl;

// Radial viscosity profile is multiplied with lateral and compositional viscosity variation
return std::max(std::min(temperature_prefactor * eta_ref * compositional_prefactor,max_eta),min_eta);
// Radial viscosity profile is multiplied by thermal and compositional prefactors
return std::max(std::min(thermal_prefactor * compositional_prefactor * eta_ref, max_eta), min_eta);
}


Expand Down Expand Up @@ -310,11 +305,9 @@ namespace aspect
eos_outputs[i].densities,
true);


if (in.requests_property(MaterialProperties::viscosity))
out.viscosities[i] = viscosity(in.temperature[i], in.pressure[i], volume_fractions[i], in.strain_rate[i], in.position[i]);


MaterialUtilities::fill_averaged_equation_of_state_outputs(eos_outputs[i], mass_fractions, volume_fractions[i], i, out);
fill_prescribed_outputs(i, volume_fractions[i], in, out);
}
Expand Down Expand Up @@ -463,10 +456,8 @@ namespace aspect
"model. Larger values will be cut off.");
prm.declare_entry ("Viscosity averaging scheme", "harmonic",
Patterns::Selection("arithmetic|harmonic|geometric|maximum composition"),
"When more than one compositional field is present at a point "
"with different viscosities, we need to come up with an average "
"viscosity at that point. Select a weighted harmonic, arithmetic, "
"geometric, or maximum composition.");
"Method to average viscosities over multiple compositional fields. "
"One of arithmetic, harmonic, geometric or maximum composition.");
prm.declare_entry ("Viscosity prefactors", "1",
Patterns::Anything(),
"List of dimensionless quantities for background mantle and compositional fields,"
Expand Down Expand Up @@ -501,9 +492,8 @@ namespace aspect
max_eta = prm.get_double ("Maximum viscosity");
max_lateral_eta_variation = prm.get_double ("Maximum lateral viscosity variation");
thermal_conductivity_value = prm.get_double ("Thermal conductivity");
viscosity_averaging = MaterialUtilities::parse_compositional_averaging_operation ("Viscosity averaging scheme",
prm);

viscosity_averaging_scheme = MaterialUtilities::parse_compositional_averaging_operation ("Viscosity averaging scheme",
prm);

// Rheological parameters
if (prm.get ("Thermal conductivity formulation") == "constant")
Expand Down Expand Up @@ -542,16 +532,15 @@ namespace aspect
equation_of_state.initialize_simulator (this->get_simulator());
equation_of_state.parse_parameters(prm);

// Check if compositional fields represent a composition
const std::vector<CompositionalFieldDescription> composition_descriptions = this->introspection().get_composition_descriptions();

// Assign background field and do some error checking
has_background_field = (equation_of_state.number_of_lookups() == this->introspection().n_chemical_composition_fields() + 1);
has_background_field = ((equation_of_state.number_of_lookups() == 1) ||
(equation_of_state.number_of_lookups() == this->introspection().n_chemical_composition_fields() + 1));
AssertThrow ((equation_of_state.number_of_lookups() == 1) ||
(equation_of_state.number_of_lookups() == this->introspection().n_chemical_composition_fields()) ||
(equation_of_state.number_of_lookups() == this->introspection().n_chemical_composition_fields() + 1),
ExcMessage("The Steinberger material model assumes that all compositional "
"fields of the type chemical composition correspond to mass fractions of "
ExcMessage("The Steinberger material model assumes that either there is a single material "
"in the simulation, or that all compositional fields of the type "
"chemical composition correspond to mass fractions of different "
"materials. There must either be one material lookup file, the same "
"number of material lookup files as compositional fields of type chemical "
"composition, or one additional file (if a background field is used). You "
Expand All @@ -561,22 +550,29 @@ namespace aspect
+ Utilities::int_to_string(this->introspection().n_chemical_composition_fields())
+ " fields of type chemical composition."));

has_background_field = (equation_of_state.number_of_lookups() == n_chemical_fields + 1);
std::vector<std::string> list_of_composition_names = this->introspection().get_composition_names();
// Parse the Viscosity prefactors parameter
std::vector<std::string> chemical_field_names = this->introspection().chemical_composition_field_names();
std::vector<std::string> compositional_field_names = this->introspection().get_composition_names();

// If there is only one lookup, the list of Viscosity prefactors should have length one.
// The following line is required for occasions when defined fields of type chemical composition
// are not intended to be used as materials.
if (equation_of_state.number_of_lookups() == 1)
list_of_composition_names.resize(1, "background");
if ((equation_of_state.number_of_lookups() == 1) && (has_background_field))
list_of_composition_names.resize(0, "background");
// if (equation_of_state.number_of_lookups() > 1)
// list_of_composition_names.resize(equation_of_state.number_of_lookups(), "background");
// if ((equation_of_state.number_of_lookups() > 1) && (has_background_field))
// list_of_composition_names.resize(equation_of_state.number_of_lookups(), "background");
if (equation_of_state.number_of_lookups() > 1)
list_of_composition_names.resize(n_chemical_fields);
viscosity_prefactors = Utilities::parse_map_to_double_array (prm.get("Viscosity prefactors"),
list_of_composition_names,
has_background_field,
"Viscosity prefactors");
{
chemical_field_names.clear();
compositional_field_names.clear();
}

if (has_background_field)
{
chemical_field_names.insert(chemical_field_names.begin(),"background");
compositional_field_names.insert(compositional_field_names.begin(),"background");
}

Utilities::MapParsing::Options options(chemical_field_names, "Viscosity prefactors");
options.list_of_allowed_keys = compositional_field_names;
viscosity_prefactors = Utilities::MapParsing::parse_map_to_double_array(prm.get("Viscosity prefactors"), options);

prm.leave_subsection();
}
prm.leave_subsection();
Expand Down
72 changes: 29 additions & 43 deletions tests/multicomponent_steinberger.prm
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
## Test the multicomponent steinberger material model. Three blobs with different viscosities move around. The reference viscosity of background is from the radial viscosity profile of Steinberger & Calderwood 2006.
# Test the multicomponent steinberger material model. Three blobs with different viscosities move around.
# The reference viscosity of background is from the radial viscosity profile of Steinberger & Calderwood 2006.

set Dimension = 2
set CFL number = 1.0
Expand All @@ -16,6 +17,9 @@ subsection Boundary temperature model
set Fixed temperature boundary indicators = top, bottom, left, right
end

subsection Boundary velocity model
set Tangential velocity boundary indicators = top, bottom, left, right
end

subsection Discretization
set Stokes velocity polynomial degree = 2
Expand All @@ -28,49 +32,56 @@ subsection Discretization
end
end


subsection Geometry model
set Model name = box
subsection Box
set X extent = 2800000
set Y extent = 2900000
set X repetitions = 4
set X extent = 2.e6
set Y extent = 2.e6
set X repetitions = 2
end
end


subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 10.0
end
end


subsection Initial temperature model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = 1600
set Function expression = 1600
end
end

subsection Compositional fields
set Number of fields = 3
end

subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = if( sqrt( (x-5.e5 )^2 + (y-5.e5)^2) < 2.e5 , 1,0); if( sqrt((x-10.e5)^2 + (y-15.e5)^2 ) < 2.e5, 1,0); if( sqrt((x-15.e5)^2 + (y-10.e5)^2) < 2.e5, 1, 0)
end
end

subsection Material model
set Model name = Steinberger
subsection Steinberger model
set Use lateral average temperature for viscosity = false
set Data directory = $ASPECT_SOURCE_DIR/data/material-model/steinberger/
set Radial viscosity file name = radial-visc.txt
set Lateral viscosity file name = temp-viscosity-prefactor.txt
set Viscosity prefactors = 1, 5, 0.1, 10
set Material file names = pyr-ringwood88.txt, pyr-ringwood88.txt, pyr-ringwood88.txt, pyr-ringwood88.txt
set Minimum viscosity = 1e19
set Maximum viscosity = 1e24
set Data directory = $ASPECT_SOURCE_DIR/data/material-model/steinberger/
set Radial viscosity file name = radial-visc.txt
set Lateral viscosity file name = temp-viscosity-prefactor.txt
set Viscosity prefactors = 1, 5, 0.1, 10
set Material file names = pyr-ringwood88.txt, pyr-ringwood88.txt, pyr-ringwood88.txt, pyr-ringwood88.txt
set Minimum viscosity = 1e19
set Maximum viscosity = 1e24
end
end


subsection Mesh refinement
set Additional refinement times =
set Initial adaptive refinement = 0
Expand All @@ -81,42 +92,17 @@ subsection Mesh refinement
set Time steps between mesh refinement = 0
end


# The parameters below this comment were created by the update script
# as replacement for the old 'Model settings' subsection. They can be
# safely merged with any existing subsections with the same name.


subsection Boundary velocity model
set Tangential velocity boundary indicators = top, bottom, left, right
end

subsection Compositional fields
set Number of fields = 3
end

subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
# set Function expression = if( x>=300000 && x<=400000 && y>=2500000 && y<=2600000, 1,0); if( x>=600000 && x<=700000 && y>=2000000 && y<=2100000, 1,0); if( x>=300000 && x<=400000 && y>=200000 && y<=300000, 1, 0)
set Function expression = if( sqrt( (x-5.e5 )^2 + (y-5.e5)^2) < 5.e4 , 1,0); if( sqrt((x-12.e5)^2 + (y-20.e5)^2 ) < 5.e4, 1,0); if( sqrt((x-17.e5)^2 + (y-25.e5)^2) < 5.e4, 1, 0)
end
end

subsection Termination criteria
set Termination criteria = end step
set End step = 2
set End step = 1
end

subsection Postprocess
set List of postprocessors = visualization
subsection Visualization
set Interpolate output = false
set List of output variables = viscosity,density
set Number of grouped files = 1
set Output format = vtu
set Time between graphical output = 60000
set Output format = gnuplot
end
end

Expand Down
13 changes: 0 additions & 13 deletions tests/multicomponent_steinberger/screen-output
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,3 @@
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
Number of active cells: 1,024 (on 5 levels)
Number of degrees of freedom: 26,647 (8,514+1,105+4,257+4,257+4,257+4,257)

Expand All @@ -29,11 +24,3 @@ Number of degrees of freedom: 26,647 (8,514+1,105+4,257+4,257+4,257+4,257)
Writing graphical output: output-multicomponent_steinberger/solution/solution-00001

Termination requested by criterion: end time


+----------------------------------------------+------------+------------+
+----------------------------------+-----------+------------+------------+
+----------------------------------+-----------+------------+------------+

-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

0 comments on commit 52d7a57

Please sign in to comment.