Merge pull request #268 from landinjm/implicit_equation_types

Support for `IMPLICIT_TIME_DEPENDENT` equation types & Advection application

applications/advection/CMakeLists.txt

@@ -0,0 +1,68 @@
+##
+#  CMake script for the PRISMS-PF applications:
+##
+
+cmake_minimum_required(VERSION 3.1.0)
+project(myapp)
+
+# Find deal.II installation
+find_package(deal.II 9.2.0 QUIET REQUIRED
+	HINTS ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR})
+
+# Check to make sure deal.II is configured with p4est
+if(NOT ${DEAL_II_WITH_P4EST})
+  message(FATAL_ERROR "deal.II must be installed with p4est.")
+endif()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+
+# Set default build type
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE Debug CACHE STRING "Build type" FORCE)
+endif()
+
+# Set up the debug, release, and run targets
+add_custom_target(debug
+  COMMAND +env ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR}
+  COMMAND +env ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
+  COMMENT "Switch CMAKE_BUILD_TYPE to Debug"
+  )
+add_custom_target(release
+  COMMAND +env ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR}
+  COMMAND +env ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
+  COMMENT "Switch CMAKE_BUILD_TYPE to Release"
+  )
+add_custom_target(run COMMAND main
+  COMMENT "Run with ${CMAKE_BUILD_TYPE} configuration"
+  )
+
+# Add postprocess.cc and nucleation.cc if they exist
+if(EXISTS "postprocess.cc")
+	add_definitions(-DPOSTPROCESS_FILE_EXISTS)
+endif()
+if(EXISTS "nucleation.cc")
+	add_definitions(-DNUCLEATION_FILE_EXISTS)
+endif()
+
+# Set location of files
+include_directories(${CMAKE_SOURCE_DIR}/../../include)
+include_directories(${CMAKE_SOURCE_DIR}/../../src)
+include_directories(${CMAKE_SOURCE_DIR})
+
+# SEt the location of the main.cc file
+set(MAIN "${CMAKE_SOURCE_DIR}/../main.cc")
+
+# Add main.cc for executable target
+add_executable(main ${MAIN})
+
+# deal.II linker
+DEAL_II_SETUP_TARGET(main)
+
+# Link libraries for the build type
+if (${CMAKE_BUILD_TYPE} STREQUAL "Release")
+	target_link_libraries(main ${CMAKE_SOURCE_DIR}/../../libprisms_pf.a)
+elseif(${CMAKE_BUILD_TYPE} STREQUAL "DebugRelease")
+  target_link_libraries(main ${CMAKE_SOURCE_DIR}/../../libprisms_pf.a)
+else()
+  target_link_libraries(main ${CMAKE_SOURCE_DIR}/../../libprisms_pf_debug.a)
+endif()

applications/advection/ICs_and_BCs.cc

@@ -0,0 +1,71 @@
+// ===========================================================================
+// FUNCTION FOR INITIAL CONDITIONS
+// ===========================================================================
+#include <deal.II/base/function_signed_distance.h>
+
+template <int dim, int degree>
+void
+customPDE<dim, degree>::setInitialCondition([[maybe_unused]] const Point<dim>  &p,
+                                            [[maybe_unused]] const unsigned int index,
+                                            [[maybe_unused]] double            &scalar_IC,
+                                            [[maybe_unused]] Vector<double>    &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<dim> center(20.0, 20.0);
+
+  dealii::Functions::SignedDistance::Sphere<dim> 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 <int dim, int degree>
+void
+customPDE<dim, degree>::setNonUniformDirichletBCs(
+  [[maybe_unused]] const Point<dim>  &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<double>    &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).
+
+  // -------------------------------------------------------------------------
+}

applications/advection/advection.md

@@ -0,0 +1,40 @@
+## PRISMS-PF: Advection Equation
+The advection equation is given by:
+
+$$
+\begin{align}
+\frac{\partial \eta}{\partial t} = -\textbf{u}\cdot\nabla\eta
+\end{align}
+$$
+
+where $\eta$ is the scalar conserved quantity and $\textbf{u}$ is a constant vector field.
+
+## Time discretization
+
+Considering backward Euler time stepping, we have the time discretized kinetics equation:
+
+$$
+\begin{align}
+\eta^{n+1} = \eta^{n} - \Delta t (\textbf{u}\cdot\nabla\eta^{n+1})
+\end{align}
+$$
+
+PRISMS-PF makes use of the Newton's method for linear and nonlinear iterations, so we can substitute $\eta^{n+1}=\Delta \eta^{n+1}+n_0^{n+1}$.
+
+$$
+\begin{align}
+\Delta \eta^{n+1} + \Delta t (\textbf{u}\cdot\nabla(\Delta \eta^{n+1}))= \eta^{n} - \eta_0^{n+1} - \Delta t (\textbf{u}\cdot\nabla\eta_0^{n+1})
+\end{align}
+$$
+
+## Weak formulation
+
+In the weak formulation, considering an arbitrary variation $w$, the above equation can be expressed as a residual equation:
+
+$$
+\begin{align}
+\int_{\Omega} w [\Delta \eta^{n+1} + \Delta t (\textbf{u}\cdot\nabla(\Delta \eta^{n+1}))] ~dV &= \int_{\Omega} w [\eta^{n} - \eta_0^{n+1} - \Delta t (\textbf{u}\cdot\nabla\eta_0^{n+1})] ~dV
+\end{align}
+$$
+
+

applications/advection/customPDE.h

@@ -0,0 +1,91 @@
+#include "matrixFreePDE.h"
+
+using namespace dealii;
+
+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,
+                      [[maybe_unused]] const unsigned int index,
+                      [[maybe_unused]] double            &scalar_IC,
+                      [[maybe_unused]] Vector<double>    &vector_IC) override;
+
+  // Function to set the non-uniform Dirichlet boundary conditions (in
+  // ICs_and_BCs.h)
+  void
+  setNonUniformDirichletBCs([[maybe_unused]] const Point<dim>  &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<double>    &vector_BC) override;
+
+private:
+#include "typeDefs.h"
+
+  const userInputParameters<dim> userInputs;
+
+  // Function to set the RHS of the governing equations for explicit time
+  // dependent equations (in equations.cc)
+  void
+  explicitEquationRHS(
+    [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>>
+                                                        &variable_list,
+    [[maybe_unused]] Point<dim, VectorizedArray<double>> q_point_loc) const override;
+
+  // Function to set the RHS of the governing equations for all other equations
+  // (in equations.cc)
+  void
+  nonExplicitEquationRHS(
+    [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>>
+                                                        &variable_list,
+    [[maybe_unused]] Point<dim, VectorizedArray<double>> q_point_loc) const override;
+
+  // Function to set the LHS of the governing equations (in equations.cc)
+  void
+  equationLHS(
+    [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>>
+                                                        &variable_list,
+    [[maybe_unused]] Point<dim, VectorizedArray<double>> q_point_loc) const override;
+
+// Function to set postprocessing expressions (in postprocess.h)
+#ifdef POSTPROCESS_FILE_EXISTS
+  void
+  postProcessedFields(
+    [[maybe_unused]] const variableContainer<dim, degree, VectorizedArray<double>>
+      &variable_list,
+    [[maybe_unused]] variableContainer<dim, degree, VectorizedArray<double>>
+                                                              &pp_variable_list,
+    [[maybe_unused]] const Point<dim, VectorizedArray<double>> q_point_loc)
+    const override;
+#endif
+
+// Function to set the nucleation probability (in nucleation.h)
+#ifdef NUCLEATION_FILE_EXISTS
+  double
+  getNucleationProbability([[maybe_unused]] variableValueContainer variable_value,
+                           [[maybe_unused]] double                 dV) const override;
+#endif
+
+  // ================================================================
+  // Methods specific to this subclass
+  // ================================================================
+
+  // ================================================================
+  // Model constants specific to this subclass
+  // ================================================================
+
+  double                 W = userInputs.get_model_constant_double("W");
+  dealii::Tensor<1, dim> velocity =
+    userInputs.get_model_constant_rank_1_tensor("velocity");
+
+  // ================================================================
+};