[WIP] Implicit solves & minor optimizations

automatic_tests/precipitateEvolution_pfunction/CMakeLists.txt → applications/advection/CMakeLists.txt

@@ -144,7 +144,7 @@ include_directories(${CMAKE_SOURCE_DIR})
 set(TARGET_SRC "${CMAKE_SOURCE_DIR}/../main.cc")
 
 # Check if there has been updates to main library
-set(dir ${PROJECT_SOURCE_DIR}/../../..)
+set(dir ${PROJECT_SOURCE_DIR}/../..)
 EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir})
 EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir})
 EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir})

automatic_tests/precipitateEvolution_pfunction/ICs_and_BCs.cc → applications/advection/ICs_and_BCs.cc

@@ -1,3 +1,5 @@
+#include <deal.II/base/function_signed_distance.h>
+
 // ===========================================================================
 // FUNCTION FOR INITIAL CONDITIONS
 // ===========================================================================
@@ -16,52 +18,28 @@ customPDE<dim, degree>::setInitialCondition([[maybe_unused]] const Point<dim>  &
   // Use "if" statements to set the initial condition for each variable
   // according to its variable index
 
-  double center[4][3] = {
-    {1.0 / 3.0, 1.0 / 3.0, 0.5},
-    {2.0 / 3.0, 2.0 / 3.0, 0.5},
-    {3.0 / 4.0, 1.0 / 4.0, 0.5},
-    {1.0 / 4.0, 3.0 / 4,   0.5}
-  };
-  double rad[4]         = {userInputs.domain_size[0] / 16.0,
-                           userInputs.domain_size[0] / 16.0,
-                           userInputs.domain_size[0] / 16.0,
-                           userInputs.domain_size[0] / 16.0};
-  double orientation[4] = {1, 1, 2, 3};
-  double dx = userInputs.domain_size[0] / ((double) userInputs.subdivisions[0]) /
-              std::pow(2.0, userInputs.refine_factor);
-  double dist;
-  scalar_IC = 0;
+  // Sphere level-set
+  double             radius = 15.0;
+  dealii::Point<dim> center(20.0, 20.0);
 
-  if (index == 0)
-    {
-      scalar_IC = 0.04;
-    }
+  dealii::Functions::SignedDistance::Sphere<dim> sphere(center, radius);
 
-  for (unsigned int i = 0; i < 4; i++)
-    {
-      dist = 0.0;
-      for (unsigned int dir = 0; dir < dim; dir++)
-        {
-          dist += (p[dir] - center[i][dir] * userInputs.domain_size[dir]) *
-                  (p[dir] - center[i][dir] * userInputs.domain_size[dir]);
-        }
-      dist = std::sqrt(dist);
+  double distance = sphere.value(p);
 
-      if (index == orientation[i])
-        {
-          scalar_IC += 0.5 * (1.0 - std::tanh((dist - rad[i]) / (dx)));
-        }
+  // tanh profile
+  double phi = 0.5 * (1.0 - std::tanh(distance / (std::sqrt(2) * W)));
+
+  if (index == 0)
+    {
+      scalar_IC = phi;
     }
 
-  if (index == 4)
+  else
     {
-      for (unsigned int d = 0; d < dim; d++)
-        {
-          vector_IC(d) = 0.0;
-        }
+      scalar_IC = 0.0;
     }
 
-  // --------------------------------------------------------------------------
+  // ---------------------------------------------------------------------
 }
 
 // ===========================================================================

applications/advection/advection.md

@@ -0,0 +1,38 @@
+## 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}
+$$

automatic_tests/precipitateEvolution_pfunction/customPDE.h → applications/advection/customPDE.h

@@ -2,25 +2,15 @@
 
 using namespace dealii;
 
-// Header files for PFunctions
-typedef VectorizedArray<double> scalarvalueType;
-#include "PLibrary/PLibrary.cc"
-#include "PLibrary/PLibrary.hh"
-
-#include <field_input/pFunction.h>
-
-// Declare the PFunctions to be used
-PFunctions::pFunction pfunct_McV("pfunct_McV"), pfunct_Mn1V("pfunct_Mn1V"),
-  pfunct_Mn2V("pfunct_Mn1V"), pfunct_Mn3V("pfunct_Mn1V"), pfunct_faV("pfunct_faV"),
-  pfunct_fbV("pfunct_fbV");
-
 template <int dim, int degree>
 class customPDE : public MatrixFreePDE<dim, degree>
 {
 public:
+  // Constructor
   customPDE(userInputParameters<dim> _userInputs)
     : MatrixFreePDE<dim, degree>(_userInputs)
     , userInputs(_userInputs) {};
+
   // Function to set the initial conditions (in ICs_and_BCs.h)
   void
   setInitialCondition([[maybe_unused]] const Point<dim>  &p,
@@ -44,7 +34,7 @@ class customPDE : public MatrixFreePDE<dim, degree>
   const userInputParameters<dim> userInputs;
 
   // Function to set the RHS of the governing equations for explicit time
-  // dependent equations (in equations.h)
+  // dependent equations (in equations.cc)
   void
   explicitEquationRHS(
     [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>>
@@ -53,15 +43,15 @@ class customPDE : public MatrixFreePDE<dim, degree>
     [[maybe_unused]] const VectorizedArray<double> element_volume) const override;
 
   // Function to set the RHS of the governing equations for all other equations
-  // (in equations.h)
+  // (in equations.cc)
   void
   nonExplicitEquationRHS(
     [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>>
                                                               &variable_list,
     [[maybe_unused]] const Point<dim, VectorizedArray<double>> q_point_loc,
     [[maybe_unused]] const VectorizedArray<double> element_volume) const override;
 
-  // Function to set the LHS of the governing equations (in equations.h)
+  // Function to set the LHS of the governing equations (in equations.cc)
   void
   equationLHS(
     [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>>
@@ -96,34 +86,9 @@ class customPDE : public MatrixFreePDE<dim, degree>
   // Model constants specific to this subclass
   // ================================================================
 
-  Tensor<2, dim> Kn1 = userInputs.get_model_constant_rank_2_tensor("Kn1");
-  Tensor<2, dim> Kn2 = userInputs.get_model_constant_rank_2_tensor("Kn2");
-  Tensor<2, dim> Kn3 = userInputs.get_model_constant_rank_2_tensor("Kn3");
-  bool           n_dependent_stiffness =
-    userInputs.get_model_constant_bool("n_dependent_stiffness");
-  Tensor<2, dim> sfts_linear1 =
-    userInputs.get_model_constant_rank_2_tensor("sfts_linear1");
-  Tensor<2, dim> sfts_const1 = userInputs.get_model_constant_rank_2_tensor("sfts_const1");
-  Tensor<2, dim> sfts_linear2 =
-    userInputs.get_model_constant_rank_2_tensor("sfts_linear2");
-  Tensor<2, dim> sfts_const2 = userInputs.get_model_constant_rank_2_tensor("sfts_const2");
-  Tensor<2, dim> sfts_linear3 =
-    userInputs.get_model_constant_rank_2_tensor("sfts_linear3");
-  Tensor<2, dim> sfts_const3 = userInputs.get_model_constant_rank_2_tensor("sfts_const3");
-  double         A4          = userInputs.get_model_constant_double("A4");
-  double         A3          = userInputs.get_model_constant_double("A3");
-  double         A2          = userInputs.get_model_constant_double("A2");
-  double         A1          = userInputs.get_model_constant_double("A1");
-  double         A0          = userInputs.get_model_constant_double("A0");
-  double         B2          = userInputs.get_model_constant_double("B2");
-  double         B1          = userInputs.get_model_constant_double("B1");
-  double         B0          = userInputs.get_model_constant_double("B0");
-
-  const static unsigned int  CIJ_tensor_size = 2 * dim - 1 + dim / 3;
-  Tensor<2, CIJ_tensor_size> CIJ_Mg =
-    userInputs.get_model_constant_elasticity_tensor("CIJ_Mg");
-  Tensor<2, CIJ_tensor_size> CIJ_Beta =
-    userInputs.get_model_constant_elasticity_tensor("CIJ_Beta");
+  double                 W = userInputs.get_model_constant_double("W");
+  dealii::Tensor<1, dim> velocity =
+    userInputs.get_model_constant_rank_1_tensor("velocity");
 
   // ================================================================
 };

applications/advection/equations.cc

@@ -0,0 +1,132 @@
+// =================================================================================
+// 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 <int dim, int degree>
+void
+customPDE<dim, degree>::explicitEquationRHS(
+  [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>> &variable_list,
+  [[maybe_unused]] const Point<dim, VectorizedArray<double>>                q_point_loc,
+  [[maybe_unused]] const VectorizedArray<double> 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 <int dim, int degree>
+void
+customPDE<dim, degree>::nonExplicitEquationRHS(
+  [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>> &variable_list,
+  [[maybe_unused]] const Point<dim, VectorizedArray<double>>                q_point_loc,
+  [[maybe_unused]] const VectorizedArray<double> 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 <int dim, int degree>
+void
+customPDE<dim, degree>::equationLHS(
+  [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>> &variable_list,
+  [[maybe_unused]] const Point<dim, VectorizedArray<double>>                q_point_loc,
+  [[maybe_unused]] const VectorizedArray<double> 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);
+}