diff --git a/include/aspect/melt.h b/include/aspect/melt.h index 2294b3946b3..146f4c961ec 100644 --- a/include/aspect/melt.h +++ b/include/aspect/melt.h @@ -222,29 +222,14 @@ namespace aspect void execute(internal::Assembly::Scratch::ScratchBase &scratch, internal::Assembly::CopyData::CopyDataBase &data) const; - }; - - /** - * Compute the residual of the advection system on a single cell in - * the case of melt migration. - */ - template - class MeltAdvectionSystemResidual : public Assemblers::ComputeResidualBase, - public aspect::SimulatorAccess - { - public: - virtual ~MeltAdvectionSystemResidual () {}; - - virtual - std::vector - execute(internal::Assembly::Scratch::ScratchBase &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 &outputs) const; + std::vector + compute_residual(internal::Assembly::Scratch::ScratchBase &scratch) const; }; /** diff --git a/include/aspect/simulator/assemblers/advection.h b/include/aspect/simulator/assemblers/advection.h index 95e21eee43f..644d732daec 100644 --- a/include/aspect/simulator/assemblers/advection.h +++ b/include/aspect/simulator/assemblers/advection.h @@ -44,6 +44,10 @@ namespace aspect void execute(internal::Assembly::Scratch::ScratchBase &scratch, internal::Assembly::CopyData::CopyDataBase &data) const; + + virtual + std::vector + compute_residual(internal::Assembly::Scratch::ScratchBase &scratch) const; }; /** @@ -79,22 +83,6 @@ namespace aspect execute(internal::Assembly::Scratch::ScratchBase &scratch, internal::Assembly::CopyData::CopyDataBase &data) const; }; - - /** - * This class computes the residual for the advection equation at every quadrature - * point of the current cell. - */ - template - class AdvectionSystemResidual : public Assemblers::ComputeResidualBase, - public SimulatorAccess - { - public: - virtual ~AdvectionSystemResidual () {}; - - virtual - std::vector - execute(internal::Assembly::Scratch::ScratchBase &scratch) const; - }; } } diff --git a/include/aspect/simulator/assemblers/interface.h b/include/aspect/simulator/assemblers/interface.h index 3da73f6b98a..de6fca178d0 100644 --- a/include/aspect/simulator/assemblers/interface.h +++ b/include/aspect/simulator/assemblers/interface.h @@ -540,44 +540,27 @@ namespace aspect * execute() is called. */ virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &) 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 - class ComputeResidualBase - { - public: - virtual ~ComputeResidualBase () {} - - virtual - std::vector - execute(internal::Assembly::Scratch::ScratchBase &) 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 &) const {} + virtual + std::vector + compute_residual(internal::Assembly::Scratch::ScratchBase &) 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(); + } }; /** @@ -656,14 +639,6 @@ namespace aspect */ std::vector > > 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 > > advection_system_residual; - /** * A structure that describes what information an assembler function * (listed as one of the assembler objects above) may need to operate. diff --git a/source/simulator/assemblers/advection.cc b/source/simulator/assemblers/advection.cc index a4b212b8bd6..5bb5cc26ece 100644 --- a/source/simulator/assemblers/advection.cc +++ b/source/simulator/assemblers/advection.cc @@ -190,7 +190,7 @@ namespace aspect template std::vector - AdvectionSystemResidual::execute(internal::Assembly::Scratch::ScratchBase &scratch_base) const + AdvectionSystem::compute_residual(internal::Assembly::Scratch::ScratchBase &scratch_base) const { internal::Assembly::Scratch::AdvectionSystem &scratch = dynamic_cast& > (scratch_base); @@ -1154,8 +1154,7 @@ namespace aspect template class AdvectionSystem; \ template class AdvectionSystemBoundaryFace; \ template class AdvectionSystemInteriorFace; \ - template class AdvectionSystemResidual; - + ASPECT_INSTANTIATE(INSTANTIATE) } } diff --git a/source/simulator/assembly.cc b/source/simulator/assembly.cc index 3b0e29b33ac..6e5bea30a1f 100644 --- a/source/simulator/assembly.cc +++ b/source/simulator/assembly.cc @@ -317,11 +317,6 @@ namespace aspect assemblers->advection_system.push_back( std_cxx11::unique_ptr > (melt_advection_system)); - aspect::Assemblers::MeltAdvectionSystemResidual *melt_advection_system_residual - = new aspect::Assemblers::MeltAdvectionSystemResidual(); - - assemblers->advection_system_residual.push_back( - std_cxx11::unique_ptr > (melt_advection_system_residual)); } else { @@ -331,11 +326,6 @@ namespace aspect assemblers->advection_system.push_back( std_cxx11::unique_ptr > (advection_system)); - aspect::Assemblers::AdvectionSystemResidual *advection_system_residual - = new aspect::Assemblers::AdvectionSystemResidual(); - - assemblers->advection_system_residual.push_back( - std_cxx11::unique_ptr > (advection_system_residual)); } if (parameters.use_discontinuous_temperature_discretization || @@ -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); } diff --git a/source/simulator/entropy_viscosity.cc b/source/simulator/entropy_viscosity.cc index c5d00daaf71..ed152e161c2 100644 --- a/source/simulator/entropy_viscosity.cc +++ b/source/simulator/entropy_viscosity.cc @@ -128,11 +128,11 @@ namespace aspect if (advection_field.is_discontinuous(introspection)) return 0.; - std::vector residual = assemblers->advection_system_residual[0]->execute(scratch); + std::vector residual = assemblers->advection_system[0]->compute_residual(scratch); - for (unsigned int i=1; iadvection_system_residual.size(); ++i) + for (unsigned int i=1; iadvection_system.size(); ++i) { - std::vector new_residual = assemblers->advection_system_residual[i]->execute(scratch); + std::vector new_residual = assemblers->advection_system[i]->compute_residual(scratch); for (unsigned int i=0; iadvection_system_residual.size(); ++i) - assemblers->advection_system_residual[i]->create_additional_material_model_outputs(scratch.material_model_outputs); + for (unsigned int i=0; iadvection_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); diff --git a/source/simulator/melt.cc b/source/simulator/melt.cc index 1e767ce9e5b..69aac525685 100644 --- a/source/simulator/melt.cc +++ b/source/simulator/melt.cc @@ -770,10 +770,11 @@ namespace aspect template std::vector - MeltAdvectionSystemResidual:: - execute(internal::Assembly::Scratch::ScratchBase &scratch_base) const + MeltAdvectionSystem:: + compute_residual(internal::Assembly::Scratch::ScratchBase &scratch_base) const { - internal::Assembly::Scratch::AdvectionSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::Scratch::AdvectionSystem &scratch = + dynamic_cast& > (scratch_base); const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; std::vector residuals(n_q_points); @@ -854,30 +855,6 @@ namespace aspect - template - void - MeltAdvectionSystemResidual::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const - { - MeltHandler::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 >() == NULL) - { - outputs.additional_outputs.push_back( - std_cxx11::shared_ptr > - (new MaterialModel::AdditionalMaterialOutputsStokesRHS (n_points))); - } - - Assert(!this->get_parameters().enable_additional_stokes_rhs - || - outputs.template get_additional_output >()->rhs_u.size() - == n_points, ExcInternalError()); - } - - - template void MeltPressureRHSCompatibilityModification:: @@ -1395,7 +1372,6 @@ namespace aspect template class MeltStokesSystem; \ template class MeltStokesSystemBoundary; \ template class MeltAdvectionSystem; \ - template class MeltAdvectionSystemResidual; \ template class MeltPressureRHSCompatibilityModification; \ template class MeltBoundaryTraction; \ } \