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.

Sign up for GitHub

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

Jump to bottom

allow a different list of assemblers for each advection field #5233

Merged
merged 1 commit into from
Jul 9, 2023
Merged

allow a different list of assemblers for each advection field #5233

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
15 changes: 10 additions & 5 deletions benchmarks/entropy_adiabat/plugins/entropy_advection.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -269,11 +269,16 @@ namespace aspect
ExcMessage("The entropy advection assembler requires "
"that adiabatic heating is disabled."));

// Replace all existing assemblers by the one for the entropy equation.
assemblers.advection_system.resize(1);
assemblers.advection_system[0] = std::make_unique<Assemblers::EntropyAdvectionSystem<dim>>();

assemblers.advection_system_assembler_properties[0].needed_update_flags = update_hessians;
// 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_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;
}
} // namespace aspect

Expand Down
4 changes: 4 additions & 0 deletions doc/modules/changes/20230709_dannberg
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
New: ASPECT now allows it to select a different list of
assemblers for each advection field.
<br>
(Juliane Dannberg, 2023/07/09)
29 changes: 15 additions & 14 deletions include/aspect/simulator/assemblers/interface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -716,28 +716,29 @@ namespace aspect
std::vector<std::unique_ptr<Assemblers::Interface<dim>>> stokes_system_on_boundary_face;

/**
* A vector of pointers containing all assemblers for the advection systems.
* A vector of vectors of pointers containing a list of all assemblers
* for each individual advection system.
* These assemblers are called once per cell.
*/
std::vector<std::unique_ptr<Assemblers::Interface<dim>>> advection_system;
std::vector<std::vector<std::unique_ptr<Assemblers::Interface<dim>>>> advection_system;

/**
* A vector of pointers containing all assemblers for the Advection
* systems that compute face contributions at boundaries. These assemblers
* are called once per boundary face with the properly initialized inputs,
* therefore they allow terms that only exist on boundary faces (e.g.
* flux boundary conditions).
* A vector of vectors of pointers containing a list of all assemblers
* for the individual advection systems that compute face contributions
* at boundaries. These assemblers are called once per boundary face with
* the properly initialized inputs, therefore they allow terms that only
* exist on boundary faces (e.g. flux boundary conditions).
*/
std::vector<std::unique_ptr<Assemblers::Interface<dim>>> advection_system_on_boundary_face;
std::vector<std::vector<std::unique_ptr<Assemblers::Interface<dim>>>> advection_system_on_boundary_face;

/**
* A vector of pointers containing all assemblers for the Advection
* systems that compute face contributions on faces between cells.
* These assemblers are called once per interior face with the properly
* initialized inputs, therefore they allow terms that only exist on
* interior faces (e.g. DG penalty terms).
* A vector of vectors of pointers containing a list of all assemblers
* for the individual advection systems that compute face contributions
* on faces between cells. These assemblers are called once per interior
* face with the properly initialized inputs, therefore they allow terms
* that only exist on interior faces (e.g. DG penalty terms).
*/
std::vector<std::unique_ptr<Assemblers::Interface<dim>>> advection_system_on_interior_face;
std::vector<std::vector<std::unique_ptr<Assemblers::Interface<dim>>>> advection_system_on_interior_face;

/**
* A structure that describes what information an assembler function
Expand Down
160 changes: 86 additions & 74 deletions source/simulator/assembly.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,68 +146,72 @@ namespace aspect
Simulator<dim>::
set_advection_assemblers()
{
assemblers->advection_system.push_back(
std::make_unique<aspect::Assemblers::AdvectionSystem<dim>>());

// add the diffusion assemblers if we have fields that use this method
if (std::find(parameters.compositional_field_methods.begin(), parameters.compositional_field_methods.end(),
Parameters<dim>::AdvectionFieldMethod::prescribed_field_with_diffusion) != parameters.compositional_field_methods.end()
|| parameters.temperature_method == Parameters<dim>::AdvectionFieldMethod::prescribed_field_with_diffusion)
assemblers->advection_system.push_back(
std::make_unique<aspect::Assemblers::DiffusionSystem<dim>>());

// add the darcy assemblers if we have fields that use this method
if (std::find(parameters.compositional_field_methods.begin(), parameters.compositional_field_methods.end(),
Parameters<dim>::AdvectionFieldMethod::fem_darcy_field) != parameters.compositional_field_methods.end())
assemblers->advection_system.push_back(
std::make_unique<aspect::Assemblers::DarcySystem<dim>>());

if (parameters.use_discontinuous_temperature_discretization ||
parameters.use_discontinuous_composition_discretization)
// Loop over all advection fields.
for (unsigned int i=0; i<1+introspection.n_compositional_fields; ++i)
{
const bool no_field_method = std::find(parameters.compositional_field_methods.begin(),
parameters.compositional_field_methods.end(),
Parameters<dim>::AdvectionFieldMethod::fem_field)
== parameters.compositional_field_methods.end();

// TODO: This currently does not work in parallel, because the sparsity
// pattern of the matrix does not seem to know about flux terms
// across periodic faces of different levels. Fix this.
AssertThrow(geometry_model->get_periodic_boundary_pairs().size() == 0 ||
Utilities::MPI::n_mpi_processes(mpi_communicator) == 1 ||
no_field_method ||
(parameters.initial_adaptive_refinement == 0 &&
parameters.adaptive_refinement_interval == 0),
ExcMessage("Combining discontinuous elements with periodic boundaries and "
"adaptive mesh refinement in parallel models is currently not supported. "
"Please switch off any of those options or run on a single process."));

assemblers->advection_system_on_boundary_face.push_back(
std::make_unique<aspect::Assemblers::AdvectionSystemBoundaryFace<dim>>());

assemblers->advection_system_on_interior_face.push_back(
std::make_unique<aspect::Assemblers::AdvectionSystemInteriorFace<dim>>());
}
assemblers->advection_system[i].push_back(
std::make_unique<aspect::Assemblers::AdvectionSystem<dim>>());

// add the diffusion assemblers if advection field i uses this method
if (std::find(parameters.compositional_field_methods.begin(), parameters.compositional_field_methods.end(),
Parameters<dim>::AdvectionFieldMethod::prescribed_field_with_diffusion) != parameters.compositional_field_methods.end()
|| parameters.temperature_method == Parameters<dim>::AdvectionFieldMethod::prescribed_field_with_diffusion)
assemblers->advection_system[i].push_back(
std::make_unique<aspect::Assemblers::DiffusionSystem<dim>>());

// add the darcy assemblers if we have fields that use this method
if (std::find(parameters.compositional_field_methods.begin(), parameters.compositional_field_methods.end(),
Parameters<dim>::AdvectionFieldMethod::fem_darcy_field) != parameters.compositional_field_methods.end())
assemblers->advection_system[i].push_back(
std::make_unique<aspect::Assemblers::DarcySystem<dim>>());

if (parameters.use_discontinuous_temperature_discretization ||
parameters.use_discontinuous_composition_discretization)
{
const bool no_field_method = std::find(parameters.compositional_field_methods.begin(),
parameters.compositional_field_methods.end(),
Parameters<dim>::AdvectionFieldMethod::fem_field)
== parameters.compositional_field_methods.end();

// TODO: This currently does not work in parallel, because the sparsity
// pattern of the matrix does not seem to know about flux terms
// across periodic faces of different levels. Fix this.
AssertThrow(geometry_model->get_periodic_boundary_pairs().size() == 0 ||
Utilities::MPI::n_mpi_processes(mpi_communicator) == 1 ||
no_field_method ||
(parameters.initial_adaptive_refinement == 0 &&
parameters.adaptive_refinement_interval == 0),
ExcMessage("Combining discontinuous elements with periodic boundaries and "
"adaptive mesh refinement in parallel models is currently not supported. "
"Please switch off any of those options or run on a single process."));

assemblers->advection_system_on_boundary_face[i].push_back(
std::make_unique<aspect::Assemblers::AdvectionSystemBoundaryFace<dim>>());

assemblers->advection_system_on_interior_face[i].push_back(
std::make_unique<aspect::Assemblers::AdvectionSystemInteriorFace<dim>>());
}

if (parameters.fixed_heat_flux_boundary_indicators.size() != 0)
{
assemblers->advection_system_on_boundary_face.push_back(
std::make_unique<aspect::Assemblers::AdvectionSystemBoundaryHeatFlux<dim>>());
}
if (parameters.fixed_heat_flux_boundary_indicators.size() != 0)
{
assemblers->advection_system_on_boundary_face[i].push_back(
std::make_unique<aspect::Assemblers::AdvectionSystemBoundaryHeatFlux<dim>>());
}

if (parameters.use_discontinuous_temperature_discretization
|| parameters.fixed_heat_flux_boundary_indicators.size() != 0)
{
assemblers->advection_system_assembler_on_face_properties[0].need_face_material_model_data = true;
assemblers->advection_system_assembler_on_face_properties[0].need_face_finite_element_evaluation = true;
}
if (parameters.use_discontinuous_temperature_discretization
|| parameters.fixed_heat_flux_boundary_indicators.size() != 0)
{
assemblers->advection_system_assembler_on_face_properties[0].need_face_material_model_data = true;
assemblers->advection_system_assembler_on_face_properties[0].need_face_finite_element_evaluation = true;
}

if (parameters.use_discontinuous_composition_discretization)
{
for (unsigned int i = 1; i<=introspection.n_compositional_fields; ++i)
if (parameters.use_discontinuous_composition_discretization)
{
assemblers->advection_system_assembler_on_face_properties[i].need_face_material_model_data = true;
assemblers->advection_system_assembler_on_face_properties[i].need_face_finite_element_evaluation = true;
for (unsigned int c = 1; c<=introspection.n_compositional_fields; ++c)
{
assemblers->advection_system_assembler_on_face_properties[c].need_face_material_model_data = true;
assemblers->advection_system_assembler_on_face_properties[c].need_face_finite_element_evaluation = true;
}
}
}
}
Expand All @@ -220,6 +224,9 @@ namespace aspect
// first let the manager delete all existing assemblers:
assemblers->reset();

assemblers->advection_system.resize(1+introspection.n_compositional_fields);
assemblers->advection_system_on_boundary_face.resize(1+introspection.n_compositional_fields);
assemblers->advection_system_on_interior_face.resize(1+introspection.n_compositional_fields);
assemblers->advection_system_assembler_properties.resize(1+introspection.n_compositional_fields);
assemblers->advection_system_assembler_on_face_properties.resize(1+introspection.n_compositional_fields);

Expand Down Expand Up @@ -258,9 +265,13 @@ namespace aspect
initialize_simulator(*this,assemblers->stokes_preconditioner);
initialize_simulator(*this,assemblers->stokes_system);
initialize_simulator(*this,assemblers->stokes_system_on_boundary_face);
initialize_simulator(*this,assemblers->advection_system);
initialize_simulator(*this,assemblers->advection_system_on_boundary_face);
initialize_simulator(*this,assemblers->advection_system_on_interior_face);

for (unsigned int i=0; i<1+introspection.n_compositional_fields; ++i)
{
initialize_simulator(*this,assemblers->advection_system[i]);
initialize_simulator(*this,assemblers->advection_system_on_boundary_face[i]);
initialize_simulator(*this,assemblers->advection_system_on_interior_face[i]);
}
}


Expand Down Expand Up @@ -881,8 +892,9 @@ namespace aspect
current_linearization_point,
true);

for (unsigned int i=0; i<assemblers->advection_system.size(); ++i)
assemblers->advection_system[i]->create_additional_material_model_outputs(scratch.material_model_outputs);
for (unsigned int i=0; i<1+introspection.n_compositional_fields; ++i)
for (unsigned int j=0; j<assemblers->advection_system[i].size(); ++j)
assemblers->advection_system[i][j]->create_additional_material_model_outputs(scratch.material_model_outputs);

heating_model_manager.create_additional_material_model_inputs_and_outputs(scratch.material_model_inputs,
scratch.material_model_outputs);
Expand Down Expand Up @@ -941,14 +953,14 @@ namespace aspect

// trigger the invocation of the various functions that actually do
// all of the assembling
for (unsigned int i=0; i<assemblers->advection_system.size(); ++i)
assemblers->advection_system[i]->execute(scratch,data);
for (unsigned int i=0; i<assemblers->advection_system[advection_field.field_index()].size(); ++i)
assemblers->advection_system[advection_field.field_index()][i]->execute(scratch,data);

// then also work on possible face terms. if necessary, initialize
// the material model data on faces
const bool has_boundary_face_assemblers = !assemblers->advection_system_on_boundary_face.empty()
const bool has_boundary_face_assemblers = !assemblers->advection_system_on_boundary_face[advection_field.field_index()].empty()
&& assemblers->advection_system_assembler_on_face_properties[advection_field.field_index()].need_face_finite_element_evaluation;
const bool has_interior_face_assemblers = !assemblers->advection_system_on_interior_face.empty()
const bool has_interior_face_assemblers = !assemblers->advection_system_on_interior_face[advection_field.field_index()].empty()
&& assemblers->advection_system_assembler_on_face_properties[advection_field.field_index()].need_face_finite_element_evaluation;

// skip the remainder if no work needs to be done on faces
Expand Down Expand Up @@ -995,11 +1007,11 @@ namespace aspect
current_linearization_point,
true);

for (unsigned int i=0; i<assemblers->advection_system_on_boundary_face.size(); ++i)
assemblers->advection_system_on_boundary_face[i]->create_additional_material_model_outputs(scratch.face_material_model_outputs);
for (unsigned int i=0; i<assemblers->advection_system_on_boundary_face[advection_field.field_index()].size(); ++i)
assemblers->advection_system_on_boundary_face[advection_field.field_index()][i]->create_additional_material_model_outputs(scratch.face_material_model_outputs);

for (unsigned int i=0; i<assemblers->advection_system_on_interior_face.size(); ++i)
assemblers->advection_system_on_interior_face[i]->create_additional_material_model_outputs(scratch.face_material_model_outputs);
for (unsigned int i=0; i<assemblers->advection_system_on_interior_face[advection_field.field_index()].size(); ++i)
assemblers->advection_system_on_interior_face[advection_field.field_index()][i]->create_additional_material_model_outputs(scratch.face_material_model_outputs);

heating_model_manager.create_additional_material_model_inputs_and_outputs(scratch.face_material_model_inputs,
scratch.face_material_model_outputs);
Expand Down Expand Up @@ -1035,11 +1047,11 @@ namespace aspect

// handle faces at periodic boundaries like interior faces
if (face->at_boundary() && !cell->has_periodic_neighbor (scratch.face_number))
for (unsigned int i=0; i<assemblers->advection_system_on_boundary_face.size(); ++i)
assemblers->advection_system_on_boundary_face[i]->execute(scratch,data);
for (unsigned int i=0; i<assemblers->advection_system_on_boundary_face[advection_field.field_index()].size(); ++i)
assemblers->advection_system_on_boundary_face[advection_field.field_index()][i]->execute(scratch,data);
else
for (unsigned int i=0; i<assemblers->advection_system_on_interior_face.size(); ++i)
assemblers->advection_system_on_interior_face[i]->execute(scratch,data);
for (unsigned int i=0; i<assemblers->advection_system_on_interior_face[advection_field.field_index()].size(); ++i)
assemblers->advection_system_on_interior_face[advection_field.field_index()][i]->execute(scratch,data);
}
}
}
Expand Down