Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Require fields of type entropy in entropy assembler and adiabatic con… #5297

Merged
merged 2 commits into from
Jul 13, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions benchmarks/entropy_adiabat/entropy_conduction.prm
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@ end
subsection Compositional fields
set Number of fields = 2
set Names of fields = entropy, density_field
set Types of fields = entropy, density
set Compositional field methods = field, prescribed field
end

Expand Down
1 change: 1 addition & 0 deletions benchmarks/entropy_adiabat/entropy_half_space.prm
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ end
subsection Compositional fields
set Number of fields = 3
set Names of fields = entropy, density_field, half_space
set Types of fields = entropy, density, generic
set Compositional field methods = field, prescribed field, static
end

Expand Down
10 changes: 7 additions & 3 deletions benchmarks/entropy_adiabat/plugins/compute_entropy_profile.cc
Original file line number Diff line number Diff line change
Expand Up @@ -63,14 +63,18 @@ namespace aspect
"PrescribedTemperatureOutputs, which is required "
"for this adiabatic conditions plugin."));

const unsigned int entropy_index = this->introspection().compositional_index_for_name("entropy");
const std::vector<unsigned int> entropy_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::entropy);

AssertThrow(entropy_indices.size() == 1,
ExcMessage("The 'compute entropy' adiabatic conditions plugin "
"requires exactly 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)
in.composition[0][entropy_index] = surface_entropy;
in.composition[0][entropy_indices[0]] = surface_entropy;

// 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
Expand Down Expand Up @@ -300,7 +304,7 @@ namespace aspect
"Of course the entropy along an adiabat is constant. "
"This plugin requires the material model to provide an "
"additional output object of type PrescribedTemperatureOutputs. "
"It also requires that there is a compositional field named "
"It also requires that there is a compositional field of type "
"'entropy' that represents the entropy of the material.")
}
}
25 changes: 16 additions & 9 deletions benchmarks/entropy_adiabat/plugins/entropy_advection.cc
Original file line number Diff line number Diff line change
Expand Up @@ -41,8 +41,9 @@ namespace aspect
const FiniteElement<dim> &fe = this->get_fe();

const typename Simulator<dim>::AdvectionField advection_field = *scratch.advection_field;
const std::vector<CompositionalFieldDescription> composition_descriptions = this->introspection().get_composition_descriptions();
if (!advection_field.is_temperature()
&& introspection.name_for_compositional_index(advection_field.compositional_variable) != "entropy")
&& composition_descriptions[advection_field.compositional_variable].type != CompositionalFieldDescription::entropy)
return;

const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
Expand Down Expand Up @@ -183,8 +184,9 @@ namespace aspect
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
std::vector<double> residuals(n_q_points,0.0);

if (advection_field.is_temperature() ||
this->introspection().name_for_compositional_index(advection_field.compositional_variable) != "entropy")
const std::vector<CompositionalFieldDescription> composition_descriptions = this->introspection().get_composition_descriptions();
if (!advection_field.is_temperature()
&& composition_descriptions[advection_field.compositional_variable].type != CompositionalFieldDescription::entropy)
return residuals;

this->get_heating_model_manager().evaluate(scratch.material_model_inputs,
Expand Down Expand Up @@ -271,14 +273,19 @@ namespace aspect

// Replace all existing assemblers for the temperature and entropy fields by the one for the entropy equation.
const unsigned int temperature_index = 0;
assemblers.advection_system[temperature_index].resize(1);
assemblers.advection_system[temperature_index][0] = std::make_unique<Assemblers::EntropyAdvectionSystem<dim>>();
assemblers.advection_system[temperature_index].clear();
assemblers.advection_system[temperature_index].emplace_back (std::make_unique<Assemblers::EntropyAdvectionSystem<dim>>());
assemblers.advection_system_assembler_properties[temperature_index].needed_update_flags = update_hessians;

const unsigned int entropy_index = 1 + simulator_access.introspection().compositional_index_for_name("entropy");
assemblers.advection_system[entropy_index].resize(1);
assemblers.advection_system[entropy_index][0] = std::make_unique<Assemblers::EntropyAdvectionSystem<dim>>();
assemblers.advection_system_assembler_properties[entropy_index].needed_update_flags = update_hessians;
const std::vector<CompositionalFieldDescription> composition_descriptions = simulator_access.introspection().get_composition_descriptions();
for (unsigned int c=0; c<composition_descriptions.size(); ++c)
if (composition_descriptions[c].type == CompositionalFieldDescription::entropy)
{
const unsigned int entropy_index = c + 1;
assemblers.advection_system[entropy_index].clear();
assemblers.advection_system[entropy_index].emplace_back (std::make_unique<Assemblers::EntropyAdvectionSystem<dim>>());
assemblers.advection_system_assembler_properties[entropy_index].needed_update_flags = update_hessians;
}
}
} // namespace aspect

Expand Down