Skip to content

Commit

Permalink
Modified entropy plugin to use entropy reader
Browse files Browse the repository at this point in the history
  • Loading branch information
@RanpengLi
RanpengLi committed Jul 15, 2023
1 parent 7de8d76 commit 0fb4bd6
Show file tree
Hide file tree
Showing 4 changed files with 92 additions and 103 deletions.
35 changes: 15 additions & 20 deletions benchmarks/entropy_adiabat/plugins/entropy_model.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,11 +42,8 @@ namespace aspect
ExcMessage("The 'entropy model' material model requires the existence of a compositional field "
"named 'entropy'. This field does not exist."));

material_lookup = std::make_unique<Utilities::StructuredDataLookup<2>>(7, 1.0);
material_lookup->load_file(data_directory + material_file_name,
this->get_mpi_communicator());

// object.initialize(this->get_mpi_communicator());
entropy_reader = std::make_unique<MaterialUtilities::Lookup::EntropyReader>();
entropy_reader->initialize(this->get_mpi_communicator(), data_directory, material_file_name);
}

template <int dim>
Expand All @@ -72,16 +69,14 @@ namespace aspect
// This is a requirement of the projected density approximation for
// the Stokes equation and not related to the entropy formulation.
// Also convert pressure from Pa to bar, bar is used in the table.
Point<2> entropy_pressure(in.composition[i][entropy_index],
this->get_adiabatic_conditions().pressure(in.position[i]) / 1.e5);
// entropy = in.composition[i][entropy_index]
// pressure = this->get_adiabatic_conditions().pressure(in.position[i]) / 1.e5
const double entropy = in.composition[i][entropy_index];
const double pressure = this->get_adiabatic_conditions().pressure(in.position[i]) / 1.e5;

out.densities[i] = material_lookup->get_data(entropy_pressure, 1);
out.thermal_expansion_coefficients[i] = material_lookup->get_data(entropy_pressure, 2);
out.specific_heat[i] = material_lookup->get_data(entropy_pressure, 3);
out.densities[i] = entropy_reader->density(entropy,pressure);
out.thermal_expansion_coefficients[i] = entropy_reader->thermal_expansivity(entropy,pressure);
out.specific_heat[i] = entropy_reader->specific_heat(entropy,pressure);

const Tensor<1, 2> density_gradient = material_lookup->get_gradients(entropy_pressure, 1);
const Tensor<1, 2> density_gradient = entropy_reader->calc_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];

Expand All @@ -102,7 +97,7 @@ namespace aspect
// set up variable to interpolate prescribed field outputs onto temperature field
if (PrescribedTemperatureOutputs<dim> *prescribed_temperature_out = out.template get_additional_output<PrescribedTemperatureOutputs<dim>>())
{
prescribed_temperature_out->prescribed_temperature_outputs[i] = material_lookup->get_data(entropy_pressure, 0);
prescribed_temperature_out->prescribed_temperature_outputs[i] = entropy_reader->temperature(entropy,pressure);
}

// Calculate Viscosity
Expand All @@ -115,15 +110,15 @@ namespace aspect
const double reference_temperature = this->get_adiabatic_conditions().is_initialized()
? this->get_adiabatic_conditions().temperature(in.position[i])
: this->get_parameters().adiabatic_surface_temperature;
const double temperature = material_lookup->get_data(entropy_pressure, 0);
const double delta_temperature = temperature - reference_temperature;
const double temperature_lookup = entropy_reader -> temperature(entropy,pressure);
const double delta_temperature = temperature_lookup - reference_temperature;

// Steinberger & Calderwood viscosity
if (temperature * reference_temperature == 0)
if (temperature_lookup * reference_temperature == 0)
out.viscosities[i] = min_eta;
else
{
double vis_lateral = std::exp(-lateral_viscosity_prefactor * delta_temperature / (temperature * reference_temperature));
double vis_lateral = std::exp(-lateral_viscosity_prefactor * delta_temperature / (temperature_lookup * reference_temperature));
if (std::isnan(vis_lateral))
vis_lateral = 1.0;

Expand All @@ -134,8 +129,8 @@ 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] = material_lookup->get_data(entropy_pressure, 4);
seismic_out->vs[i] = material_lookup->get_data(entropy_pressure, 5);
seismic_out->vp[i] = entropy_reader-> seismic_vp(entropy, pressure);
seismic_out->vs[i] = entropy_reader->seismic_vs(entropy, pressure);
}
}
}
Expand Down
6 changes: 4 additions & 2 deletions benchmarks/entropy_adiabat/plugins/entropy_model.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include <aspect/utilities.h>
#include <aspect/simulator_access.h>
#include <aspect/material_model/rheology/ascii_depth_profile.h>
#include <aspect/material_model/utilities.h>

namespace aspect
{
Expand Down Expand Up @@ -130,9 +131,10 @@ namespace aspect
std::string material_file_name;

/**
* Pointer to the StructuredDataLookup object that holds the material data.
* Pointer to the EntropyReader that read in material data for
* given entropy and pressure.
*/
std::unique_ptr<Utilities::StructuredDataLookup<2>> material_lookup;
std::unique_ptr<MaterialUtilities::Lookup::EntropyReader> entropy_reader;

/**
* Pointer to the rheology model used for depth-dependence from an
Expand Down
65 changes: 49 additions & 16 deletions include/aspect/material_model/utilities.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,6 @@
#include <deal.II/base/signaling_nan.h>
#include <deal.II/base/parameter_handler.h>

// #include <aspect/utilities.h>

namespace aspect
{
Expand Down Expand Up @@ -268,48 +267,82 @@ namespace aspect
const MPI_Comm &comm);
};

//template <int dim>
// class EntropyReader: public MaterialModel::Interface<dim> , public ::aspect::SimulatorAccess<dim>
//{

//template <int dim>
class EntropyReader // : public ::aspect::SimulatorAccess<dim>
/**
* This class reads in entropy-pressure material table and looks up for material
* proporties when entropy and pressure are given
*/
class EntropyReader
{
public:
// This function is responsible for reading the table into material_lookup

/**
* Read the table into material_lookup
*/
void
initialize(const MPI_Comm comm);
initialize(const MPI_Comm comm,
const std::string data_directory,
const std::string material_file_name);

// This function reads the specific heat for a given entropy and pressure from material_lookup
/**
* Reads the specific heat for a given entropy and pressure from material_lookup
*/
double
specific_heat(const double entropy,
const double pressure) const;

// This function reads the density for a given entropy and pressure from material_lookup
/**
* Reads the density for a given entropy and pressure from material_lookup
*/
double
density(const double entropy,
const double pressure) const;

// This function reads the thermal expansivity for a given entropy and pressure from material_lookup
/**
* Reads the thermal_expansivity for a given entropy and pressure
* from material_lookup
*/
double
thermal_expansivity(const double entropy,
const double pressure) const;

/**
* Reads the compressibility for a given entropy and pressure from material_lookup
*/
double
compressibility(const double entropy,
const double pressure) const;

/**
* Reads the temperature for a given entropy and pressure from material_lookup
*/
double
temperature(const double entropy,
const double pressure) const;

/**
* Reads the seismic p wave velocity for a given entropy and pressure
* from material_lookup
*/
double
seismic_vp(const double entropy,
const double pressure) const;

/**
* Reads the seismic s wave velocity for a given entropy and pressure
* from material_lookup
*/
double
seismic_vs(const double entropy,
const double pressure) const;

/**
* Gets density gradient for a given entropy and pressure
*/
Tensor<1, 2>
density_gradient(const double entropy,
const double pressure) const;
calc_density_gradient(const double entropy,
const double pressure) const;

private:
// This object contains the actual data (just like in entropy_model.cc)
// Utilities::StructuredDataLookup<2> material_lookup;
std::unique_ptr<Utilities::StructuredDataLookup<2>> material_lookup;

std::string data_directory;
Expand Down
89 changes: 24 additions & 65 deletions source/material_model/utilities.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -790,9 +790,11 @@ namespace aspect
}


//template <int dim>

void
EntropyReader::initialize(const MPI_Comm comm)
EntropyReader::initialize(const MPI_Comm comm,
const std::string data_directory,
const std::string material_file_name)
{
material_lookup = std::make_unique<Utilities::StructuredDataLookup<2>>(7,1.0);
material_lookup->load_file(data_directory+material_file_name,
Expand Down Expand Up @@ -844,9 +846,27 @@ namespace aspect
return temperature;
}

double
EntropyReader::seismic_vp(const double entropy,
const double pressure) const
{
Point<2> entropy_pressure(entropy, pressure);
const double seismic_vp = material_lookup->get_data(entropy_pressure, 4);
return seismic_vp;
}

double
EntropyReader::seismic_vs(const double entropy,
const double pressure) const
{
Point<2> entropy_pressure(entropy, pressure);
const double seismic_vs = material_lookup->get_data(entropy_pressure, 5);
return seismic_vs;
}

Tensor<1, 2>
EntropyReader::density_gradient(const double entropy,
const double pressure) const
EntropyReader::calc_density_gradient(const double entropy,
const double pressure) const
{
Point<2> entropy_pressure(entropy, pressure);
const Tensor<1, 2> density_gradient= material_lookup->get_gradients(entropy_pressure, 1);
Expand All @@ -856,67 +876,6 @@ namespace aspect
}


////template <int dim>
// void
// EntropyReader::initialize()
// {
// material_lookup = std::make_unique<Utilities::StructuredDataLookup<2>>(7,1.0);
// material_lookup->load_file(data_directory+material_file_name,
// this->get_mpi_communicator());
// }

// double
// EntropyReader::specific_heat(const double entropy,
// const double pressure) const
// {
// Point<2> entropy_pressure(entropy, pressure);
// return specific_heat;
// }

// double
// EntropyReader::density(const double entropy,
// const double pressure) const
// {
// Point<2> entropy_pressure(entropy, pressure);
// density = material_lookup->get_data(entropy_pressure, 1);
// return density;
// }

// double
// EntropyReader::thermal_expansivity(const double entropy,
// const double pressure) const
// {
// Point<2> entropy_pressure(entropy, pressure);
// thermal_expansivity = material_lookup->get_data(entropy_pressure, 2)
// return thermal_expansivity;
// }

// double
// EntropyReader::compressibility(const double entropy,
// const double pressure) const
// {
// Point<2> entropy_pressure(entropy, pressure);
// compressibility = material_lookup->get_data(entropy_pressure, 2)
// return compressibility;
// }

// double
// EntropyReader::temperature(const double entropy,
// const double pressure) const
// {
// Point<2> entropy_pressure(entropy, pressure);
// temperature = material_lookup->get_data(entropy_pressure, 0)
// return temperature;
// }

// Tensor<1, 2>
// EntropyReader::density_gradient(const double entropy,
// const double pressure) const
// {
// Point<2> entropy_pressure(entropy, pressure);
// density_gradient= material_lookup->get_gradients(entropy_pressure, 1);
// }


std::vector<double>
compute_only_composition_fractions(const std::vector<double> &compositional_fields,
Expand Down

0 comments on commit 0fb4bd6

Please sign in to comment.