Skip to content

Commit

Permalink
Remove Residual base class
Browse files Browse the repository at this point in the history
  • Loading branch information
@gassmoeller
gassmoeller committed Nov 20, 2017
1 parent 81b0eef commit 11b1953
Show file tree
Hide file tree
Showing 7 changed files with 37 additions and 125 deletions.
23 changes: 4 additions & 19 deletions include/aspect/melt.h
View file
Expand Up @@ -222,29 +222,14 @@ namespace aspect
void
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
internal::Assembly::CopyData::CopyDataBase<dim> &data) const;
};

/**
* Compute the residual of the advection system on a single cell in
* the case of melt migration.
*/
template <int dim>
class MeltAdvectionSystemResidual : public Assemblers::ComputeResidualBase<dim>,
public aspect::SimulatorAccess<dim>
{
public:
virtual ~MeltAdvectionSystemResidual () {};

virtual
std::vector<double>
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch) const;

/**
* Attach melt outputs.
* Compute the residual of the advection system on a single cell in
* the case of melt migration.
*/
virtual
void
create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const;
std::vector<double>
compute_residual(internal::Assembly::Scratch::ScratchBase<dim> &scratch) const;
};

/**
Expand Down
20 changes: 4 additions & 16 deletions include/aspect/simulator/assemblers/advection.h
View file
Expand Up @@ -44,6 +44,10 @@ namespace aspect
void
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
internal::Assembly::CopyData::CopyDataBase<dim> &data) const;

virtual
std::vector<double>
compute_residual(internal::Assembly::Scratch::ScratchBase<dim> &scratch) const;
};

/**
Expand Down Expand Up @@ -79,22 +83,6 @@ namespace aspect
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
internal::Assembly::CopyData::CopyDataBase<dim> &data) const;
};

/**
* This class computes the residual for the advection equation at every quadrature
* point of the current cell.
*/
template <int dim>
class AdvectionSystemResidual : public Assemblers::ComputeResidualBase<dim>,
public SimulatorAccess<dim>
{
public:
virtual ~AdvectionSystemResidual () {};

virtual
std::vector<double>
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch) const;
};
}
}

Expand Down
61 changes: 18 additions & 43 deletions include/aspect/simulator/assemblers/interface.h
View file
Expand Up @@ -540,44 +540,27 @@ namespace aspect
* execute() is called.
*/
virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &) const {}
};

/**
* A base class for objects that implement the computation of residuals
* for the advection equation. Just like the assemblers itself, the residual
* that we use to compute the necessary entropy viscosity depend on the
* equation (i.e. which terms are actually included in the
* equation). Thus different objects compute different residuals (i.e.
* the residual for a melt advection equation looks different from the
* residual for a passive compositional field).
* The interface of this class differs from the AssemblerBase class, because
* it currently returns a vector of residuals, instead of filling a copy
* data object. TODO: This could be easily changed. Is it useful?
*
* The point of this class is primarily to be able to store
* pointers to base objects in a vector. The objects are usually created
* in Simulator::set_assemblers() and destroyed in the destructor of
* the Simulator object.
*/
template <int dim>
class ComputeResidualBase
{
public:
virtual ~ComputeResidualBase () {}

virtual
std::vector<double>
execute(internal::Assembly::Scratch::ScratchBase<dim> &) const = 0;

/**
* This function gets called if a MaterialModelOutputs is created
* and allows the assembler to attach AdditionalOutputs. The
* function might be called more than once for a
* MaterialModelOutput, so it is recommended to check if
* get_additional_output() returns an instance before adding a new
* one to the additional_outputs vector.
* A required function for objects that implement the assembly of terms
* in an equation that requires the computation of residuals
* (in particular the advection equation in ASPECT).
* Just like the assemblers itself, the residual
* that we use to compute the necessary entropy viscosity depend on the
* equation (i.e. which terms are actually included in the
* equation). Thus different objects compute different residuals (i.e.
* the residual for a melt advection equation looks different from the
* residual for a passive compositional field).
*/
virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &) const {}
virtual
std::vector<double>
compute_residual(internal::Assembly::Scratch::ScratchBase<dim> &) const
{
AssertThrow(false, ExcMessage("This function should either be implemented "
"(if this is an assembler that has to compute a residual, or not be called "
"(if this is an assembler that has not to compute a residual)."));
return std::vector<double>();
}
};

/**
Expand Down Expand Up @@ -656,14 +639,6 @@ namespace aspect
*/
std::vector<std_cxx11::unique_ptr<Assemblers::Interface<dim> > > advection_system_on_interior_face;

/**
* A vector of pointers containing all assemblers that compute the
* advection system residual.
* These assemblers are called once per cell, and are currently only
* used to compute the artificial viscosity stabilization.
*/
std::vector<std_cxx11::unique_ptr<Assemblers::ComputeResidualBase<dim> > > advection_system_residual;

/**
* A structure that describes what information an assembler function
* (listed as one of the assembler objects above) may need to operate.
Expand Down
5 changes: 2 additions & 3 deletions source/simulator/assemblers/advection.cc
View file
Expand Up @@ -190,7 +190,7 @@ namespace aspect

template <int dim>
std::vector<double>
AdvectionSystemResidual<dim>::execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base) const
AdvectionSystem<dim>::compute_residual(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base) const
{
internal::Assembly::Scratch::AdvectionSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::AdvectionSystem<dim>& > (scratch_base);

Expand Down Expand Up @@ -1154,8 +1154,7 @@ namespace aspect
template class AdvectionSystem<dim>; \
template class AdvectionSystemBoundaryFace<dim>; \
template class AdvectionSystemInteriorFace<dim>; \
template class AdvectionSystemResidual<dim>;


ASPECT_INSTANTIATE(INSTANTIATE)
}
}
11 changes: 0 additions & 11 deletions source/simulator/assembly.cc
View file
Expand Up @@ -317,11 +317,6 @@ namespace aspect
assemblers->advection_system.push_back(
std_cxx11::unique_ptr<aspect::Assemblers::MeltAdvectionSystem<dim> > (melt_advection_system));

aspect::Assemblers::MeltAdvectionSystemResidual<dim> *melt_advection_system_residual
= new aspect::Assemblers::MeltAdvectionSystemResidual<dim>();

assemblers->advection_system_residual.push_back(
std_cxx11::unique_ptr<aspect::Assemblers::MeltAdvectionSystemResidual<dim> > (melt_advection_system_residual));
}
else
{
Expand All @@ -331,11 +326,6 @@ namespace aspect
assemblers->advection_system.push_back(
std_cxx11::unique_ptr<aspect::Assemblers::AdvectionSystem<dim> > (advection_system));

aspect::Assemblers::AdvectionSystemResidual<dim> *advection_system_residual
= new aspect::Assemblers::AdvectionSystemResidual<dim>();

assemblers->advection_system_residual.push_back(
std_cxx11::unique_ptr<aspect::Assemblers::AdvectionSystemResidual<dim> > (advection_system_residual));
}

if (parameters.use_discontinuous_temperature_discretization ||
Expand Down Expand Up @@ -380,7 +370,6 @@ namespace aspect
initialize_simulator(*this,assemblers->advection_system);
initialize_simulator(*this,assemblers->advection_system_on_boundary_face);
initialize_simulator(*this,assemblers->advection_system_on_interior_face);
initialize_simulator(*this,assemblers->advection_system_residual);
}


Expand Down
10 changes: 5 additions & 5 deletions source/simulator/entropy_viscosity.cc
View file
Expand Up @@ -128,11 +128,11 @@ namespace aspect
if (advection_field.is_discontinuous(introspection))
return 0.;

std::vector<double> residual = assemblers->advection_system_residual[0]->execute(scratch);
std::vector<double> residual = assemblers->advection_system[0]->compute_residual(scratch);

for (unsigned int i=1; i<assemblers->advection_system_residual.size(); ++i)
for (unsigned int i=1; i<assemblers->advection_system.size(); ++i)
{
std::vector<double> new_residual = assemblers->advection_system_residual[i]->execute(scratch);
std::vector<double> new_residual = assemblers->advection_system[i]->compute_residual(scratch);
for (unsigned int i=0; i<residual.size(); ++i)
residual[i] += new_residual[i];
}
Expand Down Expand Up @@ -395,8 +395,8 @@ namespace aspect
}
scratch.material_model_inputs.current_cell = cell;

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

material_model->evaluate(scratch.material_model_inputs,scratch.material_model_outputs);

Expand Down
32 changes: 4 additions & 28 deletions source/simulator/melt.cc
View file
Expand Up @@ -770,10 +770,11 @@ namespace aspect

template <int dim>
std::vector<double>
MeltAdvectionSystemResidual<dim>::
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base) const
MeltAdvectionSystem<dim>::
compute_residual(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base) const
{
internal::Assembly::Scratch::AdvectionSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::AdvectionSystem<dim>& > (scratch_base);
internal::Assembly::Scratch::AdvectionSystem<dim> &scratch =
dynamic_cast<internal::Assembly::Scratch::AdvectionSystem<dim>& > (scratch_base);

const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
std::vector<double> residuals(n_q_points);
Expand Down Expand Up @@ -854,30 +855,6 @@ namespace aspect



template <int dim>
void
MeltAdvectionSystemResidual<dim>::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const
{
MeltHandler<dim>::create_material_model_outputs(outputs);

const unsigned int n_points = outputs.viscosities.size();

if (this->get_parameters().enable_additional_stokes_rhs
&& outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >() == NULL)
{
outputs.additional_outputs.push_back(
std_cxx11::shared_ptr<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >
(new MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> (n_points)));
}

Assert(!this->get_parameters().enable_additional_stokes_rhs
||
outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >()->rhs_u.size()
== n_points, ExcInternalError());
}



template <int dim>
void
MeltPressureRHSCompatibilityModification<dim>::
Expand Down Expand Up @@ -1395,7 +1372,6 @@ namespace aspect
template class MeltStokesSystem<dim>; \
template class MeltStokesSystemBoundary<dim>; \
template class MeltAdvectionSystem<dim>; \
template class MeltAdvectionSystemResidual<dim>; \
template class MeltPressureRHSCompatibilityModification<dim>; \
template class MeltBoundaryTraction<dim>; \
} \
Expand Down

0 comments on commit 11b1953

Please sign in to comment.