Skip to content

Commit

Permalink
Merge pull request #6218 from jdannberg/fix_latent_heat_pressure_deri…
Browse files Browse the repository at this point in the history 
…vative

Fix pressure depth derivative in latent heat material model
  • Loading branch information
@gassmoeller
gassmoeller authored Feb 2, 2025
2 parents e148692 + 7404912 commit 6ce36dc
Show file tree
Hide file tree
Showing 4 changed files with 26 additions and 39 deletions.
52 changes: 21 additions & 31 deletions source/material_model/latent_heat.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ namespace aspect
{
const double temperature = in.temperature[i];
const double pressure = in.pressure[i];
const std::vector<double> composition = in.composition[i];
const std::vector<double> &composition = in.composition[i];
const Point<dim> position = in.position[i];

// Assign constant material properties
Expand Down Expand Up @@ -80,6 +80,20 @@ namespace aspect
out.viscosities[i] = visc_composition_dependence * visc_temperature_dependence * eta;
}

// Prepare some values for phase dependent properties
const double depth = this->get_geometry_model().depth(in.position[i]);
const double adiabatic_pressure = this->get_adiabatic_conditions().pressure(in.position[i]);
const double pressure_depth_derivative = (depth > 0.0)
?
adiabatic_pressure / depth
:
this->get_gravity_model().gravity_vector(in.position[i]).norm() * reference_rho;
MaterialUtilities::PhaseFunctionInputs<dim> phase_in(temperature,
pressure,
depth,
pressure_depth_derivative,
0);

// Calculate density
// and phase dependence of viscosity
{
Expand All @@ -90,7 +104,7 @@ namespace aspect

// second, calculate composition dependence of density
// constant density difference between peridotite and eclogite
const double density_composition_dependence = composition.size()>0
const double density_composition_dependence = composition.size() > 0
?
compositional_delta_rho * composition[0]
:
Expand All @@ -115,28 +129,16 @@ namespace aspect
// Loop through phase transitions
for (unsigned int transition_index=0; transition_index<phase_function.n_phase_transitions(); ++transition_index)
{
const double depth = this->get_geometry_model().depth(position);
const double pressure_depth_derivative = (depth > 0.0)
?
this->get_adiabatic_conditions().pressure(position) / depth
:
this->get_gravity_model().gravity_vector(position).norm() * reference_rho;

const MaterialUtilities::PhaseFunctionInputs<dim> phase_in(temperature,
pressure,
depth,
pressure_depth_derivative,
transition_index);

phase_in.phase_transition_index = transition_index;
const double phaseFunction = phase_function.compute_value(phase_in);

// For the densities we have a list of jumps, so the index used
// in the loop corresponds to the index of the phase transition.
if (composition.size()==0)
if (composition.size() == 0)
{
phase_dependence += phaseFunction * density_jumps[transition_index];
}
else if (composition.size()>0)
else if (composition.size() > 0)
{
if (transition_phases[transition_index] == 0) // 1st compositional field
phase_dependence += phaseFunction * density_jumps[transition_index] * (1.0 - composition[0]);
Expand Down Expand Up @@ -195,24 +197,12 @@ namespace aspect
if (this->get_adiabatic_conditions().is_initialized() && this->include_latent_heat())
for (unsigned int phase=0; phase<phase_function.n_phase_transitions(); ++phase)
{
const double depth = this->get_geometry_model().depth(in.position[i]);
const double pressure_depth_derivative = (this->get_adiabatic_conditions().pressure(position) > 0)
?
depth / this->get_adiabatic_conditions().pressure(position)
:
this->get_gravity_model().gravity_vector(in.position[i]).norm() * reference_rho;

const MaterialUtilities::PhaseFunctionInputs<dim> phase_in(temperature,
pressure,
depth,
pressure_depth_derivative,
phase);

phase_in.phase_transition_index = phase;
const double PhaseFunctionDerivative = phase_function.compute_derivative(phase_in);
const double clapeyron_slope = phase_function.get_transition_slope(phase);

double entropy_change = 0.0;
if (composition.size()==0) // only one compositional field
if (composition.size() == 0) // only one compositional field
entropy_change = clapeyron_slope * density_jumps[phase] / (rho * rho);
else
{
Expand Down
2 changes: 1 addition & 1 deletion source/material_model/utilities.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1216,7 +1216,7 @@ namespace aspect
start_phase_transition_index += n_phase_transitions_per_chemical_composition[n_relevant_fields];
}

const double pressure_in_bar = in.pressure/1e5;
const double pressure_in_bar = in.pressure/1.e5;

Assert (in.temperature >= minimum_temperature[n_comp] && in.temperature < maximum_temperature[n_comp], ExcInternalError());
Assert (pressure_in_bar >= minimum_pressure[n_comp] && pressure_in_bar < maximum_pressure[n_comp], ExcInternalError());
Expand Down
9 changes: 3 additions & 6 deletions tests/latent_heat_enthalpy.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ namespace aspect
{
public:
void
initialize()
initialize() override
{
const std::string datadirectory = Utilities::expand_ASPECT_SOURCE_DIR("$ASPECT_SOURCE_DIR/data/material-model/latent-heat-enthalpy-test/");
const std::string material_file_names = "testdata.txt";
Expand All @@ -49,21 +49,18 @@ namespace aspect
}

bool
is_compressible () const
is_compressible () const override
{
return false;
}

void
evaluate(const typename Interface<dim>::MaterialModelInputs &in, typename Interface<dim>::MaterialModelOutputs &out) const
evaluate(const typename Interface<dim>::MaterialModelInputs &in, typename Interface<dim>::MaterialModelOutputs &out) const override
{
const double reference_rho = 3400;
const double density_jump = 115.6;
const double reference_T = 1000;
const double eta = 8.44e21;
const double k_value = 2.38;
const double reference_specific_heat = 1000;
const double thermal_alpha = 0;

double dHdT = 0.0;
double dHdp = 0.0;
Expand Down
2 changes: 1 addition & 1 deletion tests/material_model_dependencies.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ namespace aspect
public:
virtual
std::pair<std::string,std::string>
execute (TableHandler &statistics)
execute (TableHandler & /*statistics*/) override
{
using namespace MaterialModel;
using namespace NonlinearDependence;
Expand Down

0 comments on commit 6ce36dc

Please sign in to comment.