Skip to content

Commit

Permalink
use SUNDIALS to compute reactions
Browse files Browse the repository at this point in the history
  • Loading branch information
@jdannberg
jdannberg committed Jun 14, 2024
1 parent 088b010 commit 314f518
Showing 1 changed file with 114 additions and 37 deletions.
151 changes: 114 additions & 37 deletions source/simulator/helper_functions.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,8 @@
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_values.h>

#include <deal.II/sundials/arkode.h>

#include <fstream>
#include <iostream>
#include <string>
Expand Down Expand Up @@ -1679,16 +1681,6 @@ namespace aspect
LinearAlgebra::BlockVector distributed_reaction_vector (introspection.index_sets.system_partitioning,
mpi_communicator);

// we use a different (potentially smaller) time step than in the advection scheme,
// and we want all of our reaction time steps (within one advection step) to have the same size
const unsigned int number_of_reaction_steps = std::max(static_cast<unsigned int>(time_step / parameters.reaction_time_step),
std::max(parameters.reaction_steps_per_advection_step,1U));

const double reaction_time_step_size = time_step / static_cast<double>(number_of_reaction_steps);

Assert (reaction_time_step_size > 0,
ExcMessage("Reaction time step must be greater than 0."));

pcout << " Solving composition reactions... " << std::flush;


Expand All @@ -1707,7 +1699,8 @@ namespace aspect
// For each field (T or composition) store where each support point is found in the list of indices of support points
// unique_support_points. This means index=support_point_index_by_field[0][3] is the index of the third support point
// of the temperature field and unique_support_points[index] is its location.
std::vector<std::vector<unsigned int>> support_point_index_by_field(1+introspection.n_compositional_fields);
const unsigned int n_fields = 1 + introspection.n_compositional_fields;
std::vector<std::vector<unsigned int>> support_point_index_by_field(n_fields);

// Fill in the support_point_index_by_field data structure. This is a bit complicated...
{
Expand Down Expand Up @@ -1746,10 +1739,11 @@ namespace aspect
combined_support_points,
update_quadrature_points | update_values | update_gradients);

const unsigned int n_q_points = combined_support_points.size();
std::vector<types::global_dof_index> local_dof_indices (dof_handler.get_fe().dofs_per_cell);
MaterialModel::MaterialModelInputs<dim> in(combined_support_points.size(), introspection.n_compositional_fields);
MaterialModel::MaterialModelOutputs<dim> out(combined_support_points.size(), introspection.n_compositional_fields);
HeatingModel::HeatingModelOutputs heating_model_outputs(combined_support_points.size(), introspection.n_compositional_fields);
MaterialModel::MaterialModelInputs<dim> in(n_q_points, introspection.n_compositional_fields);
MaterialModel::MaterialModelOutputs<dim> out(n_q_points, introspection.n_compositional_fields);
HeatingModel::HeatingModelOutputs heating_model_outputs(n_q_points, introspection.n_compositional_fields);

// add reaction rate outputs
material_model->create_additional_named_outputs(out);
Expand All @@ -1766,6 +1760,70 @@ namespace aspect
// some heating models require the additional outputs
heating_model_manager.create_additional_material_model_inputs_and_outputs(in, out);

// We use SUNDIALs ARKode to compute the reactions. Set up the required parameters.
// TODO: Should we change some of these based on the Reaction time step input parameter?
using VectorType = Vector<double>;
SUNDIALS::ARKode<VectorType>::AdditionalData data;
data.initial_time = time;
data.final_time = time + time_step;
data.initial_step_size = 0.001 * time_step;
data.output_period = time_step;
data.minimum_step_size = 1.e-6 * time_step;

// Both tolerances are added, but the composition might become 0.
// We therefore set the absolute tolerance to a very small value.
data.relative_tolerance = 1e-6;
data.absolute_tolerance = 1e-10;

SUNDIALS::ARKode<VectorType> ode(data);

std::vector<std::vector<double>> initial_values_C (n_q_points, std::vector<double> (introspection.n_compositional_fields));
std::vector<double> initial_values_T (n_q_points);

// We have to store all values of the temperature and composition fields in one long vector.
// Create the vector and functions for transferring values between this vector and the material model inputs object.
VectorType fields (n_q_points * n_fields);

auto copy_fields_into_one_vector = [n_q_points,n_fields](const MaterialModel::MaterialModelInputs<dim> &in,
VectorType &fields)
{
for (unsigned int j=0; j<n_q_points; ++j)
for (unsigned int f=0; f<n_fields; ++f)
if (f==0)
fields[j*n_fields+f] = in.temperature[j];
else
fields[j*n_fields+f] = in.composition[j][f-1];
return;
};

auto copy_fields_into_material_model_inputs = [n_q_points,n_fields](const VectorType &fields,
MaterialModel::MaterialModelInputs<dim> &in)
{
for (unsigned int j=0; j<n_q_points; ++j)
for (unsigned int f=0; f<n_fields; ++f)
if (f==0)
in.temperature[j] = fields[j*n_fields+f];
else
in.composition[j][f-1] = fields[j*n_fields+f];
return;
};

auto copy_rates_into_one_vector = [n_q_points,n_fields](const MaterialModel::ReactionRateOutputs<dim> *reaction_out,
const HeatingModel::HeatingModelOutputs &heating_out,
VectorType &rates)
{
for (unsigned int j=0; j<n_q_points; ++j)
for (unsigned int f=0; f<n_fields; ++f)
if (f==0)
rates[j*n_fields+f] = heating_out.rates_of_temperature_change[j];
else
rates[j*n_fields+f] = reaction_out->reaction_rates[j][f-1];
return;
};

unsigned int total_iteration_count = 0;
unsigned int number_of_solves = 0;

// Make a loop first over all cells, than over all reaction time steps, and then over
// all degrees of freedom in each element to compute the reactions. This is possible
// because the reactions only depend on the temperature and composition values at a given
Expand All @@ -1789,36 +1847,49 @@ namespace aspect
fe_values.reinit (cell);
in.reinit(fe_values, cell, introspection, solution);

std::vector<std::vector<double>> accumulated_reactions_C (combined_support_points.size(), std::vector<double> (introspection.n_compositional_fields));
std::vector<double> accumulated_reactions_T (combined_support_points.size());
std::vector<std::vector<double>> accumulated_reactions_C (n_q_points, std::vector<double> (introspection.n_compositional_fields));
std::vector<double> accumulated_reactions_T (n_q_points);

copy_fields_into_one_vector (in, fields);

initial_values_C = in.composition;
initial_values_T = in.temperature;

ode.explicit_function = [&] (const double /*time*/,
const VectorType &y,
VectorType &ydot)
{
copy_fields_into_material_model_inputs (y, in);
material_model->fill_additional_material_model_inputs(in, solution, fe_values, introspection);
material_model->evaluate(in, out);
heating_model_manager.evaluate(in, out, heating_model_outputs);
copy_rates_into_one_vector (reaction_rate_outputs, heating_model_outputs, ydot);
};

// Make the reaction time steps: We have to update the values of compositional fields and the temperature.
// We can reuse the same material model inputs and outputs structure for each reaction time step.
// We store the computed updates to temperature and composition in a separate (accumulated_reactions) vector,
// so that we can later copy it over to the solution vector.
for (unsigned int i=0; i<number_of_reaction_steps; ++i)
{
// Loop over composition element
material_model->fill_additional_material_model_inputs(in, solution, fe_values, introspection);
const unsigned int iteration_count = ode.solve_ode(fields);

material_model->evaluate(in, out);
heating_model_manager.evaluate(in, out, heating_model_outputs);
total_iteration_count += iteration_count;
number_of_solves += 1;

for (unsigned int j=0; j<combined_support_points.size(); ++j)
{
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
// simple forward euler
in.composition[j][c] = in.composition[j][c]
+ reaction_time_step_size * reaction_rate_outputs->reaction_rates[j][c];
accumulated_reactions_C[j][c] += reaction_time_step_size * reaction_rate_outputs->reaction_rates[j][c];
}
in.temperature[j] = in.temperature[j]
+ reaction_time_step_size * heating_model_outputs.rates_of_temperature_change[j];
for (unsigned int j=0; j<n_q_points; ++j)
for (unsigned int f=0; f<n_fields; ++f)
{
if (f==0)
{
in.temperature[j] = fields[j*n_fields+f];
accumulated_reactions_T[j] = in.temperature[j] - initial_values_T[j];
}
else
{
in.composition[j][f-1] = fields[j*n_fields+f];
accumulated_reactions_C[j][f-1] = in.composition[j][f-1] - initial_values_C[j][f-1];
}
}

accumulated_reactions_T[j] += reaction_time_step_size * heating_model_outputs.rates_of_temperature_change[j];
}
}

cell->get_dof_indices (local_dof_indices);
const unsigned int component_idx_T = introspection.component_indices.temperature;
Expand Down Expand Up @@ -1877,8 +1948,14 @@ namespace aspect

initialize_current_linearization_point();

const double average_iteration_count = number_of_solves > 0
?
total_iteration_count / number_of_solves
:
total_iteration_count;

pcout << "in "
<< number_of_reaction_steps
<< average_iteration_count
<< " substep(s)."
<< std::endl;
}
Expand Down

0 comments on commit 314f518

Please sign in to comment.