diff --git a/applications/advection/CMakeLists.txt b/applications/advection/CMakeLists.txt new file mode 100644 index 000000000..31e2527fe --- /dev/null +++ b/applications/advection/CMakeLists.txt @@ -0,0 +1,68 @@ +## +# CMake script for the PRISMS-PF applications: +## + +cmake_minimum_required(VERSION 3.1.0) +project(myapp) + +# Find deal.II installation +find_package(deal.II 9.2.0 QUIET REQUIRED + HINTS ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}) + +# Check to make sure deal.II is configured with p4est +if(NOT ${DEAL_II_WITH_P4EST}) + message(FATAL_ERROR "deal.II must be installed with p4est.") +endif() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() + +# Set default build type +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Debug CACHE STRING "Build type" FORCE) +endif() + +# Set up the debug, release, and run targets +add_custom_target(debug + COMMAND +env ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR} + COMMAND +env ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all + COMMENT "Switch CMAKE_BUILD_TYPE to Debug" + ) +add_custom_target(release + COMMAND +env ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR} + COMMAND +env ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all + COMMENT "Switch CMAKE_BUILD_TYPE to Release" + ) +add_custom_target(run COMMAND main + COMMENT "Run with ${CMAKE_BUILD_TYPE} configuration" + ) + +# Add postprocess.cc and nucleation.cc if they exist +if(EXISTS "postprocess.cc") + add_definitions(-DPOSTPROCESS_FILE_EXISTS) +endif() +if(EXISTS "nucleation.cc") + add_definitions(-DNUCLEATION_FILE_EXISTS) +endif() + +# Set location of files +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) +include_directories(${CMAKE_SOURCE_DIR}) + +# SEt the location of the main.cc file +set(MAIN "${CMAKE_SOURCE_DIR}/../main.cc") + +# Add main.cc for executable target +add_executable(main ${MAIN}) + +# deal.II linker +DEAL_II_SETUP_TARGET(main) + +# Link libraries for the build type +if (${CMAKE_BUILD_TYPE} STREQUAL "Release") + target_link_libraries(main ${CMAKE_SOURCE_DIR}/../../libprisms_pf.a) +elseif(${CMAKE_BUILD_TYPE} STREQUAL "DebugRelease") + target_link_libraries(main ${CMAKE_SOURCE_DIR}/../../libprisms_pf.a) +else() + target_link_libraries(main ${CMAKE_SOURCE_DIR}/../../libprisms_pf_debug.a) +endif() diff --git a/applications/advection/ICs_and_BCs.cc b/applications/advection/ICs_and_BCs.cc new file mode 100644 index 000000000..cf986f563 --- /dev/null +++ b/applications/advection/ICs_and_BCs.cc @@ -0,0 +1,71 @@ +// =========================================================================== +// FUNCTION FOR INITIAL CONDITIONS +// =========================================================================== +#include + +template +void +customPDE::setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) +{ + // --------------------------------------------------------------------- + // ENTER THE INITIAL CONDITIONS HERE + // --------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the initial condition for each variable + // according to its variable index + + // Sphere level-set + double radius = 15.0; + dealii::Point center(20.0, 20.0); + + dealii::Functions::SignedDistance::Sphere sphere(center, radius); + + double distance = sphere.value(p); + + // tanh profile + double phi = 0.5 * (1.0 - std::tanh(distance / (std::sqrt(2) * W))); + + if (index == 0) + { + scalar_IC = phi; + } + + else + { + scalar_IC = 0.0; + } + + // --------------------------------------------------------------------- +} + +// =========================================================================== +// FUNCTION FOR NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS +// =========================================================================== + +template +void +customPDE::setNonUniformDirichletBCs( + [[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) +{ + // -------------------------------------------------------------------------- + // ENTER THE NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS HERE + // -------------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the boundary condition for each variable + // according to its variable index. This function can be left blank if there + // are no non-uniform Dirichlet boundary conditions. For BCs that change in + // time, you can access the current time through the variable "time". The + // boundary index can be accessed via the variable "direction", which starts + // at zero and uses the same order as the BC specification in parameters.in + // (i.e. left = 0, right = 1, bottom = 2, top = 3, front = 4, back = 5). + + // ------------------------------------------------------------------------- +} diff --git a/applications/advection/advection.md b/applications/advection/advection.md new file mode 100644 index 000000000..5f411c4ac --- /dev/null +++ b/applications/advection/advection.md @@ -0,0 +1,40 @@ +## PRISMS-PF: Advection Equation +The advection equation is given by: + +$$ +\begin{align} +\frac{\partial \eta}{\partial t} = -\textbf{u}\cdot\nabla\eta +\end{align} +$$ + +where $\eta$ is the scalar conserved quantity and $\textbf{u}$ is a constant vector field. + +## Time discretization + +Considering backward Euler time stepping, we have the time discretized kinetics equation: + +$$ +\begin{align} +\eta^{n+1} = \eta^{n} - \Delta t (\textbf{u}\cdot\nabla\eta^{n+1}) +\end{align} +$$ + +PRISMS-PF makes use of the Newton's method for linear and nonlinear iterations, so we can substitute $\eta^{n+1}=\Delta \eta^{n+1}+n_0^{n+1}$. + +$$ +\begin{align} +\Delta \eta^{n+1} + \Delta t (\textbf{u}\cdot\nabla(\Delta \eta^{n+1}))= \eta^{n} - \eta_0^{n+1} - \Delta t (\textbf{u}\cdot\nabla\eta_0^{n+1}) +\end{align} +$$ + +## Weak formulation + +In the weak formulation, considering an arbitrary variation $w$, the above equation can be expressed as a residual equation: + +$$ +\begin{align} +\int_{\Omega} w [\Delta \eta^{n+1} + \Delta t (\textbf{u}\cdot\nabla(\Delta \eta^{n+1}))] ~dV &= \int_{\Omega} w [\eta^{n} - \eta_0^{n+1} - \Delta t (\textbf{u}\cdot\nabla\eta_0^{n+1})] ~dV +\end{align} +$$ + + diff --git a/applications/advection/customPDE.h b/applications/advection/customPDE.h new file mode 100644 index 000000000..440eba387 --- /dev/null +++ b/applications/advection/customPDE.h @@ -0,0 +1,91 @@ +#include "matrixFreePDE.h" + +using namespace dealii; + +template +class customPDE : public MatrixFreePDE +{ +public: + // Constructor + customPDE(userInputParameters _userInputs) + : MatrixFreePDE(_userInputs) + , userInputs(_userInputs) {}; + + // Function to set the initial conditions (in ICs_and_BCs.h) + void + setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) override; + + // Function to set the non-uniform Dirichlet boundary conditions (in + // ICs_and_BCs.h) + void + setNonUniformDirichletBCs([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) override; + +private: +#include "typeDefs.h" + + const userInputParameters userInputs; + + // Function to set the RHS of the governing equations for explicit time + // dependent equations (in equations.cc) + void + explicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] Point> q_point_loc) const override; + + // Function to set the RHS of the governing equations for all other equations + // (in equations.cc) + void + nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] Point> q_point_loc) const override; + + // Function to set the LHS of the governing equations (in equations.cc) + void + equationLHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] Point> q_point_loc) const override; + +// Function to set postprocessing expressions (in postprocess.h) +#ifdef POSTPROCESS_FILE_EXISTS + void + postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc) + const override; +#endif + +// Function to set the nucleation probability (in nucleation.h) +#ifdef NUCLEATION_FILE_EXISTS + double + getNucleationProbability([[maybe_unused]] variableValueContainer variable_value, + [[maybe_unused]] double dV) const override; +#endif + + // ================================================================ + // Methods specific to this subclass + // ================================================================ + + // ================================================================ + // Model constants specific to this subclass + // ================================================================ + + double W = userInputs.get_model_constant_double("W"); + dealii::Tensor<1, dim> velocity = + userInputs.get_model_constant_rank_1_tensor("velocity"); + + // ================================================================ +}; diff --git a/applications/advection/equations.cc b/applications/advection/equations.cc new file mode 100644 index 000000000..826c094ce --- /dev/null +++ b/applications/advection/equations.cc @@ -0,0 +1,129 @@ +// ================================================================================= +// Set the attributes of the primary field variables +// ================================================================================= +// This function sets attributes for each variable/equation in the app. The +// attributes are set via standardized function calls. The first parameter for +// each function call is the variable index (starting at zero). The first set of +// variable/equation attributes are the variable name (any string), the variable +// type (SCALAR/VECTOR), and the equation type (EXPLICIT_TIME_DEPENDENT/ +// TIME_INDEPENDENT/AUXILIARY). The next set of attributes describe the +// dependencies for the governing equation on the values and derivatives of the +// other variables for the value term and gradient term of the RHS and the LHS. +// The final pair of attributes determine whether a variable represents a field +// that can nucleate and whether the value of the field is needed for nucleation +// rate calculations. + +void +variableAttributeLoader::loadVariableAttributes() +{ + // Variable 0 + set_variable_name(0, "n"); + set_variable_type(0, SCALAR); + set_variable_equation_type(0, IMPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(0, "n, old(n), grad(n)"); + set_dependencies_value_term_LHS(0, "change(n), grad(change(n))"); +} + +// ============================================================================================= +// explicitEquationRHS (needed only if one or more equation is explict time +// dependent) +// ============================================================================================= +// This function calculates the right-hand-side of the explicit time-dependent +// equations for each variable. It takes "variable_list" as an input, which is a +// list of the value and derivatives of each of the variables at a specific +// quadrature point. The (x,y,z) location of that quadrature point is given by +// "q_point_loc". The function outputs two terms to variable_list -- one +// proportional to the test function and one proportional to the gradient of the +// test function. The index for each variable in this list corresponds to the +// index given at the top of this file. + +template +void +customPDE::explicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] Point> q_point_loc) const +{} + +// ============================================================================================= +// nonExplicitEquationRHS (needed only if one or more equation is time +// independent or auxiliary) +// ============================================================================================= +// This function calculates the right-hand-side of all of the equations that are +// not explicit time-dependent equations. It takes "variable_list" as an input, +// which is a list of the value and derivatives of each of the variables at a +// specific quadrature point. The (x,y,z) location of that quadrature point is +// given by "q_point_loc". The function outputs two terms to variable_list -- +// one proportional to the test function and one proportional to the gradient of +// the test function. The index for each variable in this list corresponds to +// the index given at the top of this file. + +template +void +customPDE::nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] Point> q_point_loc) const +{ + // Getting necessary variables + scalarvalueType n = variable_list.get_scalar_value(0); + scalargradType nx = variable_list.get_scalar_gradient(0); + scalarvalueType n_old = variable_list.get_old_scalar_value(0); + + // Converting prescibed velocity to a vectorvalueType + vectorvalueType vel; + for (unsigned int i = 0; i < dim; i++) + { + vel[i] = constV(velocity[i]); + } + + // Norm of the local velocity + scalarvalueType u_l2norm = 1.0e-12 + vel.norm_square(); + + // Submission terms + scalarvalueType residual = (n_old - n - constV(userInputs.dtValue) * vel * nx); + scalarvalueType eq_n = residual; + + variable_list.set_scalar_value_term_RHS(0, eq_n); +} + +// ============================================================================================= +// equationLHS (needed only if at least one equation is time independent) +// ============================================================================================= +// This function calculates the left-hand-side of time-independent equations. It +// takes "variable_list" as an input, which is a list of the value and +// derivatives of each of the variables at a specific quadrature point. The +// (x,y,z) location of that quadrature point is given by "q_point_loc". The +// function outputs two terms to variable_list -- one proportional to the test +// function and one proportional to the gradient of the test function -- for the +// left-hand-side of the equation. The index for each variable in this list +// corresponds to the index given at the top of this file. If there are multiple +// elliptic equations, conditional statements should be sed to ensure that the +// correct residual is being submitted. The index of the field being solved can +// be accessed by "this->currentFieldIndex". + +template +void +customPDE::equationLHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] Point> q_point_loc) const +{ + // Getting necessary variables + scalarvalueType change_n = variable_list.get_change_in_scalar_value(0); + scalargradType change_nx = variable_list.get_change_in_scalar_gradient(0); + + // Converting prescibed velocity to a vectorvalueType + vectorvalueType vel; + for (unsigned int i = 0; i < dim; i++) + { + vel[i] = constV(velocity[i]); + } + + // Norm of the local velocity + scalarvalueType u_l2norm = 1.0e-12 + vel.norm_square(); + + // Submission terms + scalarvalueType residual = (change_n + constV(userInputs.dtValue) * vel * change_nx); + scalarvalueType eq_n = residual; + + variable_list.set_scalar_value_term_LHS(0, eq_n); +} diff --git a/applications/advection/parameters.prm b/applications/advection/parameters.prm new file mode 100644 index 000000000..983d05247 --- /dev/null +++ b/applications/advection/parameters.prm @@ -0,0 +1,134 @@ +# ================================================================================= +# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) +# ================================================================================= +set Number of dimensions = 2 + +# ================================================================================= +# Set the length of the domain in all three dimensions +# (Domain size Z ignored in 2D) +# ================================================================================= +# Each axes spans from zero to the specified length +set Domain size X = 100 +set Domain size Y = 100 +set Domain size Z = 100 + +# ================================================================================= +# Set the element parameters +# ================================================================================= +# The number of elements in each direction is 2^(refineFactor) * subdivisions +# Subdivisions Z ignored in 2D +# For optimal performance, use refineFactor primarily to determine the element size +set Subdivisions X = 1 +set Subdivisions Y = 1 +set Subdivisions Z = 1 + +set Refine factor = 7 + +# Set the polynomial degree of the element (allowed values: 1, 2, or 3) +set Element degree = 2 + +# ================================================================================= +# Set the adaptive mesh refinement parameters +# ================================================================================= +# Set the flag determining if adaptive meshing is activated +set Mesh adaptivity = false + +# Set the maximum and minimum level of refinement +# When adaptive meshing is enabled, the refine factor set in the block above is +# only used to generate the first pass of the mesh as the initial conditions are +# applied. It should be set somewhere between the max and min levels below. +set Max refinement level = 7 +set Min refinement level = 4 + +# Set the number of time steps between remeshing operations +set Steps between remeshing operations = 100 + +# Set the criteria for adapting the mesh +subsection Refinement criterion: n + # Select whether the mesh is refined based on the variable value (VALUE), + # its gradient (GRADIENT), or both (VALUE_AND_GRADIENT) + set Criterion type = VALUE + # Set the lower and upper bounds for the value-based refinement window + set Value lower bound = 0.01 + set Value upper bound = 1.1 +end + +# ================================================================================= +# Set the time step parameters +# ================================================================================= +# The size of the time step +set Time step = 1.0e-4 + +# The simulation ends when either the number of time steps is reached or the +# simulation time is reached. +set Number of time steps = 100000 + +# ================================================================================= +# Set the output parameters +# ================================================================================= +# Type of spacing between outputs ("EQUAL_SPACING", "LOG_SPACING", "N_PER_DECADE", +# or "LIST") +set Output condition = EQUAL_SPACING + +# Number of times the program outputs the fields (total number for "EQUAL_SPACING" +# and "LOG_SPACING", number per decade for "N_PER_DECADE", ignored for "LIST") +set Number of outputs = 5 + +# Whether to print timing information every time the code outputs +set Print timing information with output = false + +# The number of time steps between updates being printed to the screen +set Skip print steps = 1000 + +# ================================================================================= +# Set the boundary conditions +# ================================================================================= +# Set the boundary condition for each variable, where each variable is given by +# its name, as defined in equations.cc. The four boundary condition +# types are NATURAL, DIRICHLET, NON_UNIFORM_DIRICHLET and PERIODIC. If all +# of the boundaries have the same boundary condition, only one boundary condition +# type needs to be given. If multiple boundary condition types are needed, give a +# comma-separated list of the types. The order is the miniumum of x, maximum of x, +# minimum of y, maximum of y, minimum of z, maximum of z (i.e left, right, bottom, +# top in 2D and left, right, bottom, top, front, back in 3D). The value of a +# Dirichlet BC is specfied in the following way -- DIRCHILET: val -- where 'val' is +# the desired value. If the boundary condition is NON_UNIFORM_DIRICHLET, the +# boundary condition should be specified in the appropriate function in 'ICs_and_BCs.h'. +# Example 1: All periodic BCs for variable 'c' +# set Boundary condition for variable c = PERIODIC +# Example 2: Zero-derivative BCs on the left and right, Dirichlet BCs with value +# 1.5 on the top and bottom for variable 'n' in 2D +# set Boundary condition for variable n = NATURAL, NATURAL, DIRICHLET: 1.5, DIRICHLET: 1.5 + +set Boundary condition for variable n = PERIODIC + +# ================================================================================= +# Set the model constants +# ================================================================================= +# Set the user-defined model constants, which must have a counter-part given in +# customPDE.h. These are most often used in the residual equations in equations.cc, +# but may also be used for initial conditions and nucleation calculations. The type +# options currently are DOUBLE, INT, BOOL, TENSOR, and [symmetry] ELASTIC CONSTANTS +# where [symmetry] is ISOTROPIC, TRANSVERSE, ORTHOTROPIC, or ANISOTROPIC. + +# The interfacial width +set Model constant W = 0.5, DOUBLE + +# The prescribed velocity +set Model constant velocity = (1.0, 1.0, 0.0), TENSOR + +# ================================================================================= +# Set the linear solver parameters +# ================================================================================= + +subsection Linear solver parameters: n + # Whether the tolerance value is compared to the residual (ABSOLUTE_RESIDUAL) + # or the change in the residual (RELATIVE_RESIDUAL_CHANGE) + set Tolerance type = RELATIVE_RESIDUAL_CHANGE + + # The tolerance for convergence (L2 norm) + set Tolerance value = 1e-3 + + # The maximum number of linear solver iterations per solve + set Maximum linear solver iterations = 1000 +end diff --git a/include/EquationDependencyParser.h b/include/EquationDependencyParser.h index 4a8dd1b03..4f784ed39 100644 --- a/include/EquationDependencyParser.h +++ b/include/EquationDependencyParser.h @@ -38,6 +38,8 @@ class EquationDependencyParser std::vector eval_flags_nonexplicit_RHS; std::vector eval_flags_nonexplicit_LHS; std::vector eval_flags_change_nonexplicit_LHS; + std::vector eval_flags_old_RHS; + std::vector eval_flags_old_LHS; std::vector eval_flags_residual_explicit_RHS; std::vector @@ -69,6 +71,7 @@ class EquationDependencyParser std::string &value_dependencies, std::string &gradient_dependencies, std::vector &evaluation_flags, + std::vector &old_flags, std::vector &residual_flags, bool &is_nonlinear); @@ -84,6 +87,7 @@ class EquationDependencyParser std::string &value_dependencies, std::string &gradient_dependencies, std::vector &evaluation_flags, + std::vector &old_flags, std::vector &change_flags, std::vector &residual_flags, bool &is_nonlinear); diff --git a/include/matrixFreePDE.h b/include/matrixFreePDE.h index b157da3ab..818f62ca4 100644 --- a/include/matrixFreePDE.h +++ b/include/matrixFreePDE.h @@ -226,6 +226,12 @@ class MatrixFreePDE : public Subscriptor /*Vector all the solution vectors in the problem. In a multi-field problem, each primal * field has a solution vector associated with it.*/ std::vector solutionSet; + + /** + * \brief Map of old solution vectors and their corresponding global field index. + */ + boost::unordered_map> solutionSet_previous; + /*Vector all the residual (RHS) vectors in the problem. In a multi-field problem, each * primal field has a residual vector associated with it.*/ std::vector residualSet; @@ -270,6 +276,14 @@ class MatrixFreePDE : public Subscriptor void applyBCs(unsigned int fieldIndex); + /** + * \brief Copy the solutions from the current timestep to previous time states. + */ + void + copy_solution_vectors(const std::vector &solutionSet, + boost::unordered_map> + &solutionSet_previous); + /*Method to compute the right hand side (RHS) residual vectors*/ void computeExplicitRHS(); diff --git a/include/model_variables.h b/include/model_variables.h index 6f3da0c2e..4d78b8af9 100644 --- a/include/model_variables.h +++ b/include/model_variables.h @@ -1,5 +1,3 @@ -// Model Variables Class - #ifndef INCLUDE_MODELVARIABLE_H_ #define INCLUDE_MODELVARIABLE_H_ @@ -33,11 +31,12 @@ class modelResidual struct variable_info { - bool is_scalar; unsigned int global_var_index; dealii::EvaluationFlags::EvaluationFlags evaluation_flags; dealii::EvaluationFlags::EvaluationFlags residual_flags; - bool var_needed; + bool var_needed = false; + bool is_scalar = true; + bool is_implicit = false; }; -#endif /* INCLUDE_MODELVARIABLE_H_ */ +#endif diff --git a/include/userInputParameters.h b/include/userInputParameters.h index cb42e886d..b372e6b8e 100644 --- a/include/userInputParameters.h +++ b/include/userInputParameters.h @@ -254,13 +254,19 @@ class userInputParameters std::vector var_nonlinear; // Variables needed to calculate the RHS - unsigned int num_var_explicit_RHS, num_var_nonexplicit_RHS; - std::vector varInfoListExplicitRHS, varInfoListNonexplicitRHS; + unsigned int num_var_explicit_RHS; + unsigned int num_var_nonexplicit_RHS; + unsigned int num_var_old_RHS; + std::vector varInfoListExplicitRHS; + std::vector varInfoListNonexplicitRHS; + std::vector varInfoList_old_RHS; // Variables needed to calculate the LHS unsigned int num_var_LHS; + unsigned int num_var_old_LHS; std::vector varInfoListLHS; std::vector varChangeInfoListLHS; + std::vector varInfoList_old_LHS; // Variables for loading in initial conditions std::vector load_ICs; diff --git a/include/variableContainer.h b/include/variableContainer.h index e5a0fb9be..9974ecd0d 100644 --- a/include/variableContainer.h +++ b/include/variableContainer.h @@ -1,5 +1,3 @@ -// This class permits the access of a subset of indexed fields and gives an -// error if any non-allowed fields are requested #ifndef VARIBLECONTAINER_H #define VARIBLECONTAINER_H @@ -11,132 +9,407 @@ #include +#include "model_variables.h" #include "userInputParameters.h" +/** + * \brief This class handles the access and submission of variables and residuals. + * + * \tparam dim The number of dimensions in the problem. + * \tparam degree The polynomial degree of the shape functions. + * \tparam T The type of the scalar variable. + */ template class variableContainer { public: #include "typeDefs.h" - // Constructors - - // Standard contructor, used for most situations + /** + * \brief Standard constructor for nonexplicit LHS + */ variableContainer(const dealii::MatrixFree &data, const std::vector &_varInfoList, - const std::vector &_varChangeInfoList); + const std::vector &_varChangeInfoList, + const std::vector &_varOldInfoList); + /** + * \brief Standard constructor for explicit & nonexplicit RHS + */ variableContainer(const dealii::MatrixFree &data, - const std::vector &_varInfoList); - // Nonstandard constructor, used when only one index of "data" should be used, - // use with care! + const std::vector &_varInfoList, + const std::vector &_varOldInfoList); + + /** + * \brief Nonstandard constructor for postprocessing when only one index of "data" + * should be used. + */ variableContainer(const dealii::MatrixFree &data, const std::vector &_varInfoList, const unsigned int &fixed_index); - // Methods to get the value/grad/hess in the residual method (this is how the - // user gets these values in equations.h) + /** + * \brief Return the value of the specified scalar field. + * + * \param global_variable_index The global index of the field. + */ T get_scalar_value(unsigned int global_variable_index) const; + + /** + * \brief Return the gradient of the specified scalar field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<1, dim, T> get_scalar_gradient(unsigned int global_variable_index) const; + + /** + * \brief Return the hessian of the specified scalar field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<2, dim, T> get_scalar_hessian(unsigned int global_variable_index) const; + + /** + * \brief Return the value of the specified vector field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<1, dim, T> get_vector_value(unsigned int global_variable_index) const; + + /** + * \brief Return the gradient of the specified vector field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<2, dim, T> get_vector_gradient(unsigned int global_variable_index) const; + + /** + * \brief Return the hessian of the specified vector field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<3, dim, T> get_vector_hessian(unsigned int global_variable_index) const; + /** + * \brief Return the change in value of the specified scalar field. + * + * \param global_variable_index The global index of the field. + */ T get_change_in_scalar_value(unsigned int global_variable_index) const; + + /** + * \brief Return the change in gradient of the specified scalar field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<1, dim, T> get_change_in_scalar_gradient(unsigned int global_variable_index) const; + + /** + * \brief Return the change in hessian of the specified scalar field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<2, dim, T> get_change_in_scalar_hessian(unsigned int global_variable_index) const; + + /** + * \brief Return the change in value of the specified vector field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<1, dim, T> get_change_in_vector_value(unsigned int global_variable_index) const; + + /** + * \brief Return the change in gradient of the specified vector field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<2, dim, T> get_change_in_vector_gradient(unsigned int global_variable_index) const; + + /** + * \brief Return the change in hessian of the specified vector field. + * + * \param global_variable_index The global index of the field. + */ dealii::Tensor<3, dim, T> get_change_in_vector_hessian(unsigned int global_variable_index) const; - // Methods to set the value residual and the gradient residual (this is how - // the user sets these values in equations.h) + /** + * \brief Return the old value of the specified scalar field. + * + * \param global_variable_index The global index of the field. + */ + T + get_old_scalar_value(unsigned int global_variable_index) const; + + /** + * \brief Return the old gradient of the specified scalar field. + * + * \param global_variable_index The global index of the field. + */ + dealii::Tensor<1, dim, T> + get_old_scalar_gradient(unsigned int global_variable_index) const; + + /** + * \brief Return the old hessian of the specified scalar field. + * + * \param global_variable_index The global index of the field. + */ + dealii::Tensor<2, dim, T> + get_old_scalar_hessian(unsigned int global_variable_index) const; + + /** + * \brief Return the old value of the specified vector field. + * + * \param global_variable_index The global index of the field. + */ + dealii::Tensor<1, dim, T> + get_old_vector_value(unsigned int global_variable_index) const; + + /** + * \brief Return the old gradient of the specified vector field. + * + * \param global_variable_index The global index of the field. + */ + dealii::Tensor<2, dim, T> + get_old_vector_gradient(unsigned int global_variable_index) const; + + /** + * \brief Return the old hessian of the specified vector field. + * + * \param global_variable_index The global index of the field. + */ + dealii::Tensor<3, dim, T> + get_old_vector_hessian(unsigned int global_variable_index) const; + + /** + * \brief Set the value residual of the specified scalar field. + * + * \param global_variable_index The global index of the field. + * \param val The value term. + */ void set_scalar_value_term_RHS(unsigned int global_variable_index, T val); + + /** + * \brief Set the gradient residual of the specified scalar field. + * + * \param global_variable_index The global index of the field. + * \param grad The gradient term. + */ void set_scalar_gradient_term_RHS(unsigned int global_variable_index, dealii::Tensor<1, dim, T> grad); + + /** + * \brief Set the value residual of the specified vector field. + * + * \param global_variable_index The global index of the field. + * \param val The value term. + */ void set_vector_value_term_RHS(unsigned int global_variable_index, dealii::Tensor<1, dim, T> val); + + /** + * \brief Set the gradient residual of the specified vector field. + * + * \param global_variable_index The global index of the field. + * \param grad The gradient term. + */ void set_vector_gradient_term_RHS(unsigned int global_variable_index, dealii::Tensor<2, dim, T> grad); + /** + * \brief Set the value residual of the specified scalar field. + * + * \param global_variable_index The global index of the field. + * \param val The value term. + */ void set_scalar_value_term_LHS(unsigned int global_variable_index, T val); + + /** + * \brief Set the gradient residual of the specified scalar field. + * + * \param global_variable_index The global index of the field. + * \param grad The gradient term. + */ void set_scalar_gradient_term_LHS(unsigned int global_variable_index, dealii::Tensor<1, dim, T> grad); + + /** + * \brief Set the value residual of the specified vector field. + * + * \param global_variable_index The global index of the field. + * \param val The value term. + */ void set_vector_value_term_LHS(unsigned int global_variable_index, dealii::Tensor<1, dim, T> val); + + /** + * \brief Set the gradient residual of the specified vector field. + * + * \param global_variable_index The global index of the field. + * \param grad The gradient term. + */ void set_vector_gradient_term_LHS(unsigned int global_variable_index, dealii::Tensor<2, dim, T> grad); - // Initialize, read DOFs, and set evaulation flags for each variable + /** + * \brief Initialize, read DOFs, and set evaulation flags for each variable. + * + * \param src The source vector. + * \param cell The cell where this is done. + */ void reinit_and_eval(const std::vector &src, unsigned int cell); + + /** + * \brief Initialize, read DOFs, and set evaulation flags for the change in each + * variable. + * + * \param src The source vector. + * \param cell The cell where this is done. + * \param var_being_solved The variable that we are evaluating. + */ void reinit_and_eval_change_in_solution(const vectorType &src, unsigned int cell, unsigned int var_being_solved); - // Only initialize the FEEvaluation object for each variable (used for - // post-processing) + /** + * \brief Initialize, read DOFs, and set evaulation flags for the previous timestep of + * each variable. + * + * \param src The source vector. + * \param cell The cell where this is done. + */ + void + reinit_and_eval_old_solution(const std::vector &src, unsigned int cell); + + /** + * Initialize the FEEvaluation object for each variable for post-processing since + * evaluation is uneccessary. + * + * \param cell The cell where this is done. + */ void reinit(unsigned int cell); - // Integrate the residuals and distribute from local to global + /** + * \brief Integrate the residuals and distribute from local to global for the RHS. + * + * \param dst The destination vector. + */ void integrate_and_distribute(std::vector &dst); + + /** + * \brief Integrate the residuals and distribute from local to global for the LHS. + * + * \param dst The destination vector. + * \param var_being_solved The variable that we are evaluating. + */ void integrate_and_distribute_change_in_solution_LHS(vectorType &dst, const unsigned int var_being_solved); - // The quadrature point index, a method to get the number of quadrature points - // per cell, and a method to get the xyz coordinates for the quadrature point + /** + * \brief The quadrature point index. + */ unsigned int q_point; - unsigned int + /** + * \brief Return the number of quadrature points per cell. + */ + [[nodiscard]] unsigned int get_num_q_points() const; + /** + * \brief Return the xyz coordinates of the quadrature point. + */ dealii::Point get_q_point_location() const; private: - // Vectors of the actual FEEvaluation objects for each active variable, split - // into scalar variables and vector variables for type reasons using scalar_FEEval = dealii::FEEvaluation; using vector_FEEval = dealii::FEEvaluation; + /** + * \brief A map of FEEvaluation objects for each active scalar variable at the current + * timestep, specified by global variable index. + */ boost::unordered_map> scalar_vars_map; + + /** + * \brief A map of FEEvaluation objects for each active vector variable at the current + * timestep, specified by global variable index. + */ boost::unordered_map> vector_vars_map; + /** + * \brief A map of FEEvaluation objects for the change in each active scalar variable at + * the current timestep, specified by global variable index. + */ boost::unordered_map> scalar_change_in_vars_map; + + /** + * \brief A map of FEEvaluation objects for the change in each active vector variable at + * the current timestep, specified by global variable index. + */ boost::unordered_map> vector_change_in_vars_map; - // Object containing some information about each variable (indices, whether - // the val/grad/hess is needed, etc) + /** + * \brief A map of FEEvaluation objects for each active scalar variable from the + * previous timestep, specified by global variable index. + */ + boost::unordered_map> scalar_old_vars_map; + + /** + * \brief A map of FEEvaluation objects for each active vector variable from the + * previous timestep, specified by global variable index. + */ + boost::unordered_map> vector_old_vars_map; + + // Object containing some information about each variable + + /** + * \brief Vector container for information about each variable (e.g., indices, whether + * the val/grad/hess is needed, etc). + */ std::vector varInfoList; + + /** + * \brief Vector container for information about the change in each variable (e.g., + * indices, whether the val/grad/hess is needed, etc). + */ std::vector varChangeInfoList; - // The number of variables + /** + * \brief Vector container for information about the previous timestep of each variable + * (e.g., indices, whether the val/grad/hess is needed, etc). + */ + std::vector varOldInfoList; + + /** + * \brief The number of fields. + */ unsigned int num_var; }; diff --git a/src/EquationDependencyParser/EquationDependencyParser.cc b/src/EquationDependencyParser/EquationDependencyParser.cc index ba199ba84..906b28e57 100644 --- a/src/EquationDependencyParser/EquationDependencyParser.cc +++ b/src/EquationDependencyParser/EquationDependencyParser.cc @@ -29,6 +29,8 @@ EquationDependencyParser::parse(std::vector &var_name, eval_flags_nonexplicit_RHS.resize(n_variables, dealii::EvaluationFlags::nothing); eval_flags_nonexplicit_LHS.resize(n_variables, dealii::EvaluationFlags::nothing); eval_flags_change_nonexplicit_LHS.resize(n_variables, dealii::EvaluationFlags::nothing); + eval_flags_old_RHS.resize(n_variables, dealii::EvaluationFlags::nothing); + eval_flags_old_LHS.resize(n_variables, dealii::EvaluationFlags::nothing); // Resize the residual evaluation flag vectors eval_flags_residual_explicit_RHS.resize(n_variables, dealii::EvaluationFlags::nothing); @@ -55,6 +57,7 @@ EquationDependencyParser::parse(std::vector &var_name, sorted_dependencies_value_RHS[i], sorted_dependencies_gradient_RHS[i], eval_flags_explicit_RHS, + eval_flags_old_RHS, eval_flags_residual_explicit_RHS, single_var_nonlinear); @@ -70,6 +73,7 @@ EquationDependencyParser::parse(std::vector &var_name, sorted_dependencies_value_RHS[i], sorted_dependencies_gradient_RHS[i], eval_flags_nonexplicit_RHS, + eval_flags_old_RHS, eval_flags_residual_nonexplicit_RHS, single_var_nonlinear); @@ -86,6 +90,7 @@ EquationDependencyParser::parse(std::vector &var_name, sorted_dependencies_value_RHS[i], sorted_dependencies_gradient_RHS[i], eval_flags_nonexplicit_RHS, + eval_flags_old_RHS, eval_flags_residual_nonexplicit_RHS, single_var_nonlinear_RHS); @@ -95,6 +100,7 @@ EquationDependencyParser::parse(std::vector &var_name, sorted_dependencies_value_LHS[i], sorted_dependencies_gradient_LHS[i], eval_flags_nonexplicit_LHS, + eval_flags_old_LHS, eval_flags_change_nonexplicit_LHS, eval_flags_residual_nonexplicit_LHS, single_var_nonlinear_LHS); @@ -112,6 +118,7 @@ EquationDependencyParser::parseDependencyListRHS( std::string &value_dependencies, std::string &gradient_dependencies, std::vector &evaluation_flags, + std::vector &old_flags, std::vector &residual_flags, bool &is_nonlinear) { @@ -165,6 +172,21 @@ EquationDependencyParser::parseDependencyListRHS( variable.begin(), variable.end()); + std::string old_value_variable = {"old()"}; + old_value_variable.insert(--old_value_variable.end(), + variable.begin(), + variable.end()); + + std::string old_gradient_variable = {"grad(old())"}; + old_gradient_variable.insert(--(--old_gradient_variable.end()), + variable.begin(), + variable.end()); + + std::string old_hessian_variable = {"hess(old())"}; + old_hessian_variable.insert(--(--old_hessian_variable.end()), + variable.begin(), + variable.end()); + // Is the variable we are finding the dependencies for explicit bool variable_is_explicit = variable_eq_type[variable_index] == EXPLICIT_TIME_DEPENDENT; @@ -197,6 +219,36 @@ EquationDependencyParser::parseDependencyListRHS( dealii::EvaluationFlags::hessians; dependency_entry_assigned = true; } + // Case if the dependency is old(x) + else if (dependency == old_value_variable) + { + AssertThrow(!dependency_variable_is_explicit, + dealii::ExcMessage( + "PRISMS-PF Error: Currently, old variables are only " + "available for implicit equations.")); + old_flags[dependency_variable_index] |= dealii::EvaluationFlags::values; + dependency_entry_assigned = true; + } + // Case if the dependency is grad(old(x)) + else if (dependency == old_gradient_variable) + { + AssertThrow(!dependency_variable_is_explicit, + dealii::ExcMessage( + "PRISMS-PF Error: Currently, old variables are only " + "available for implicit equations.")); + old_flags[dependency_variable_index] |= dealii::EvaluationFlags::gradients; + dependency_entry_assigned = true; + } + // Case if the dependency is hess(old(x)) + else if (dependency == old_hessian_variable) + { + AssertThrow(!dependency_variable_is_explicit, + dealii::ExcMessage( + "PRISMS-PF Error: Currently, old variables are only " + "available for implicit equations.")); + old_flags[dependency_variable_index] |= dealii::EvaluationFlags::hessians; + dependency_entry_assigned = true; + } // Check for nonlinearity is_nonlinear = @@ -220,6 +272,7 @@ EquationDependencyParser::parseDependencyListLHS( std::string &value_dependencies, std::string &gradient_dependencies, std::vector &evaluation_flags, + std::vector &old_flags, std::vector &change_flags, std::vector &residual_flags, bool &is_nonlinear) @@ -290,10 +343,29 @@ EquationDependencyParser::parseDependencyListLHS( variable.begin(), variable.end()); + std::string old_value_variable = {"old()"}; + old_value_variable.insert(--old_value_variable.end(), + variable.begin(), + variable.end()); + + std::string old_gradient_variable = {"grad(old())"}; + old_gradient_variable.insert(--(--old_gradient_variable.end()), + variable.begin(), + variable.end()); + + std::string old_hessian_variable = {"hess(old())"}; + old_hessian_variable.insert(--(--old_hessian_variable.end()), + variable.begin(), + variable.end()); + // Is the variable we are finding the dependencies for explicit bool dependency_variable_is_explicit = variable_eq_type[dependency_variable_index] == EXPLICIT_TIME_DEPENDENT; + // Is the variable we are finding the dependencies for implicit + bool dependency_variable_is_implicit = + variable_eq_type[dependency_variable_index] == IMPLICIT_TIME_DEPENDENT; + // Case if the dependency is x if (dependency == variable) { @@ -321,6 +393,45 @@ EquationDependencyParser::parseDependencyListLHS( dealii::EvaluationFlags::hessians; dependency_entry_assigned = true; + // Check for nonlinearity + is_nonlinear = !dependency_variable_is_explicit; + } + // Case if the dependency is old(x) + else if (dependency == old_value_variable) + { + AssertThrow(dependency_variable_is_implicit, + dealii::ExcMessage( + "PRISMS-PF Error: Currently, old variables are only " + "available for implicit equations.")); + old_flags[dependency_variable_index] |= dealii::EvaluationFlags::values; + dependency_entry_assigned = true; + + // Check for nonlinearity + is_nonlinear = !dependency_variable_is_explicit; + } + // Case if the dependency is grad(old(x)) + else if (dependency == old_gradient_variable) + { + AssertThrow(dependency_variable_is_implicit, + dealii::ExcMessage( + "PRISMS-PF Error: Currently, old variables are only " + "available for implicit equations.")); + old_flags[dependency_variable_index] |= dealii::EvaluationFlags::gradients; + dependency_entry_assigned = true; + + // Check for nonlinearity + is_nonlinear = !dependency_variable_is_explicit; + } + // Case if the dependency is hess(old(x)) + else if (dependency == old_hessian_variable) + { + AssertThrow(dependency_variable_is_implicit, + dealii::ExcMessage( + "PRISMS-PF Error: Currently, old variables are only " + "available for implicit equations.")); + old_flags[dependency_variable_index] |= dealii::EvaluationFlags::hessians; + dependency_entry_assigned = true; + // Check for nonlinearity is_nonlinear = !dependency_variable_is_explicit; } diff --git a/src/matrixfree/computeLHS.cc b/src/matrixfree/computeLHS.cc index 322e51f67..334d02ee4 100644 --- a/src/matrixfree/computeLHS.cc +++ b/src/matrixfree/computeLHS.cc @@ -57,8 +57,11 @@ MatrixFreePDE::getLHS( const vectorType &src, const std::pair &cell_range) const { - variableContainer> - variable_list(data, userInputs.varInfoListLHS, userInputs.varChangeInfoListLHS); + variableContainer> variable_list( + data, + userInputs.varInfoListLHS, + userInputs.varChangeInfoListLHS, + userInputs.varInfoList_old_LHS); // loop over cells for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) @@ -66,6 +69,7 @@ MatrixFreePDE::getLHS( // Initialize, read DOFs, and set evaulation flags for each variable variable_list.reinit_and_eval(solutionSet, cell); variable_list.reinit_and_eval_change_in_solution(src, cell, currentFieldIndex); + variable_list.reinit_and_eval_old_solution(solutionSet, cell); unsigned int num_q_points = variable_list.get_num_q_points(); diff --git a/src/matrixfree/computeRHS.cc b/src/matrixfree/computeRHS.cc index 2e1a48d3d..2876b01ac 100644 --- a/src/matrixfree/computeRHS.cc +++ b/src/matrixfree/computeRHS.cc @@ -32,13 +32,15 @@ MatrixFreePDE::getExplicitRHS( { variableContainer> variable_list( data, - userInputs.varInfoListExplicitRHS); + userInputs.varInfoListExplicitRHS, + userInputs.varInfoList_old_RHS); // loop over cells for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) { // Initialize, read DOFs, and set evaulation flags for each variable variable_list.reinit_and_eval(src, cell); + variable_list.reinit_and_eval_old_solution(src, cell); unsigned int num_q_points = variable_list.get_num_q_points(); @@ -87,13 +89,15 @@ MatrixFreePDE::getNonexplicitRHS( { variableContainer> variable_list( data, - userInputs.varInfoListNonexplicitRHS); + userInputs.varInfoListNonexplicitRHS, + userInputs.varInfoList_old_RHS); // loop over cells for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) { // Initialize, read DOFs, and set evaulation flags for each variable variable_list.reinit_and_eval(src, cell); + variable_list.reinit_and_eval_old_solution(src, cell); unsigned int num_q_points = variable_list.get_num_q_points(); diff --git a/src/matrixfree/init.cc b/src/matrixfree/init.cc index fb9bb16ef..10cb38eae 100644 --- a/src/matrixfree/init.cc +++ b/src/matrixfree/init.cc @@ -95,10 +95,6 @@ MatrixFreePDE::init() { isTimeDependentBVP = true; hasNonExplicitEquation = true; - std::cerr << "PRISMS-PF Error: IMPLICIT_TIME_DEPENDENT equation " - "types are not currently supported" - << std::endl; - abort(); } else if (field.pdetype == AUXILIARY) { @@ -261,18 +257,29 @@ MatrixFreePDE::init() pcout << "initializing parallel::distributed residual and solution vectors\n"; for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) { - vectorType *U, *R; + auto *U = new vectorType; + auto *R = new vectorType; - U = new vectorType; - R = new vectorType; solutionSet.push_back(U); residualSet.push_back(R); + matrixFreeObject.initialize_dof_vector(*R, fieldIndex); *R = 0; matrixFreeObject.initialize_dof_vector(*U, fieldIndex); *U = 0; + // If the equation type is implicit initialize a solution vector for the previous + // timestep. In the future, we should support arbitrary saving of old timesteps. + if (fields[fieldIndex].pdetype == TIME_INDEPENDENT) + { + std::unique_ptr U_previous = std::make_unique(); + solutionSet_previous.emplace(fieldIndex, std::move(U_previous)); + matrixFreeObject.initialize_dof_vector(*solutionSet_previous[fieldIndex], + fieldIndex); + *solutionSet_previous[fieldIndex] = 0; + } + // Initializing temporary dU vector required for implicit solves of the // elliptic equation. if (fields[fieldIndex].pdetype == TIME_INDEPENDENT || diff --git a/src/matrixfree/invM.cc b/src/matrixfree/invM.cc index bd39ab70a..158f7a4a4 100644 --- a/src/matrixfree/invM.cc +++ b/src/matrixfree/invM.cc @@ -44,20 +44,6 @@ MatrixFreePDE::computeInvM() } } - // Check if invM has been found. If not, print an "error" message - if (!invMscalarFound) - { - pcout << "matrixFreePDE.h: no PARABOLIC scalar field... hence setting " - "parabolicScalarFieldIndex to 0 and marching ahead withn invM " - "computation\n"; - } - else if (!invMvectorFound) - { - pcout << "matrixFreePDE.h: no PARABOLIC vector field... hence setting " - "parabolicVectorFieldIndex to 0 and marching ahead withn invM " - "computation\n"; - } - // Initialize invM and clear its values matrixFreeObject.initialize_dof_vector(invMscalar, parabolicScalarFieldIndex); invMscalar = 0.0; diff --git a/src/matrixfree/postprocessor.cc b/src/matrixfree/postprocessor.cc index b92162a9b..2c8cd01cf 100644 --- a/src/matrixfree/postprocessor.cc +++ b/src/matrixfree/postprocessor.cc @@ -39,9 +39,8 @@ MatrixFreePDE::getPostProcessedFields( const std::pair &cell_range) { // initialize FEEvaulation objects - variableContainer> variable_list( - data, - userInputs.pp_baseVarInfoList); + variableContainer> + variable_list(data, userInputs.pp_baseVarInfoList, userInputs.varInfoList_old_RHS); variableContainer> pp_variable_list(data, userInputs.pp_varInfoList, 0); @@ -50,6 +49,7 @@ MatrixFreePDE::getPostProcessedFields( { // Initialize, read DOFs, and set evaulation flags for each variable variable_list.reinit_and_eval(src, cell); + variable_list.reinit_and_eval_old_solution(src, cell); pp_variable_list.reinit(cell); unsigned int num_q_points = variable_list.get_num_q_points(); diff --git a/src/matrixfree/reinit.cc b/src/matrixfree/reinit.cc index 097b173a5..f06c8c5c0 100644 --- a/src/matrixfree/reinit.cc +++ b/src/matrixfree/reinit.cc @@ -122,13 +122,20 @@ MatrixFreePDE::reinit() pcout << "initializing parallel::distributed residual and solution vectors\n"; for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) { - vectorType *U; - - U = solutionSet.at(fieldIndex); + vectorType *U = solutionSet.at(fieldIndex); matrixFreeObject.initialize_dof_vector(*U, fieldIndex); *U = 0; + // If the equation type is implicit reinitialize a solution vector for the previous + // timestep. In the future, we should support arbitrary saving of old timesteps. + if (fields[fieldIndex].pdetype == TIME_INDEPENDENT) + { + matrixFreeObject.initialize_dof_vector(*solutionSet_previous[fieldIndex], + fieldIndex); + *solutionSet_previous[fieldIndex] = 0; + } + // Initializing temporary dU vector required for implicit solves of the // elliptic equation. if (fields[fieldIndex].pdetype == TIME_INDEPENDENT || diff --git a/src/matrixfree/solveIncrement.cc b/src/matrixfree/solveIncrement.cc index 3199cf0ff..37778098f 100644 --- a/src/matrixfree/solveIncrement.cc +++ b/src/matrixfree/solveIncrement.cc @@ -14,6 +14,9 @@ MatrixFreePDE::solveIncrement(bool skip_time_dependent) Timer time; char buffer[200]; + // Copy solution vectors + copy_solution_vectors(solutionSet, solutionSet_previous); + // Get the RHS of the explicit equations if (hasExplicitEquation && !skip_time_dependent) { @@ -547,4 +550,28 @@ MatrixFreePDE::updateImplicitSolution(unsigned int fieldIndex, return nonlinear_it_converged; } +template +void +MatrixFreePDE::copy_solution_vectors( + const std::vector &solutionSet, + boost::unordered_map> &solutionSet_previous) +{ + for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) + { + // Skip the copy if not tracking prior solutions for the variable index. + if (!userInputs.varInfoListNonexplicitRHS[fieldIndex].is_implicit) + { + continue; + } + + // Copy the solution vector. + for (unsigned int dof = 0; dof < solutionSet[fieldIndex]->locally_owned_size(); + ++dof) + { + solutionSet_previous[fieldIndex]->local_element(dof) = + solutionSet[fieldIndex]->local_element(dof); + } + } +} + #include "../../include/matrixFreePDE_template_instantiations.h" diff --git a/src/userInputParameters/loadVariableAttributes.cc b/src/userInputParameters/loadVariableAttributes.cc index 942936973..e85491088 100644 --- a/src/userInputParameters/loadVariableAttributes.cc +++ b/src/userInputParameters/loadVariableAttributes.cc @@ -69,7 +69,7 @@ userInputParameters::loadVariableAttributes( varInfoListExplicitRHS.push_back(varInfo); } - // Load variable information for calculating the RHS for nonexplicit equations + // Load the variable information for calculating the RHS for nonexplicit equations num_var_nonexplicit_RHS = 0; for (unsigned int i = 0; i < number_of_variables; i++) { @@ -96,10 +96,41 @@ userInputParameters::loadVariableAttributes( varInfo.is_scalar = var_type[i] == SCALAR; + varInfo.is_implicit = var_eq_type[i] == IMPLICIT_TIME_DEPENDENT; + varInfoListNonexplicitRHS.push_back(varInfo); } + // Load the old variable information for calculating the RHS for nonexplicit equations + num_var_old_RHS = 0; + for (unsigned int i = 0; i < number_of_variables; i++) + { + if (!(variable_attributes.equation_dependency_parser.eval_flags_old_RHS[i] & + dealii::EvaluationFlags::nothing)) + { + num_var_old_RHS++; + } + } + varInfoList_old_RHS.reserve(num_var_old_RHS); + for (unsigned int i = 0; i < number_of_variables; i++) + { + variable_info varInfo; + + varInfo.evaluation_flags = + variable_attributes.equation_dependency_parser.eval_flags_old_RHS[i]; + + varInfo.residual_flags = variable_attributes.equation_dependency_parser + .eval_flags_residual_nonexplicit_RHS[i]; + + varInfo.global_var_index = i; + + varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing); + + varInfo.is_scalar = var_type[i] == SCALAR; + + varInfoList_old_RHS.push_back(varInfo); + } - // Load variable information for calculating the LHS + // Load variable information for calculating the LHS for nonexplicit equations num_var_LHS = 0; for (unsigned int i = 0; i < number_of_variables; i++) { @@ -127,9 +158,12 @@ userInputParameters::loadVariableAttributes( varInfo.is_scalar = var_type[i] == SCALAR; + varInfo.is_implicit = var_eq_type[i] == IMPLICIT_TIME_DEPENDENT; + varInfoListLHS.push_back(varInfo); } - + // Load the change variable information for calculating the LHS for nonexplicit + // equations varChangeInfoListLHS.reserve(num_var_LHS); for (unsigned int i = 0; i < number_of_variables; i++) { @@ -138,7 +172,6 @@ userInputParameters::loadVariableAttributes( varInfo.evaluation_flags = variable_attributes.equation_dependency_parser .eval_flags_change_nonexplicit_LHS[i]; - // FOR NOW, TAKING THESE FROM THE VARIABLE ITSELF!! varInfo.residual_flags = variable_attributes.equation_dependency_parser .eval_flags_residual_nonexplicit_LHS[i]; @@ -151,6 +184,36 @@ userInputParameters::loadVariableAttributes( varChangeInfoListLHS.push_back(varInfo); } + // Load the old variable information for calculating the LHS for nonexplicit equations + num_var_old_LHS = 0; + for (unsigned int i = 0; i < number_of_variables; i++) + { + if (!(variable_attributes.equation_dependency_parser.eval_flags_old_LHS[i] & + dealii::EvaluationFlags::nothing)) + { + num_var_old_LHS++; + } + } + varInfoList_old_LHS.reserve(num_var_old_LHS); + for (unsigned int i = 0; i < number_of_variables; i++) + { + variable_info varInfo; + + varInfo.evaluation_flags = + variable_attributes.equation_dependency_parser.eval_flags_old_LHS[i]; + + varInfo.residual_flags = variable_attributes.equation_dependency_parser + .eval_flags_residual_nonexplicit_LHS[i]; + + varInfo.global_var_index = i; + + varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing); + + varInfo.is_scalar = var_type[i] == SCALAR; + + varInfoList_old_LHS.push_back(varInfo); + } + // Load variable information for postprocessing // First, the info list for the base field variables pp_baseVarInfoList.reserve(number_of_variables); diff --git a/src/variableContainer/variableContainer.cc b/src/variableContainer/variableContainer.cc index 35ca18ffe..682c315d8 100644 --- a/src/variableContainer/variableContainer.cc +++ b/src/variableContainer/variableContainer.cc @@ -4,15 +4,18 @@ template variableContainer::variableContainer( const dealii::MatrixFree &data, const std::vector &_varInfoList, - const std::vector &_varChangeInfoList) + const std::vector &_varChangeInfoList, + const std::vector &_varOldInfoList) : varInfoList(_varInfoList) , varChangeInfoList(_varChangeInfoList) + , varOldInfoList(_varOldInfoList) , num_var(varInfoList.size()) { for (unsigned int i = 0; i < num_var; i++) { const auto &var_info = varInfoList[i]; const auto &var_change_info = varChangeInfoList[i]; + const auto &var_old_info = varOldInfoList[i]; if (var_info.var_needed) { @@ -45,34 +48,69 @@ variableContainer::variableContainer( std::make_unique(data, i)); } } + + if (var_old_info.var_needed) + { + const unsigned int var_index = var_old_info.global_var_index; + + if (var_old_info.is_scalar) + { + scalar_old_vars_map.emplace(var_index, + std::make_unique(data, i)); + } + else + { + vector_old_vars_map.emplace(var_index, + std::make_unique(data, i)); + } + } } } template variableContainer::variableContainer( const dealii::MatrixFree &data, - const std::vector &_varInfoList) + const std::vector &_varInfoList, + const std::vector &_varOldInfoList) : varInfoList(_varInfoList) + , varOldInfoList(_varOldInfoList) , num_var(varInfoList.size()) { for (unsigned int i = 0; i < num_var; i++) { - const auto &var_info = varInfoList[i]; + const auto &var_info = varInfoList[i]; + const auto &var_old_info = varOldInfoList[i]; - if (!var_info.var_needed) + if (var_info.var_needed) { - continue; - } - - const unsigned int var_index = var_info.global_var_index; + const unsigned int var_index = var_info.global_var_index; - if (var_info.is_scalar) - { - scalar_vars_map.emplace(var_index, std::make_unique(data, i)); + if (var_info.is_scalar) + { + scalar_vars_map.emplace(var_index, + std::make_unique(data, i)); + } + else + { + vector_vars_map.emplace(var_index, + std::make_unique(data, i)); + } } - else + + if (var_old_info.var_needed) { - vector_vars_map.emplace(var_index, std::make_unique(data, i)); + const unsigned int var_index = var_old_info.global_var_index; + + if (var_old_info.is_scalar) + { + scalar_old_vars_map.emplace(var_index, + std::make_unique(data, i)); + } + else + { + vector_old_vars_map.emplace(var_index, + std::make_unique(data, i)); + } } } } @@ -229,6 +267,40 @@ variableContainer::reinit_and_eval_change_in_solution( } } +template +void +variableContainer::reinit_and_eval_old_solution( + const std::vector &src, + unsigned int cell) +{ + for (unsigned int i = 0; i < num_var; i++) + { + const auto &var_old_info = varOldInfoList[i]; + + if (!var_old_info.var_needed) + { + continue; + } + + const unsigned int var_index = var_old_info.global_var_index; + + if (var_old_info.is_scalar) + { + auto *scalar_FEEval_ptr = scalar_old_vars_map[var_index].get(); + scalar_FEEval_ptr->reinit(cell); + scalar_FEEval_ptr->read_dof_values(*src[i]); + scalar_FEEval_ptr->evaluate(var_old_info.evaluation_flags); + } + else + { + auto *vector_FEEval_ptr = vector_old_vars_map[var_index].get(); + vector_FEEval_ptr->reinit(cell); + vector_FEEval_ptr->read_dof_values(*src[i]); + vector_FEEval_ptr->evaluate(var_old_info.evaluation_flags); + } + } +} + template void variableContainer::reinit(unsigned int cell) @@ -307,26 +379,20 @@ variableContainer::integrate_and_distribute_change_in_solution_L } } -// Need to add index checking to these functions so that an error is thrown if -// the index wasn't set template T variableContainer::get_scalar_value( unsigned int global_variable_index) const { - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::values) - { - return scalar_vars_map.at(global_variable_index)->get_value(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable value that " - "was not marked as needed in 'equations.cc'. The attempted " - "access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::values, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a scalar variable value " + "that was not marked as needed in 'equations.cc'. The attempted " + "access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return scalar_vars_map.at(global_variable_index)->get_value(q_point); } template @@ -334,19 +400,15 @@ dealii::Tensor<1, dim, T> variableContainer::get_scalar_gradient( unsigned int global_variable_index) const { - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::gradients) - { - return scalar_vars_map.at(global_variable_index)->get_gradient(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable gradient " - "that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::gradients, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a scalar variable gradient " + "that was not marked as needed in 'equations.cc'. The attempted " + "access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return scalar_vars_map.at(global_variable_index)->get_gradient(q_point); } template @@ -354,19 +416,15 @@ dealii::Tensor<2, dim, T> variableContainer::get_scalar_hessian( unsigned int global_variable_index) const { - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::hessians) - { - return scalar_vars_map.at(global_variable_index)->get_hessian(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable hessian " - "that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::hessians, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a scalar variable hessian " + "that was not marked as needed in 'equations.cc'. The attempted " + "access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return scalar_vars_map.at(global_variable_index)->get_hessian(q_point); } template @@ -374,19 +432,15 @@ dealii::Tensor<1, dim, T> variableContainer::get_vector_value( unsigned int global_variable_index) const { - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::values) - { - return vector_vars_map.at(global_variable_index)->get_value(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable value that " - "was not marked as needed in 'equations.cc'. The attempted " - "access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::values, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a vector variable value " + "that was not marked as needed in 'equations.cc'. The attempted " + "access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return vector_vars_map.at(global_variable_index)->get_value(q_point); } template @@ -394,19 +448,15 @@ dealii::Tensor<2, dim, T> variableContainer::get_vector_gradient( unsigned int global_variable_index) const { - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::gradients) - { - return vector_vars_map.at(global_variable_index)->get_gradient(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable gradient " - "that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::gradients, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a vector variable gradient " + "that was not marked as needed in 'equations.cc'. The attempted " + "access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return vector_vars_map.at(global_variable_index)->get_gradient(q_point); } template @@ -414,41 +464,31 @@ dealii::Tensor<3, dim, T> variableContainer::get_vector_hessian( unsigned int global_variable_index) const { - if (varInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::hessians) - { - return vector_vars_map.at(global_variable_index)->get_hessian(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a variable hessian " - "that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::hessians, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a vector variable hessian " + "that was not marked as needed in 'equations.cc'. The attempted " + "access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return vector_vars_map.at(global_variable_index)->get_hessian(q_point); } -// Need to add index checking to these functions so that an error is thrown if -// the index wasn't set template T variableContainer::get_change_in_scalar_value( unsigned int global_variable_index) const { - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::values) - { - return scalar_change_in_vars_map.at(global_variable_index)->get_value(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "value that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varChangeInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::values, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a change in scalar variable " + "value that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return scalar_change_in_vars_map.at(global_variable_index)->get_value(q_point); } template @@ -456,19 +496,15 @@ dealii::Tensor<1, dim, T> variableContainer::get_change_in_scalar_gradient( unsigned int global_variable_index) const { - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::gradients) - { - return scalar_change_in_vars_map.at(global_variable_index)->get_gradient(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "gradient that was not marked as needed in 'equations.cc'. " - "The attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varChangeInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::gradients, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a change in scalar variable " + "gradient that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return scalar_change_in_vars_map.at(global_variable_index)->get_gradient(q_point); } template @@ -476,19 +512,15 @@ dealii::Tensor<2, dim, T> variableContainer::get_change_in_scalar_hessian( unsigned int global_variable_index) const { - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::hessians) - { - return scalar_change_in_vars_map.at(global_variable_index)->get_hessian(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "hessian that was not marked as needed in 'equations.cc'. " - "The attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varChangeInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::hessians, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a change in scalar variable " + "hessian that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return scalar_change_in_vars_map.at(global_variable_index)->get_hessian(q_point); } template @@ -496,19 +528,15 @@ dealii::Tensor<1, dim, T> variableContainer::get_change_in_vector_value( unsigned int global_variable_index) const { - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::values) - { - return vector_change_in_vars_map.at(global_variable_index)->get_value(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "value that was not marked as needed in 'equations.cc'. The " - "attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varChangeInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::values, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a change in vector variable " + "value that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return vector_change_in_vars_map.at(global_variable_index)->get_value(q_point); } template @@ -516,19 +544,15 @@ dealii::Tensor<2, dim, T> variableContainer::get_change_in_vector_gradient( unsigned int global_variable_index) const { - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::gradients) - { - return vector_change_in_vars_map.at(global_variable_index)->get_gradient(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "gradient that was not marked as needed in 'equations.cc'. " - "The attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varChangeInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::gradients, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a change in vector variable " + "gradient that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return vector_change_in_vars_map.at(global_variable_index)->get_gradient(q_point); } template @@ -536,19 +560,111 @@ dealii::Tensor<3, dim, T> variableContainer::get_change_in_vector_hessian( unsigned int global_variable_index) const { - if (varChangeInfoList[global_variable_index].evaluation_flags & - dealii::EvaluationFlags::hessians) - { - return vector_change_in_vars_map.at(global_variable_index)->get_hessian(q_point); - } - else - { - std::cerr << "PRISMS-PF Error: Attempted access of a change in variable " - "hessian that was not marked as needed in 'equations.cc'. " - "The attempted access was for variable with index " - << global_variable_index << " ." << std::endl; - abort(); - } + Assert(varChangeInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::hessians, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a change in vector variable " + "hessian that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return vector_change_in_vars_map.at(global_variable_index)->get_hessian(q_point); +} + +template +T +variableContainer::get_old_scalar_value( + unsigned int global_variable_index) const +{ + Assert(varOldInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::values, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a previous scalar variable " + "value that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return scalar_old_vars_map.at(global_variable_index)->get_value(q_point); +} + +template +dealii::Tensor<1, dim, T> +variableContainer::get_old_scalar_gradient( + unsigned int global_variable_index) const +{ + Assert(varOldInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::gradients, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a previous scalar variable " + "gradient that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return scalar_old_vars_map.at(global_variable_index)->get_gradient(q_point); +} + +template +dealii::Tensor<2, dim, T> +variableContainer::get_old_scalar_hessian( + unsigned int global_variable_index) const +{ + Assert(varOldInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::hessians, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a previous scalar variable " + "hessian that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + ".")); + + return scalar_old_vars_map.at(global_variable_index)->get_hessian(q_point); +} + +template +dealii::Tensor<1, dim, T> +variableContainer::get_old_vector_value( + unsigned int global_variable_index) const +{ + Assert(varOldInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::values, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a previous vector variable " + "value that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + " .")); + + return vector_old_vars_map.at(global_variable_index)->get_value(q_point); +} + +template +dealii::Tensor<2, dim, T> +variableContainer::get_old_vector_gradient( + unsigned int global_variable_index) const +{ + Assert(varOldInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::gradients, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a previous vector variable " + "gradient that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + " .")); + + return vector_old_vars_map.at(global_variable_index)->get_gradient(q_point); +} + +template +dealii::Tensor<3, dim, T> +variableContainer::get_old_vector_hessian( + unsigned int global_variable_index) const +{ + Assert(varOldInfoList[global_variable_index].evaluation_flags & + dealii::EvaluationFlags::hessians, + dealii::ExcMessage( + std::string("PRISMS-PF Error: Attempted access of a previous vector variable " + "hessian that was not marked as needed in 'equations.cc'. The " + "attempted access was for variable with index ") + + std::to_string(global_variable_index) + " .")); + + return vector_old_vars_map.at(global_variable_index)->get_hessian(q_point); } // The methods to set the residual terms