include density_gradient

include/aspect/material_model/equation_of_state/interface.h

@@ -99,6 +99,7 @@ namespace aspect
       std::vector<double> entropy_derivative_temperature;
       std::vector<double> vp;
       std::vector<double> vs;
+      std::vector<Tensor<1, 2>> density_gradient;
     };
 
 

source/material_model/entropy_model.cc

@@ -86,14 +86,15 @@ namespace aspect
                              "iterates over the advection equations but a non iterating solver scheme was selected. "
                              "Please check the consistency of your solver scheme."));
 
-    //for (unsigned i = 0; i < n_material_lookups; ++i)
-      for (unsigned i = 0; i < material_file_names.size(); ++i)
+      AssertThrow(material_file_names.size() == 1,
+                  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>());
           entropy_reader[i]->initialize(this->get_mpi_communicator(), data_directory, material_file_names[i]);
           }
-     // entropy_reader = std::make_unique<MaterialUtilities::Lookup::EntropyReader>();
-     // entropy_reader->initialize(this->get_mpi_communicator(), data_directory, material_file_name);
 
       lateral_viscosity_prefactor_lookup = std::make_unique<internal::LateralViscosityLookup>(data_directory+lateral_viscosity_file_name,
                                            this->get_mpi_communicator());
@@ -179,19 +180,17 @@ namespace aspect
           const double entropy = in.composition[i][entropy_index];
           const double pressure = this->get_adiabatic_conditions().pressure(in.position[i]) / 1.e5;
 
-/*
-          eos_outputs[i].densities[0] = entropy_reader[0]->density(entropy,pressure);
-          eos_outputs[i].thermal_expansion_coefficients[0] = entropy_reader[0]->thermal_expansivity(entropy,pressure);
-          eos_outputs[i].specific_heat_capacities[0] = entropy_reader[0]->specific_heat(entropy,pressure);
- */ 
-
 
-       //   for (unsigned int j=0; j<eos_outputs[i].densities.size(); ++j)
+
        for (unsigned int j=0; j<material_file_names.size(); ++j)
             {
               eos_outputs[i].densities[j] = entropy_reader[j]->density(entropy, pressure);
               eos_outputs[i].thermal_expansion_coefficients[j] = entropy_reader[j]->thermal_expansivity(entropy,pressure);
               eos_outputs[i].specific_heat_capacities[j] = entropy_reader[j]->specific_heat(entropy,pressure);
+              eos_outputs[i].density_gradient[j] = entropy_reader[j]->density_gradient(entropy,pressure);
+
+              const Tensor<1, 2> pressure_unit_vector({0.0, 1.0});
+              eos_outputs[i].compressibilities[j] = (eos_outputs[i].density_gradient[j] * pressure_unit_vector) / eos_outputs[i].densities[j];
             }
 
 
@@ -218,23 +217,18 @@ namespace aspect
                                                                                true);
 
           out.densities[i] = MaterialUtilities::average_value (volume_fractions[i], eos_outputs[i].densities, MaterialUtilities::arithmetic);
-
-     //     out.compressibilities[i] = MaterialUtilities::average_value (volume_fractions[i], eos_outputs.compressibilities, MaterialUtilities::arithmetic);
-
+
           out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value (volume_fractions[i], eos_outputs[i].thermal_expansion_coefficients, MaterialUtilities::arithmetic);
 
           out.specific_heat[i] = MaterialUtilities::average_value (mass_fractions, eos_outputs[i].specific_heat_capacities, MaterialUtilities::arithmetic);
 
+          out.compressibilities[i] = MaterialUtilities::average_value (mass_fractions, eos_outputs[i].compressibilities, MaterialUtilities::arithmetic);
 
-
-//          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);
-
+/*
           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];
-
+*/
 
           // Thermal conductivity can be pressure temperature dependent
           const double temperature_lookup =  entropy_reader[0]->temperature(entropy,pressure);
@@ -472,8 +466,7 @@ namespace aspect
         {
           data_directory              = Utilities::expand_ASPECT_SOURCE_DIR(prm.get ("Data directory"));
           material_file_names          = Utilities::split_string_list(prm.get ("Material file name"));
-      //  n_material_lookups           = material_file_names.size();          
-      //  material_file_name          = prm.get ("Material file name");
+
           lateral_viscosity_file_name  = prm.get ("Lateral viscosity file name");
           min_eta                     = prm.get_double ("Minimum viscosity");
           max_eta                     = prm.get_double ("Maximum viscosity");

source/material_model/equation_of_state/interface.cc

@@ -37,7 +37,8 @@ namespace aspect
       entropy_derivative_pressure(n_individual_compositions_and_phases, numbers::signaling_nan<double>()),
       entropy_derivative_temperature(n_individual_compositions_and_phases, numbers::signaling_nan<double>()),
       vp(n_individual_compositions_and_phases, numbers::signaling_nan<double>()),
-      vs(n_individual_compositions_and_phases, numbers::signaling_nan<double>())
+      vs(n_individual_compositions_and_phases, numbers::signaling_nan<double>()),
+      density_gradient(n_individual_compositions_and_phases, numbers::signaling_nan<Tensor<1, 2>>())
     {}