Multicomponent entropy method - average the looked-up material properties

source/adiabatic_conditions/compute_entropy_profile.cc

@@ -64,18 +64,29 @@ namespace aspect
  "for this adiabatic conditions plugin."));
 
  const std::vector<unsigned int> entropy_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::entropy);
-
- AssertThrow(entropy_indices.size() == 1,
+ // TODO : need to make it work for more than one field
+ AssertThrow(entropy_indices.size() >= 1,
  ExcMessage("The 'compute entropy' adiabatic conditions plugin "
- "requires exactly one field of type 'entropy'."));
+ "requires at least one field of type 'entropy'."));
 
  // Constant properties on the reference profile
  // We only need the material model to compute the density
  in.requested_properties = MaterialModel::MaterialProperties::density | MaterialModel::MaterialProperties::additional_outputs;
  in.velocity[0] = Tensor <1,dim> ();
  // The entropy along an adiabat is constant (equals the surface entropy)
+ // TODO : need to make it work for more than one entropy field
+ // TODO : provide more ways to specify compositional fields like in compute_profile.cc
+ for (unsigned int i=0; i < this->n_compositional_fields(); ++i)
+ in.composition[0][i] = 0;
+
  in.composition[0][entropy_indices[0]] = surface_entropy;
 
+ if (this->introspection().composition_type_exists(CompositionalFieldDescription::Type::chemical_composition))
+ {
+ const std::vector<unsigned int> &chemical_indices = this->introspection().chemical_composition_field_indices();
+ in.composition[0][chemical_indices[0]] = 1;
+ }
+
  // Check whether gravity is pointing up / out or down / in. In the normal case it should
  // point down / in and therefore gravity should be positive, leading to increasing
  // adiabatic pressures and temperatures with depth. In some cases it will point up / out

source/material_model/entropy_model.cc

@@ -28,6 +28,7 @@
 #include <iostream>
 #include <aspect/material_model/rheology/visco_plastic.h>
 #include <aspect/material_model/steinberger.h>
+#include <aspect/material_model/equation_of_state/interface.h>
 
 namespace aspect
 {
@@ -76,7 +77,7 @@ namespace aspect
  "'projected density field' approximation "
  "for the mass conservation equation, which is not selected."));
 
- AssertThrow (this->introspection().compositional_name_exists("entropy"),
+ AssertThrow (this->introspection().composition_type_exists(CompositionalFieldDescription::Type::entropy),
  ExcMessage("The 'entropy model' material model requires the existence of a compositional field "
  "named 'entropy'. This field does not exist."));
 
@@ -85,10 +86,12 @@ namespace aspect
  "iterates over the advection equations but a non iterating solver scheme was selected. "
  "Please check the consistency of your solver scheme."));
 
- AssertThrow(material_file_names.size() == 1,
+ AssertThrow(material_file_names.size() == 1 || SimulatorAccess<dim>::get_end_time () == 0,
  ExcMessage("The 'entropy model' material model can only handle one composition, "
  "and can therefore only read one material lookup table."));
 
+
+
  for (unsigned int i = 0; i < material_file_names.size(); ++i)
  {
  entropy_reader.push_back(std::make_unique<MaterialUtilities::Lookup::EntropyReader>());
@@ -161,7 +164,18 @@ namespace aspect
  MaterialModel::MaterialModelOutputs<dim> &out) const
  {
  const unsigned int projected_density_index = this->introspection().compositional_index_for_name("density_field");
- const unsigned int entropy_index = this->introspection().compositional_index_for_name("entropy");
+ //TODO : need to make it work for more than one field
+ const std::vector<unsigned int> &entropy_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::entropy);
+ const unsigned int entropy_index = entropy_indices[0];
+ const std::vector<unsigned int> &composition_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::chemical_composition);
+
+ AssertThrow(composition_indices.size() == material_file_names.size() - 1,
+ ExcMessage("The 'entropy model' material model assumes that there exists a background field in addition to the compositional fields, "
+ "and therefore it requires one more lookup table than there are chemical compositional fields."));
+
+ EquationOfStateOutputs<dim> eos_outputs (material_file_names.size());
+ std::vector<double> volume_fractions (material_file_names.size());
+ std::vector<double> mass_fractions (material_file_names.size());
 
  for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
  {
@@ -173,14 +187,38 @@ namespace aspect
  const double entropy = in.composition[i][entropy_index];
  const double pressure = this->get_adiabatic_conditions().pressure(in.position[i]) / 1.e5;
 
- out.densities[i] = entropy_reader[0]->density(entropy,pressure);
- out.thermal_expansion_coefficients[i] = entropy_reader[0]->thermal_expansivity(entropy,pressure);
- out.specific_heat[i] = entropy_reader[0]->specific_heat(entropy,pressure);
+ // Loop over all material files, and store the looked-up values for all compositions.
+ for (unsigned int j=0; j<material_file_names.size(); ++j)
+ {
+ eos_outputs.densities[j] = entropy_reader[j]->density(entropy, pressure);
+ eos_outputs.thermal_expansion_coefficients[j] = entropy_reader[j]->thermal_expansivity(entropy,pressure);
+ eos_outputs.specific_heat_capacities[j] = entropy_reader[j]->specific_heat(entropy,pressure);
 
- const Tensor<1, 2> density_gradient = entropy_reader[0]->density_gradient(entropy,pressure);
- const Tensor<1, 2> pressure_unit_vector({0.0, 1.0});
- out.compressibilities[i] = (density_gradient * pressure_unit_vector) / out.densities[i];
+  const Tensor<1, 2> pressure_unit_vector({0.0, 1.0});
+  eos_outputs.compressibilities[j] = ((entropy_reader[j]->density_gradient(entropy,pressure)) * pressure_unit_vector) / eos_outputs.densities[j];
+  }
 
+ // Calculate volume fractions from mass fractions
+ // If there is only one lookup table, set the mass and volume fractions to 1
+ if (material_file_names.size() == 1)
+ mass_fractions [0] = 1.0;
+
+ else
+ {
+ // We only want to compute mass/volume fractions for fields that are chemical compositions.
+ mass_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i], this->introspection().chemical_composition_field_indices());
+ }
+
+ volume_fractions = MaterialUtilities::compute_volumes_from_masses(mass_fractions,
+ eos_outputs.densities,
+ true);
+
+ out.densities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.densities, MaterialUtilities::arithmetic);
+ out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.thermal_expansion_coefficients, MaterialUtilities::arithmetic);
+
+ out.specific_heat[i] = MaterialUtilities::average_value (mass_fractions, eos_outputs.specific_heat_capacities, MaterialUtilities::arithmetic);
+
+ out.compressibilities[i] = MaterialUtilities::average_value (mass_fractions, eos_outputs.compressibilities, MaterialUtilities::arithmetic);
 
  // Thermal conductivity can be pressure temperature dependent
  const double temperature_lookup = entropy_reader[0]->temperature(entropy,pressure);
@@ -261,8 +299,16 @@ namespace aspect
  // fill seismic velocities outputs if they exist
  if (SeismicAdditionalOutputs<dim> *seismic_out = out.template get_additional_output<SeismicAdditionalOutputs<dim>>())
  {
- seismic_out->vp[i] = entropy_reader[0]->seismic_vp(entropy, pressure);
- seismic_out->vs[i] = entropy_reader[0]->seismic_vs(entropy, pressure);
+
+ std::vector<double> vp (material_file_names.size());
+ std::vector<double> vs (material_file_names.size());
+ for (unsigned int j=0; j<material_file_names.size(); ++j)
+ {
+ vp[j] = entropy_reader[j]->seismic_vp(entropy,pressure);
+ vs[j] = entropy_reader[j]->seismic_vs(entropy,pressure);
+ }
+ seismic_out->vp[i] = MaterialUtilities::average_value (volume_fractions, vp, MaterialUtilities::arithmetic);
+ seismic_out->vs[i] = MaterialUtilities::average_value (volume_fractions, vs, MaterialUtilities::arithmetic);
  }
  }
  }
@@ -286,7 +332,7 @@ namespace aspect
  "files located in the `data/' subdirectory of ASPECT.");
  prm.declare_entry ("Material file name", "material_table.txt",
  Patterns::List (Patterns::Anything()),
- "The file name of the material data.");
+ "The file name of the material data. The first material data file is intended for the background composition. ");
  prm.declare_entry ("Reference viscosity", "1e22",
  Patterns::Double(0),
  "The viscosity that is used in this model. "

tests/entropy_plasticity.prm

@@ -1,4 +1,4 @@
-# A box model for convection under early-Earth conditions, with plasticity,
+# A 2-D annulus model for convection under early-Earth conditions, with plasticity,
 # reference viscosity profile, and p-T dependent conductivity
 
 set Dimension = 2

tests/entropy_read_multi_table.cc

@@ -0,0 +1,21 @@
+/*
+ Copyright (C) 2022 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file LICENSE. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+
+#include "../benchmarks/entropy_adiabat/plugins/entropy_advection.cc"

tests/entropy_read_multi_table.prm

@@ -0,0 +1,120 @@
+# A box model with two compositions. It reads in two look-up tables
+# and computes the initial property based on the chemical compositions.
+
+set Dimension = 2
+set Use years in output instead of seconds = true
+set End time = 0
+set Maximum time step = 1e6
+set Nonlinear solver scheme = iterated Advection and Stokes
+set Max nonlinear iterations = 50
+set Surface pressure = 0
+set Adiabatic surface temperature = 1560.29
+
+subsection Formulation
+ set Mass conservation = projected density field
+end
+
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X extent = 13000
+ set Y extent = 13000
+ set Y repetitions = 1
+ end
+end
+
+
+subsection Adiabatic conditions model
+ # 'compute entropy profile' computes arbitrary adiabats
+ # internally, based on the data table.
+ set Model name = compute entropy profile
+
+ subsection Compute entropy profile
+ # Entropy equivalent to T=1560.29 K according to the used table
+ set Surface entropy = 2500
+ end
+end
+
+# We prescribe temperatures according to the data table.
+# This output is computed in the material model.
+subsection Temperature field
+ set Temperature method = prescribed field with diffusion
+end
+
+# We solve the entropy equation for the compositional field with type
+# 'entropy'. Temperature and density are then computed based on entropy and
+# pressure.
+subsection Compositional fields
+ set Number of fields = 4
+ set Names of fields = percent_bas, entropy_bas, entropy_haz, density_field
+ set Types of fields = chemical composition, entropy, entropy, density
+ set Compositional field methods = field, field, field, prescribed field
+end
+
+
+# Prescribing closed boundaries.
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = top, left, right, bottom
+end
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 0
+ end
+end
+
+# We start at the same entropy as the one used for the
+# adiabatic profile. The last compositional field gives
+# the analytical solution for the temperature at the
+# end of the model run based on the half-space cooling
+# model.
+subsection Initial composition model
+ set List of model names = function
+
+ subsection Function
+ # Entropy equivalent to T=1560.29 K according to table
+ set Function expression = 0.2; 2500; 2500; 0.0
+
+ end
+end
+
+# We use a data table for orthopyroxene computed using the thermodynamic
+# modeling software HeFESTo.
+# We use a very large viscosity to make sure the material does not move.
+subsection Material model
+ set Model name = entropy model
+
+ subsection Entropy model
+ set Data directory = $ASPECT_SOURCE_DIR/data/material-model/entropy-table/opxtable/
+ set Material file name = material_table.txt, multi_composition_test.txt
+ set Maximum viscosity = 1e30
+ set Minimum viscosity = 1e30
+ set Lateral viscosity file name = constant_lateral_vis_prefactor.txt
+ set Thermal conductivity formulation = constant
+ set Thermal conductivity = 4.7
+ end
+end
+
+subsection Mesh refinement
+ set Initial global refinement = 0
+ set Initial adaptive refinement = 0
+ set Time steps between mesh refinement = 0
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization, velocity statistics, temperature statistics, mass flux statistics, composition statistics
+
+ subsection Visualization
+ set Time between graphical output = 1e8
+ set List of output variables = material properties, adiabat, nonadiabatic pressure, nonadiabatic temperature
+
+ subsection Material properties
+ set List of material properties = density,thermal expansivity,specific heat,viscosity, compressibility
+ end
+
+ set Output format = gnuplot
+ end
+end

tests/entropy_read_multi_table/statistics

@@ -0,0 +1,39 @@
+# 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 composition solver 3
+# 13: Iterations for composition solver 4
+# 14: Iterations for Stokes solver
+# 15: Velocity iterations in Stokes preconditioner
+# 16: Schur complement iterations in Stokes preconditioner
+# 17: Visualization file name
+# 18: RMS velocity (m/year)
+# 19: Max. velocity (m/year)
+# 20: Minimal temperature (K)
+# 21: Average temperature (K)
+# 22: Maximal temperature (K)
+# 23: Outward mass flux through boundary with indicator 0 ("left") (kg/yr)
+# 24: Outward mass flux through boundary with indicator 1 ("right") (kg/yr)
+# 25: Outward mass flux through boundary with indicator 2 ("bottom") (kg/yr)
+# 26: Outward mass flux through boundary with indicator 3 ("top") (kg/yr)
+# 27: Minimal value for composition percent_bas
+# 28: Maximal value for composition percent_bas
+# 29: Global mass for composition percent_bas
+# 30: Minimal value for composition entropy_bas
+# 31: Maximal value for composition entropy_bas
+# 32: Global mass for composition entropy_bas
+# 33: Minimal value for composition entropy_haz
+# 34: Maximal value for composition entropy_haz
+# 35: Global mass for composition entropy_haz
+# 36: Minimal value for composition density_field
+# 37: Maximal value for composition density_field
+# 38: Global mass for composition density_field
+0 0.000000000000e+00 0.000000000000e+00 1 22 9 36 2 1 0 0 0 4294967294 4294967294 0 0 output-entropy_read_multi_table/solution/solution-00000 0.00000000e+00 0.00000000e+00 1.56029000e+03 1.56029000e+03 1.56029000e+03 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.00000000e-01 2.00000000e-01 3.38000000e+07 2.50000000e+03 2.50000000e+03 4.22500000e+11 2.50000000e+03 2.50000000e+03 4.22500000e+11 3.74396667e+03 3.74396667e+03 6.32730367e+11