diff --git a/cookbooks/inner_core_convection/inner_core_assembly.cc b/cookbooks/inner_core_convection/inner_core_assembly.cc index 26a5a280f7b..ed5d1b41489 100644 --- a/cookbooks/inner_core_convection/inner_core_assembly.cc +++ b/cookbooks/inner_core_convection/inner_core_assembly.cc @@ -1,7 +1,7 @@ #include #include #include -#include +#include #include #include @@ -40,7 +40,7 @@ namespace aspect */ template class PhaseBoundaryAssembler : - public aspect::internal::Assembly::Assemblers::AssemblerBase, + public aspect::Assemblers::Interface, public SimulatorAccess { @@ -48,20 +48,21 @@ namespace aspect virtual void - phase_change_boundary_conditions (const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); const unsigned int n_q_points = scratch.face_finite_element_values.n_quadrature_points; //assemble force terms for the matrix for all boundary faces - if (cell->face(face_no)->at_boundary()) + if (scratch.cell->face(scratch.face_number)->at_boundary()) { - scratch.face_finite_element_values.reinit (cell, face_no); + scratch.face_finite_element_values.reinit (scratch.cell, scratch.face_number); for (unsigned int q=0; q void set_assemblers_phase_boundary(const SimulatorAccess &simulator_access, - internal::Assembly::AssemblerLists &assemblers, - std::vector > > &assembler_objects) + Assemblers::Manager &assemblers) { AssertThrow (dynamic_cast*> (&simulator_access.get_material_model()) != 0, ExcMessage ("The phase boundary assembler can only be used with the " "material model 'inner core material'!")); - std_cxx11::shared_ptr > phase_boundary_assembler - (new PhaseBoundaryAssembler()); - assembler_objects.push_back (phase_boundary_assembler); - - // add the terms for phase change boundary conditions - assemblers.local_assemble_stokes_system_on_boundary_face - .connect (std_cxx11::bind(&PhaseBoundaryAssembler::phase_change_boundary_conditions, - std_cxx11::cref (*phase_boundary_assembler), - std_cxx11::_1, - std_cxx11::_2, - // discard pressure_scaling, - // discard rebuild_stokes_matrix, - std_cxx11::_5, - std_cxx11::_6)); - - + PhaseBoundaryAssembler *phase_boundary_assembler = new PhaseBoundaryAssembler(); + assemblers.stokes_system_on_boundary_face.push_back (std_cxx11::unique_ptr > (phase_boundary_assembler)); } } diff --git a/doc/update_source_files_to_2.0.0.sed b/doc/update_source_files_to_2.0.0.sed index 6d547f083a3..d947268b7ef 100644 --- a/doc/update_source_files_to_2.0.0.sed +++ b/doc/update_source_files_to_2.0.0.sed @@ -62,3 +62,11 @@ s/tracer particle/particle/g s/tracer/particle/g s/Tracer/Particle/g s/TRACER/PARTICLE/g + +# Rename assembler base class +s/internal::Assembly::AssemblerLists/Assemblers::AssemblerLists/g +s/internal::Assembly::Assemblers/Assemblers/g +s/struct AssemblerLists/class Manager/g +s/AssemblerLists/Manager/g +s/AssemblerBase/Interface/g +s:assembly.h:simulator/assemblers/interface.h:g diff --git a/include/aspect/assembly.h b/include/aspect/assembly.h deleted file mode 100644 index e08f0b41358..00000000000 --- a/include/aspect/assembly.h +++ /dev/null @@ -1,811 +0,0 @@ -/* - Copyright (C) 2016 - 2017 by the authors of the ASPECT code. - - This file is part of ASPECT. - - ASPECT is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2, or (at your option) - any later version. - - ASPECT is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with ASPECT; see the file LICENSE. If not see - . -*/ - - -#ifndef _aspect_assembly_h -#define _aspect_assembly_h - -#include -#include -#include - -#include - - -namespace aspect -{ - using namespace dealii; - - template - class Simulator; - - struct AdvectionField; - - namespace internal - { - namespace Assembly - { - namespace Scratch - { - template - struct StokesPreconditioner - { - StokesPreconditioner (const FiniteElement &finite_element, - const Quadrature &quadrature, - const Mapping &mapping, - const UpdateFlags update_flags, - const unsigned int n_compositional_fields, - const unsigned int stokes_dofs_per_cell, - const bool add_compaction_pressure); - StokesPreconditioner (const StokesPreconditioner &data); - - virtual ~StokesPreconditioner (); - - FEValues finite_element_values; - - std::vector local_dof_indices; - std::vector dof_component_indices; - std::vector > grads_phi_u; - std::vector div_phi_u; - std::vector phi_p; - std::vector phi_p_c; - std::vector > grad_phi_p; - - /** - * Material model inputs and outputs computed at the current - * linearization point. - */ - MaterialModel::MaterialModelInputs material_model_inputs; - MaterialModel::MaterialModelOutputs material_model_outputs; - }; - - - - // We derive the StokesSystem scratch array from the - // StokesPreconditioner array. We do this because all the objects that - // are necessary for the assembly of the preconditioner are also - // needed for the actual system matrix and right hand side, plus some - // extra data that we need for the time stepping and traction boundaries - // on the right hand side. - template - struct StokesSystem : public StokesPreconditioner - { - StokesSystem (const FiniteElement &finite_element, - const Mapping &mapping, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const UpdateFlags update_flags, - const UpdateFlags face_update_flags, - const unsigned int n_compositional_fields, - const unsigned int stokes_dofs_per_cell, - const bool add_compaction_pressure, - const bool use_reference_profile); - - StokesSystem (const StokesSystem &data); - - FEFaceValues face_finite_element_values; - - std::vector > phi_u; - std::vector > grads_phi_u; - std::vector > velocity_values; - std::vector velocity_divergence; - - /** - * Material model inputs and outputs computed at the current - * linearization point. - * - * In contrast to the variables above, the following two - * variables are used in the assembly at quadrature points - * on faces, not on cells. - */ - MaterialModel::MaterialModelInputs face_material_model_inputs; - MaterialModel::MaterialModelOutputs face_material_model_outputs; - - /** - * In some approximations of the Stokes equations the density used - * for the mass conservation (/continuity) equation is some form - * of reference density, while the density used for calculating - * the buoyancy force is the full density. In case such a formulation - * is used the reference density (and its derivative in depth - * direction) is queried from the adiabatic conditions plugin - * and is stored in these variables. - */ - std::vector reference_densities; - std::vector reference_densities_depth_derivative; - }; - - - - template - struct AdvectionSystem - { - AdvectionSystem (const FiniteElement &finite_element, - const FiniteElement &advection_element, - const Mapping &mapping, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const UpdateFlags update_flags, - const UpdateFlags face_update_flags, - const unsigned int n_compositional_fields); - AdvectionSystem (const AdvectionSystem &data); - - FEValues finite_element_values; - - std_cxx11::shared_ptr > face_finite_element_values; - std_cxx11::shared_ptr > neighbor_face_finite_element_values; - std_cxx11::shared_ptr > subface_finite_element_values; - - std::vector local_dof_indices; - - /** - * Variables describing the values and gradients of the - * shape functions at the quadrature points, as they are - * used in the advection assembly function. note that the sizes - * of these arrays are equal to the number of shape functions - * corresponding to the advected field (and not of the overall - * field!), and that they are also correspondingly indexed. - */ - std::vector phi_field; - std::vector > grad_phi_field; - std::vector face_phi_field; - std::vector > face_grad_phi_field; - std::vector neighbor_face_phi_field; - std::vector > neighbor_face_grad_phi_field; - - std::vector > old_velocity_values; - std::vector > old_old_velocity_values; - - std::vector old_pressure; - std::vector old_old_pressure; - std::vector > old_pressure_gradients; - std::vector > old_old_pressure_gradients; - - std::vector > old_strain_rates; - std::vector > old_old_strain_rates; - - std::vector old_temperature_values; - std::vector old_old_temperature_values; - - std::vector old_field_values; - std::vector old_old_field_values; - std::vector > old_field_grads; - std::vector > old_old_field_grads; - std::vector old_field_laplacians; - std::vector old_old_field_laplacians; - - std::vector > old_composition_values; - std::vector > old_old_composition_values; - - std::vector current_temperature_values; - std::vector > current_velocity_values; - std::vector > face_current_velocity_values; - std::vector > mesh_velocity_values; - std::vector > face_mesh_velocity_values; - - std::vector > current_strain_rates; - std::vector > current_composition_values; - std::vector current_velocity_divergences; - - /** - * Material model inputs and outputs computed at the current - * linearization point. - */ - MaterialModel::MaterialModelInputs material_model_inputs; - MaterialModel::MaterialModelOutputs material_model_outputs; - - MaterialModel::MaterialModelInputs face_material_model_inputs; - MaterialModel::MaterialModelOutputs face_material_model_outputs; - - MaterialModel::MaterialModelInputs neighbor_face_material_model_inputs; - MaterialModel::MaterialModelOutputs neighbor_face_material_model_outputs; - - /** - * Heating model outputs computed at the quadrature points of the - * current cell at the time of the current linearization point. - * As explained in the class documentation of - * HeatingModel::HeatingModelOutputs each term contains the sum of all - * enabled heating mechanism contributions. - */ - HeatingModel::HeatingModelOutputs heating_model_outputs; - HeatingModel::HeatingModelOutputs face_heating_model_outputs; - HeatingModel::HeatingModelOutputs neighbor_face_heating_model_outputs; - }; - - - - } - - - // The CopyData arrays are similar to the - // Scratch arrays. They provide a - // constructor, a copy operation, and - // some arrays for local matrix, local - // vectors and the relation between local - // and global degrees of freedom (a.k.a. - // local_dof_indices). - namespace CopyData - { - template - struct StokesPreconditioner - { - StokesPreconditioner (const unsigned int stokes_dofs_per_cell); - StokesPreconditioner (const StokesPreconditioner &data); - - virtual ~StokesPreconditioner (); - - FullMatrix local_matrix; - std::vector local_dof_indices; - - /** - * Extract the values listed in @p all_dof_indices only if - * it corresponds to the Stokes component and copy it to the variable - * local_dof_indices declared above in the same class as this function - */ - void extract_stokes_dof_indices(const std::vector &all_dof_indices, - const Introspection &introspection, - const dealii::FiniteElement &finite_element); - }; - - - - template - struct StokesSystem : public StokesPreconditioner - { - StokesSystem (const unsigned int stokes_dofs_per_cell, - const bool do_pressure_rhs_compatibility_modification); - StokesSystem (const StokesSystem &data); - - Vector local_rhs; - Vector local_pressure_shape_function_integrals; - }; - - - - template - struct AdvectionSystem - { - /** - * Constructor. - * - * @param finite_element The element that describes the field for - * which we are trying to assemble a linear system. Not - * the global finite element. - * @param field_is_discontinuous If true, the field is a DG element. - */ - AdvectionSystem (const FiniteElement &finite_element, - const bool field_is_discontinuous); - - /** - * Copy constructor. - */ - AdvectionSystem (const AdvectionSystem &data); - - /** - * Local contributions to the global matrix and right hand side - * that correspond only to the variables listed in local_dof_indices - */ - FullMatrix local_matrix; - /** Local contributions to the global matrix from the face terms in the - * discontinuous Galerkin method. The vectors are of length - * GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell - * so as to hold one matrix for each possible face or subface of the cell. - * The discontinuous Galerkin bilinear form contains terms arising from internal - * (to the cell) values and external (to the cell) values. - * _int_ext and ext_int hold the terms arising from the pairing between a cell - * and its neighbor, while _ext_ext is the pairing of the neighbor's dofs with - * themselves. In the continuous Galerkin case, these are unused, and set to size zero. - **/ - std::vector > local_matrices_int_ext; - std::vector > local_matrices_ext_int; - std::vector > local_matrices_ext_ext; - Vector local_rhs; - - /** Denotes which face matrices have actually been assembled in the DG field - * assembly. Entries for matrices not used (for example, those corresponding - * to non-existent subfaces; or faces being assembled by the neighboring cell) - * are set to false. - **/ - std::vector assembled_matrices; - - /** - * Indices of those degrees of freedom that actually correspond - * to the temperature or compositional field. since this structure - * is used to represent just contributions to the advection - * systems, there will be no contributions to other parts of the - * system and consequently, we do not need to list here indices - * that correspond to velocity or pressure degrees (or, in fact - * any other variable outside the block we are currently considering) - */ - std::vector local_dof_indices; - /** Indices of the degrees of freedom corresponding to the temperature - * or composition field on all possible neighboring cells. This is used - * in the discontinuous Galerkin method. The outer std::vector has - * length GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell, - * and has size zero if in the continuous Galerkin case. - **/ - std::vector > neighbor_dof_indices; - }; - } - - /** - * A namespace for the definition of sub-namespaces in which we - * implement assemblers for the various terms in the linear systems - * ASPECT solves. - */ - namespace Assemblers - { - /** - * A base class for objects that implement assembly - * operations. - * - * The point of this class is primarily so that we can store - * pointers to such objects in a list. The objects are created - * in Simulator::set_assemblers() (which is the only place where - * we know their actual type) and destroyed in the destructor of - * class Simulation. - */ - template - class AssemblerBase - { - public: - virtual ~AssemblerBase () {} - - /** - * 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. - */ - virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &) const {} - }; - } - - /** - * A structure that consists of member variables representing - * each one or more functions that need to be called when - * assembling right hand side vectors, matrices, or complete linear - * systems. We use this approach in order to support the following - * cases: - * - Assembling different formulations: When assembling either the - * full equations or only the Boussinesq approximation (to give just - * two examples), one needs different terms. This could be achieved - * using a large number of switch or if - * statements in the code, or one could encapsulate each equation - * or approximation into a collection of functions for this particular - * purpose. The approach chosen here in essence allows the - * implementation of each set of equations in its own scope, and we - * then just need to store a pointer to the function that - * assembles the Stokes system (for example) for the selected - * approximation. The pointer to this function is stored in the - * appropriate member variable of this class. - * - Sometimes, we want to assemble a number of terms that build on - * each other. An example is the addition of free boundary terms - * to the Stokes matrix. Rather than having to "know" in one - * place about all of the terms that need to be assembled, - * we simply add the function that computes these terms as - * another "slot" to the appropriate "signal" declared in this - * class. - */ - template - struct AssemblerLists - { - /** - * A structure used to accumulate the results of the - * compute_advection_system_residual slot functions below. - * It takes an iterator range of vectors and returns the element-wise - * sum of the values. - */ - template - struct ResidualWeightSum - { - typedef VectorType result_type; - - template - VectorType operator()(InputIterator first, InputIterator last) const - { - // If there are no slots to call, just return the - // default-constructed value - if (first == last) - return VectorType(); - - VectorType sum = *first++; - while (first != last) - { - Assert(sum.size() == first->size(), - ExcMessage("Invalid vector length for residual summation.")); - - for (unsigned int i = 0; i < first->size(); ++i) - sum[i] += (*first)[i]; - - ++first; - } - - return sum; - } - }; - - - /** - * A signal that is called from Simulator::local_assemble_stokes_preconditioner() - * and whose slots are supposed to assemble terms that together form the - * Stokes preconditioner matrix. - * - * The arguments to the slots are as follows: - * - The Simulator::pressure_scaling value used to scale velocity - * and pressure components against each other. - * - The scratch object in which temporary data is stored that - * assemblers may need. - * - The copy object into which assemblers add up their contributions. - */ - boost::signals2::signal &, - internal::Assembly::CopyData::StokesPreconditioner &)> local_assemble_stokes_preconditioner; - - /** - * A signal that is called from Simulator::local_assemble_stokes_system() - * and whose slots are supposed to assemble terms that together form the - * Stokes system matrix and right hand side. - * - * The arguments to the slots are as follows: - * - The cell on which we currently assemble. - * - The Simulator::pressure_scaling value used to scale velocity - * and pressure components against each other. - * - Whether or not to actually assemble the matrix. If @p false, - * then only assemble the right hand side of the Stokes system. - * - The scratch object in which temporary data is stored that - * assemblers may need. - * - The copy object into which assemblers add up their contributions. - */ - boost::signals2::signal::active_cell_iterator &, - const double, - const bool, - internal::Assembly::Scratch::StokesSystem &, - internal::Assembly::CopyData::StokesSystem &)> local_assemble_stokes_system; - - /** - * A signal that is called from Simulator::local_assemble_stokes_system() - * and whose slots are supposed to assemble terms that together form the - * Stokes system matrix and right hand side. This signal is called - * once for each boundary face. - * - * The arguments to the slots are as follows: - * - The cell on which we currently assemble. - * - The number of the face on which we intend to assemble. This - * face (of the current cell) will be at the boundary of the - * domain. - * - The Simulator::pressure_scaling value used to scale velocity - * and pressure components against each other. - * - Whether or not to actually assemble the matrix. If @p false, - * then only assemble the right hand side of the Stokes system. - * - The scratch object in which temporary data is stored that - * assemblers may need. - * - The copy object into which assemblers add up their contributions. - */ - boost::signals2::signal::active_cell_iterator &, - const unsigned int, - const double, - const bool, - internal::Assembly::Scratch::StokesSystem &, - internal::Assembly::CopyData::StokesSystem &)> local_assemble_stokes_system_on_boundary_face; - - /** - * A signal that is called from Simulator::local_assemble_advection_system() - * and whose slots are supposed to assemble terms that together form the - * Advection system matrix and right hand side. - * - * The arguments to the slots are as follows: - * - The cell on which we currently assemble. - * - The advection field that is currently assembled. - * - The artificial viscosity computed for the current cell - * - The scratch object in which temporary data is stored that - * assemblers may need. - * - The copy object into which assemblers add up their contributions. - */ - boost::signals2::signal::active_cell_iterator &, - const typename Simulator::AdvectionField &, - const double, - internal::Assembly::Scratch::AdvectionSystem &, - internal::Assembly::CopyData::AdvectionSystem &)> local_assemble_advection_system; - - /** - * A signal that is called from Simulator::local_assemble_advection_system() - * and whose slots are supposed to compute an advection system residual - * according to the terms their assemblers implement. The residuals are - * computed and returned in a std::vector with one entry per - * quadrature point of this cell. If more than one - * assembler is connected the residuals are simply added up by the - * ResidualWeightSum structure. - * - * The arguments to the slots are as follows: - * - The cell on which we currently assemble. - * - The advection field that is currently assembled. - * - The scratch object in which temporary data is stored that - * assemblers may need. - * The return value of each slot is: - * - A std::vector containing one residual entry per quadrature - * point. - */ - boost::signals2::signal (const typename DoFHandler::active_cell_iterator &, - const typename Simulator::AdvectionField &, - internal::Assembly::Scratch::AdvectionSystem &), - ResidualWeightSum > > compute_advection_system_residual; - - /** - * A signal that is called from Simulator::local_assemble_advection_system() - * and whose slots are supposed to assemble terms that together form the - * boundary contribution of the Advection system matrix and right hand side. This signal is called - * once for each boundary face. - * - * The arguments to the slots are as follows: - * - The cell on which we currently assemble. - * - The number of the face on which we intend to assemble. This - * face (of the current cell) will be at the boundary of the - * domain. - * - The advection field that is currently assembled. - * - The scratch object in which temporary data is stored that - * assemblers may need. - * - The copy object into which assemblers add up their contributions. - */ - boost::signals2::signal::active_cell_iterator &, - const unsigned int, - const typename Simulator::AdvectionField &, - internal::Assembly::Scratch::AdvectionSystem &, - internal::Assembly::CopyData::AdvectionSystem &)> local_assemble_advection_system_on_boundary_face; - - /** - * A signal that is called from Simulator::local_assemble_advection_system() - * and whose slots are supposed to assemble terms that together form the - * interior face contribution of the Advection system matrix and right - * hand side. This signal is called once for each interior face. - * - * The arguments to the slots are as follows: - * - The cell on which we currently assemble. - * - The number of the face (of the current cell) on which we intend to - * assemble. - * - The advection field that is currently assembled. - * - The scratch object in which temporary data is stored that - * assemblers may need. - * - The copy object into which assemblers add up their contributions. - */ - boost::signals2::signal::active_cell_iterator &, - const unsigned int, - const typename Simulator::AdvectionField &, - internal::Assembly::Scratch::AdvectionSystem &, - internal::Assembly::CopyData::AdvectionSystem &)> local_assemble_advection_system_on_interior_face; - - /** - * A structure that describes what information an assembler function - * (listed as one of the signals/slots above) may need to operate. - * - * There are a number of pieces of information that are always - * assumed to be needed. For example, the Stokes and advection - * assemblers will always need to have access to the material - * model outputs. But the Stokes assembler may or may not need - * access to material model output for quadrature points on faces. - * - * These properties are all preset in a conservative way - * (i.e., disabled) in the constructor of this class, but can - * be enabled in Simulator::set_assemblers() when adding - * individual assemblers. Functions such as - * Simulator::local_assemble_stokes_preconditioner(), - * Simulator::local_assemble_stokes_system() will then query - * these flags to determine whether something has to be - * initialized for at least one of the assemblers they call. - */ - struct Properties - { - /** - * Constructor. Disable all properties as described in the - * class documentation. - */ - Properties (); - - /** - * Whether or not at least one of the assembler slots in - * a signal requires the initialization and re-computation of - * a MaterialModelOutputs object for each face. This - * property is only relevant to assemblers that operate on - * boundary faces. - */ - bool need_face_material_model_data; - - /** - * Whether or not at least one of the assembler slots in a - * signal requires the evaluation of the FEFaceValues object. - */ - bool need_face_finite_element_evaluation; - - /** - * Whether or not at least one of the assembler slots in a - * signal requires the computation of the viscosity. - */ - bool need_viscosity; - - /** - * A list of FEValues UpdateFlags that are necessary for - * a given operation. Assembler objects may add to this list - * as necessary; it will be initialized with a set of - * "default" flags that will always be set. - */ - UpdateFlags needed_update_flags; - }; - - /** - * A list of properties of the various types of assemblers. - * These property lists are set in Simulator::set_assemblers() - * where we add individual functions to the signals above. - */ - Properties stokes_preconditioner_assembler_properties; - Properties stokes_system_assembler_properties; - Properties stokes_system_assembler_on_boundary_face_properties; - std::vector advection_system_assembler_properties; - std::vector advection_system_assembler_on_face_properties; - }; - } - } - - namespace Assemblers - { - /** - * A class containing the functions to assemble the different terms of the advection equation. - */ - template - class AdvectionAssembler : public aspect::internal::Assembly::Assemblers::AssemblerBase, - public SimulatorAccess - { - public: - /** - * This function assembles the terms for the matrix and right-hand-side of the advection - * equation for the current cell. - */ - void - local_assemble_advection_system (const typename Simulator::AdvectionField &advection_field, - const double artificial_viscosity, - internal::Assembly::Scratch::AdvectionSystem &scratch, - internal::Assembly::CopyData::AdvectionSystem &data) const; - - /** - * This function computes the residual for the advection equation at every quadrature - * point of the current cell and returns the residual in a vector of doubles. - */ - std::vector - compute_advection_system_residual(const typename Simulator::AdvectionField &advection_field, - internal::Assembly::Scratch::AdvectionSystem &scratch) const; - - /** - * This function assembles the face terms for the matrix and right-hand-side of the advection - * equation for a face at the boundary of the domain. - */ - void - local_assemble_discontinuous_advection_boundary_face_terms(const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - const typename Simulator::AdvectionField &advection_field, - internal::Assembly::Scratch::AdvectionSystem &scratch, - internal::Assembly::CopyData::AdvectionSystem &data) const; - - /** - * This function assembles the face terms for the matrix and right-hand-side of the advection - * equation for a face in the interior of the domain. - */ - void - local_assemble_discontinuous_advection_interior_face_terms(const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - const typename Simulator::AdvectionField &advection_field, - internal::Assembly::Scratch::AdvectionSystem &scratch, - internal::Assembly::CopyData::AdvectionSystem &data) const; - }; - - /** - * A class containing the functions to assemble the different terms of the Stokes equations. - */ - template - class StokesAssembler : public aspect::internal::Assembly::Assemblers::AssemblerBase, - public SimulatorAccess - { - public: - - /** - * Create AdditionalMaterialOutputsStokesRHS if we need to do so. - */ - virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const; - - /** - * This function assembles the terms of the Stokes preconditioner matrix for the current cell. - */ - void - preconditioner (const double pressure_scaling, - internal::Assembly::Scratch::StokesPreconditioner &scratch, - internal::Assembly::CopyData::StokesPreconditioner &data) const; - - /** - * This function assembles the terms of the Stokes preconditioner matrix for the current cell. - */ - void - compressible_preconditioner (const double pressure_scaling, - internal::Assembly::Scratch::StokesPreconditioner &scratch, - internal::Assembly::CopyData::StokesPreconditioner &data) const; - - /** - * This function assembles the terms for the matrix and right-hand-side of the incompressible - * Stokes equation for the current cell. - */ - void - incompressible_terms (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const; - - /** - * This function assembles the term that arises in the viscosity term of Stokes matrix for - * compressible models, because the divergence of the velocity is not longer zero. - */ - void - compressible_strain_rate_viscosity_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const; - - /** - * This function assembles the right-hand-side term of the Stokes equation - * that is caused by the compressibility in the mass conservation equation. - * This function approximates this term as - * $- \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u}$ - */ - void - reference_density_compressibility_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const; - - /** - * This function assembles the compressibility term of the Stokes equation - * that is caused by the compressibility in the mass conservation equation. - * It includes this term implicitly in the matrix, - * which is therefore not longer symmetric. - * This function approximates this term as - * $ - \nabla \mathbf{u} - \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$ - */ - void - implicit_reference_density_compressibility_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const; - - /** - * This function assembles the right-hand-side term of the Stokes equation - * that is caused by the compressibility in the mass conservation equation. - * This function approximates this term as - * $ - \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$ - */ - void - isothermal_compression_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const; - }; - } -} - -#endif diff --git a/include/aspect/free_surface.h b/include/aspect/free_surface.h index 96fccab4e2c..0e87ac85c34 100644 --- a/include/aspect/free_surface.h +++ b/include/aspect/free_surface.h @@ -23,12 +23,30 @@ #define _aspect_free_surface_h #include +#include namespace aspect { using namespace dealii; - template + namespace Assemblers + { + /** + * Apply stabilization to a cell of the system matrix. The + * stabilization is only added to cells on a free surface. The + * scheme is based on that of Kaus et. al., 2010. Called during + * assembly of the system matrix. + */ + template + class ApplyStabilization: public Assemblers::Interface, + public SimulatorAccess + { + void execute (internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + } + + template class FreeSurfaceHandler { public: @@ -60,16 +78,6 @@ namespace aspect */ void setup_dofs(); - /** - * Apply stabilization to a cell of the system matrix. The - * stabilization is only added to cells on a free surface. The - * scheme is based on that of Kaus et. al., 2010. Called during - * assembly of the system matrix. - */ - void apply_stabilization (const typename DoFHandler::active_cell_iterator &cell, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data); - /** * Declare parameters for the free surface handling. */ @@ -81,6 +89,14 @@ namespace aspect */ void parse_parameters (ParameterHandler &prm); + + /** + * Return the chosen stabilization term. See + * Kaus et al 2010 for details on the meaning of + * this term. + */ + double get_stabilization_term () const; + private: /** * Set the boundary conditions for the solution of the elliptic @@ -207,19 +223,6 @@ namespace aspect */ std::set tangential_mesh_boundary_indicators; - /** - * A handle on the connection that connects the Stokes assembler - * signal of the main simulator object to the apply_stabilization() - * function. We keep track of this connection because we need to - * break it once the current free surface handler object goes out - * of scope. - * - * With the current variable, the connection is broken once the - * scoped_connection goes out of scope, i.e., when the surrounding - * class is destroyed. - */ - boost::signals2::scoped_connection assembler_connection; - friend class Simulator; friend class SimulatorAccess; }; diff --git a/include/aspect/melt.h b/include/aspect/melt.h index c5a9b867f43..e5706202752 100644 --- a/include/aspect/melt.h +++ b/include/aspect/melt.h @@ -24,7 +24,7 @@ #include #include -#include +#include #include #include @@ -137,132 +137,147 @@ namespace aspect namespace Assemblers { /** - * A class for the definition of functions that implement the - * linear system terms for the *melt* migration compressible or - * incompressible equations. - */ + * A base class for the definition of assemblers that implement the + * linear system terms for the *melt* migration compressible or + * incompressible equations. + */ template - class MeltEquations : public aspect::internal::Assembly::Assemblers::AssemblerBase, + class MeltInterface : public aspect::Assemblers::Interface, public SimulatorAccess { public: + virtual ~MeltInterface () {}; + /** - * Attach melt outputs. + * Attach melt outputs. Since most melt assemblers require the + * melt material model properties they are created in this base class + * already. */ virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const; + }; + /** + * Compute the integrals for the preconditioner for the Stokes system in + * the case of melt migration on a single cell. + */ + template + class MeltStokesPreconditioner : public MeltInterface + { + public: + virtual ~MeltStokesPreconditioner () {}; - /** - * Compute the integrals for the preconditioner for the Stokes system in - * the case of melt migration on a single cell. - */ + virtual void - local_assemble_stokes_preconditioner_melt (const double pressure_scaling, - internal::Assembly::Scratch::StokesPreconditioner &scratch, - internal::Assembly::CopyData::StokesPreconditioner &data) const; + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; - /** - * Compute the integrals for the Stokes matrix and right hand side in - * the case of melt migration on a single cell. - */ + /** + * Compute the integrals for the Stokes matrix and right hand side in + * the case of melt migration on a single cell. + */ + template + class MeltStokesSystem : public MeltInterface + { + public: + virtual ~MeltStokesSystem () {}; + + virtual void - local_assemble_stokes_system_melt (const typename DoFHandler::active_cell_iterator &cell, - const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const; + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; - /** - * Compute the boundary integrals for the Stokes right hand side in - * the case of melt migration on a single cell. These boundary terms - * are used to describe Neumann boundary conditions for the fluid pressure. - */ + + /** + * Compute the boundary integrals for the Stokes right hand side in + * the case of melt migration on a single cell. These boundary terms + * are used to describe Neumann boundary conditions for the fluid pressure. + */ + template + class MeltStokesSystemBoundary : public MeltInterface + { + public: + virtual ~MeltStokesSystemBoundary () {}; + + virtual void - local_assemble_stokes_system_melt_boundary (const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - const double pressure_scaling, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const; + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; - /** - * Compute the integrals for the Advection system matrix and right hand side in - * the case of melt migration on a single cell. - */ + /** + * Compute the integrals for the Advection system matrix and right hand side in + * the case of melt migration on a single cell. + */ + template + class MeltAdvectionSystem : public MeltInterface + { + public: + virtual ~MeltAdvectionSystem () {}; + + virtual void - local_assemble_advection_system_melt (const typename Simulator::AdvectionField &advection_field, - const double artificial_viscosity, - internal::Assembly::Scratch::AdvectionSystem &scratch, - internal::Assembly::CopyData::AdvectionSystem &data) const; + 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. - */ - std::vector - compute_advection_system_residual_melt(const typename Simulator::AdvectionField &advection_field, - internal::Assembly::Scratch::AdvectionSystem &scratch) 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; - private: /** - * Compute the right-hand side of the fluid pressure equation of the Stokes - * system in case the simulation uses melt transport. This includes a term - * derived from Darcy's law, a term including the melting rate and a term dependent - * on the densities and velocities of fluid and solid. + * Attach melt outputs. */ - double - compute_fluid_pressure_rhs(const internal::Assembly::Scratch::StokesSystem &scratch, - MaterialModel::MaterialModelInputs &material_model_inputs, - MaterialModel::MaterialModelOutputs &material_model_outputs, - const unsigned int q_point) const; + virtual + void + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const; + }; - /** - * Compute the right-hand side for the porosity system. This function - * returns 0 for the temperature and all of the compositional fields, except - * for the porosity. It includes the melting rate and a term dependent - * on the density and velocity. - * - * This function is implemented in - * source/simulator/melt.cc. - */ - double - compute_melting_RHS(const internal::Assembly::Scratch::AdvectionSystem &scratch, - const typename MaterialModel::Interface::MaterialModelInputs &material_model_inputs, - const typename MaterialModel::Interface::MaterialModelOutputs &material_model_outputs, - const typename Simulator::AdvectionField &advection_field, - const unsigned int q_point) const; + /** + * Integrate the local fluid pressure shape functions on a single cell + * for models with melt migration, so that they can later be used to do + * the pressure right-hand side compatibility modification. + */ + template + class MeltPressureRHSCompatibilityModification : public MeltInterface + { + public: + virtual ~MeltPressureRHSCompatibilityModification () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; }; /** - * A namespace for the definition of functions that implement various - * other terms that need to occasionally or always be assembled. + * Assemble traction boundary condition terms for models with melt. */ - namespace OtherTerms + template + class MeltBoundaryTraction : public MeltInterface { - /** - * Integrate the local fluid pressure shape functions on a single cell - * for models with melt migration, so that they can later be used to do - * the pressure right-hand side compatibility modification. - */ - template - void - pressure_rhs_compatibility_modification_melt (const SimulatorAccess &simulator_access, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data); + public: + virtual ~MeltBoundaryTraction () {}; - /** - * Assemble traction boundary condition terms for models with melt. - */ - template - void - boundary_traction_melt (const SimulatorAccess &simulator_access, - const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data); - } + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; } diff --git a/include/aspect/newton.h b/include/aspect/newton.h index 437a24677c9..6f5e625ee9b 100644 --- a/include/aspect/newton.h +++ b/include/aspect/newton.h @@ -22,7 +22,7 @@ #ifndef _aspect__newton_h #define _aspect__newton_h -#include +#include #include #include #include @@ -113,88 +113,129 @@ namespace aspect namespace Assemblers { /** - * A class containing the functions to assemble the different terms of the Newton Stokes system. - */ + * A base class for the definition of assemblers that implement the + * linear system terms for the NewtonStokes solver scheme. + */ template - class NewtonStokesAssembler : public aspect::internal::Assembly::Assemblers::AssemblerBase, + class NewtonInterface : public aspect::Assemblers::Interface, public SimulatorAccess { public: + virtual ~NewtonInterface () {}; + /** - * This function assembles the terms of the Newton Stokes preconditioner matrix for the current cell. + * Attach Newton outputs. Since most Newton assemblers require the + * material model derivatives they are created in this base class + * already. */ + virtual void - preconditioner (const double pressure_scaling, - internal::Assembly::Scratch::StokesPreconditioner &scratch, - internal::Assembly::CopyData::StokesPreconditioner &data) const; + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const; + }; + + /** + * This class assembles the terms of the Newton Stokes preconditioner matrix for the current cell. + */ + template + class NewtonStokesPreconditioner : public NewtonInterface + { + public: + virtual ~NewtonStokesPreconditioner () {}; - /** - * This function assembles the terms for the matrix and right-hand-side of the incompressible - * Newton Stokes system for the current cell. - */ void - incompressible_terms (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const; + execute (internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the terms for the matrix and right-hand-side of the incompressible + * Newton Stokes system for the current cell. + */ + template + class NewtonStokesIncompressibleTerms : public NewtonInterface + { + public: + virtual ~NewtonStokesIncompressibleTerms () {}; - /** - * This function assembles the term that arises in the viscosity term of Newton Stokes matrix for - * compressible models, because the divergence of the velocity is no longer zero. - */ void - compressible_strain_rate_viscosity_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const; + execute (internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; - /** - * This function assembles the right-hand-side term of the Newton Stokes system - * that is caused by the compressibility in the mass conservation equation. - * This function approximates this term as - * $- \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u}$ - */ + /** + * This class assembles the term that arises in the viscosity term of the Newton Stokes matrix for + * compressible models, because the divergence of the velocity is not longer zero. + */ + template + class NewtonStokesCompressibleStrainRateViscosityTerm : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~NewtonStokesCompressibleStrainRateViscosityTerm () {}; + + virtual void - reference_density_compressibility_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const; + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; - /** - * This function assembles the compressibility term of the Newton Stokes system - * that is caused by the compressibility in the mass conservation equation. - * It includes this term implicitly in the matrix, - * which is therefore not longer symmetric. - * This function approximates this term as - * $ - \nabla \mathbf{u} - \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$ - */ + /** + * This class assembles the right-hand-side term of the Newton Stokes system + * that is caused by the compressibility in the mass conservation equation. + * This function approximates this term as + * $- \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u}$ + */ + template + class NewtonStokesReferenceDensityCompressibilityTerm : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~NewtonStokesReferenceDensityCompressibilityTerm () {}; + + virtual void - implicit_reference_density_compressibility_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const; + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; - /** - * This function assembles the right-hand-side term of the Newton Stokes system - * that is caused by the compressibility in the mass conservation equation. - * This function approximates this term as - * $ - \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$ - */ + /** + * This class assembles the compressibility term of the Newton Stokes system + * that is caused by the compressibility in the mass conservation equation. + * It includes this term implicitly in the matrix, + * which is therefore not longer symmetric. + * This function approximates this term as + * $ - \nabla \mathbf{u} - \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$ + */ + template + class NewtonStokesImplicitReferenceDensityCompressibilityTerm : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~NewtonStokesImplicitReferenceDensityCompressibilityTerm () {}; + + virtual void - isothermal_compression_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const; + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the right-hand-side term of the Newton Stokes system + * that is caused by the compressibility in the mass conservation equation. + * This function approximates this term as + * $ - \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$ + */ + template + class NewtonStokesIsothermalCompressionTerm : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~NewtonStokesIsothermalCompressionTerm () {}; - /** - * Attach derivatives outputs. - */ virtual void - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const; + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; }; } } diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h index ada5d7c1698..e39db3a8cad 100644 --- a/include/aspect/simulator.h +++ b/include/aspect/simulator.h @@ -99,15 +99,15 @@ namespace aspect template struct StokesSystem; template struct AdvectionSystem; } - - namespace Assemblers - { - template class AssemblerBase; - } - template struct AssemblerLists; } } + namespace Assemblers + { + template class Interface; + template class Manager; + } + /** * This is the main class of ASPECT. It implements the overall simulation * algorithm using the numerical methods discussed in the papers and manuals @@ -740,51 +740,13 @@ namespace aspect * matrices, and right hand side vectors. * * One would probably want this variable to just be a member of type - * internal::Assembly::AssemblerLists, but this requires that + * Assemblers::Manager, but this requires that * this type is declared in the current scope, and that would require - * including which we don't want because it's big. + * including which we don't want because it's big. * Consequently, we just store a pointer to such an object, and create * the object pointed to at the top of set_assemblers(). */ - std_cxx11::unique_ptr > assemblers; - - /** - * A collection of objects that implement member functions that may - * appear in the assembler signal lists. What the objects do is not - * actually important, but individual assembler objects may encapsulate - * data that is used by concrete assemblers. - * - * The objects pointed to by this vector are created in - * set_assemblers(), and are later destroyed by the destructor - * of the current class. - */ - std::vector > > assembler_objects; - - /** - * Material models, through functions derived from - * MaterialModel::Interface::evaluate(), put their computed material - * parameters into a structure of type MaterialModel::MaterialModelOutputs. - * By default, material models will compute those parameters that - * correspond to the member variables of that structure. However, - * there are situations where parts of the simulator need additional - * pieces of information; a typical example would be the use of a - * Newton scheme that also requires the computation of derivatives - * of material parameters with respect to pressure, temperature, and - * possibly other variables. - * - * The computation of such additional information is controlled by - * the presence of a collection of pointers in - * MaterialModel::MaterialModelOutputs that point to additional - * objects. Whether or not one needs these additional objects depends - * on what linear system is being assembled, or what postprocessing - * wants to compute. For the purpose of assembly, the current - * function creates the additional objects (such as the one that stores - * derivatives) and adds pointers to them to the collection, based on - * what assemblers are selected. It does so by calling the - * internal::Assemblers::AssemblerBase::create_additional_material_model_outputs() - * functions from each object in Simulator::assembler_objects. - */ - void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &) const; + std_cxx11::unique_ptr > assemblers; /** * Determine, based on the run-time parameters of the current simulation, diff --git a/include/aspect/simulator/assemblers/advection.h b/include/aspect/simulator/assemblers/advection.h new file mode 100644 index 00000000000..95e21eee43f --- /dev/null +++ b/include/aspect/simulator/assemblers/advection.h @@ -0,0 +1,102 @@ +/* + Copyright (C) 2017 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file doc/COPYING. If not see + . +*/ + +#ifndef _aspect_simulator_assemblers_advection_h +#define _aspect_simulator_assemblers_advection_h + + +#include +#include + +namespace aspect +{ + namespace Assemblers + { + /** + * This class assembles the terms for the matrix and right-hand-side of the advection + * equation for the current cell. + */ + template + class AdvectionSystem : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~AdvectionSystem () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the face terms for the matrix and right-hand-side of + * the discontinuous advection equation for a face at the boundary of the domain. + */ + template + class AdvectionSystemBoundaryFace : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~AdvectionSystemBoundaryFace () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the face terms for the matrix and right-hand-side of + * the discontinuous advection equation for a face in the interior of the domain. + */ + template + class AdvectionSystemInteriorFace : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~AdvectionSystemInteriorFace () {}; + + virtual + void + 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; + }; + } +} + + +#endif diff --git a/include/aspect/simulator/assemblers/interface.h b/include/aspect/simulator/assemblers/interface.h new file mode 100644 index 00000000000..7487f3d1b33 --- /dev/null +++ b/include/aspect/simulator/assemblers/interface.h @@ -0,0 +1,629 @@ +/* + Copyright (C) 2017 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + +#ifndef _aspect_simulator_assemblers_interface_h +#define _aspect_simulator_assemblers_interface_h + +#include +#include +#include + +#include + +namespace aspect +{ + using namespace dealii; + + template + class Simulator; + + struct AdvectionField; + + /** + * A namespace that is used for internal scratch objects, i.e. objects that define + * the inputs and outputs for the individual assembler objects. The different classes + * are provided for the different systems of equations, and the inputs are filled + * by the local_assemble_... functions in assembly.cc, + * before being handed over to the assemblers. + */ + namespace internal + { + namespace Assembly + { + namespace Scratch + { + template + struct ScratchBase + { + ScratchBase() {}; + + virtual ~ScratchBase () {}; + + /** + * Cell object on which we currently operate. + */ + typename DoFHandler::active_cell_iterator cell; + + /** + * The number of the face object with respect to the current + * cell on which we operate. If we currently + * operate on a cell, this is not meaningful. + */ + unsigned face_number; + }; + + template + struct StokesPreconditioner: public ScratchBase + { + StokesPreconditioner (const FiniteElement &finite_element, + const Quadrature &quadrature, + const Mapping &mapping, + const UpdateFlags update_flags, + const unsigned int n_compositional_fields, + const unsigned int stokes_dofs_per_cell, + const bool add_compaction_pressure, + const bool rebuild_stokes_matrix); + StokesPreconditioner (const StokesPreconditioner &data); + + virtual ~StokesPreconditioner (); + + FEValues finite_element_values; + + std::vector local_dof_indices; + std::vector dof_component_indices; + std::vector > grads_phi_u; + std::vector div_phi_u; + std::vector phi_p; + std::vector phi_p_c; + std::vector > grad_phi_p; + + double pressure_scaling; + + /** + * Material model inputs and outputs computed at the current + * linearization point. + */ + MaterialModel::MaterialModelInputs material_model_inputs; + MaterialModel::MaterialModelOutputs material_model_outputs; + + /** + * Whether the Stokes matrix should be rebuild during this + * assembly. If the matrix does not change, assembling the right + * hand side is sufficient. + */ + const bool rebuild_stokes_matrix; + }; + + + + // We derive the StokesSystem scratch array from the + // StokesPreconditioner array. We do this because all the objects that + // are necessary for the assembly of the preconditioner are also + // needed for the actual system matrix and right hand side, plus some + // extra data that we need for the time stepping and traction boundaries + // on the right hand side. + template + struct StokesSystem : public StokesPreconditioner + { + StokesSystem (const FiniteElement &finite_element, + const Mapping &mapping, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const UpdateFlags update_flags, + const UpdateFlags face_update_flags, + const unsigned int n_compositional_fields, + const unsigned int stokes_dofs_per_cell, + const bool add_compaction_pressure, + const bool use_reference_profile, + const bool rebuild_stokes_matrix, + const bool rebuild_stokes_newton_matrix); + + StokesSystem (const StokesSystem &data); + + FEFaceValues face_finite_element_values; + + std::vector > phi_u; + std::vector > velocity_values; + std::vector velocity_divergence; + + /** + * Material model inputs and outputs computed at the current + * linearization point. + * + * In contrast to the variables above, the following two + * variables are used in the assembly at quadrature points + * on faces, not on cells. + */ + MaterialModel::MaterialModelInputs face_material_model_inputs; + MaterialModel::MaterialModelOutputs face_material_model_outputs; + + /** + * In some approximations of the Stokes equations the density used + * for the mass conservation (/continuity) equation is some form + * of reference density, while the density used for calculating + * the buoyancy force is the full density. In case such a formulation + * is used the reference density (and its derivative in depth + * direction) is queried from the adiabatic conditions plugin + * and is stored in these variables. + */ + std::vector reference_densities; + std::vector reference_densities_depth_derivative; + + /** + * Whether the Newton solver Stokes matrix should be rebuild during + * this assembly. If the matrix does not change, assembling the right + * hand side is sufficient. + */ + const bool rebuild_newton_stokes_matrix; + }; + + + + template + struct AdvectionSystem: public ScratchBase + { + AdvectionSystem (const FiniteElement &finite_element, + const FiniteElement &advection_element, + const Mapping &mapping, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const UpdateFlags update_flags, + const UpdateFlags face_update_flags, + const unsigned int n_compositional_fields, + const typename Simulator::AdvectionField &advection_field); + AdvectionSystem (const AdvectionSystem &data); + + FEValues finite_element_values; + + std_cxx11::shared_ptr > face_finite_element_values; + std_cxx11::shared_ptr > neighbor_face_finite_element_values; + std_cxx11::shared_ptr > subface_finite_element_values; + + std::vector local_dof_indices; + + /** + * Variables describing the values and gradients of the + * shape functions at the quadrature points, as they are + * used in the advection assembly function. note that the sizes + * of these arrays are equal to the number of shape functions + * corresponding to the advected field (and not of the overall + * field!), and that they are also correspondingly indexed. + */ + std::vector phi_field; + std::vector > grad_phi_field; + std::vector face_phi_field; + std::vector > face_grad_phi_field; + std::vector neighbor_face_phi_field; + std::vector > neighbor_face_grad_phi_field; + + std::vector > old_velocity_values; + std::vector > old_old_velocity_values; + + std::vector old_pressure; + std::vector old_old_pressure; + std::vector > old_pressure_gradients; + std::vector > old_old_pressure_gradients; + + std::vector > old_strain_rates; + std::vector > old_old_strain_rates; + + std::vector old_temperature_values; + std::vector old_old_temperature_values; + + std::vector old_field_values; + std::vector old_old_field_values; + std::vector > old_field_grads; + std::vector > old_old_field_grads; + std::vector old_field_laplacians; + std::vector old_old_field_laplacians; + + std::vector > old_composition_values; + std::vector > old_old_composition_values; + + std::vector current_temperature_values; + std::vector > current_velocity_values; + std::vector > face_current_velocity_values; + std::vector > mesh_velocity_values; + std::vector > face_mesh_velocity_values; + + std::vector > current_strain_rates; + std::vector > current_composition_values; + std::vector current_velocity_divergences; + + /** + * Material model inputs and outputs computed at the current + * linearization point. + */ + MaterialModel::MaterialModelInputs material_model_inputs; + MaterialModel::MaterialModelOutputs material_model_outputs; + + MaterialModel::MaterialModelInputs face_material_model_inputs; + MaterialModel::MaterialModelOutputs face_material_model_outputs; + + MaterialModel::MaterialModelInputs neighbor_face_material_model_inputs; + MaterialModel::MaterialModelOutputs neighbor_face_material_model_outputs; + + /** + * Heating model outputs computed at the quadrature points of the + * current cell at the time of the current linearization point. + * As explained in the class documentation of + * HeatingModel::HeatingModelOutputs each term contains the sum of all + * enabled heating mechanism contributions. + */ + HeatingModel::HeatingModelOutputs heating_model_outputs; + HeatingModel::HeatingModelOutputs face_heating_model_outputs; + HeatingModel::HeatingModelOutputs neighbor_face_heating_model_outputs; + + const typename Simulator::AdvectionField *advection_field; + double artificial_viscosity; + }; + + + + } + + + // The CopyData arrays are similar to the + // Scratch arrays. They provide a + // constructor, a copy operation, and + // some arrays for local matrix, local + // vectors and the relation between local + // and global degrees of freedom (a.k.a. + // local_dof_indices). + namespace CopyData + { + template + struct CopyDataBase + { + CopyDataBase() {}; + + virtual ~CopyDataBase () {}; + }; + + template + struct StokesPreconditioner: public CopyDataBase + { + StokesPreconditioner (const unsigned int stokes_dofs_per_cell); + StokesPreconditioner (const StokesPreconditioner &data); + + virtual ~StokesPreconditioner (); + + FullMatrix local_matrix; + std::vector local_dof_indices; + + /** + * Extract the values listed in @p all_dof_indices only if + * it corresponds to the Stokes component and copy it to the variable + * local_dof_indices declared above in the same class as this function + */ + void extract_stokes_dof_indices(const std::vector &all_dof_indices, + const Introspection &introspection, + const FiniteElement &finite_element); + }; + + + + template + struct StokesSystem : public StokesPreconditioner + { + StokesSystem (const unsigned int stokes_dofs_per_cell, + const bool do_pressure_rhs_compatibility_modification); + StokesSystem (const StokesSystem &data); + + Vector local_rhs; + Vector local_pressure_shape_function_integrals; + }; + + + + template + struct AdvectionSystem: public CopyDataBase + { + /** + * Constructor. + * + * @param finite_element The element that describes the field for + * which we are trying to assemble a linear system. Not + * the global finite element. + * @param field_is_discontinuous If true, the field is a DG element. + */ + AdvectionSystem (const FiniteElement &finite_element, + const bool field_is_discontinuous); + + /** + * Copy constructor. + */ + AdvectionSystem (const AdvectionSystem &data); + + /** + * Local contributions to the global matrix + * that correspond only to the variables listed in local_dof_indices + */ + FullMatrix local_matrix; + + /** Local contributions to the global matrix from the face terms in the + * discontinuous Galerkin method. The vectors are of length + * GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell + * so as to hold one matrix for each possible face or subface of the cell. + * The discontinuous Galerkin bilinear form contains terms arising from internal + * (to the cell) values and external (to the cell) values. + * _int_ext and ext_int hold the terms arising from the pairing between a cell + * and its neighbor, while _ext_ext is the pairing of the neighbor's dofs with + * themselves. In the continuous Galerkin case, these are unused, and set to size zero. + **/ + std::vector > local_matrices_int_ext; + std::vector > local_matrices_ext_int; + std::vector > local_matrices_ext_ext; + + /** + * Local contributions to the right hand side + * that correspond only to the variables listed in local_dof_indices + */ + Vector local_rhs; + + /** Denotes which face matrices have actually been assembled in the DG field + * assembly. Entries for matrices not used (for example, those corresponding + * to non-existent subfaces; or faces being assembled by the neighboring cell) + * are set to false. + **/ + std::vector assembled_matrices; + + /** + * Indices of those degrees of freedom that actually correspond + * to the temperature or compositional field. since this structure + * is used to represent just contributions to the advection + * systems, there will be no contributions to other parts of the + * system and consequently, we do not need to list here indices + * that correspond to velocity or pressure degrees (or, in fact + * any other variable outside the block we are currently considering) + */ + std::vector local_dof_indices; + /** Indices of the degrees of freedom corresponding to the temperature + * or composition field on all possible neighboring cells. This is used + * in the discontinuous Galerkin method. The outer std::vector has + * length GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell, + * and has size zero if in the continuous Galerkin case. + **/ + std::vector > neighbor_dof_indices; + }; + } + } + } + + /** + * A namespace for the definition of assemblers for the various terms in the + * linear systems ASPECT solves. + */ + namespace Assemblers + { + /** + * A base class for objects that implement assembly + * operations. + * + * The point of this class is primarily so that we can store + * pointers to such objects in a list. The objects are created + * in Simulator::set_assemblers() (which is the only place where + * we know their actual type) and destroyed in the destructor of + * the Simulator object. + */ + template + class Interface + { + public: + virtual ~Interface () {} + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &, + internal::Assembly::CopyData::CopyDataBase &) 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. + */ + virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &) const {} + }; + + /** + * A base class for objects that implement assembly + * operations. + * + * The point of this class is primarily so that we can store + * pointers to such objects in a list. The objects are created + * in Simulator::set_assemblers() (which is the only place where + * we know their actual type) 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. + */ + virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &) const {} + }; + + /** + * A class that owns member variables representing + * all assemblers that need to be called when + * assembling right hand side vectors, matrices, or complete linear + * systems. We use this approach in order to support the following + * cases: + * - Assembling different formulations: When assembling either the + * full equations or only the Boussinesq approximation (to give just + * two examples), one needs different terms. This could be achieved + * using a large number of switch or if + * statements in the code, or one could encapsulate each equation + * or approximation into a collection of assemblers for this particular + * purpose. The approach chosen here in essence allows the + * implementation of each set of equations in its own scope, and we + * then just need to store a pointer to the object that + * assembles the Stokes system (for example) for the selected + * approximation. The pointer to this function is stored in the + * appropriate member variable of this class. + * - Sometimes, we want to assemble a number of terms that build on + * each other. An example is the addition of free boundary terms + * to the Stokes matrix. Rather than having to "know" in one + * place about all of the terms that need to be assembled, + * we simply add the function that computes these terms as + * another object to the appropriate set of assemblers declared + * in this class. + */ + template + class Manager + { + public: + /** + * A vector of pointers containing all assemblers for the Stokes preconditioner. + * These assemblers are called once per cell. + */ + std::vector > > stokes_preconditioner; + + /** + * A vector of pointers containing all assemblers for the Stokes system. + * These assemblers are called once per cell. + */ + std::vector > > stokes_system; + + /** + * A vector of pointers containing all assemblers for the Stokes system. 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. traction boundary + * conditions). + */ + std::vector > > stokes_system_on_boundary_face; + + /** + * A vector of pointers containing all assemblers for the advection systems. + * These assemblers are called once per cell. + */ + std::vector > > advection_system; + + /** + * A vector of pointers containing all assemblers for the Advection system. 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 > > advection_system_on_boundary_face; + + /** + * A vector of pointers containing all assemblers for the Advection system. 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 > > advection_system_on_interior_face; + + /** + * A vector of pointers containing all assemblers for 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. + * + * There are a number of pieces of information that are always + * assumed to be needed. For example, the Stokes and advection + * assemblers will always need to have access to the material + * model outputs. But the Stokes assembler may or may not need + * access to material model output for quadrature points on faces. + * + * These properties are all preset in a conservative way + * (i.e., disabled) in the constructor of this class, but can + * be enabled in Simulator::set_assemblers() when adding + * individual assemblers. Functions such as + * Simulator::local_assemble_stokes_preconditioner(), + * Simulator::local_assemble_stokes_system() will then query + * these flags to determine whether something has to be + * initialized for at least one of the assemblers they call. + */ + struct Properties + { + /** + * Constructor. Disable all properties as described in the + * class documentation. + */ + Properties (); + + /** + * Whether or not at least one of the assembler slots in + * a signal requires the initialization and re-computation of + * a MaterialModelOutputs object for each face. This + * property is only relevant to assemblers that operate on + * boundary faces. + */ + bool need_face_material_model_data; + + /** + * Whether or not at least one of the assembler slots in a + * signal requires the evaluation of the FEFaceValues object. + */ + bool need_face_finite_element_evaluation; + + /** + * Whether or not at least one of the assembler slots in a + * signal requires the computation of the viscosity. + */ + bool need_viscosity; + + /** + * A list of FEValues UpdateFlags that are necessary for + * a given operation. Assembler objects may add to this list + * as necessary; it will be initialized with a set of + * "default" flags that will always be set. + */ + UpdateFlags needed_update_flags; + }; + + /** + * A list of properties of the various types of assemblers. + * These property lists are set in Simulator::set_assemblers() + * where we add individual functions to the signals above. + */ + Properties stokes_preconditioner_assembler_properties; + Properties stokes_system_assembler_properties; + Properties stokes_system_assembler_on_boundary_face_properties; + std::vector advection_system_assembler_properties; + std::vector advection_system_assembler_on_face_properties; + }; + } +} + + +#endif diff --git a/include/aspect/simulator/assemblers/stokes.h b/include/aspect/simulator/assemblers/stokes.h new file mode 100644 index 00000000000..3636d801ccb --- /dev/null +++ b/include/aspect/simulator/assemblers/stokes.h @@ -0,0 +1,204 @@ +/* + Copyright (C) 2017 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file doc/COPYING. If not see + . +*/ + +#ifndef _aspect_simulator_assemblers_stokes_h +#define _aspect_simulator_assemblers_stokes_h + + +#include +#include + +namespace aspect +{ + namespace Assemblers + { + /** + * A class containing the functions to assemble the Stokes preconditioner. + */ + template + class StokesPreconditioner : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~StokesPreconditioner () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * A class containing the functions to assemble the compressible adjustment + * to the Stokes preconditioner. + */ + template + class StokesCompressiblePreconditioner : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~StokesCompressiblePreconditioner () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the terms for the matrix and right-hand-side of the incompressible + * Stokes equation for the current cell. + */ + template + class StokesIncompressibleTerms : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~StokesIncompressibleTerms () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + + /** + * Create AdditionalMaterialOutputsStokesRHS if we need to do so. + */ + virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const; + }; + + /** + * This class assembles the term that arises in the viscosity term of Stokes matrix for + * compressible models, because the divergence of the velocity is not longer zero. + */ + template + class StokesCompressibleStrainRateViscosityTerm : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~StokesCompressibleStrainRateViscosityTerm () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the right-hand-side term of the Stokes equation + * that is caused by the compressibility in the mass conservation equation. + * This class approximates this term as + * $- \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u}$ + */ + template + class StokesReferenceDensityCompressibilityTerm : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~StokesReferenceDensityCompressibilityTerm () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the compressibility term of the Stokes equation + * that is caused by the compressibility in the mass conservation equation. + * It includes this term implicitly in the matrix, + * which is therefore not longer symmetric. + * This class approximates this term as + * $ - \nabla \mathbf{u} - \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$ + */ + template + class StokesImplicitReferenceDensityCompressibilityTerm : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~StokesImplicitReferenceDensityCompressibilityTerm () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the right-hand-side term of the Stokes equation + * that is caused by the compressibility in the mass conservation equation. + * This class approximates this term as + * $ - \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$ + */ + template + class StokesIsothermalCompressionTerm : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~StokesIsothermalCompressionTerm () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the right-hand-side term of the Stokes equation + * that is caused by the compressibility in the mass conservation equation. + * This class approximates this term as + * $ - \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$ + */ + template + class StokesBoundaryTraction : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~StokesBoundaryTraction () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + + /** + * This class assembles the right-hand-side term of the Stokes equation + * that is caused by the compressibility in the mass conservation equation. + * This class approximates this term as + * $ - \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$ + */ + template + class StokesPressureRHSCompatibilityModification : public Assemblers::Interface, + public SimulatorAccess + { + public: + virtual ~StokesPressureRHSCompatibilityModification () {}; + + virtual + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const; + }; + } +} + + +#endif diff --git a/include/aspect/simulator_access.h b/include/aspect/simulator_access.h index 54bac2751dd..7e5bd975839 100644 --- a/include/aspect/simulator_access.h +++ b/include/aspect/simulator_access.h @@ -99,6 +99,7 @@ namespace aspect } template class MeltHandler; + template class FreeSurfaceHandler; template class NewtonHandler; @@ -479,13 +480,6 @@ namespace aspect const bool compute_strainrate, MaterialModel::MaterialModelInputs &material_model_inputs) const; - /** - * This function simply calls Simulator::create_additional_material_model_outputs() - * with the given arguments. - */ - void - create_additional_material_model_outputs (MaterialModel::MaterialModelOutputs &) const; - /** * Return a pointer to the gravity model description. */ @@ -645,6 +639,13 @@ namespace aspect const NewtonHandler & get_newton_handler () const; + /** + * Return a set of boundary indicators that describes which of the + * boundaries have a free surface boundary condition + */ + const FreeSurfaceHandler & + get_free_surface_handler () const; + /** * Return a reference to the lateral averaging object owned * by the simulator, which can be used to query lateral averages diff --git a/include/aspect/simulator_signals.h b/include/aspect/simulator_signals.h index 72d833ce9f5..654aec5af15 100644 --- a/include/aspect/simulator_signals.h +++ b/include/aspect/simulator_signals.h @@ -33,19 +33,13 @@ namespace aspect { - namespace internal + namespace Assemblers { - namespace Assembly - { - namespace Assemblers - { - template - class AssemblerBase; - } + template + class Manager; - template - struct AssemblerLists; - } + template + class Interface; } /** @@ -263,8 +257,7 @@ namespace aspect * allows modification of the assembly objects active in this simulation. */ boost::signals2::signal &, - aspect::internal::Assembly::AssemblerLists &, - std::vector > > &)> + aspect::Assemblers::Manager &)> set_assemblers; }; diff --git a/source/simulator/assemblers/advection.cc b/source/simulator/assemblers/advection.cc index 5c1b5eab2d3..157402f5393 100644 --- a/source/simulator/assemblers/advection.cc +++ b/source/simulator/assemblers/advection.cc @@ -18,10 +18,10 @@ . */ +#include + #include #include -#include -#include namespace aspect { @@ -52,13 +52,16 @@ namespace aspect template void - AdvectionAssembler::local_assemble_advection_system (const typename Simulator::AdvectionField &advection_field, - const double artificial_viscosity, - internal::Assembly::Scratch::AdvectionSystem &scratch, - internal::Assembly::CopyData::AdvectionSystem &data) const + AdvectionSystem::execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::AdvectionSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::AdvectionSystem &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); + + const typename Simulator::AdvectionField advection_field = *scratch.advection_field; const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; const unsigned int advection_dofs_per_cell = data.local_dof_indices.size(); @@ -171,7 +174,7 @@ namespace aspect { data.local_matrix(i,j) += ( - (time_step * (conductivity + artificial_viscosity) + (time_step * (conductivity + scratch.artificial_viscosity) * (scratch.grad_phi_field[i] * scratch.grad_phi_field[j])) + ((time_step * (scratch.phi_field[i] * (current_u * scratch.grad_phi_field[j]))) + (bdf2_factor * scratch.phi_field[i] * scratch.phi_field[j])) * @@ -187,9 +190,11 @@ namespace aspect template std::vector - AdvectionAssembler::compute_advection_system_residual(const typename Simulator::AdvectionField &advection_field, - internal::Assembly::Scratch::AdvectionSystem &scratch) const + AdvectionSystemResidual::execute(internal::Assembly::Scratch::ScratchBase &scratch_base) const { + internal::Assembly::Scratch::AdvectionSystem &scratch = dynamic_cast& > (scratch_base); + + const typename Simulator::AdvectionField advection_field = *scratch.advection_field; const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; std::vector residuals(n_q_points); @@ -240,16 +245,21 @@ namespace aspect template void - AdvectionAssembler::local_assemble_discontinuous_advection_boundary_face_terms(const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - const typename Simulator::AdvectionField &advection_field, - internal::Assembly::Scratch::AdvectionSystem &scratch, - internal::Assembly::CopyData::AdvectionSystem &data) const + AdvectionSystemBoundaryFace::execute(internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::AdvectionSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::AdvectionSystem &data = dynamic_cast& > (data_base); + const Parameters ¶meters = this->get_parameters(); const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); + const typename Simulator::AdvectionField advection_field = *scratch.advection_field; + + const unsigned int face_no = scratch.face_number; + const typename DoFHandler::face_iterator face = scratch.cell->face(face_no); + const unsigned int n_q_points = scratch.face_finite_element_values->n_quadrature_points; const double time_step = this->get_timestep(); @@ -269,7 +279,6 @@ namespace aspect const FEValuesExtractors::Scalar solution_field = advection_field.scalar_extractor(introspection); - typename DoFHandler::face_iterator face = cell->face (face_no); if (((parameters.fixed_temperature_boundary_indicators.find( face->boundary_id() @@ -351,11 +360,11 @@ namespace aspect const double dirichlet_value = (advection_field.is_temperature() ? this->get_boundary_temperature_manager().boundary_temperature( - cell->face(face_no)->boundary_id(), + face->boundary_id(), scratch.face_finite_element_values->quadrature_point(q)) : this->get_boundary_composition().boundary_composition( - cell->face(face_no)->boundary_id(), + face->boundary_id(), scratch.face_finite_element_values->quadrature_point(q), advection_field.compositional_variable)); @@ -437,11 +446,6 @@ namespace aspect } } } - else if (cell->has_periodic_neighbor (face_no)) - { - // Periodic temperature/composition term: consider the corresponding periodic faces as the case of interior faces - this->local_assemble_discontinuous_advection_interior_face_terms(cell, face_no, advection_field, scratch, data); - } else { // Neumann temperature term - no non-zero contribution as only homogeneous Neumann boundary conditions are implemented elsewhere for temperature @@ -452,16 +456,22 @@ namespace aspect template void - AdvectionAssembler::local_assemble_discontinuous_advection_interior_face_terms(const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - const typename Simulator::AdvectionField &advection_field, - internal::Assembly::Scratch::AdvectionSystem &scratch, - internal::Assembly::CopyData::AdvectionSystem &data) const + AdvectionSystemInteriorFace::execute(internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::AdvectionSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::AdvectionSystem &data = dynamic_cast& > (data_base); + const Parameters ¶meters = this->get_parameters(); const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); + const typename DoFHandler::active_cell_iterator cell = scratch.cell; + const unsigned int face_no = scratch.face_number; + const typename DoFHandler::face_iterator face = cell->face(face_no); + + const typename Simulator::AdvectionField advection_field = *scratch.advection_field; + const unsigned int n_q_points = scratch.face_finite_element_values->n_quadrature_points; const double time_step = this->get_timestep(); @@ -482,8 +492,6 @@ namespace aspect const FEValuesExtractors::Scalar solution_field = advection_field.scalar_extractor(introspection); - typename DoFHandler::face_iterator face = cell->face (face_no); - // interior face or periodic face - no contribution on RHS const typename DoFHandler::cell_iterator @@ -1143,8 +1151,10 @@ namespace aspect namespace Assemblers { #define INSTANTIATE(dim) \ - template class \ - AdvectionAssembler; + template class AdvectionSystem; \ + template class AdvectionSystemBoundaryFace; \ + template class AdvectionSystemInteriorFace; \ + template class AdvectionSystemResidual; ASPECT_INSTANTIATE(INSTANTIATE) } diff --git a/source/simulator/assemblers/interface.cc b/source/simulator/assemblers/interface.cc new file mode 100644 index 00000000000..d9073156173 --- /dev/null +++ b/source/simulator/assemblers/interface.cc @@ -0,0 +1,491 @@ +/* + Copyright (C) 2016 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file doc/COPYING. If not see + . +*/ + +#include + +#include +#include + +#include + +namespace aspect +{ + namespace internal + { + namespace Assembly + { + namespace Scratch + { + template + StokesPreconditioner:: + StokesPreconditioner (const FiniteElement &finite_element, + const Quadrature &quadrature, + const Mapping &mapping, + const UpdateFlags update_flags, + const unsigned int n_compositional_fields, + const unsigned int stokes_dofs_per_cell, + const bool add_compaction_pressure, + const bool rebuild_matrix) + : + finite_element_values (mapping, finite_element, quadrature, + update_flags), + local_dof_indices (finite_element.dofs_per_cell), + dof_component_indices(stokes_dofs_per_cell), + grads_phi_u (stokes_dofs_per_cell, numbers::signaling_nan >()), + div_phi_u (stokes_dofs_per_cell, numbers::signaling_nan()), + phi_p (stokes_dofs_per_cell, numbers::signaling_nan()), + phi_p_c (add_compaction_pressure ? stokes_dofs_per_cell : 0, numbers::signaling_nan()), + grad_phi_p (add_compaction_pressure ? stokes_dofs_per_cell : 0, numbers::signaling_nan >()), + pressure_scaling(numbers::signaling_nan()), + material_model_inputs(quadrature.size(), n_compositional_fields), + material_model_outputs(quadrature.size(), n_compositional_fields), + rebuild_stokes_matrix(rebuild_matrix) + {} + + + + template + StokesPreconditioner:: + StokesPreconditioner (const StokesPreconditioner &scratch) + : + finite_element_values (scratch.finite_element_values.get_mapping(), + scratch.finite_element_values.get_fe(), + scratch.finite_element_values.get_quadrature(), + scratch.finite_element_values.get_update_flags()), + local_dof_indices (scratch.local_dof_indices), + dof_component_indices( scratch.dof_component_indices), + grads_phi_u (scratch.grads_phi_u), + div_phi_u (scratch.div_phi_u), + phi_p (scratch.phi_p), + phi_p_c (scratch.phi_p_c), + grad_phi_p(scratch.grad_phi_p), + pressure_scaling(scratch.pressure_scaling), + material_model_inputs(scratch.material_model_inputs), + material_model_outputs(scratch.material_model_outputs), + rebuild_stokes_matrix(scratch.rebuild_stokes_matrix) + {} + + + template + StokesPreconditioner:: + ~StokesPreconditioner () + {} + + + + + template + StokesSystem:: + StokesSystem (const FiniteElement &finite_element, + const Mapping &mapping, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const UpdateFlags update_flags, + const UpdateFlags face_update_flags, + const unsigned int n_compositional_fields, + const unsigned int stokes_dofs_per_cell, + const bool add_compaction_pressure, + const bool use_reference_density_profile, + const bool rebuild_stokes_matrix, + const bool rebuild_newton_stokes_matrix) + : + StokesPreconditioner (finite_element, quadrature, + mapping, + update_flags, + n_compositional_fields, + stokes_dofs_per_cell, + add_compaction_pressure, + rebuild_stokes_matrix), + + face_finite_element_values (mapping, + finite_element, + face_quadrature, + face_update_flags), + + phi_u (stokes_dofs_per_cell, numbers::signaling_nan >()), + velocity_values (quadrature.size(), numbers::signaling_nan >()), + velocity_divergence(quadrature.size(), numbers::signaling_nan()), + face_material_model_inputs(face_quadrature.size(), n_compositional_fields), + face_material_model_outputs(face_quadrature.size(), n_compositional_fields), + reference_densities(use_reference_density_profile ? quadrature.size() : 0, numbers::signaling_nan()), + reference_densities_depth_derivative(use_reference_density_profile ? quadrature.size() : 0, numbers::signaling_nan()), + rebuild_newton_stokes_matrix(rebuild_newton_stokes_matrix) + {} + + + + template + StokesSystem:: + StokesSystem (const StokesSystem &scratch) + : + StokesPreconditioner (scratch), + + face_finite_element_values (scratch.face_finite_element_values.get_mapping(), + scratch.face_finite_element_values.get_fe(), + scratch.face_finite_element_values.get_quadrature(), + scratch.face_finite_element_values.get_update_flags()), + + phi_u (scratch.phi_u), + velocity_values (scratch.velocity_values), + velocity_divergence (scratch.velocity_divergence), + face_material_model_inputs(scratch.face_material_model_inputs), + face_material_model_outputs(scratch.face_material_model_outputs), + reference_densities(scratch.reference_densities), + reference_densities_depth_derivative(scratch.reference_densities_depth_derivative), + rebuild_newton_stokes_matrix(scratch.rebuild_newton_stokes_matrix) + {} + + + + template + AdvectionSystem:: + AdvectionSystem (const FiniteElement &finite_element, + const FiniteElement &advection_element, + const Mapping &mapping, + const Quadrature &quadrature, + const Quadrature &face_quadrature, + const UpdateFlags update_flags, + const UpdateFlags face_update_flags, + const unsigned int n_compositional_fields, + const typename Simulator::AdvectionField &field) + : + finite_element_values (mapping, + finite_element, quadrature, + update_flags), + face_finite_element_values ((face_quadrature.size() > 0 + ? + new FEFaceValues (mapping, + finite_element, face_quadrature, + face_update_flags) + : + NULL)), + neighbor_face_finite_element_values ((face_quadrature.size() > 0 + ? + new FEFaceValues (mapping, + finite_element, face_quadrature, + face_update_flags) + : + NULL)), + subface_finite_element_values ((face_quadrature.size() > 0 + ? + new FESubfaceValues (mapping, + finite_element, face_quadrature, + face_update_flags) + : + NULL)), + local_dof_indices (finite_element.dofs_per_cell), + + phi_field (advection_element.dofs_per_cell, numbers::signaling_nan()), + grad_phi_field (advection_element.dofs_per_cell, numbers::signaling_nan >()), + face_phi_field ((face_quadrature.size() > 0 ? advection_element.dofs_per_cell : 0), + numbers::signaling_nan()), + face_grad_phi_field ((face_quadrature.size() > 0 ? advection_element.dofs_per_cell : 0), + numbers::signaling_nan >()), + neighbor_face_phi_field ((face_quadrature.size() > 0 ? advection_element.dofs_per_cell : 0), + numbers::signaling_nan()), + neighbor_face_grad_phi_field ((face_quadrature.size() > 0 ? advection_element.dofs_per_cell : 0), + numbers::signaling_nan >()), + old_velocity_values (quadrature.size(), numbers::signaling_nan >()), + old_old_velocity_values (quadrature.size(), numbers::signaling_nan >()), + old_pressure (quadrature.size(), numbers::signaling_nan()), + old_old_pressure (quadrature.size(), numbers::signaling_nan()), + old_pressure_gradients (quadrature.size(), numbers::signaling_nan >()), + old_old_pressure_gradients (quadrature.size(), numbers::signaling_nan >()), + old_strain_rates (quadrature.size(), numbers::signaling_nan >()), + old_old_strain_rates (quadrature.size(), numbers::signaling_nan >()), + old_temperature_values (quadrature.size(), numbers::signaling_nan()), + old_old_temperature_values(quadrature.size(), numbers::signaling_nan()), + old_field_values (quadrature.size(), numbers::signaling_nan()), + old_old_field_values(quadrature.size(), numbers::signaling_nan()), + old_field_grads(quadrature.size(), numbers::signaling_nan >()), + old_old_field_grads(quadrature.size(), numbers::signaling_nan >()), + old_field_laplacians(quadrature.size(), numbers::signaling_nan()), + old_old_field_laplacians(quadrature.size(), numbers::signaling_nan()), + old_composition_values(n_compositional_fields, + std::vector(quadrature.size(), numbers::signaling_nan())), + old_old_composition_values(n_compositional_fields, + std::vector(quadrature.size(), numbers::signaling_nan())), + current_temperature_values(quadrature.size(), numbers::signaling_nan()), + current_velocity_values(quadrature.size(), numbers::signaling_nan >()), + face_current_velocity_values(face_quadrature.size(), numbers::signaling_nan >()), + mesh_velocity_values(quadrature.size(), numbers::signaling_nan >()), + face_mesh_velocity_values(face_quadrature.size(), numbers::signaling_nan >()), + current_strain_rates(quadrature.size(), numbers::signaling_nan >()), + current_composition_values(n_compositional_fields, + std::vector(quadrature.size(), numbers::signaling_nan())), + current_velocity_divergences(quadrature.size(), numbers::signaling_nan()), + material_model_inputs(quadrature.size(), n_compositional_fields), + material_model_outputs(quadrature.size(), n_compositional_fields), + face_material_model_inputs(face_quadrature.size(), n_compositional_fields), + face_material_model_outputs(face_quadrature.size(), n_compositional_fields), + neighbor_face_material_model_inputs(face_quadrature.size(), n_compositional_fields), + neighbor_face_material_model_outputs(face_quadrature.size(), n_compositional_fields), + heating_model_outputs(quadrature.size(), n_compositional_fields), + face_heating_model_outputs(face_quadrature.size(), n_compositional_fields), + neighbor_face_heating_model_outputs(face_quadrature.size(), n_compositional_fields), + advection_field(&field), + artificial_viscosity(numbers::signaling_nan()) + {} + + + + template + AdvectionSystem:: + AdvectionSystem (const AdvectionSystem &scratch) + : + finite_element_values (scratch.finite_element_values.get_mapping(), + scratch.finite_element_values.get_fe(), + scratch.finite_element_values.get_quadrature(), + scratch.finite_element_values.get_update_flags()), + face_finite_element_values (scratch.face_finite_element_values), + neighbor_face_finite_element_values (scratch.neighbor_face_finite_element_values), + subface_finite_element_values (scratch.subface_finite_element_values), + local_dof_indices (scratch.finite_element_values.get_fe().dofs_per_cell), + + phi_field (scratch.phi_field), + grad_phi_field (scratch.grad_phi_field), + face_phi_field (scratch.face_phi_field), + face_grad_phi_field (scratch.face_grad_phi_field), + neighbor_face_phi_field (scratch.neighbor_face_phi_field), + neighbor_face_grad_phi_field (scratch.neighbor_face_grad_phi_field), + old_velocity_values (scratch.old_velocity_values), + old_old_velocity_values (scratch.old_old_velocity_values), + old_pressure (scratch.old_pressure), + old_old_pressure (scratch.old_old_pressure), + old_pressure_gradients (scratch.old_pressure_gradients), + old_old_pressure_gradients (scratch.old_old_pressure_gradients), + old_strain_rates (scratch.old_strain_rates), + old_old_strain_rates (scratch.old_old_strain_rates), + old_temperature_values (scratch.old_temperature_values), + old_old_temperature_values (scratch.old_old_temperature_values), + old_field_values (scratch.old_field_values), + old_old_field_values(scratch.old_old_field_values), + old_field_grads (scratch.old_field_grads), + old_old_field_grads (scratch.old_old_field_grads), + old_field_laplacians (scratch.old_field_laplacians), + old_old_field_laplacians (scratch.old_old_field_laplacians), + old_composition_values(scratch.old_composition_values), + old_old_composition_values(scratch.old_old_composition_values), + current_temperature_values(scratch.current_temperature_values), + current_velocity_values(scratch.current_velocity_values), + face_current_velocity_values(scratch.face_current_velocity_values), + mesh_velocity_values(scratch.mesh_velocity_values), + face_mesh_velocity_values(scratch.face_mesh_velocity_values), + current_strain_rates(scratch.current_strain_rates), + current_composition_values(scratch.current_composition_values), + current_velocity_divergences(scratch.current_velocity_divergences), + material_model_inputs(scratch.material_model_inputs), + material_model_outputs(scratch.material_model_outputs), + face_material_model_inputs(scratch.face_material_model_inputs), + face_material_model_outputs(scratch.face_material_model_outputs), + neighbor_face_material_model_inputs(scratch.neighbor_face_material_model_inputs), + neighbor_face_material_model_outputs(scratch.neighbor_face_material_model_outputs), + heating_model_outputs(scratch.heating_model_outputs), + face_heating_model_outputs(scratch.face_heating_model_outputs), + neighbor_face_heating_model_outputs(scratch.neighbor_face_heating_model_outputs), + advection_field(scratch.advection_field), + artificial_viscosity(scratch.artificial_viscosity) + {} + } + + + namespace CopyData + { + + template + StokesPreconditioner:: + StokesPreconditioner (const unsigned int stokes_dofs_per_cell) + : + local_matrix (stokes_dofs_per_cell, + stokes_dofs_per_cell), + local_dof_indices (stokes_dofs_per_cell) + {} + + + + template + StokesPreconditioner:: + StokesPreconditioner (const StokesPreconditioner &data) + : + local_matrix (data.local_matrix), + local_dof_indices (data.local_dof_indices) + {} + + + template + void StokesPreconditioner:: + extract_stokes_dof_indices(const std::vector &all_dof_indices, + const Introspection &introspection, + const FiniteElement &finite_element) + { + const unsigned int dofs_per_cell = finite_element.dofs_per_cell; + + for (unsigned int i=0, i_stokes=0; ilocal_dof_indices[i_stokes] = all_dof_indices[i]; + ++i_stokes; + } + ++i; + } + } + + + template + StokesPreconditioner:: + ~StokesPreconditioner () + {} + + + + template + StokesSystem:: + StokesSystem (const unsigned int stokes_dofs_per_cell, + const bool do_pressure_rhs_compatibility_modification) + : + StokesPreconditioner (stokes_dofs_per_cell), + local_rhs (stokes_dofs_per_cell), + local_pressure_shape_function_integrals (do_pressure_rhs_compatibility_modification ? + stokes_dofs_per_cell + : + 0) + {} + + + + template + StokesSystem:: + StokesSystem (const StokesSystem &data) + : + StokesPreconditioner (data), + local_rhs (data.local_rhs), + local_pressure_shape_function_integrals (data.local_pressure_shape_function_integrals.size()) + {} + + + + + template + AdvectionSystem:: + AdvectionSystem (const FiniteElement &finite_element, + const bool field_is_discontinuous) + : + local_matrix (finite_element.dofs_per_cell, + finite_element.dofs_per_cell), + local_matrices_int_ext ((field_is_discontinuous + ? + GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell + : + 0), + FullMatrix(finite_element.dofs_per_cell, + finite_element.dofs_per_cell)), + local_matrices_ext_int ((field_is_discontinuous + ? + GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell + : + 0), + FullMatrix(finite_element.dofs_per_cell, + finite_element.dofs_per_cell)), + local_matrices_ext_ext ((field_is_discontinuous + ? + GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell + : + 0), + FullMatrix(finite_element.dofs_per_cell, + finite_element.dofs_per_cell)), + local_rhs (finite_element.dofs_per_cell), + + assembled_matrices ((field_is_discontinuous + ? + GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell + : + 0), false), + + local_dof_indices (finite_element.dofs_per_cell), + neighbor_dof_indices ((field_is_discontinuous + ? + GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell + : + 0), + std::vector(finite_element.dofs_per_cell)) + {} + + + + template + AdvectionSystem:: + AdvectionSystem (const AdvectionSystem &data) + : + local_matrix (data.local_matrix), + local_matrices_int_ext (data.local_matrices_int_ext), + local_matrices_ext_int (data.local_matrices_ext_int), + local_matrices_ext_ext (data.local_matrices_ext_ext), + local_rhs (data.local_rhs), + + assembled_matrices (data.assembled_matrices), + + local_dof_indices (data.local_dof_indices), + neighbor_dof_indices (data.neighbor_dof_indices) + {} + + } + } + } + + + + namespace Assemblers + { + template + Manager::Properties::Properties () + : + need_face_material_model_data (false), + need_face_finite_element_evaluation(false), + need_viscosity(false), + needed_update_flags () + {} + } +} // namespace aspect + +// explicit instantiation of the functions we implement in this file +namespace aspect +{ + +#define INSTANTIATE(dim) \ + namespace internal { \ + namespace Assembly { \ + namespace Scratch { \ + template struct StokesPreconditioner; \ + template struct StokesSystem; \ + template struct AdvectionSystem; \ + } \ + namespace CopyData { \ + template struct StokesPreconditioner; \ + template struct StokesSystem; \ + template struct AdvectionSystem; \ + } \ + } \ + } \ + namespace Assemblers { \ + template struct Manager::Properties; \ + } + ASPECT_INSTANTIATE(INSTANTIATE) +} diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc index 9f8558b495a..66c69cf1dfc 100644 --- a/source/simulator/assemblers/newton_stokes.cc +++ b/source/simulator/assemblers/newton_stokes.cc @@ -29,11 +29,20 @@ namespace aspect { template void - NewtonStokesAssembler:: - preconditioner (const double pressure_scaling, - internal::Assembly::Scratch::StokesPreconditioner &scratch, - internal::Assembly::CopyData::StokesPreconditioner &data) const + NewtonInterface::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const { + NewtonHandler::create_material_model_outputs(outputs); + } + + template + void + NewtonStokesPreconditioner:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const + { + internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); @@ -96,7 +105,7 @@ namespace aspect if (scratch.dof_component_indices[i] == scratch.dof_component_indices[j]) data.local_matrix(i, j) += ((2.0 * eta * (scratch.grads_phi_u[i] * scratch.grads_phi_u[j])) - + one_over_eta * pressure_scaling * pressure_scaling + + one_over_eta * scratch.pressure_scaling * scratch.pressure_scaling * (scratch.phi_p[i] * scratch.phi_p[j])) * JxW; } @@ -121,10 +130,9 @@ namespace aspect data.local_matrix(i, j) += ((2.0 * eta * (scratch.grads_phi_u[i] * scratch.grads_phi_u[j])) + derivative_scaling_factor * alpha * (scratch.grads_phi_u[i] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[j]) * strain_rate + scratch.grads_phi_u[j] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[i]) * strain_rate) - + one_over_eta * pressure_scaling - * pressure_scaling - * (scratch.phi_p[i] * scratch - .phi_p[j])) + + one_over_eta * scratch.pressure_scaling + * scratch.pressure_scaling + * (scratch.phi_p[i] * scratch.phi_p[j])) * JxW; } @@ -168,12 +176,13 @@ namespace aspect template void - NewtonStokesAssembler:: - incompressible_terms (const double pressure_scaling, - const bool assemble_newton_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const + NewtonStokesIncompressibleTerms:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); @@ -227,12 +236,12 @@ namespace aspect for (unsigned int i=0; i void - NewtonStokesAssembler:: - compressible_strain_rate_viscosity_term (const double /*pressure_scaling*/, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const + NewtonStokesCompressibleStrainRateViscosityTerm:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { - if (!rebuild_stokes_matrix) + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + + if (!scratch.rebuild_stokes_matrix) return; + + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); @@ -409,17 +421,16 @@ namespace aspect template void - NewtonStokesAssembler:: - reference_density_compressibility_term (const double pressure_scaling, - const bool /*rebuild_stokes_matrix*/, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const + NewtonStokesReferenceDensityCompressibilityTerm:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + // assemble RHS of: // - div u = 1/rho * drho/dz g/||g||* u - (void)parameters; - Assert(parameters.formulation_mass_conservation == + Assert(this->get_parameters().formulation_mass_conservation == Parameters::Formulation::MassConservation::reference_density_profile, ExcInternalError()); @@ -448,7 +459,7 @@ namespace aspect const double JxW = scratch.finite_element_values.JxW(q); for (unsigned int i=0; i void - NewtonStokesAssembler:: - implicit_reference_density_compressibility_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const + NewtonStokesImplicitReferenceDensityCompressibilityTerm:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + // assemble compressibility term of: // - div u - 1/rho * drho/dz g/||g||* u = 0 - (void)parameters; - Assert(parameters.formulation_mass_conservation == + Assert(this->get_parameters().formulation_mass_conservation == Parameters::Formulation::MassConservation::implicit_reference_density_profile, ExcInternalError()); - if (!rebuild_stokes_matrix) + if (!scratch.rebuild_stokes_matrix) return; const Introspection &introspection = this->introspection(); @@ -502,7 +512,7 @@ namespace aspect for (unsigned int i=0; i void - NewtonStokesAssembler:: - isothermal_compression_term (const double pressure_scaling, - const bool /*rebuild_stokes_matrix*/, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const + NewtonStokesIsothermalCompressionTerm:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + // assemble RHS of: // - div u = 1/rho * drho/dp rho * g * u - (void)parameters; - Assert(parameters.formulation_mass_conservation == + Assert(this->get_parameters().formulation_mass_conservation == Parameters::Formulation::MassConservation::isothermal_compression, ExcInternalError()); @@ -558,7 +567,7 @@ namespace aspect // to the manual, this term seems to have the wrong sign, but this // is because we negate the entire equation to make sure we get // -div(u) as the adjoint operator of grad(p) - (pressure_scaling * + (scratch.pressure_scaling * compressibility * density * (scratch.velocity_values[q] * gravity) * scratch.phi_p[i]) @@ -566,15 +575,6 @@ namespace aspect * JxW; } } - - - template - void - NewtonStokesAssembler::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const - { - NewtonHandler::create_material_model_outputs(outputs); - } - } } // namespace aspect @@ -584,8 +584,12 @@ namespace aspect namespace Assemblers { #define INSTANTIATE(dim) \ - template class \ - NewtonStokesAssembler; + template class NewtonStokesPreconditioner; \ + template class NewtonStokesIncompressibleTerms; \ + template class NewtonStokesCompressibleStrainRateViscosityTerm; \ + template class NewtonStokesReferenceDensityCompressibilityTerm; \ + template class NewtonStokesImplicitReferenceDensityCompressibilityTerm; \ + template class NewtonStokesIsothermalCompressionTerm; ASPECT_INSTANTIATE(INSTANTIATE) } diff --git a/source/simulator/assemblers/stokes.cc b/source/simulator/assemblers/stokes.cc index 60c6edc736b..d7d9f077dd7 100644 --- a/source/simulator/assemblers/stokes.cc +++ b/source/simulator/assemblers/stokes.cc @@ -18,10 +18,9 @@ . */ +#include #include #include -#include -#include #include @@ -31,35 +30,18 @@ namespace aspect { template void - StokesAssembler:: - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const + StokesPreconditioner:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { - const unsigned int n_points = outputs.viscosities.size(); + internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast& > (data_base); - 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 - StokesAssembler:: - preconditioner (const double pressure_scaling, - internal::Assembly::Scratch::StokesPreconditioner &scratch, - internal::Assembly::CopyData::StokesPreconditioner &data) const - { const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + const double pressure_scaling = scratch.pressure_scaling; // First loop over all dofs and find those that are in the Stokes system // save the component (pressure and dim velocities) each belongs to. @@ -114,8 +96,8 @@ namespace aspect * scratch.grads_phi_u[j])) + one_over_eta * pressure_scaling * pressure_scaling - * (scratch.phi_p[i] * scratch - .phi_p[j])) + * (scratch.phi_p[i] + * scratch.phi_p[j])) * JxW; } } @@ -124,11 +106,13 @@ namespace aspect template void - StokesAssembler:: - compressible_preconditioner (const double /*pressure_scaling*/, - internal::Assembly::Scratch::StokesPreconditioner &scratch, - internal::Assembly::CopyData::StokesPreconditioner &data) const + StokesCompressiblePreconditioner:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); @@ -188,16 +172,18 @@ namespace aspect template void - StokesAssembler:: - incompressible_terms (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const + StokesIncompressibleTerms:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + const double pressure_scaling = scratch.pressure_scaling; const MaterialModel::AdditionalMaterialOutputsStokesRHS *force = scratch.material_model_outputs.template get_additional_output >(); @@ -210,7 +196,7 @@ namespace aspect { scratch.phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].value (i,q); scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q); - if (rebuild_stokes_matrix) + if (scratch.rebuild_stokes_matrix) { scratch.grads_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].symmetric_gradient(i,q); scratch.div_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].divergence (i, q); @@ -222,7 +208,7 @@ namespace aspect // Viscosity scalar - const double eta = (rebuild_stokes_matrix + const double eta = (scratch.rebuild_stokes_matrix ? scratch.material_model_outputs.viscosities[q] : @@ -248,7 +234,7 @@ namespace aspect + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i]) * JxW; - if (rebuild_stokes_matrix) + if (scratch.rebuild_stokes_matrix) for (unsigned int j=0; j void - StokesAssembler:: - compressible_strain_rate_viscosity_term (const double /*pressure_scaling*/, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const + StokesIncompressibleTerms:: + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const + { + 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 + StokesCompressibleStrainRateViscosityTerm:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { - if (!rebuild_stokes_matrix) + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + + if (!scratch.rebuild_stokes_matrix) return; const Introspection &introspection = this->introspection(); @@ -327,17 +336,16 @@ namespace aspect template void - StokesAssembler:: - reference_density_compressibility_term (const double pressure_scaling, - const bool /*rebuild_stokes_matrix*/, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const + StokesReferenceDensityCompressibilityTerm:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + // assemble RHS of: // - div u = 1/rho * drho/dz g/||g||* u - (void)parameters; - Assert(parameters.formulation_mass_conservation == + Assert(this->get_parameters().formulation_mass_conservation == Parameters::Formulation::MassConservation::reference_density_profile, ExcInternalError()); @@ -345,6 +353,7 @@ namespace aspect const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + const double pressure_scaling = scratch.pressure_scaling; for (unsigned int q=0; q void - StokesAssembler:: - implicit_reference_density_compressibility_term (const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const + StokesImplicitReferenceDensityCompressibilityTerm:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + // assemble compressibility term of: // - div u - 1/rho * drho/dz g/||g||* u = 0 - (void)parameters; - Assert(parameters.formulation_mass_conservation == + Assert(this->get_parameters().formulation_mass_conservation == Parameters::Formulation::MassConservation::implicit_reference_density_profile, ExcInternalError()); - if (!rebuild_stokes_matrix) + if (!scratch.rebuild_stokes_matrix) return; const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + const double pressure_scaling = scratch.pressure_scaling; for (unsigned int q=0; q void - StokesAssembler:: - isothermal_compression_term (const double pressure_scaling, - const bool /*rebuild_stokes_matrix*/, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data, - const Parameters ¶meters) const + StokesIsothermalCompressionTerm:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + // assemble RHS of: // - div u = 1/rho * drho/dp rho * g * u - (void)parameters; - Assert(parameters.formulation_mass_conservation == + Assert(this->get_parameters().formulation_mass_conservation == Parameters::Formulation::MassConservation::isothermal_compression, ExcInternalError()); @@ -448,6 +456,7 @@ namespace aspect const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + const double pressure_scaling = scratch.pressure_scaling; for (unsigned int q=0; q + void + StokesPressureRHSCompatibilityModification::execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const + { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + + const Introspection &introspection = this->introspection(); + const FiniteElement &fe = scratch.finite_element_values.get_fe(); + + const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); + const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + + for (unsigned int q=0; q + void + StokesBoundaryTraction::execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const + { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + + const Introspection &introspection = this->introspection(); + const FiniteElement &fe = scratch.finite_element_values.get_fe(); + + // see if any of the faces are traction boundaries for which + // we need to assemble force terms for the right hand side + const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); + + const typename DoFHandler::face_iterator face = scratch.cell->face(scratch.face_number); + + if (this->get_boundary_traction() + .find (face->boundary_id()) + != + this->get_boundary_traction().end()) + { + scratch.face_finite_element_values.reinit (scratch.cell, scratch.face_number); + + for (unsigned int q=0; q traction + = this->get_boundary_traction().find( + face->boundary_id() + )->second + ->boundary_traction (face->boundary_id(), + scratch.face_finite_element_values.quadrature_point(q), + scratch.face_finite_element_values.normal_vector(q)); + + for (unsigned int i=0, i_stokes=0; i_stokes; + template class StokesPreconditioner; \ + template class StokesCompressiblePreconditioner; \ + template class StokesIncompressibleTerms; \ + template class StokesCompressibleStrainRateViscosityTerm; \ + template class StokesReferenceDensityCompressibilityTerm; \ + template class StokesImplicitReferenceDensityCompressibilityTerm; \ + template class StokesIsothermalCompressionTerm; \ + template class StokesPressureRHSCompatibilityModification; \ + template class StokesBoundaryTraction; ASPECT_INSTANTIATE(INSTANTIATE) } diff --git a/source/simulator/assembly.cc b/source/simulator/assembly.cc index 6173591da5f..d98d921e6c8 100644 --- a/source/simulator/assembly.cc +++ b/source/simulator/assembly.cc @@ -21,11 +21,14 @@ #include #include -#include #include + +#include #include #include #include +#include +#include #include #include @@ -43,432 +46,6 @@ namespace aspect { - namespace internal - { - namespace Assembly - { - namespace Scratch - { - template - StokesPreconditioner:: - StokesPreconditioner (const FiniteElement &finite_element, - const Quadrature &quadrature, - const Mapping &mapping, - const UpdateFlags update_flags, - const unsigned int n_compositional_fields, - const unsigned int stokes_dofs_per_cell, - const bool add_compaction_pressure) - : - finite_element_values (mapping, finite_element, quadrature, - update_flags), - local_dof_indices (finite_element.dofs_per_cell), - dof_component_indices(stokes_dofs_per_cell), - grads_phi_u (stokes_dofs_per_cell, numbers::signaling_nan >()), - div_phi_u (stokes_dofs_per_cell, numbers::signaling_nan()), - phi_p (stokes_dofs_per_cell, numbers::signaling_nan()), - phi_p_c (add_compaction_pressure ? stokes_dofs_per_cell : 0, numbers::signaling_nan()), - grad_phi_p (add_compaction_pressure ? stokes_dofs_per_cell : 0, numbers::signaling_nan >()), - material_model_inputs(quadrature.size(), n_compositional_fields), - material_model_outputs(quadrature.size(), n_compositional_fields) - {} - - - - template - StokesPreconditioner:: - StokesPreconditioner (const StokesPreconditioner &scratch) - : - finite_element_values (scratch.finite_element_values.get_mapping(), - scratch.finite_element_values.get_fe(), - scratch.finite_element_values.get_quadrature(), - scratch.finite_element_values.get_update_flags()), - local_dof_indices (scratch.local_dof_indices), - dof_component_indices( scratch.dof_component_indices), - grads_phi_u (scratch.grads_phi_u), - div_phi_u (scratch.div_phi_u), - phi_p (scratch.phi_p), - phi_p_c (scratch.phi_p_c), - grad_phi_p(scratch.grad_phi_p), - material_model_inputs(scratch.material_model_inputs), - material_model_outputs(scratch.material_model_outputs) - {} - - - template - StokesPreconditioner:: - ~StokesPreconditioner () - {} - - - - - template - StokesSystem:: - StokesSystem (const FiniteElement &finite_element, - const Mapping &mapping, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const UpdateFlags update_flags, - const UpdateFlags face_update_flags, - const unsigned int n_compositional_fields, - const unsigned int stokes_dofs_per_cell, - const bool add_compaction_pressure, - const bool use_reference_density_profile) - : - StokesPreconditioner (finite_element, quadrature, - mapping, - update_flags, - n_compositional_fields, - stokes_dofs_per_cell, - add_compaction_pressure), - - face_finite_element_values (mapping, - finite_element, - face_quadrature, - face_update_flags), - - phi_u (stokes_dofs_per_cell, numbers::signaling_nan >()), - grads_phi_u (stokes_dofs_per_cell, numbers::signaling_nan >()), - velocity_values (quadrature.size(), numbers::signaling_nan >()), - velocity_divergence (quadrature.size(), numbers::signaling_nan()), - face_material_model_inputs(face_quadrature.size(), n_compositional_fields), - face_material_model_outputs(face_quadrature.size(), n_compositional_fields), - reference_densities(use_reference_density_profile ? quadrature.size() : 0, numbers::signaling_nan()), - reference_densities_depth_derivative(use_reference_density_profile ? quadrature.size() : 0, numbers::signaling_nan()) - {} - - - - template - StokesSystem:: - StokesSystem (const StokesSystem &scratch) - : - StokesPreconditioner (scratch), - - face_finite_element_values (scratch.face_finite_element_values.get_mapping(), - scratch.face_finite_element_values.get_fe(), - scratch.face_finite_element_values.get_quadrature(), - scratch.face_finite_element_values.get_update_flags()), - - phi_u (scratch.phi_u), - grads_phi_u (scratch.grads_phi_u), - velocity_values (scratch.velocity_values), - velocity_divergence (scratch.velocity_divergence), - face_material_model_inputs(scratch.face_material_model_inputs), - face_material_model_outputs(scratch.face_material_model_outputs), - reference_densities(scratch.reference_densities), - reference_densities_depth_derivative(scratch.reference_densities_depth_derivative) - {} - - - - template - AdvectionSystem:: - AdvectionSystem (const FiniteElement &finite_element, - const FiniteElement &advection_element, - const Mapping &mapping, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const UpdateFlags update_flags, - const UpdateFlags face_update_flags, - const unsigned int n_compositional_fields) - : - finite_element_values (mapping, - finite_element, quadrature, - update_flags), - face_finite_element_values ((face_quadrature.size() > 0 - ? - new FEFaceValues (mapping, - finite_element, face_quadrature, - face_update_flags) - : - NULL)), - neighbor_face_finite_element_values ((face_quadrature.size() > 0 - ? - new FEFaceValues (mapping, - finite_element, face_quadrature, - face_update_flags) - : - NULL)), - subface_finite_element_values ((face_quadrature.size() > 0 - ? - new FESubfaceValues (mapping, - finite_element, face_quadrature, - face_update_flags) - : - NULL)), - local_dof_indices (finite_element.dofs_per_cell), - - phi_field (advection_element.dofs_per_cell, numbers::signaling_nan()), - grad_phi_field (advection_element.dofs_per_cell, numbers::signaling_nan >()), - face_phi_field ((face_quadrature.size() > 0 ? advection_element.dofs_per_cell : 0), - numbers::signaling_nan()), - face_grad_phi_field ((face_quadrature.size() > 0 ? advection_element.dofs_per_cell : 0), - numbers::signaling_nan >()), - neighbor_face_phi_field ((face_quadrature.size() > 0 ? advection_element.dofs_per_cell : 0), - numbers::signaling_nan()), - neighbor_face_grad_phi_field ((face_quadrature.size() > 0 ? advection_element.dofs_per_cell : 0), - numbers::signaling_nan >()), - old_velocity_values (quadrature.size(), numbers::signaling_nan >()), - old_old_velocity_values (quadrature.size(), numbers::signaling_nan >()), - old_pressure (quadrature.size(), numbers::signaling_nan()), - old_old_pressure (quadrature.size(), numbers::signaling_nan()), - old_pressure_gradients (quadrature.size(), numbers::signaling_nan >()), - old_old_pressure_gradients (quadrature.size(), numbers::signaling_nan >()), - old_strain_rates (quadrature.size(), numbers::signaling_nan >()), - old_old_strain_rates (quadrature.size(), numbers::signaling_nan >()), - old_temperature_values (quadrature.size(), numbers::signaling_nan()), - old_old_temperature_values(quadrature.size(), numbers::signaling_nan()), - old_field_values (quadrature.size(), numbers::signaling_nan()), - old_old_field_values(quadrature.size(), numbers::signaling_nan()), - old_field_grads(quadrature.size(), numbers::signaling_nan >()), - old_old_field_grads(quadrature.size(), numbers::signaling_nan >()), - old_field_laplacians(quadrature.size(), numbers::signaling_nan()), - old_old_field_laplacians(quadrature.size(), numbers::signaling_nan()), - old_composition_values(n_compositional_fields, - std::vector(quadrature.size(), numbers::signaling_nan())), - old_old_composition_values(n_compositional_fields, - std::vector(quadrature.size(), numbers::signaling_nan())), - current_temperature_values(quadrature.size(), numbers::signaling_nan()), - current_velocity_values(quadrature.size(), numbers::signaling_nan >()), - face_current_velocity_values(face_quadrature.size(), numbers::signaling_nan >()), - mesh_velocity_values(quadrature.size(), numbers::signaling_nan >()), - face_mesh_velocity_values(face_quadrature.size(), numbers::signaling_nan >()), - current_strain_rates(quadrature.size(), numbers::signaling_nan >()), - current_composition_values(n_compositional_fields, - std::vector(quadrature.size(), numbers::signaling_nan())), - current_velocity_divergences(quadrature.size(), numbers::signaling_nan()), - material_model_inputs(quadrature.size(), n_compositional_fields), - material_model_outputs(quadrature.size(), n_compositional_fields), - face_material_model_inputs(face_quadrature.size(), n_compositional_fields), - face_material_model_outputs(face_quadrature.size(), n_compositional_fields), - neighbor_face_material_model_inputs(face_quadrature.size(), n_compositional_fields), - neighbor_face_material_model_outputs(face_quadrature.size(), n_compositional_fields), - heating_model_outputs(quadrature.size(), n_compositional_fields), - face_heating_model_outputs(face_quadrature.size(), n_compositional_fields), - neighbor_face_heating_model_outputs(face_quadrature.size(), n_compositional_fields) - {} - - - - template - AdvectionSystem:: - AdvectionSystem (const AdvectionSystem &scratch) - : - finite_element_values (scratch.finite_element_values.get_mapping(), - scratch.finite_element_values.get_fe(), - scratch.finite_element_values.get_quadrature(), - scratch.finite_element_values.get_update_flags()), - face_finite_element_values (scratch.face_finite_element_values), - neighbor_face_finite_element_values (scratch.neighbor_face_finite_element_values), - subface_finite_element_values (scratch.subface_finite_element_values), - local_dof_indices (scratch.finite_element_values.get_fe().dofs_per_cell), - - phi_field (scratch.phi_field), - grad_phi_field (scratch.grad_phi_field), - face_phi_field (scratch.face_phi_field), - face_grad_phi_field (scratch.face_grad_phi_field), - neighbor_face_phi_field (scratch.neighbor_face_phi_field), - neighbor_face_grad_phi_field (scratch.neighbor_face_grad_phi_field), - old_velocity_values (scratch.old_velocity_values), - old_old_velocity_values (scratch.old_old_velocity_values), - old_pressure (scratch.old_pressure), - old_old_pressure (scratch.old_old_pressure), - old_pressure_gradients (scratch.old_pressure_gradients), - old_old_pressure_gradients (scratch.old_old_pressure_gradients), - old_strain_rates (scratch.old_strain_rates), - old_old_strain_rates (scratch.old_old_strain_rates), - old_temperature_values (scratch.old_temperature_values), - old_old_temperature_values (scratch.old_old_temperature_values), - old_field_values (scratch.old_field_values), - old_old_field_values(scratch.old_old_field_values), - old_field_grads (scratch.old_field_grads), - old_old_field_grads (scratch.old_old_field_grads), - old_field_laplacians (scratch.old_field_laplacians), - old_old_field_laplacians (scratch.old_old_field_laplacians), - old_composition_values(scratch.old_composition_values), - old_old_composition_values(scratch.old_old_composition_values), - current_temperature_values(scratch.current_temperature_values), - current_velocity_values(scratch.current_velocity_values), - face_current_velocity_values(scratch.face_current_velocity_values), - mesh_velocity_values(scratch.mesh_velocity_values), - face_mesh_velocity_values(scratch.face_mesh_velocity_values), - current_strain_rates(scratch.current_strain_rates), - current_composition_values(scratch.current_composition_values), - current_velocity_divergences(scratch.current_velocity_divergences), - material_model_inputs(scratch.material_model_inputs), - material_model_outputs(scratch.material_model_outputs), - face_material_model_inputs(scratch.face_material_model_inputs), - face_material_model_outputs(scratch.face_material_model_outputs), - neighbor_face_material_model_inputs(scratch.neighbor_face_material_model_inputs), - neighbor_face_material_model_outputs(scratch.neighbor_face_material_model_outputs), - heating_model_outputs(scratch.heating_model_outputs), - face_heating_model_outputs(scratch.face_heating_model_outputs), - neighbor_face_heating_model_outputs(scratch.neighbor_face_heating_model_outputs) - {} - - } - - - namespace CopyData - { - - template - StokesPreconditioner:: - StokesPreconditioner (const unsigned int stokes_dofs_per_cell) - : - local_matrix (stokes_dofs_per_cell, - stokes_dofs_per_cell), - local_dof_indices (stokes_dofs_per_cell) - {} - - - - template - StokesPreconditioner:: - StokesPreconditioner (const StokesPreconditioner &data) - : - local_matrix (data.local_matrix), - local_dof_indices (data.local_dof_indices) - {} - - - template - void StokesPreconditioner:: - extract_stokes_dof_indices(const std::vector &all_dof_indices, - const Introspection &introspection, - const dealii::FiniteElement &finite_element) - { - const unsigned int dofs_per_cell = finite_element.dofs_per_cell; - - for (unsigned int i=0, i_stokes=0; ilocal_dof_indices[i_stokes] = all_dof_indices[i]; - ++i_stokes; - } - ++i; - } - } - - - template - StokesPreconditioner:: - ~StokesPreconditioner () - {} - - - - template - StokesSystem:: - StokesSystem (const unsigned int stokes_dofs_per_cell, - const bool do_pressure_rhs_compatibility_modification) - : - StokesPreconditioner (stokes_dofs_per_cell), - local_rhs (stokes_dofs_per_cell), - local_pressure_shape_function_integrals (do_pressure_rhs_compatibility_modification ? - stokes_dofs_per_cell - : - 0) - {} - - - - template - StokesSystem:: - StokesSystem (const StokesSystem &data) - : - StokesPreconditioner (data), - local_rhs (data.local_rhs), - local_pressure_shape_function_integrals (data.local_pressure_shape_function_integrals.size()) - {} - - - - - template - AdvectionSystem:: - AdvectionSystem (const FiniteElement &finite_element, - const bool field_is_discontinuous) - : - local_matrix (finite_element.dofs_per_cell, - finite_element.dofs_per_cell), - local_matrices_int_ext ((field_is_discontinuous - ? - GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell - : - 0), - FullMatrix(finite_element.dofs_per_cell, - finite_element.dofs_per_cell)), - local_matrices_ext_int ((field_is_discontinuous - ? - GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell - : - 0), - FullMatrix(finite_element.dofs_per_cell, - finite_element.dofs_per_cell)), - local_matrices_ext_ext ((field_is_discontinuous - ? - GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell - : - 0), - FullMatrix(finite_element.dofs_per_cell, - finite_element.dofs_per_cell)), - local_rhs (finite_element.dofs_per_cell), - - assembled_matrices ((field_is_discontinuous - ? - GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell - : - 0), false), - - local_dof_indices (finite_element.dofs_per_cell), - neighbor_dof_indices ((field_is_discontinuous - ? - GeometryInfo::max_children_per_face * GeometryInfo::faces_per_cell - : - 0), - std::vector(finite_element.dofs_per_cell)) - {} - - - - template - AdvectionSystem:: - AdvectionSystem (const AdvectionSystem &data) - : - local_matrix (data.local_matrix), - local_matrices_int_ext (data.local_matrices_int_ext), - local_matrices_ext_int (data.local_matrices_ext_int), - local_matrices_ext_ext (data.local_matrices_ext_ext), - local_rhs (data.local_rhs), - - assembled_matrices (data.assembled_matrices), - - local_dof_indices (data.local_dof_indices), - neighbor_dof_indices (data.neighbor_dof_indices) - {} - - } - - - - template - AssemblerLists::Properties::Properties () - : - need_face_material_model_data (false), - need_face_finite_element_evaluation(false), - need_viscosity(false), - needed_update_flags () - {} - - } - } - - - - template void Simulator:: @@ -525,216 +102,93 @@ namespace aspect - namespace Assemblers + namespace { - /** - * A namespace for the definition of functions that implement various - * other terms that need to occasionally or always be assembled. - */ - namespace OtherTerms + template + void + initialize_simulator(Simulator &simulator, + std::vector > &assemblers) { - template - void - pressure_rhs_compatibility_modification (const SimulatorAccess &simulator_access, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) - { - const Introspection &introspection = simulator_access.introspection(); - const FiniteElement &fe = simulator_access.get_fe(); - - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; - - for (unsigned int q=0; q - void - boundary_traction (const SimulatorAccess &simulator_access, - const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) - { - const Introspection &introspection = simulator_access.introspection(); - const FiniteElement &fe = scratch.finite_element_values.get_fe(); - - // see if any of the faces are traction boundaries for which - // we need to assemble force terms for the right hand side - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - - if (simulator_access.get_boundary_traction() - .find (cell->face(face_no)->boundary_id()) - != - simulator_access.get_boundary_traction().end()) - { - scratch.face_finite_element_values.reinit (cell, face_no); - - for (unsigned int q=0; q traction - = simulator_access.get_boundary_traction().find( - cell->face(face_no)->boundary_id() - )->second - ->boundary_traction (cell->face(face_no)->boundary_id(), - scratch.face_finite_element_values.quadrature_point(q), - scratch.face_finite_element_values.normal_vector(q)); - - for (unsigned int i=0, i_stokes=0; i_stokes *p = dynamic_cast* >(assemblers[i].get())) + p->initialize_simulator(simulator); } } - template - void - Simulator:: - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const - { - typedef typename std::vector > > - assembler_vector_t; - - for (typename assembler_vector_t::const_iterator it = assembler_objects.begin(); - it != assembler_objects.end(); - ++it) - { - (*it)->create_additional_material_model_outputs(outputs); - } - - } - - - template void Simulator:: set_assemblers () { - // create an object for the complete equations assembly; add its - // member functions to the signals and add the object the list - // of assembler objects - aspect::Assemblers::StokesAssembler *stokes_assembler = NULL; - aspect::Assemblers::NewtonStokesAssembler *newton_stokes_assembler = NULL; - - if (assemble_newton_stokes_system) - newton_stokes_assembler = new aspect::Assemblers::NewtonStokesAssembler(); - else - stokes_assembler = new aspect::Assemblers::StokesAssembler(); - - aspect::Assemblers::AdvectionAssembler *adv_assembler - = new aspect::Assemblers::AdvectionAssembler(); - assemblers->advection_system_assembler_properties.resize(1+introspection.n_compositional_fields); assemblers->advection_system_assembler_on_face_properties.resize(1+introspection.n_compositional_fields); - aspect::Assemblers::MeltEquations *melt_equation_assembler = NULL; - if (parameters.include_melt_transport) - melt_equation_assembler = new aspect::Assemblers::MeltEquations(); if (parameters.include_melt_transport) - assemblers->local_assemble_stokes_preconditioner - .connect (std_cxx11::bind(&aspect::Assemblers::MeltEquations::local_assemble_stokes_preconditioner_melt, - std_cxx11::cref (*melt_equation_assembler), - std_cxx11::_1, std_cxx11::_2, std_cxx11::_3)); + { + aspect::Assemblers::MeltStokesPreconditioner *preconditioner = new aspect::Assemblers::MeltStokesPreconditioner(); + assemblers->stokes_preconditioner.push_back( + std_cxx11::unique_ptr > (preconditioner)); + } else if (assemble_newton_stokes_system) - assemblers->local_assemble_stokes_preconditioner - .connect (std_cxx11::bind(&aspect::Assemblers::NewtonStokesAssembler::preconditioner, - std_cxx11::cref (*newton_stokes_assembler), - std_cxx11::_1, - std_cxx11::_2, - std_cxx11::_3)); + { + aspect::Assemblers::NewtonStokesPreconditioner *preconditioner = new aspect::Assemblers::NewtonStokesPreconditioner(); + assemblers->stokes_preconditioner.push_back( + std_cxx11::unique_ptr > (preconditioner)); + } else { - assemblers->local_assemble_stokes_preconditioner - .connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler::preconditioner, - std_cxx11::cref (*stokes_assembler), - std_cxx11::_1, std_cxx11::_2, std_cxx11::_3)); + aspect::Assemblers::StokesPreconditioner *preconditioner = new aspect::Assemblers::StokesPreconditioner(); + assemblers->stokes_preconditioner.push_back( + std_cxx11::unique_ptr > (preconditioner)); if (material_model->is_compressible()) - assemblers->local_assemble_stokes_preconditioner - .connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler::compressible_preconditioner, - std_cxx11::cref (*stokes_assembler), - std_cxx11::_1, std_cxx11::_2, std_cxx11::_3)); + { + aspect::Assemblers::StokesCompressiblePreconditioner *preconditioner = new aspect::Assemblers::StokesCompressiblePreconditioner(); + assemblers->stokes_preconditioner.push_back( + std_cxx11::unique_ptr > (preconditioner)); + } } if (parameters.include_melt_transport) - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::MeltEquations::local_assemble_stokes_system_melt, - std_cxx11::cref (*melt_equation_assembler), - std_cxx11::_1, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5)); - + { + aspect::Assemblers::MeltStokesSystem *melt_stokes_system = new aspect::Assemblers::MeltStokesSystem(); + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (melt_stokes_system)); + } else if (assemble_newton_stokes_system) { - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::NewtonStokesAssembler::incompressible_terms, - std_cxx11::cref (*newton_stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5)); + aspect::Assemblers::NewtonStokesIncompressibleTerms *incompressible_terms = new aspect::Assemblers::NewtonStokesIncompressibleTerms(); + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (incompressible_terms)); if (material_model->is_compressible()) - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::NewtonStokesAssembler::compressible_strain_rate_viscosity_term, - std_cxx11::cref (*newton_stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5)); + { + aspect::Assemblers::NewtonStokesCompressibleStrainRateViscosityTerm *compressible_terms + = new aspect::Assemblers::NewtonStokesCompressibleStrainRateViscosityTerm(); + + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (compressible_terms)); + } if (parameters.formulation_mass_conservation == Parameters::Formulation::MassConservation::implicit_reference_density_profile) - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::NewtonStokesAssembler::implicit_reference_density_compressibility_term, - std_cxx11::cref (*newton_stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5, - std_cxx11::cref (this->parameters))); + { + aspect::Assemblers::NewtonStokesImplicitReferenceDensityCompressibilityTerm *compressible_terms + = new aspect::Assemblers::NewtonStokesImplicitReferenceDensityCompressibilityTerm(); + + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (compressible_terms)); + } else if (parameters.formulation_mass_conservation == Parameters::Formulation::MassConservation::reference_density_profile) { - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::NewtonStokesAssembler::reference_density_compressibility_term, - std_cxx11::cref (*newton_stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5, - std_cxx11::cref (this->parameters))); + aspect::Assemblers::NewtonStokesReferenceDensityCompressibilityTerm *compressible_terms + = new aspect::Assemblers::NewtonStokesReferenceDensityCompressibilityTerm(); + + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (compressible_terms)); } else if (parameters.formulation_mass_conservation == Parameters::Formulation::MassConservation::incompressible) @@ -742,61 +196,47 @@ namespace aspect // do nothing, because we assembled div u =0 above already } else - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::NewtonStokesAssembler::isothermal_compression_term, - std_cxx11::cref (*newton_stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5, - std_cxx11::cref (this->parameters))); + { + aspect::Assemblers::NewtonStokesIsothermalCompressionTerm *compressible_terms + = new aspect::Assemblers::NewtonStokesIsothermalCompressionTerm(); + + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (compressible_terms)); + } } else { - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler::incompressible_terms, - std_cxx11::cref (*stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5)); + aspect::Assemblers::StokesIncompressibleTerms *incompressible_terms = new aspect::Assemblers::StokesIncompressibleTerms(); + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (incompressible_terms)); if (material_model->is_compressible()) - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler::compressible_strain_rate_viscosity_term, - std_cxx11::cref (*stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5)); + { + aspect::Assemblers::StokesCompressibleStrainRateViscosityTerm *compressible_terms + = new aspect::Assemblers::StokesCompressibleStrainRateViscosityTerm(); + + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (compressible_terms)); + } if (parameters.formulation_mass_conservation == Parameters::Formulation::MassConservation::implicit_reference_density_profile) - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler::implicit_reference_density_compressibility_term, - std_cxx11::cref (*stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5, - std_cxx11::cref (this->parameters))); + { + aspect::Assemblers::StokesImplicitReferenceDensityCompressibilityTerm *compressible_terms + = new aspect::Assemblers::StokesImplicitReferenceDensityCompressibilityTerm(); + + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (compressible_terms)); + } else if (parameters.formulation_mass_conservation == Parameters::Formulation::MassConservation::reference_density_profile) { - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler::reference_density_compressibility_term, - std_cxx11::cref (*stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5, - std_cxx11::cref (this->parameters))); + aspect::Assemblers::StokesReferenceDensityCompressibilityTerm *compressible_terms + = new aspect::Assemblers::StokesReferenceDensityCompressibilityTerm(); + + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (compressible_terms)); } else if (parameters.formulation_mass_conservation == Parameters::Formulation::MassConservation::incompressible) @@ -804,32 +244,15 @@ namespace aspect // do nothing, because we assembled div u =0 above already } else - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler::isothermal_compression_term, - std_cxx11::cref (*stokes_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5, - std_cxx11::cref (this->parameters))); - - } - - if (assemble_newton_stokes_system) - assembler_objects.push_back (std_cxx11::shared_ptr > - (newton_stokes_assembler)); - else - assembler_objects.push_back (std_cxx11::shared_ptr > - (stokes_assembler)); - - assembler_objects.push_back (std_cxx11::shared_ptr > - (adv_assembler)); + { + aspect::Assemblers::StokesIsothermalCompressionTerm *compressible_terms + = new aspect::Assemblers::StokesIsothermalCompressionTerm(); + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (compressible_terms)); + } - if (parameters.include_melt_transport) - assembler_objects.push_back (std_cxx11::shared_ptr > - (melt_equation_assembler)); + } // add the boundary integral for melt migration if (parameters.include_melt_transport) @@ -838,123 +261,94 @@ namespace aspect assemblers->stokes_system_assembler_on_boundary_face_properties.needed_update_flags = (update_values | update_quadrature_points | update_normal_vectors | update_gradients | update_JxW_values); - assemblers->local_assemble_stokes_system_on_boundary_face - .connect (std_cxx11::bind(&aspect::Assemblers::MeltEquations::local_assemble_stokes_system_melt_boundary, - std_cxx11::cref (*melt_equation_assembler), - std_cxx11::_1, - std_cxx11::_2, - std_cxx11::_3, - // discard rebuild_stokes_matrix, - std_cxx11::_5, - std_cxx11::_6)); + + aspect::Assemblers::MeltStokesSystemBoundary *melt_stokes_system_boundary = new aspect::Assemblers::MeltStokesSystemBoundary(); + assemblers->stokes_system_on_boundary_face.push_back( + std_cxx11::unique_ptr > (melt_stokes_system_boundary)); } // add the terms for traction boundary conditions if (parameters.include_melt_transport) { - assemblers->local_assemble_stokes_system_on_boundary_face - .connect (std_cxx11::bind(&aspect::Assemblers::OtherTerms::boundary_traction_melt, - SimulatorAccess(*this), - std_cxx11::_1, - std_cxx11::_2, - // discard pressure_scaling, - // discard rebuild_stokes_matrix, - std_cxx11::_5, - std_cxx11::_6)); + aspect::Assemblers::MeltBoundaryTraction *stokes_boundary_traction + = new aspect::Assemblers::MeltBoundaryTraction(); + + assemblers->stokes_system_on_boundary_face.push_back( + std_cxx11::unique_ptr > (stokes_boundary_traction)); } else { - assemblers->local_assemble_stokes_system_on_boundary_face - .connect (std_cxx11::bind(&aspect::Assemblers::OtherTerms::boundary_traction, - SimulatorAccess(*this), - std_cxx11::_1, - std_cxx11::_2, - // discard pressure_scaling, - // discard rebuild_stokes_matrix, - std_cxx11::_5, - std_cxx11::_6)); + aspect::Assemblers::StokesBoundaryTraction *stokes_boundary_traction + = new aspect::Assemblers::StokesBoundaryTraction(); + + assemblers->stokes_system_on_boundary_face.push_back( + std_cxx11::unique_ptr > (stokes_boundary_traction)); } // add the terms necessary to normalize the pressure if (do_pressure_rhs_compatibility_modification) { if (parameters.include_melt_transport) - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::OtherTerms::pressure_rhs_compatibility_modification_melt, - SimulatorAccess(*this), - // discard cell, - // discard pressure_scaling, - // discard rebuild_stokes_matrix, - std_cxx11::_4, - std_cxx11::_5)); + { + aspect::Assemblers::MeltPressureRHSCompatibilityModification *melt_pressure_RHS_modification + = new aspect::Assemblers::MeltPressureRHSCompatibilityModification(); + + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (melt_pressure_RHS_modification)); + } else - assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&aspect::Assemblers::OtherTerms::pressure_rhs_compatibility_modification, - SimulatorAccess(*this), - // discard cell, - // discard pressure_scaling, - // discard rebuild_stokes_matrix, - std_cxx11::_4, - std_cxx11::_5)); + { + aspect::Assemblers::StokesPressureRHSCompatibilityModification *stokes_pressure_RHS_modification + = new aspect::Assemblers::StokesPressureRHSCompatibilityModification(); + + assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (stokes_pressure_RHS_modification)); + } } if (parameters.include_melt_transport) { - assemblers->local_assemble_advection_system - .connect (std_cxx11::bind(&aspect::Assemblers::MeltEquations::local_assemble_advection_system_melt, - std_cxx11::cref (*melt_equation_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5)); - - assemblers->compute_advection_system_residual - .connect (std_cxx11::bind(&aspect::Assemblers::MeltEquations::compute_advection_system_residual_melt, - std_cxx11::cref (*melt_equation_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3)); + aspect::Assemblers::MeltAdvectionSystem *melt_advection_system + = new aspect::Assemblers::MeltAdvectionSystem(); + + 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 { - assemblers->local_assemble_advection_system - .connect (std_cxx11::bind(&aspect::Assemblers::AdvectionAssembler::local_assemble_advection_system, - std_cxx11::cref (*adv_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5)); - - assemblers->compute_advection_system_residual - .connect (std_cxx11::bind(&aspect::Assemblers::AdvectionAssembler::compute_advection_system_residual, - std_cxx11::cref (*adv_assembler), - // discard cell, - std_cxx11::_2, - std_cxx11::_3)); + aspect::Assemblers::AdvectionSystem *advection_system + = new aspect::Assemblers::AdvectionSystem(); + + 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 || parameters.use_discontinuous_composition_discretization) { - assemblers->local_assemble_advection_system_on_interior_face - .connect(std_cxx11::bind(&aspect::Assemblers::AdvectionAssembler::local_assemble_discontinuous_advection_interior_face_terms, - std_cxx11::cref (*adv_assembler), - std_cxx11::_1, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5)); - - assemblers->local_assemble_advection_system_on_boundary_face - .connect(std_cxx11::bind(&aspect::Assemblers::AdvectionAssembler::local_assemble_discontinuous_advection_boundary_face_terms, - std_cxx11::cref (*adv_assembler), - std_cxx11::_1, - std_cxx11::_2, - std_cxx11::_3, - std_cxx11::_4, - std_cxx11::_5)); + aspect::Assemblers::AdvectionSystemBoundaryFace *advection_system_boundary_face + = new aspect::Assemblers::AdvectionSystemBoundaryFace(); + + assemblers->advection_system_on_boundary_face.push_back( + std_cxx11::unique_ptr > (advection_system_boundary_face)); + + aspect::Assemblers::AdvectionSystemInteriorFace *advection_system_interior_face + = new aspect::Assemblers::AdvectionSystemInteriorFace(); + + assemblers->advection_system_on_interior_face.push_back( + std_cxx11::unique_ptr > (advection_system_interior_face)); if (parameters.use_discontinuous_temperature_discretization) { @@ -973,13 +367,17 @@ namespace aspect } // allow other assemblers to add themselves or modify the existing ones by firing the signal - this->signals.set_assemblers(*this, *assemblers, assembler_objects); + this->signals.set_assemblers(*this, *assemblers); // ensure that all assembler objects have access to the SimulatorAccess // base class - for (unsigned int i=0; i *p = dynamic_cast*>(assembler_objects[i].get())) - p->initialize_simulator(*this); + 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); + initialize_simulator(*this,assemblers->advection_system_residual); } @@ -1001,6 +399,9 @@ namespace aspect // Prepare the data structures for assembly scratch.finite_element_values.reinit (cell); + scratch.cell = cell; + scratch.pressure_scaling = material_model->reference_viscosity() + / geometry_model->length_scale(); data.local_matrix = 0; @@ -1009,7 +410,9 @@ namespace aspect cell, true, scratch.material_model_inputs); - create_additional_material_model_outputs(scratch.material_model_outputs); + + for (unsigned int i=0; istokes_preconditioner.size(); ++i) + assemblers->stokes_preconditioner[i]->create_additional_material_model_outputs(scratch.material_model_outputs); material_model->evaluate(scratch.material_model_inputs, scratch.material_model_outputs); @@ -1019,9 +422,8 @@ namespace aspect scratch.finite_element_values.get_mapping(), scratch.material_model_outputs); - // trigger the invocation of the various functions that actually do - // all of the assembling - assemblers->local_assemble_stokes_preconditioner(pressure_scaling, scratch, data); + for (unsigned int i=0; istokes_preconditioner.size(); ++i) + assemblers->stokes_preconditioner[i]->execute(scratch,data); } @@ -1086,7 +488,8 @@ namespace aspect cell_update_flags, introspection.n_compositional_fields, stokes_dofs_per_cell, - parameters.include_melt_transport), + parameters.include_melt_transport, + rebuild_stokes_matrix), internal::Assembly::CopyData:: StokesPreconditioner (stokes_dofs_per_cell)); @@ -1189,6 +592,9 @@ namespace aspect // Prepare the data structures for assembly scratch.finite_element_values.reinit (cell); + scratch.cell = cell; + scratch.pressure_scaling = material_model->reference_viscosity() + / geometry_model->length_scale(); if (rebuild_stokes_matrix) data.local_matrix = 0; @@ -1202,7 +608,9 @@ namespace aspect cell, assemble_newton_stokes_system ? true : rebuild_stokes_matrix, scratch.material_model_inputs); - create_additional_material_model_outputs(scratch.material_model_outputs); + + for (unsigned int i=0; istokes_system.size(); ++i) + assemblers->stokes_system[i]->create_additional_material_model_outputs(scratch.material_model_outputs); material_model->evaluate(scratch.material_model_inputs, scratch.material_model_outputs); @@ -1231,16 +639,15 @@ namespace aspect // trigger the invocation of the various functions that actually do // all of the assembling - assemblers->local_assemble_stokes_system(cell, pressure_scaling, - (!assemble_newton_stokes_system && rebuild_stokes_matrix) || (assemble_newton_stokes_system && assemble_newton_stokes_matrix), - scratch, data); + for (unsigned int i=0; istokes_system.size(); ++i) + assemblers->stokes_system[i]->execute(scratch,data); // then also work on possible face terms. if necessary, initialize // the material model data on faces - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - if (cell->at_boundary(face_no)) + for (scratch.face_number=0; scratch.face_number::faces_per_cell; ++scratch.face_number) + if (cell->at_boundary(scratch.face_number)) { - scratch.face_finite_element_values.reinit (cell, face_no); + scratch.face_finite_element_values.reinit (cell, scratch.face_number); if (assemblers->stokes_system_assembler_on_boundary_face_properties.need_face_material_model_data) { @@ -1252,7 +659,9 @@ namespace aspect cell, need_viscosity, scratch.face_material_model_inputs); - create_additional_material_model_outputs(scratch.face_material_model_outputs); + + for (unsigned int i=0; istokes_system_on_boundary_face.size(); ++i) + assemblers->stokes_system_on_boundary_face[i]->create_additional_material_model_outputs(scratch.face_material_model_outputs); material_model->evaluate(scratch.face_material_model_inputs, scratch.face_material_model_outputs); @@ -1272,10 +681,8 @@ namespace aspect // scratch.face_material_model_outputs); } - assemblers->local_assemble_stokes_system_on_boundary_face(cell, face_no, - pressure_scaling, - (!assemble_newton_stokes_system && rebuild_stokes_matrix) || (assemble_newton_stokes_system && assemble_newton_stokes_matrix), - scratch, data); + for (unsigned int i=0; istokes_system_on_boundary_face.size(); ++i) + assemblers->stokes_system_on_boundary_face[i]->execute(scratch,data); } } @@ -1397,7 +804,9 @@ namespace aspect introspection.n_compositional_fields, stokes_dofs_per_cell, parameters.include_melt_transport, - use_reference_density_profile), + use_reference_density_profile, + rebuild_stokes_matrix, + assemble_newton_stokes_matrix), internal::Assembly::CopyData:: StokesSystem (stokes_dofs_per_cell, do_pressure_rhs_compatibility_modification)); @@ -1472,6 +881,7 @@ namespace aspect const unsigned int solution_component = advection_field.component_index(introspection); scratch.finite_element_values.reinit (cell); + scratch.cell = cell; // get all dof indices on the current cell, then extract those // that correspond to the solution_field we are interested in @@ -1513,7 +923,9 @@ namespace aspect cell, true, scratch.material_model_inputs); - 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); @@ -1560,18 +972,19 @@ namespace aspect // TODO: Compute artificial viscosity once per timestep instead of each time // temperature system is assembled (as this might happen more than once per // timestep for iterative solvers) - const double nu = viscosity_per_cell[cell->active_cell_index()]; - Assert (nu >= 0, ExcMessage ("The artificial viscosity needs to be a non-negative quantity.")); + scratch.artificial_viscosity = viscosity_per_cell[cell->active_cell_index()]; + Assert (scratch.artificial_viscosity >= 0, ExcMessage ("The artificial viscosity needs to be a non-negative quantity.")); // trigger the invocation of the various functions that actually do // all of the assembling - assemblers->local_assemble_advection_system(cell, advection_field, nu, scratch, data); + for (unsigned int i=0; iadvection_system.size(); ++i) + assemblers->advection_system[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->local_assemble_advection_system_on_boundary_face.empty() + const bool has_boundary_face_assemblers = !assemblers->advection_system_on_boundary_face.empty() && assemblers->advection_system_assembler_on_face_properties[advection_field.field_index()].need_face_finite_element_evaluation; - const bool has_interior_face_assemblers = !assemblers->local_assemble_advection_system_on_interior_face.empty() + const bool has_interior_face_assemblers = !assemblers->advection_system_on_interior_face.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 @@ -1591,14 +1004,14 @@ namespace aspect } } - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + for (scratch.face_number=0; scratch.face_number::faces_per_cell; ++scratch.face_number) { - const typename DoFHandler::face_iterator face = cell->face (face_no); + typename DoFHandler::face_iterator face = cell->face (scratch.face_number); if ((has_boundary_face_assemblers && face->at_boundary()) || (has_interior_face_assemblers && !face->at_boundary())) { - (*scratch.face_finite_element_values).reinit (cell, face_no); + (*scratch.face_finite_element_values).reinit (cell, scratch.face_number); (*scratch.face_finite_element_values)[introspection.extractors.velocities].get_function_values(current_linearization_point, scratch.face_current_velocity_values); @@ -1616,7 +1029,11 @@ namespace aspect true, scratch.face_material_model_inputs); - create_additional_material_model_outputs(scratch.face_material_model_outputs); + for (unsigned int i=0; iadvection_system_on_boundary_face.size(); ++i) + assemblers->advection_system_on_boundary_face[i]->create_additional_material_model_outputs(scratch.material_model_outputs); + + for (unsigned int i=0; iadvection_system_on_interior_face.size(); ++i) + assemblers->advection_system_on_interior_face[i]->create_additional_material_model_outputs(scratch.material_model_outputs); material_model->evaluate(scratch.face_material_model_inputs, scratch.face_material_model_outputs); @@ -1633,14 +1050,13 @@ namespace aspect // scratch.face_material_model_outputs); } - if (face->at_boundary()) - assemblers->local_assemble_advection_system_on_boundary_face(cell, face_no, - advection_field, - scratch, data); + // handle faces at periodic boundaries like interior faces + if (face->at_boundary() && !cell->has_periodic_neighbor (scratch.face_number)) + for (unsigned int i=0; iadvection_system_on_boundary_face.size(); ++i) + assemblers->advection_system_on_boundary_face[i]->execute(scratch,data); else - assemblers->local_assemble_advection_system_on_interior_face(cell, face_no, - advection_field, - scratch, data); + for (unsigned int i=0; iadvection_system_on_interior_face.size(); ++i) + assemblers->advection_system_on_interior_face[i]->execute(scratch,data); } } } @@ -1662,7 +1078,7 @@ namespace aspect /* In the following, we copy DG contributions element by element. This * is allowed since there are no constraints imposed on discontinuous fields. */ - if (!assemblers->local_assemble_advection_system_on_interior_face.empty() && + if (!assemblers->advection_system_on_interior_face.empty() && assemblers->advection_system_assembler_on_face_properties[advection_field.field_index()].need_face_finite_element_evaluation) { for (unsigned int f=0; f::max_children_per_face @@ -1729,10 +1145,10 @@ namespace aspect + (parameters.stokes_velocity_degree+1)/2; - const bool allocate_face_quadrature = (!assemblers->local_assemble_advection_system_on_boundary_face.empty() || - !assemblers->local_assemble_advection_system_on_interior_face.empty()) && + const bool allocate_face_quadrature = (!assemblers->advection_system_on_boundary_face.empty() || + !assemblers->advection_system_on_interior_face.empty()) && assemblers->advection_system_assembler_on_face_properties[advection_field.field_index()].need_face_finite_element_evaluation; - const bool allocate_neighbor_contributions = !assemblers->local_assemble_advection_system_on_interior_face.empty() && + const bool allocate_neighbor_contributions = !assemblers->advection_system_on_interior_face.empty() && assemblers->advection_system_assembler_on_face_properties[advection_field.field_index()].need_face_finite_element_evaluation;; const UpdateFlags update_flags = update_values | @@ -1780,7 +1196,8 @@ namespace aspect Quadrature ()), update_flags, face_update_flags, - introspection.n_compositional_fields), + introspection.n_compositional_fields, + advection_field), internal::Assembly::CopyData:: AdvectionSystem (finite_element.base_element(advection_field.base_element(introspection)), allocate_neighbor_contributions)); @@ -1798,15 +1215,6 @@ namespace aspect namespace aspect { #define INSTANTIATE(dim) \ - namespace internal { \ - namespace Assembly { \ - namespace Scratch { \ - template struct AdvectionSystem; \ - template struct StokesSystem; \ - } \ - template struct AssemblerLists::Properties; \ - } \ - } \ template void Simulator::set_assemblers (); \ template void Simulator::local_assemble_stokes_preconditioner ( \ const DoFHandler::active_cell_iterator &cell, \ @@ -1840,9 +1248,8 @@ namespace aspect const FEValuesBase &input_finite_element_values, \ const DoFHandler::active_cell_iterator &cell, \ const bool compute_strainrate, \ - MaterialModel::MaterialModelInputs &material_model_inputs) const; \ - template void Simulator::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const; \ - + MaterialModel::MaterialModelInputs &material_model_inputs) const; + ASPECT_INSTANTIATE(INSTANTIATE) diff --git a/source/simulator/core.cc b/source/simulator/core.cc index d2cc7713163..ce2e743fdfd 100644 --- a/source/simulator/core.cc +++ b/source/simulator/core.cc @@ -21,12 +21,12 @@ #include #include -#include #include #include #include #include +#include #include #include @@ -148,7 +148,7 @@ namespace aspect Simulator::Simulator (const MPI_Comm mpi_communicator_, ParameterHandler &prm) : - assemblers (new internal::Assembly::AssemblerLists()), + assemblers (new Assemblers::Manager()), parameters (prm, mpi_communicator_), melt_handler (parameters.include_melt_transport ? new MeltHandler (prm) : NULL), newton_handler (parameters.nonlinear_solver == NonlinearSolver::Newton_Stokes ? new NewtonHandler () : NULL), diff --git a/source/simulator/entropy_viscosity.cc b/source/simulator/entropy_viscosity.cc index 5f1923bdec0..c5d00daaf71 100644 --- a/source/simulator/entropy_viscosity.cc +++ b/source/simulator/entropy_viscosity.cc @@ -19,7 +19,7 @@ */ #include -#include +#include #include #include @@ -128,9 +128,15 @@ namespace aspect if (advection_field.is_discontinuous(introspection)) return 0.; - std::vector residual = assemblers->compute_advection_system_residual(scratch.material_model_inputs.current_cell, - advection_field, - scratch); + std::vector residual = assemblers->advection_system_residual[0]->execute(scratch); + + for (unsigned int i=1; iadvection_system_residual.size(); ++i) + { + std::vector new_residual = assemblers->advection_system_residual[i]->execute(scratch); + for (unsigned int i=0; i (), update_flags, face_update_flags, - introspection.n_compositional_fields); + introspection.n_compositional_fields, + advection_field); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); for (; celladvection_system_residual.size(); ++i) + assemblers->advection_system_residual[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/free_surface.cc b/source/simulator/free_surface.cc index 5cca37b8797..903819e2172 100644 --- a/source/simulator/free_surface.cc +++ b/source/simulator/free_surface.cc @@ -22,7 +22,6 @@ #include #include #include -#include #include #include @@ -37,9 +36,96 @@ - namespace aspect { + namespace Assemblers + { + template + void + ApplyStabilization:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const + { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + + if (!this->get_parameters().free_surface_enabled) + return; + + const Introspection &introspection = this->introspection(); + const FiniteElement &fe = this->get_fe(); + + const typename DoFHandler::active_cell_iterator cell (&this->get_triangulation(), + scratch.finite_element_values.get_cell()->level(), + scratch.finite_element_values.get_cell()->index(), + &this->get_dof_handler()); + + const unsigned int n_face_q_points = scratch.face_finite_element_values.n_quadrature_points; + const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); + + // only apply on free surface faces + if (cell->at_boundary() && cell->is_locally_owned()) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + if (cell->face(face_no)->at_boundary()) + { + const types::boundary_id boundary_indicator + = cell->face(face_no)->boundary_id(); + + if (this->get_parameters().free_surface_boundary_indicators.find(boundary_indicator) + == this->get_parameters().free_surface_boundary_indicators.end()) + continue; + + scratch.face_finite_element_values.reinit(cell, face_no); + + this->compute_material_model_input_values (this->get_solution(), + scratch.face_finite_element_values, + cell, + false, + scratch.face_material_model_inputs); + + this->get_material_model().evaluate(scratch.face_material_model_inputs, scratch.face_material_model_outputs); + + for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) + { + for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) + { + if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) + { + scratch.phi_u[i_stokes] = scratch.face_finite_element_values[introspection.extractors.velocities].value(i, q_point); + ++i_stokes; + } + ++i; + } + + const Tensor<1,dim> + gravity = this->get_gravity_model().gravity_vector(scratch.face_finite_element_values.quadrature_point(q_point)); + double g_norm = gravity.norm(); + + // construct the relevant vectors + const Tensor<1,dim> n_hat = scratch.face_finite_element_values.normal_vector(q_point); + const Tensor<1,dim> g_hat = (g_norm == 0.0 ? Tensor<1,dim>() : gravity/g_norm); + + const double pressure_perturbation = scratch.face_material_model_outputs.densities[q_point] * + this->get_timestep() * + this->get_free_surface_handler().get_stabilization_term() * + g_norm; + + // see Kaus et al 2010 for details of the stabilization term + for (unsigned int i=0; i< stokes_dofs_per_cell; ++i) + for (unsigned int j=0; j< stokes_dofs_per_cell; ++j) + { + // The fictive stabilization stress is (phi_u[i].g)*(phi_u[j].n) + const double stress_value = -pressure_perturbation* + (scratch.phi_u[i]*g_hat) * (scratch.phi_u[j]*n_hat) + *scratch.face_finite_element_values.JxW(q_point); + + data.local_matrix(i,j) += stress_value; + } + } + } + } + } + template FreeSurfaceHandler::FreeSurfaceHandler (Simulator &simulator, ParameterHandler &prm) @@ -49,15 +135,11 @@ namespace aspect { parse_parameters(prm); - assembler_connection = - sim.assemblers->local_assemble_stokes_system - .connect (std_cxx11::bind(&FreeSurfaceHandler::apply_stabilization, - std_cxx11::ref(*this), - std_cxx11::_1, - // discard pressure_scaling, - // discard rebuild_stokes_matrix, - std_cxx11::_4, - std_cxx11::_5)); + aspect::Assemblers::ApplyStabilization *surface_stabilization + = new aspect::Assemblers::ApplyStabilization(); + + sim.assemblers->stokes_system.push_back( + std_cxx11::unique_ptr > (surface_stabilization)); // Note that we do not want face_material_model_data, because we do not // connect to a face assembler. We instead connect to a normal assembler, @@ -657,81 +739,9 @@ namespace aspect } template - void - FreeSurfaceHandler:: - apply_stabilization (const typename DoFHandler::active_cell_iterator &cell, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) + double FreeSurfaceHandler::get_stabilization_term() const { - if (!sim.parameters.free_surface_enabled) - return; - - const Introspection &introspection = sim.introspection; - const FiniteElement &fe = sim.finite_element; - - const unsigned int n_face_q_points = scratch.face_finite_element_values.n_quadrature_points; - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - - // only apply on free surface faces - if (cell->at_boundary() && cell->is_locally_owned()) - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - if (cell->face(face_no)->at_boundary()) - { - const types::boundary_id boundary_indicator - = cell->face(face_no)->boundary_id(); - - if (sim.parameters.free_surface_boundary_indicators.find(boundary_indicator) - == sim.parameters.free_surface_boundary_indicators.end()) - continue; - - scratch.face_finite_element_values.reinit(cell, face_no); - - sim.compute_material_model_input_values (sim.solution, - scratch.face_finite_element_values, - cell, - false, - scratch.face_material_model_inputs); - - sim.material_model->evaluate(scratch.face_material_model_inputs, scratch.face_material_model_outputs); - - for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) - { - for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) - { - if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) - { - scratch.phi_u[i_stokes] = scratch.face_finite_element_values[introspection.extractors.velocities].value(i, q_point); - ++i_stokes; - } - ++i; - } - - const Tensor<1,dim> - gravity = sim.gravity_model->gravity_vector(scratch.face_finite_element_values.quadrature_point(q_point)); - double g_norm = gravity.norm(); - - // construct the relevant vectors - const Tensor<1,dim> n_hat = scratch.face_finite_element_values.normal_vector(q_point); - const Tensor<1,dim> g_hat = (g_norm == 0.0 ? Tensor<1,dim>() : gravity/g_norm); - - double pressure_perturbation = scratch.face_material_model_outputs.densities[q_point] * - sim.time_step * free_surface_theta * g_norm; - - // see Kaus et al 2010 for details of the stabilization term - for (unsigned int i=0; i< stokes_dofs_per_cell; ++i) - for (unsigned int j=0; j< stokes_dofs_per_cell; ++j) - { - // The fictive stabilization stress is (phi_u[i].g)*(phi_u[j].n) - const double stress_value = -pressure_perturbation* - (scratch.phi_u[i]*g_hat) * (scratch.phi_u[j]*n_hat) - *scratch.face_finite_element_values.JxW(q_point); - - data.local_matrix(i,j) += stress_value; - } - } - } - - + return free_surface_theta; } } diff --git a/source/simulator/melt.cc b/source/simulator/melt.cc index 6417b926181..823e86664be 100644 --- a/source/simulator/melt.cc +++ b/source/simulator/melt.cc @@ -86,7 +86,7 @@ namespace aspect template void - MeltEquations::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const + MeltInterface::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const { MeltHandler::create_material_model_outputs(outputs); @@ -109,15 +109,19 @@ namespace aspect template void - MeltEquations:: - local_assemble_stokes_preconditioner_melt (const double pressure_scaling, - internal::Assembly::Scratch::StokesPreconditioner &scratch, - internal::Assembly::CopyData::StokesPreconditioner &data) const + MeltStokesPreconditioner:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + const double pressure_scaling = scratch.pressure_scaling; + MaterialModel::MeltOutputs *melt_outputs = scratch.material_model_outputs.template get_additional_output >(); @@ -218,17 +222,81 @@ namespace aspect } } + namespace + { + /** + * Compute the right-hand side of the fluid pressure equation of the Stokes + * system in case the simulation uses melt transport. This includes a term + * derived from Darcy's law, a term including the melting rate and a term dependent + * on the densities and velocities of fluid and solid. + */ + template + double + compute_fluid_pressure_rhs(const SimulatorAccess *simulator_access, + const internal::Assembly::Scratch::StokesSystem &scratch, + const unsigned int q_point) + { + if (!simulator_access->get_parameters().include_melt_transport) + return 0.0; + + const Introspection &introspection = simulator_access->introspection(); + const MaterialModel::MeltOutputs *melt_out = scratch.material_model_outputs.template get_additional_output >(); + + Assert(melt_out != NULL, ExcInternalError()); + + const unsigned int porosity_index = introspection.compositional_index_for_name("porosity"); + const unsigned int is_compressible = simulator_access->get_material_model().is_compressible(); + + const double melting_rate = scratch.material_model_outputs.reaction_terms[q_point][porosity_index]; + const double solid_density = scratch.material_model_outputs.densities[q_point]; + const double fluid_density = melt_out->fluid_densities[q_point]; + const double solid_compressibility = scratch.material_model_outputs.compressibilities[q_point]; + const Tensor<1,dim> fluid_density_gradient = melt_out->fluid_density_gradients[q_point]; + const Tensor<1,dim> current_u = scratch.velocity_values[q_point]; + const double porosity = std::max(scratch.material_model_inputs.composition[q_point][porosity_index],0.0); + const double K_D = (porosity > simulator_access->get_melt_handler().melt_transport_threshold + ? + melt_out->permeabilities[q_point] / melt_out->fluid_viscosities[q_point] + : + 0.0); + + const Tensor<1,dim> + gravity = simulator_access->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q_point)); + + double fluid_pressure_RHS = 0.0; + + // melting term + fluid_pressure_RHS -= melting_rate * (1.0/fluid_density - 1.0/solid_density); + + // compression term + // The whole expression for the first term on the RHS would be + // (u_s \cdot g) (\phi \rho_f \kappa_f + (1 - \phi) \rho_s \kappa_s). + // However, we already have the term (u_s \cdot g) \rho_s \kappa_s in the + // assembly of the stokes system without melt. Because of that, we only + // need to have -\phi \rho_s \kappa_s here. + fluid_pressure_RHS += is_compressible + ? + (current_u * fluid_density_gradient) * porosity / fluid_density + - (current_u * gravity) * porosity * solid_density * solid_compressibility + + K_D * (fluid_density_gradient * gravity) + : + 0.0; + + return fluid_pressure_RHS; + } + } + template void - MeltEquations:: - local_assemble_stokes_system_melt (const typename DoFHandler::active_cell_iterator &/*cell*/, - const double pressure_scaling, - const bool rebuild_stokes_matrix, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const + MeltStokesSystem:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); @@ -247,6 +315,8 @@ namespace aspect Assert(dynamic_cast*>(&this->get_material_model()) != NULL, ExcMessage("Error: The current material model needs to be derived from MeltInterface to use melt transport.")); + const double pressure_scaling = scratch.pressure_scaling; + for (unsigned int i=0, i_stokes=0; i_stokes()); - const double eta_two_thirds = (rebuild_stokes_matrix + const double eta_two_thirds = (scratch.rebuild_stokes_matrix ? scratch.material_model_outputs.viscosities[q] * 2.0 / 3.0 : @@ -318,9 +388,8 @@ namespace aspect const double viscosity_c = melt_outputs->compaction_viscosities[q]; const Tensor<1,dim> density_gradient_f = melt_outputs->fluid_density_gradients[q]; const double density_f = melt_outputs->fluid_densities[q]; - const double p_f_RHS = compute_fluid_pressure_rhs(scratch, - scratch.material_model_inputs, - scratch.material_model_outputs, + const double p_f_RHS = compute_fluid_pressure_rhs(this, + scratch, q); const double bulk_density = (1.0 - porosity) * density_s + porosity * density_f; @@ -359,7 +428,7 @@ namespace aspect ) * JxW; - if (rebuild_stokes_matrix) + if (scratch.rebuild_stokes_matrix) for (unsigned int j=0; j void - MeltEquations:: - local_assemble_stokes_system_melt_boundary (const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - const double pressure_scaling, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const + MeltStokesSystemBoundary:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); @@ -417,11 +486,13 @@ namespace aspect const unsigned int p_f_component_index = introspection.variable("fluid pressure").first_component_index; const unsigned int p_c_component_index = introspection.variable("compaction pressure").first_component_index; + const typename DoFHandler::face_iterator face = scratch.cell->face(scratch.face_number); + MaterialModel::MeltOutputs *melt_outputs = scratch.face_material_model_outputs.template get_additional_output >(); std::vector grad_p_f(n_face_q_points); this->get_melt_handler().boundary_fluid_pressure->fluid_pressure_gradient( - cell->face(face_no)->boundary_id(), + face->boundary_id(), scratch.face_material_model_inputs, scratch.face_material_model_outputs, #if DEAL_II_VERSION_GTE(9,0,0) @@ -454,7 +525,7 @@ namespace aspect { // apply the fluid pressure boundary condition data.local_rhs(i_stokes) += (scratch.face_finite_element_values[ex_p_f].value(i, q) - * pressure_scaling * K_D * + * scratch.pressure_scaling * K_D * (density_f * (scratch.face_finite_element_values.normal_vector(q) * gravity) - grad_p_f[q]) @@ -466,15 +537,60 @@ namespace aspect } } + namespace + { + /** + * Compute the right-hand side for the porosity system. + * It includes the melting rate and a term dependent + * on the density and velocity. + * + * This function is implemented in + * source/simulator/melt.cc. + */ + template + double + compute_melting_RHS(const SimulatorAccess *simulator_access, + const internal::Assembly::Scratch::AdvectionSystem &scratch, + const unsigned int q_point) + { + Assert (scratch.material_model_outputs.densities[q_point] > 0, + ExcMessage ("The density needs to be a positive quantity " + "when melt transport is included in the simulation.")); + + const double melting_rate = scratch.material_model_outputs.reaction_terms[q_point][scratch.advection_field->compositional_variable]; + const double density = scratch.material_model_outputs.densities[q_point]; + const double current_phi = scratch.material_model_inputs.composition[q_point][scratch.advection_field->compositional_variable]; + const double divergence_u = scratch.current_velocity_divergences[q_point]; + const double compressibility = (simulator_access->get_material_model().is_compressible() + ? + scratch.material_model_outputs.compressibilities[q_point] + : + 0.0); + const Tensor<1,dim> current_u = scratch.current_velocity_values[q_point]; + const Tensor<1,dim> + gravity = simulator_access->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q_point)); + + double melt_transport_RHS = melting_rate / density + + divergence_u + compressibility * density * (current_u * gravity); + + if (current_phi < simulator_access->get_melt_handler().melt_transport_threshold + && melting_rate < simulator_access->get_melt_handler().melt_transport_threshold) + melt_transport_RHS = melting_rate / density; + + return melt_transport_RHS; + } + } + template void - MeltEquations:: - local_assemble_advection_system_melt (const typename Simulator::AdvectionField &advection_field, - const double artificial_viscosity, - internal::Assembly::Scratch::AdvectionSystem &scratch, - internal::Assembly::CopyData::AdvectionSystem &data) const + MeltAdvectionSystem:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { + internal::Assembly::Scratch::AdvectionSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::AdvectionSystem &data = dynamic_cast& > (data_base); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); @@ -485,9 +601,9 @@ namespace aspect const double time_step = this->get_timestep(); const double old_time_step = this->get_old_timestep(); - const unsigned int solution_component = advection_field.component_index(introspection); + const unsigned int solution_component = scratch.advection_field->component_index(introspection); - const FEValuesExtractors::Scalar solution_field = advection_field.scalar_extractor(introspection); + const FEValuesExtractors::Scalar solution_field = scratch.advection_field->scalar_extractor(introspection); MaterialModel::MeltOutputs *melt_outputs = scratch.material_model_outputs.template get_additional_output >(); @@ -522,7 +638,7 @@ namespace aspect const double bulk_density = (1.0 - porosity) * scratch.material_model_outputs.densities[q] + porosity * melt_outputs->fluid_densities[q]; const double density_c_P = - ((advection_field.is_temperature()) + ((scratch.advection_field->is_temperature()) ? bulk_density * scratch.material_model_outputs.specific_heat[q] @@ -534,13 +650,13 @@ namespace aspect "non-negative quantity.")); const double conductivity = - (advection_field.is_temperature() + (scratch.advection_field->is_temperature() ? scratch.material_model_outputs.thermal_conductivities[q] : 0.0); const double latent_heat_LHS = - ((advection_field.is_temperature()) + ((scratch.advection_field->is_temperature()) ? scratch.heating_model_outputs.lhs_latent_heat_terms[q] : @@ -550,24 +666,25 @@ namespace aspect "to the left hand side needs to be a non-negative quantity.")); const double gamma = - ((advection_field.is_temperature()) + ((scratch.advection_field->is_temperature()) ? scratch.heating_model_outputs.heating_source_terms[q] : 0.0); const double reaction_term = - ((advection_field.is_temperature() || this->get_melt_handler().is_porosity(advection_field)) + ((scratch.advection_field->is_temperature() || this->get_melt_handler().is_porosity(*scratch.advection_field)) ? 0.0 : - scratch.material_model_outputs.reaction_terms[q][advection_field.compositional_variable]); + scratch.material_model_outputs.reaction_terms[q][scratch.advection_field->compositional_variable]); - const double melt_transport_RHS = compute_melting_RHS (scratch, - scratch.material_model_inputs, - scratch.material_model_outputs, - advection_field, - q); + const double melt_transport_RHS = (this->get_melt_handler().is_porosity(*scratch.advection_field) ? + compute_melting_RHS (this, + scratch, + q) + : + 0.0); const double field_term_for_rhs = (use_bdf2_scheme ? @@ -588,7 +705,7 @@ namespace aspect current_u -= scratch.mesh_velocity_values[q]; const double melt_transport_LHS = - (this->get_melt_handler().is_porosity(advection_field) + (this->get_melt_handler().is_porosity(*scratch.advection_field) ? scratch.current_velocity_divergences[q] + (this->get_material_model().is_compressible() @@ -610,7 +727,7 @@ namespace aspect double density_c_P_solid = density_c_P; double density_c_P_melt = 0.0; - if (advection_field.is_temperature() && porosity >= this->get_melt_handler().melt_transport_threshold + if (scratch.advection_field->is_temperature() && porosity >= this->get_melt_handler().melt_transport_threshold && this->get_melt_handler().heat_advection_by_melt) { density_c_P_solid = (1.0 - porosity) * scratch.material_model_outputs.densities[q] * scratch.material_model_outputs.specific_heat[q]; @@ -635,7 +752,7 @@ namespace aspect { data.local_matrix(i,j) += ( - (time_step * (conductivity + artificial_viscosity) + (time_step * (conductivity + scratch.artificial_viscosity) * (scratch.grad_phi_field[i] * scratch.grad_phi_field[j])) + ((time_step * (scratch.phi_field[i] * (current_u * scratch.grad_phi_field[j]))) + (factor * scratch.phi_field[i] * scratch.phi_field[j])) * @@ -653,10 +770,11 @@ namespace aspect template std::vector - MeltEquations:: - compute_advection_system_residual_melt(const typename Simulator::AdvectionField &advection_field, - internal::Assembly::Scratch::AdvectionSystem &scratch) const + MeltAdvectionSystemResidual:: + execute(internal::Assembly::Scratch::ScratchBase &scratch_base) const { + 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); @@ -676,7 +794,7 @@ namespace aspect const double u_grad_field = u * (scratch.old_field_grads[q] + scratch.old_old_field_grads[q]) / 2; - if (advection_field.is_temperature()) + if (scratch.advection_field->is_temperature()) { const double density = scratch.material_model_outputs.densities[q]; const double conductivity = scratch.material_model_outputs.thermal_conductivities[q]; @@ -695,22 +813,24 @@ namespace aspect const double field = ((scratch.old_field_values)[q] + (scratch.old_old_field_values)[q]) / 2; const double dreaction_term_dt = - (this->get_old_timestep() == 0 || (this->get_melt_handler().is_porosity(advection_field) + (this->get_old_timestep() == 0 || (this->get_melt_handler().is_porosity(*scratch.advection_field) && this->include_melt_transport())) ? 0.0 : - (scratch.material_model_outputs.reaction_terms[q][advection_field.compositional_variable] + (scratch.material_model_outputs.reaction_terms[q][scratch.advection_field->compositional_variable] / this->get_old_timestep()); - const double melt_transport_RHS = compute_melting_RHS (scratch, - scratch.material_model_inputs, - scratch.material_model_outputs, - advection_field, - q); + + const double melt_transport_RHS = (this->get_melt_handler().is_porosity(*scratch.advection_field) ? + compute_melting_RHS (this, + scratch, + q) + : + 0.0); const double melt_transport_LHS = - (this->get_melt_handler().is_porosity(advection_field) + (this->get_melt_handler().is_porosity(*scratch.advection_field) ? scratch.current_velocity_divergences[q] + (this->get_material_model().is_compressible() @@ -735,187 +855,127 @@ namespace aspect template - double - MeltEquations:: - compute_fluid_pressure_rhs(const internal::Assembly::Scratch::StokesSystem &scratch, - MaterialModel::MaterialModelInputs &material_model_inputs, - MaterialModel::MaterialModelOutputs &material_model_outputs, - const unsigned int q_point) const + void + MeltAdvectionSystemResidual::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const { - if (!this->get_parameters().include_melt_transport) - return 0.0; + MeltHandler::create_material_model_outputs(outputs); - const Introspection &introspection = this->introspection(); - MaterialModel::MeltOutputs *melt_out = material_model_outputs.template get_additional_output >(); - - Assert(melt_out != NULL, ExcInternalError()); - - const unsigned int porosity_index = introspection.compositional_index_for_name("porosity"); - const unsigned int is_compressible = this->get_material_model().is_compressible(); - - const double melting_rate = material_model_outputs.reaction_terms[q_point][porosity_index]; - const double solid_density = material_model_outputs.densities[q_point]; - const double fluid_density = melt_out->fluid_densities[q_point]; - const double solid_compressibility = material_model_outputs.compressibilities[q_point]; - const Tensor<1,dim> fluid_density_gradient = melt_out->fluid_density_gradients[q_point]; - const Tensor<1,dim> current_u = scratch.velocity_values[q_point]; - const double porosity = std::max(material_model_inputs.composition[q_point][porosity_index],0.0); - const double K_D = (porosity > this->get_melt_handler().melt_transport_threshold - ? - melt_out->permeabilities[q_point] / melt_out->fluid_viscosities[q_point] - : - 0.0); - - const Tensor<1,dim> - gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q_point)); - - double fluid_pressure_RHS = 0.0; - - // melting term - fluid_pressure_RHS -= melting_rate * (1.0/fluid_density - 1.0/solid_density); - - // compression term - // The whole expression for the first term on the RHS would be - // (u_s \cdot g) (\phi \rho_f \kappa_f + (1 - \phi) \rho_s \kappa_s). - // However, we already have the term (u_s \cdot g) \rho_s \kappa_s in the - // assembly of the stokes system without melt. Because of that, we only - // need to have -\phi \rho_s \kappa_s here. - fluid_pressure_RHS += is_compressible - ? - (current_u * fluid_density_gradient) * porosity / fluid_density - - (current_u * gravity) * porosity * solid_density * solid_compressibility - + K_D * (fluid_density_gradient * gravity) - : - 0.0; + const unsigned int n_points = outputs.viscosities.size(); - return fluid_pressure_RHS; + 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 - double - MeltEquations:: - compute_melting_RHS(const internal::Assembly::Scratch::AdvectionSystem &scratch, - const typename MaterialModel::Interface::MaterialModelInputs &material_model_inputs, - const typename MaterialModel::Interface::MaterialModelOutputs &material_model_outputs, - const typename Simulator::AdvectionField &advection_field, - const unsigned int q_point) const + void + MeltPressureRHSCompatibilityModification:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { - if ((!this->get_melt_handler().is_porosity(advection_field)) || (!this->include_melt_transport())) - return 0.0; - - Assert (material_model_outputs.densities[q_point] > 0, - ExcMessage ("The density needs to be a positive quantity " - "when melt transport is included in the simulation.")); - - const double melting_rate = material_model_outputs.reaction_terms[q_point][advection_field.compositional_variable]; - const double density = material_model_outputs.densities[q_point]; - const double current_phi = material_model_inputs.composition[q_point][advection_field.compositional_variable]; - const double divergence_u = scratch.current_velocity_divergences[q_point]; - const double compressibility = (this->get_material_model().is_compressible() - ? - material_model_outputs.compressibilities[q_point] - : - 0.0); - const Tensor<1,dim> current_u = scratch.current_velocity_values[q_point]; - const Tensor<1,dim> - gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q_point)); - - double melt_transport_RHS = melting_rate / density - + divergence_u + compressibility * density * (current_u * gravity); - - if (current_phi < this->get_melt_handler().melt_transport_threshold - && melting_rate < this->get_melt_handler().melt_transport_threshold) - melt_transport_RHS = melting_rate / density; - - return melt_transport_RHS; - } + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); + + const Introspection &introspection = this->introspection(); + const FiniteElement &fe = scratch.finite_element_values.get_fe(); + + const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); + const unsigned int p_f_component_index = introspection.variable("fluid pressure").first_component_index; + const unsigned int p_c_component_index = introspection.variable("compaction pressure").first_component_index; + const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; - namespace OtherTerms + const FEValuesExtractors::Scalar ex_p_f = introspection.variable("fluid pressure").extractor_scalar(); + + for (unsigned int q=0; q + void + MeltBoundaryTraction:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const { - template - void - pressure_rhs_compatibility_modification_melt (const SimulatorAccess &simulator_access, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) - { - const Introspection &introspection = simulator_access.introspection(); - const FiniteElement &fe = scratch.finite_element_values.get_fe(); + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast& > (data_base); - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int p_f_component_index = introspection.variable("fluid pressure").first_component_index; - const unsigned int p_c_component_index = introspection.variable("compaction pressure").first_component_index; + const Introspection &introspection = this->introspection(); + const FiniteElement &fe = this->get_fe(); - const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + // see if any of the faces are traction boundaries for which + // we need to assemble force terms for the right hand side + const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); + const unsigned int p_f_component_index = introspection.variable("fluid pressure").first_component_index; + const unsigned int p_c_component_index = introspection.variable("compaction pressure").first_component_index; - const FEValuesExtractors::Scalar ex_p_f = introspection.variable("fluid pressure").extractor_scalar(); + const typename DoFHandler::active_cell_iterator cell (&this->get_triangulation(), + scratch.face_finite_element_values.get_cell()->level(), + scratch.face_finite_element_values.get_cell()->index(), + &this->get_dof_handler()); - for (unsigned int q=0; q::faces_per_cell; ++face_no) + if (scratch.face_finite_element_values.get_face_index() == cell->face_index(face_no)) + break; + Assert(face_no != numbers::invalid_unsigned_int,ExcInternalError()); - if (is_velocity_or_pressures(introspection,p_c_component_index,p_f_component_index,component_index_i)) - { - scratch.phi_p[i_stokes] = scratch.finite_element_values[ex_p_f].value (i, q); - data.local_pressure_shape_function_integrals(i_stokes) += scratch.phi_p[i_stokes] * scratch.finite_element_values.JxW(q); - ++i_stokes; - } - ++i; - } - } + const typename DoFHandler::face_iterator face = cell->face(face_no); - template - void - boundary_traction_melt (const SimulatorAccess &simulator_access, - const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) - { - const Introspection &introspection = simulator_access.introspection(); - const FiniteElement &fe = scratch.finite_element_values.get_fe(); - - // see if any of the faces are traction boundaries for which - // we need to assemble force terms for the right hand side - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int p_f_component_index = introspection.variable("fluid pressure").first_component_index; - const unsigned int p_c_component_index = introspection.variable("compaction pressure").first_component_index; - - if (simulator_access.get_boundary_traction() - .find (cell->face(face_no)->boundary_id()) - != - simulator_access.get_boundary_traction().end()) - { - scratch.face_finite_element_values.reinit (cell, face_no); + if (this->get_boundary_traction() + .find (face->boundary_id()) + != + this->get_boundary_traction().end()) + { + scratch.face_finite_element_values.reinit (cell, face_no); - for (unsigned int q=0; q traction - = simulator_access.get_boundary_traction().find( - cell->face(face_no)->boundary_id() - )->second - ->boundary_traction (cell->face(face_no)->boundary_id(), - scratch.face_finite_element_values.quadrature_point(q), - scratch.face_finite_element_values.normal_vector(q)); - - for (unsigned int i=0, i_stokes=0; i_stokes traction + = this->get_boundary_traction().find( + face->boundary_id() + )->second + ->boundary_traction (face->boundary_id(), + scratch.face_finite_element_values.quadrature_point(q), + scratch.face_finite_element_values.normal_vector(q)); + + for (unsigned int i=0, i_stokes=0; i_stokes; \ - \ - \ - namespace OtherTerms \ - { \ - template void pressure_rhs_compatibility_modification_melt (const SimulatorAccess &simulator_access, \ - internal::Assembly::Scratch::StokesSystem &scratch, \ - internal::Assembly::CopyData::StokesSystem &data); \ - template void boundary_traction_melt (const SimulatorAccess &simulator_access, \ - const typename DoFHandler::active_cell_iterator &cell, \ - const unsigned int face_no, \ - internal::Assembly::Scratch::StokesSystem &scratch, \ - internal::Assembly::CopyData::StokesSystem &data); \ - } \ + template class MeltStokesPreconditioner; \ + template class MeltStokesSystem; \ + template class MeltStokesSystemBoundary; \ + template class MeltAdvectionSystem; \ + template class MeltAdvectionSystemResidual; \ + template class MeltPressureRHSCompatibilityModification; \ + template class MeltBoundaryTraction; \ } \ ASPECT_INSTANTIATE(INSTANTIATE) diff --git a/source/simulator/simulator_access.cc b/source/simulator/simulator_access.cc index 988dc46ad8a..b8a9b7f7786 100644 --- a/source/simulator/simulator_access.cc +++ b/source/simulator/simulator_access.cc @@ -370,13 +370,6 @@ namespace aspect } - template - void - SimulatorAccess::create_additional_material_model_outputs (MaterialModel::MaterialModelOutputs &out) const - { - simulator->create_additional_material_model_outputs(out); - } - template const std::map > > & @@ -566,6 +559,15 @@ namespace aspect return *(simulator->newton_handler); } + template + const FreeSurfaceHandler & + SimulatorAccess::get_free_surface_handler () const + { + Assert (simulator->free_surface.get() != 0, + ExcMessage("You can not call this function if the free surface is not enabled.")); + return *(simulator->free_surface); + } + template void SimulatorAccess::get_composition_values_at_q_point (const std::vector > &composition_values, diff --git a/tests/additional_outputs_02.cc b/tests/additional_outputs_02.cc index 9f1af325592..a3c5d85e2eb 100644 --- a/tests/additional_outputs_02.cc +++ b/tests/additional_outputs_02.cc @@ -2,7 +2,7 @@ #include #include #include -#include +#include #include @@ -67,7 +67,7 @@ namespace aspect template class TestAssembler : - public aspect::internal::Assembly::Assemblers::AssemblerBase + public aspect::Assemblers::Interface { public: @@ -83,9 +83,11 @@ namespace aspect } - virtual void assemble_stokes(internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem ©) + virtual void execute(internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &/*data_base*/) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + MaterialModel::AdditionalOutputs1 *additional = scratch.material_model_outputs.template get_additional_output >(); @@ -102,8 +104,7 @@ namespace aspect template void set_assemblers1(const SimulatorAccess &, - internal::Assembly::AssemblerLists &assemblers, - std::vector > > &assembler_objects) + Assemblers::Manager &assemblers) { std::cout << "* set_assemblers()" << std::endl; @@ -113,19 +114,7 @@ namespace aspect quiet = false; TestAssembler *test_assembler = new TestAssembler(); - assembler_objects.push_back(std_cxx11::shared_ptr >(test_assembler)); - - assemblers.local_assemble_stokes_system - .connect (std_cxx11::bind(&TestAssembler::assemble_stokes, - std_cxx11::ref (*test_assembler), - // discard cell, - // discard pressure_scaling, - // discard bool, - - std_cxx11::_4, - std_cxx11::_5)); - - + assemblers.stokes_system.push_back(std_cxx11::unique_ptr >(test_assembler)); } } diff --git a/tests/additional_outputs_02/screen-output b/tests/additional_outputs_02/screen-output index 53869150db7..24ae954330e 100644 --- a/tests/additional_outputs_02/screen-output +++ b/tests/additional_outputs_02/screen-output @@ -8,24 +8,14 @@ Number of active cells: 4 (on 2 levels) Number of degrees of freedom: 84 (50+9+25) *** Timestep 0: t=0 seconds -* create_additional_material_model_outputs() called - creating additional output! -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called - creating additional output! -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! Solving temperature system... 0 iterations. * create_additional_material_model_outputs() called creating additional output! @@ -44,15 +34,10 @@ Number of degrees of freedom: 84 (50+9+25) * evaluate called with additional outputs! * local_assemble_stokes call, have additional? 1 value = 42 - Rebuilding Stokes preconditioner...* create_additional_material_model_outputs() called - creating additional output! -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! + Rebuilding Stokes preconditioner...* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! Solving Stokes system... 11+0 iterations. @@ -70,4 +55,4 @@ Termination requested by criterion: end time +---------------------------------+-----------+------------+------------+ +---------------------------------+-----------+------------+------------+ -called without: 4 with: 16 +called without: 16 with: 4 diff --git a/tests/additional_outputs_03.cc b/tests/additional_outputs_03.cc index 610e84f0ab0..629333341ab 100644 --- a/tests/additional_outputs_03.cc +++ b/tests/additional_outputs_03.cc @@ -2,7 +2,7 @@ #include #include #include -#include +#include #include @@ -81,7 +81,7 @@ namespace aspect template class TestAssembler : - public aspect::internal::Assembly::Assemblers::AssemblerBase + public aspect::Assemblers::Interface { public: @@ -97,9 +97,11 @@ namespace aspect } - virtual void assemble_stokes(internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem ©) + virtual void execute(internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &/*data_base*/) const { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast& > (scratch_base); + MaterialModel::AdditionalOutputs1 *additional = scratch.material_model_outputs.template get_additional_output >(); @@ -117,8 +119,7 @@ namespace aspect template void set_assemblers1(const SimulatorAccess &, - internal::Assembly::AssemblerLists &assemblers, - std::vector > > &assembler_objects) + Assemblers::Manager &assemblers) { std::cout << "* set_assemblers()" << std::endl; std::cout << "called without: " << counter_without << " with: " << counter_with << std::endl; @@ -127,19 +128,7 @@ namespace aspect quiet = false; TestAssembler *test_assembler = new TestAssembler(); - assembler_objects.push_back(std_cxx11::shared_ptr >(test_assembler)); - - assemblers.local_assemble_stokes_system - .connect (std_cxx11::bind(&TestAssembler::assemble_stokes, - std_cxx11::ref (*test_assembler), - // discard cell, - // discard pressure_scaling, - // discard bool, - - std_cxx11::_4, - std_cxx11::_5)); - - + assemblers.stokes_system.push_back(std_cxx11::unique_ptr >(test_assembler)); } } diff --git a/tests/additional_outputs_03/screen-output b/tests/additional_outputs_03/screen-output index 487aa69e9c4..5a745fb0bf5 100644 --- a/tests/additional_outputs_03/screen-output +++ b/tests/additional_outputs_03/screen-output @@ -8,32 +8,14 @@ Number of active cells: 4 (on 2 levels) Number of degrees of freedom: 84 (50+9+25) *** Timestep 0: t=0 seconds -* create_additional_material_model_outputs() called - creating additional output! -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called - creating additional output! -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -averaging! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! Solving temperature system... 0 iterations. * create_additional_material_model_outputs() called creating additional output! @@ -56,19 +38,10 @@ averaging! averaging! * local_assemble_stokes call, have additional? 1 value = 21 21 - Rebuilding Stokes preconditioner...* create_additional_material_model_outputs() called - creating additional output! -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -averaging! -* create_additional_material_model_outputs() called -* evaluate called with additional outputs! -averaging! + Rebuilding Stokes preconditioner...* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! +* evaluate called without additional outputs! Solving Stokes system... 9+0 iterations. @@ -86,4 +59,4 @@ Termination requested by criterion: end time +---------------------------------+-----------+------------+------------+ +---------------------------------+-----------+------------+------------+ -called without: 4 with: 16 +called without: 16 with: 4 diff --git a/tests/inner_core.cc b/tests/inner_core.cc index 62b0ca713c4..d0fdca2e22f 100644 --- a/tests/inner_core.cc +++ b/tests/inner_core.cc @@ -1,426 +1 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -namespace aspect -{ - namespace MaterialModel - { - using namespace dealii; - - template - class InnerCore : public MaterialModel::Simple - { - public: - virtual void evaluate(const MaterialModel::MaterialModelInputs &in, - MaterialModel::MaterialModelOutputs &out) const; - - /** - * A function that is called at the beginning of each time step to - * indicate what the model time is for which the boundary values will - * next be evaluated. For the current class, the function passes to - * the parsed function what the current time is. - */ - virtual - void - update (); - - /** - * Declare the parameters this class takes through input files. The - * default implementation of this function does not describe any - * parameters. Consequently, derived classes do not have to overload - * this function if they do not take any runtime parameters. - */ - static - void - declare_parameters (ParameterHandler &prm); - - /** - * Read the parameters this class declares from the parameter file. - * The default implementation of this function does not read any - * parameters. Consequently, derived classes do not have to overload - * this function if they do not take any runtime parameters. - */ - virtual - void - parse_parameters (ParameterHandler &prm); - - /** - * A function object representing resistance to phase change at the - * inner core boundary as a function the position (and, optionally, - * the model time). - */ - Functions::ParsedFunction resistance_to_phase_change; - }; - - } -} - -namespace aspect -{ - namespace MaterialModel - { - - template - void - InnerCore:: - evaluate(const MaterialModel::MaterialModelInputs &in, - MaterialModel::MaterialModelOutputs &out) const - { - // First, we use the material descriptions of the 'simple' material model to fill all of the material - // model outputs. Below, we will then overwrite selected properties (the specific heat) to make the - // product of density and specific heat a constant. - Simple::evaluate(in, out); - - // We want the right-hand side of the momentum equation to be (- Ra T gravity) and - // density * cp to be 1 - for (unsigned int q=0; q < in.position.size(); ++q) - { - out.densities[q] = - out.thermal_expansion_coefficients[q] * in.temperature[q]; - if (std::abs(out.densities[q]) > 0.0) - out.specific_heat[q] /= out.densities[q]; - } - } - - - template - void - InnerCore::update() - { - // we get time passed as seconds (always) but may want - // to reinterpret it in years - if (this->convert_output_to_years()) - resistance_to_phase_change.set_time (this->get_time() / year_in_seconds); - else - resistance_to_phase_change.set_time (this->get_time()); - } - - - template - void - InnerCore::declare_parameters (ParameterHandler &prm) - { - Simple::declare_parameters (prm); - - prm.enter_subsection("Material model"); - { - prm.enter_subsection("Inner core"); - { - prm.enter_subsection("Phase change resistance function"); - { - Functions::ParsedFunction::declare_parameters (prm, 1); - } - prm.leave_subsection(); - } - prm.leave_subsection(); - } - prm.leave_subsection(); - } - - - template - void - InnerCore::parse_parameters (ParameterHandler &prm) - { - Simple::parse_parameters (prm); - - prm.enter_subsection("Material model"); - { - prm.enter_subsection("Inner core"); - { - prm.enter_subsection("Phase change resistance function"); - try - { - resistance_to_phase_change.parse_parameters (prm); - } - catch (...) - { - std::cerr << "ERROR: FunctionParser failed to parse\n" - << "\t'Phase boundary model.Function'\n" - << "with expression\n" - << "\t'" << prm.get("Function expression") << "'" - << "More information about the cause of the parse error \n" - << "is shown below.\n"; - throw; - } - prm.leave_subsection(); - } - prm.leave_subsection(); - } - prm.leave_subsection(); - } - - } -} - - - -namespace aspect -{ - namespace HeatingModel - { - using namespace dealii; - - /** - * A class that implements a constant radiogenic heating rate. - * - * @ingroup HeatingModels - */ - template - class ConstantCoreHeating : public Interface - { - public: - /** - * Return the heating terms. For the current class, this - * function obviously simply returns a constant value. - */ - virtual - void - evaluate (const MaterialModel::MaterialModelInputs &material_model_inputs, - const MaterialModel::MaterialModelOutputs &material_model_outputs, - HeatingModel::HeatingModelOutputs &heating_model_outputs) const; - - /** - * @name Functions used in dealing with run-time parameters - * @{ - */ - /** - * Declare the parameters this class takes through input files. - */ - static - void - declare_parameters (ParameterHandler &prm); - - /** - * Read the parameters this class declares from the parameter file. - */ - virtual - void - parse_parameters (ParameterHandler &prm); - /** - * @} - */ - - private: - double radiogenic_heating_rate; - }; - } -} - -namespace aspect -{ - namespace HeatingModel - { - template - void - ConstantCoreHeating:: - evaluate (const MaterialModel::MaterialModelInputs &/*material_model_inputs*/, - const MaterialModel::MaterialModelOutputs &/*material_model_outputs*/, - HeatingModel::HeatingModelOutputs &heating_model_outputs) const - { - for (unsigned int q=0; q - void - ConstantCoreHeating::declare_parameters (ParameterHandler &prm) - { - prm.enter_subsection("Heating model"); - { - prm.enter_subsection("Constant core heating"); - { - prm.declare_entry ("Radiogenic heating rate", "0e0", - Patterns::Double (0), - "The specific rate of heating due to radioactive decay (or other bulk sources " - "you may want to describe). This parameter corresponds to the variable " - "$H$ in the temperature equation stated in the manual, and the heating " - "term is $\rho H$. Units: W/kg."); - } - prm.leave_subsection(); - } - prm.leave_subsection(); - } - - - - template - void - ConstantCoreHeating::parse_parameters (ParameterHandler &prm) - { - prm.enter_subsection("Heating model"); - { - prm.enter_subsection("Constant core heating"); - { - radiogenic_heating_rate = prm.get_double ("Radiogenic heating rate"); - } - prm.leave_subsection(); - } - prm.leave_subsection(); - } - } -} - - -namespace aspect -{ - /** - * A new assembler class that implements boundary conditions for the - * normal stress and the normal velocity that take into account the - * rate of phase change (melting/freezing) at the inner-outer core - * boundary. The model is based on Deguen, Alboussiere, and Cardin - * (2013), Thermal convection in Earth’s inner core with phase change - * at its boundary. GJI, 194, 1310-133. - * - * The mechanical boundary conditions for the inner core are - * tangential stress-free and continuity of the normal stress at the - * inner-outer core boundary. For the non-dimensional equations, that - * means that we can define a 'phase change number' $\mathcal{P}$ so - * that the normal stress at the boundary is $-\mathcal{P} u_r$ with - * the radial velocity $u_r$. This number characterizes the resistance - * to phase change at the boundary, with $\mathcal{P}\rightarrow\infty$ - * corresponding to infinitely slow melting/freezing (free slip - * boundary), and $\mathcal{P}\rightarrow0$ corresponding to - * instantaneous melting/freezing (zero normal stress, open boundary). - * - * In the weak form, this results in boundary conditions of the form - * of a surface integral: - * $$\int_S \mathcal{P} (\mathbf u \cdot \mathbf n) (\mathbf v \cdot \mathbf n) dS,$$, - * with the normal vector $\mathbf n$. - * - * The function value of $\mathcal{P}$ is taken from the inner core - * material model. - */ - template - class PhaseBoundaryAssembler : - public aspect::internal::Assembly::Assemblers::AssemblerBase, - public SimulatorAccess - { - - public: - - virtual - void - phase_change_boundary_conditions (const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face_no, - internal::Assembly::Scratch::StokesSystem &scratch, - internal::Assembly::CopyData::StokesSystem &data) const - { - const Introspection &introspection = this->introspection(); - const FiniteElement &fe = this->get_fe(); - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int n_q_points = scratch.face_finite_element_values.n_quadrature_points; - - // assemble force terms for the matrix for all boundary faces - if (cell->face(face_no)->at_boundary()) - { - scratch.face_finite_element_values.reinit (cell, face_no); - - for (unsigned int q=0; q&> - (this->get_material_model()).resistance_to_phase_change - .value(scratch.material_model_inputs.position[q]); - - for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) - { - if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) - { - scratch.phi_u[i_stokes] = scratch.face_finite_element_values[introspection - .extractors.velocities].value(i, q); - ++i_stokes; - } - ++i; - } - - const Tensor<1,dim> normal_vector = scratch.face_finite_element_values.normal_vector(q); - const double JxW = scratch.face_finite_element_values.JxW(q); - - // boundary term: P*u*n*v*n*JxW(q) - for (unsigned int i=0; i - void set_assemblers_phase_boundary(const SimulatorAccess &simulator_access, - internal::Assembly::AssemblerLists &assemblers, - std::vector > > &assembler_objects) - { - AssertThrow (dynamic_cast*> - (&simulator_access.get_material_model()) != 0, - ExcMessage ("The phase boundary assembler can only be used with the " - "material model 'inner core material'!")); - - std_cxx11::shared_ptr > phase_boundary_assembler - (new PhaseBoundaryAssembler()); - assembler_objects.push_back (phase_boundary_assembler); - - // add the terms for phase change boundary conditions - assemblers.local_assemble_stokes_system_on_boundary_face - .connect (std_cxx11::bind(&PhaseBoundaryAssembler::phase_change_boundary_conditions, - std_cxx11::cref (*phase_boundary_assembler), - std_cxx11::_1, - std_cxx11::_2, - // discard pressure_scaling, - // discard rebuild_stokes_matrix, - std_cxx11::_5, - std_cxx11::_6)); - - - } -} - -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - signals.set_assemblers.connect (&aspect::set_assemblers_phase_boundary); -} - -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>) - -// explicit instantiations -namespace aspect -{ - namespace MaterialModel - { - ASPECT_REGISTER_MATERIAL_MODEL(InnerCore, - "inner core material", - "A simple material model that is like the " - "'Simple' model, but has a constant $\rho c_p$, " - "and implements a function that characterizes the " - "resistance to melting/freezing at the inner core " - "boundary.") - } - - namespace HeatingModel - { - ASPECT_REGISTER_HEATING_MODEL(ConstantCoreHeating, - "constant core heating", - "Implementation of a model in which the heating " - "rate is constant.") - } -} +#include "../cookbooks/inner_core_convection/inner_core_assembly.cc" diff --git a/tests/melt_introspection.cc b/tests/melt_introspection.cc index f627b7200aa..7d6694ac547 100644 --- a/tests/melt_introspection.cc +++ b/tests/melt_introspection.cc @@ -3,7 +3,7 @@ #include #include #include -#include +#include #include #include diff --git a/tests/signal_fem.cc b/tests/signal_fem.cc index 8b546db9871..146c4bb1776 100644 --- a/tests/signal_fem.cc +++ b/tests/signal_fem.cc @@ -2,7 +2,7 @@ #include #include #include -#include +#include #include #include