From 7cd35c977429188e6aebc1660a628be01e79623a Mon Sep 17 00:00:00 2001 From: landinjm Date: Thu, 31 Oct 2024 16:43:15 -0400 Subject: [PATCH] Revert "Merge pull request #268 from landinjm/implicit_equation_types" This reverts commit f97e3990ffbd68833c803787d4210ffd3853f375, reversing changes made to 7d3638790115517c0b83607292c6c824722b5c5f. --- applications/advection/CMakeLists.txt | 68 --- applications/advection/ICs_and_BCs.cc | 71 --- applications/advection/advection.md | 40 -- applications/advection/customPDE.h | 94 ---- applications/advection/equations.cc | 132 ----- applications/advection/parameters.prm | 134 ----- include/EquationDependencyParser.h | 4 - include/matrixFreePDE.h | 14 - include/model_variables.h | 9 +- include/userInputParameters.h | 10 +- include/variableContainer.h | 323 +----------- .../EquationDependencyParser.cc | 111 ----- src/matrixfree/computeLHS.cc | 8 +- src/matrixfree/computeRHS.cc | 8 +- src/matrixfree/init.cc | 21 +- src/matrixfree/invM.cc | 14 + src/matrixfree/postprocessor.cc | 6 +- src/matrixfree/reinit.cc | 13 +- src/matrixfree/solveIncrement.cc | 27 - .../loadVariableAttributes.cc | 71 +-- src/variableContainer/variableContainer.cc | 462 +++++++----------- 21 files changed, 240 insertions(+), 1400 deletions(-) delete mode 100644 applications/advection/CMakeLists.txt delete mode 100644 applications/advection/ICs_and_BCs.cc delete mode 100644 applications/advection/advection.md delete mode 100644 applications/advection/customPDE.h delete mode 100644 applications/advection/equations.cc delete mode 100644 applications/advection/parameters.prm diff --git a/applications/advection/CMakeLists.txt b/applications/advection/CMakeLists.txt deleted file mode 100644 index 9522ce21e..000000000 --- a/applications/advection/CMakeLists.txt +++ /dev/null @@ -1,68 +0,0 @@ -## -# 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 deleted file mode 100644 index cf986f563..000000000 --- a/applications/advection/ICs_and_BCs.cc +++ /dev/null @@ -1,71 +0,0 @@ -// =========================================================================== -// 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 deleted file mode 100644 index 5f411c4ac..000000000 --- a/applications/advection/advection.md +++ /dev/null @@ -1,40 +0,0 @@ -## 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 deleted file mode 100644 index 718d77fcd..000000000 --- a/applications/advection/customPDE.h +++ /dev/null @@ -1,94 +0,0 @@ -#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]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) 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]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) const override; - - // Function to set the LHS of the governing equations (in equations.cc) - void - equationLHS( - [[maybe_unused]] variableContainer> - &variable_list, - [[maybe_unused]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) 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, - [[maybe_unused]] const VectorizedArray element_volume) 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 deleted file mode 100644 index aeda295e9..000000000 --- a/applications/advection/equations.cc +++ /dev/null @@ -1,132 +0,0 @@ -// ================================================================================= -// 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]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) 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]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) 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]] const Point> q_point_loc, - [[maybe_unused]] const VectorizedArray element_volume) 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 deleted file mode 100644 index 983d05247..000000000 --- a/applications/advection/parameters.prm +++ /dev/null @@ -1,134 +0,0 @@ -# ================================================================================= -# 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 4f784ed39..4a8dd1b03 100644 --- a/include/EquationDependencyParser.h +++ b/include/EquationDependencyParser.h @@ -38,8 +38,6 @@ 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 @@ -71,7 +69,6 @@ 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); @@ -87,7 +84,6 @@ 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 6c6737a86..27be76f2b 100644 --- a/include/matrixFreePDE.h +++ b/include/matrixFreePDE.h @@ -226,12 +226,6 @@ 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; @@ -287,14 +281,6 @@ class MatrixFreePDE : public Subscriptor */ dealii::AlignedVector> element_volume; - /** - * \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 4d78b8af9..6f3da0c2e 100644 --- a/include/model_variables.h +++ b/include/model_variables.h @@ -1,3 +1,5 @@ +// Model Variables Class + #ifndef INCLUDE_MODELVARIABLE_H_ #define INCLUDE_MODELVARIABLE_H_ @@ -31,12 +33,11 @@ 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 = false; - bool is_scalar = true; - bool is_implicit = false; + bool var_needed; }; -#endif +#endif /* INCLUDE_MODELVARIABLE_H_ */ diff --git a/include/userInputParameters.h b/include/userInputParameters.h index 55ea4e2ef..92af67224 100644 --- a/include/userInputParameters.h +++ b/include/userInputParameters.h @@ -259,19 +259,13 @@ class userInputParameters std::vector var_nonlinear; // Variables needed to calculate the RHS - 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; + unsigned int num_var_explicit_RHS, num_var_nonexplicit_RHS; + std::vector varInfoListExplicitRHS, varInfoListNonexplicitRHS; // 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 9974ecd0d..e5a0fb9be 100644 --- a/include/variableContainer.h +++ b/include/variableContainer.h @@ -1,3 +1,5 @@ +// 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 @@ -9,407 +11,132 @@ #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" - /** - * \brief Standard constructor for nonexplicit LHS - */ - variableContainer(const dealii::MatrixFree &data, - const std::vector &_varInfoList, - const std::vector &_varChangeInfoList, - const std::vector &_varOldInfoList); + // Constructors - /** - * \brief Standard constructor for explicit & nonexplicit RHS - */ + // Standard contructor, used for most situations variableContainer(const dealii::MatrixFree &data, const std::vector &_varInfoList, - const std::vector &_varOldInfoList); + const std::vector &_varChangeInfoList); - /** - * \brief Nonstandard constructor for postprocessing when only one index of "data" - * should be used. - */ + variableContainer(const dealii::MatrixFree &data, + const std::vector &_varInfoList); + // Nonstandard constructor, used when only one index of "data" should be used, + // use with care! variableContainer(const dealii::MatrixFree &data, const std::vector &_varInfoList, const unsigned int &fixed_index); - /** - * \brief Return the value of the specified scalar field. - * - * \param global_variable_index The global index of the field. - */ + // Methods to get the value/grad/hess in the residual method (this is how the + // user gets these values in equations.h) 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; - /** - * \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. - */ + // Methods to set the value residual and the gradient residual (this is how + // the user sets these values in equations.h) 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); - /** - * \brief Initialize, read DOFs, and set evaulation flags for each variable. - * - * \param src The source vector. - * \param cell The cell where this is done. - */ + // Initialize, read DOFs, and set evaulation flags for each variable 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); - /** - * \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. - */ + // Only initialize the FEEvaluation object for each variable (used for + // post-processing) void reinit(unsigned int cell); - /** - * \brief Integrate the residuals and distribute from local to global for the RHS. - * - * \param dst The destination vector. - */ + // Integrate the residuals and distribute from local to global 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); - /** - * \brief The quadrature point index. - */ + // 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 unsigned int q_point; - /** - * \brief Return the number of quadrature points per cell. - */ - [[nodiscard]] unsigned int + 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; - /** - * \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). - */ + // Object containing some information about each variable (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; - /** - * \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. - */ + // The number of variables unsigned int num_var; }; diff --git a/src/EquationDependencyParser/EquationDependencyParser.cc b/src/EquationDependencyParser/EquationDependencyParser.cc index 906b28e57..ba199ba84 100644 --- a/src/EquationDependencyParser/EquationDependencyParser.cc +++ b/src/EquationDependencyParser/EquationDependencyParser.cc @@ -29,8 +29,6 @@ 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); @@ -57,7 +55,6 @@ 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); @@ -73,7 +70,6 @@ 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); @@ -90,7 +86,6 @@ 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); @@ -100,7 +95,6 @@ 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); @@ -118,7 +112,6 @@ 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) { @@ -172,21 +165,6 @@ 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; @@ -219,36 +197,6 @@ 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 = @@ -272,7 +220,6 @@ 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) @@ -343,29 +290,10 @@ 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) { @@ -393,45 +321,6 @@ 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 41115a74e..29fe2aa0d 100644 --- a/src/matrixfree/computeLHS.cc +++ b/src/matrixfree/computeLHS.cc @@ -57,11 +57,8 @@ MatrixFreePDE::getLHS( const vectorType &src, const std::pair &cell_range) const { - variableContainer> variable_list( - data, - userInputs.varInfoListLHS, - userInputs.varChangeInfoListLHS, - userInputs.varInfoList_old_LHS); + variableContainer> + variable_list(data, userInputs.varInfoListLHS, userInputs.varChangeInfoListLHS); // loop over cells for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) @@ -69,7 +66,6 @@ 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 4267abecd..8b0cd96c6 100644 --- a/src/matrixfree/computeRHS.cc +++ b/src/matrixfree/computeRHS.cc @@ -32,15 +32,13 @@ MatrixFreePDE::getExplicitRHS( { variableContainer> variable_list( data, - userInputs.varInfoListExplicitRHS, - userInputs.varInfoList_old_RHS); + userInputs.varInfoListExplicitRHS); // 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(); @@ -91,15 +89,13 @@ MatrixFreePDE::getNonexplicitRHS( { variableContainer> variable_list( data, - userInputs.varInfoListNonexplicitRHS, - userInputs.varInfoList_old_RHS); + userInputs.varInfoListNonexplicitRHS); // 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 e07217b7a..58d4342d9 100644 --- a/src/matrixfree/init.cc +++ b/src/matrixfree/init.cc @@ -95,6 +95,10 @@ 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) { @@ -260,29 +264,18 @@ MatrixFreePDE::init() pcout << "initializing parallel::distributed residual and solution vectors\n"; for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) { - auto *U = new vectorType; - auto *R = new vectorType; + vectorType *U, *R; + 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 158f7a4a4..bd39ab70a 100644 --- a/src/matrixfree/invM.cc +++ b/src/matrixfree/invM.cc @@ -44,6 +44,20 @@ 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 f057c1c01..06fa05f05 100644 --- a/src/matrixfree/postprocessor.cc +++ b/src/matrixfree/postprocessor.cc @@ -39,8 +39,9 @@ MatrixFreePDE::getPostProcessedFields( const std::pair &cell_range) { // initialize FEEvaulation objects - variableContainer> - variable_list(data, userInputs.pp_baseVarInfoList, userInputs.varInfoList_old_RHS); + variableContainer> variable_list( + data, + userInputs.pp_baseVarInfoList); variableContainer> pp_variable_list(data, userInputs.pp_varInfoList, 0); @@ -49,7 +50,6 @@ 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 dff10d4f0..031598846 100644 --- a/src/matrixfree/reinit.cc +++ b/src/matrixfree/reinit.cc @@ -125,20 +125,13 @@ MatrixFreePDE::reinit() pcout << "initializing parallel::distributed residual and solution vectors\n"; for (unsigned int fieldIndex = 0; fieldIndex < fields.size(); fieldIndex++) { - vectorType *U = solutionSet.at(fieldIndex); + vectorType *U; + + 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 37778098f..3199cf0ff 100644 --- a/src/matrixfree/solveIncrement.cc +++ b/src/matrixfree/solveIncrement.cc @@ -14,9 +14,6 @@ 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) { @@ -550,28 +547,4 @@ 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 e85491088..942936973 100644 --- a/src/userInputParameters/loadVariableAttributes.cc +++ b/src/userInputParameters/loadVariableAttributes.cc @@ -69,7 +69,7 @@ userInputParameters::loadVariableAttributes( varInfoListExplicitRHS.push_back(varInfo); } - // Load the variable information for calculating the RHS for nonexplicit equations + // Load 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,41 +96,10 @@ 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 for nonexplicit equations + // Load variable information for calculating the LHS num_var_LHS = 0; for (unsigned int i = 0; i < number_of_variables; i++) { @@ -158,12 +127,9 @@ 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++) { @@ -172,6 +138,7 @@ 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]; @@ -184,36 +151,6 @@ 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 682c315d8..35ca18ffe 100644 --- a/src/variableContainer/variableContainer.cc +++ b/src/variableContainer/variableContainer.cc @@ -4,18 +4,15 @@ template variableContainer::variableContainer( const dealii::MatrixFree &data, const std::vector &_varInfoList, - const std::vector &_varChangeInfoList, - const std::vector &_varOldInfoList) + const std::vector &_varChangeInfoList) : 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) { @@ -48,69 +45,34 @@ 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 &_varOldInfoList) + const std::vector &_varInfoList) : 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_old_info = varOldInfoList[i]; + const auto &var_info = varInfoList[i]; - if (var_info.var_needed) + if (!var_info.var_needed) { - 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)); - } - else - { - vector_vars_map.emplace(var_index, - std::make_unique(data, i)); - } + continue; } - if (var_old_info.var_needed) - { - const unsigned int var_index = var_old_info.global_var_index; + const unsigned int var_index = var_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)); - } + 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)); } } } @@ -267,40 +229,6 @@ 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) @@ -379,20 +307,26 @@ 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 { - 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); + 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(); + } } template @@ -400,15 +334,19 @@ dealii::Tensor<1, dim, T> variableContainer::get_scalar_gradient( unsigned int global_variable_index) const { - 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); + 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(); + } } template @@ -416,15 +354,19 @@ dealii::Tensor<2, dim, T> variableContainer::get_scalar_hessian( unsigned int global_variable_index) const { - 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); + 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(); + } } template @@ -432,15 +374,19 @@ dealii::Tensor<1, dim, T> variableContainer::get_vector_value( unsigned int global_variable_index) const { - 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); + 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(); + } } template @@ -448,15 +394,19 @@ dealii::Tensor<2, dim, T> variableContainer::get_vector_gradient( unsigned int global_variable_index) const { - 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); + 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(); + } } template @@ -464,31 +414,41 @@ dealii::Tensor<3, dim, T> variableContainer::get_vector_hessian( unsigned int global_variable_index) const { - 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); + 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(); + } } +// 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 { - 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); + 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(); + } } template @@ -496,15 +456,19 @@ dealii::Tensor<1, dim, T> variableContainer::get_change_in_scalar_gradient( unsigned int global_variable_index) const { - 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); + 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(); + } } template @@ -512,15 +476,19 @@ dealii::Tensor<2, dim, T> variableContainer::get_change_in_scalar_hessian( unsigned int global_variable_index) const { - 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); + 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(); + } } template @@ -528,15 +496,19 @@ dealii::Tensor<1, dim, T> variableContainer::get_change_in_vector_value( unsigned int global_variable_index) const { - 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); + 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(); + } } template @@ -544,15 +516,19 @@ dealii::Tensor<2, dim, T> variableContainer::get_change_in_vector_gradient( unsigned int global_variable_index) const { - 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); + 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(); + } } template @@ -560,111 +536,19 @@ dealii::Tensor<3, dim, T> variableContainer::get_change_in_vector_hessian( unsigned int global_variable_index) const { - 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); + 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(); + } } // The methods to set the residual terms