deal.II version GIT relicensing-2202-g5ed228918a 2024-12-01 17:10:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The 'Generalized Swift-Hohenberg Equation Solve' code gallery program

Default constructor, initializes all variables and objects with default values.

This program was contributed by Samuel Scheuerman <sam.scheuerman.051@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Pictures from this code gallery program

Annotated version of README.md

Introduction

This program is used to solve the generalized Swift-Hohenberg equation

\begin{align*} \frac{\partial u}{\partial t} = ru - (k_c + \Delta)^2 u + g_1 u^2 - u^3 \end{align*}

where \(k_c\) is the wave number, \(r\) is some fixed constant, and \(g_1\) is a parameter which determines the behavior of the solutions. Note that the equation is simply called the Swift-Hoheneberg equation if \(g_1 = 0\). For this solver, we chose \(k_c = 1\) and \(r = 0.3\). Choosing \(k_c\) to be 1 will mean that our solutions have a pattern wavelength of \(2\pi\). We choose \(r = 0.3\) because solutions are reasonably well behaved for small values of \(r\) and \(g_1\), but there are interesting behaviors that occur when \(g_1\) is smaller or larger than \(r\) in magnitude, so this allows us room to vary \(g_1\) and explore these behavior. Additionally, we choose \(r = 0.3\) because this matches the parameters used by Gurevich in [1]. We chose our parameters to match so that we could compare the output of our program to the results presented in [1], which was useful for validating that our code was functioning properly during the development process. To summarize, this code solves:

\begin{align*} \frac{\partial u}{\partial t} = 0.3u - (1 + \Delta)^2 u + g_1 u^2 - u^3 \end{align*}

Discretization and Solving the Bilaplacian

The equation has two aspects which are challenging to solve: the nonlinear terms \(g_1u^2 - u^3\) and the Bilaplacian operator \((1 + \Delta)^2\), which introduces \(4^{th}\) derivatives. To deal with the Bilaplacian, we introduce a variable \(v\) and construct a system of PDEs:

\begin{align*} \frac{\partial u}{\partial t} &= 0.3u - (1 + \Delta) v + g_1 u^2 - u^3\\ (1 + \Delta)u &= v \end{align*}

We can solve these two equations simultaneously by treating our finite elements as vector valued, and interpreting our system of equations as a single vector-valued PDE. We can handle the nonlinear terms by treating them fully explicitly. If we discretize in time and rearrange terms, our system of equations becomes

\begin{align*} (1 - kr)U_n + k(1 + \Delta)V_n &= U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\ (1 + \Delta)U_n - V_n &= 0 \end{align*}

where \(k\) is the discrete timestep, \(U_n\) and \(V_n\) are the solutions for \(u\) and \(v\) at the current timestep, and \(U_{n-1}\) and \(V_{n-1}\) are the solutions for \(u\) and \(v\) at the previous timestep. We then reframe this system as a vector valued problem

\begin{align*} \left(\begin{matrix} 1 - kr & k(1 + \Delta)\\ 1 + \Delta & -1 \end{matrix}\right) \left(\begin{matrix} U_n\\ V_n \end{matrix}\right) &= \left(\begin{matrix} U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\ 0 \end{matrix}\right) \end{align*}

As usual, we multiply each side of the equation by a test function

\begin{align*} \overrightarrow\varphi_i = \left(\begin{matrix} \phi_i\\ \psi_i \end{matrix}\right) \end{align*}

and then integrate over the domain \(\Omega\) to get the equation

\begin{align*} \int_\Omega \left(\begin{matrix} \phi_i\\ \psi_i \end{matrix}\right)\cdot\left(\begin{matrix} 1 - kr & k(1 + \Delta)\\ 1 + \Delta & -1 \end{matrix}\right) \left(\begin{matrix} U_n\\ V_n \end{matrix}\right) &= \int_\Omega \left(\begin{matrix} \phi_i\\ \psi_i \end{matrix}\right)\cdot\left(\begin{matrix} U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\ 0 \end{matrix}\right)\\ \end{align*}

We can expand our solution vector in this basis

\begin{align*} \int_\Omega \sum_j u_j\left(\begin{matrix} \phi_i\\ \psi_i \end{matrix}\right)\cdot\left(\begin{matrix} 1 - kr & k(1 + \Delta)\\ 1 + \Delta & -1 \end{matrix}\right) \left(\begin{matrix} \phi_j\\ \psi_j \end{matrix}\right) &= \int_\Omega\left(\begin{matrix} \phi_i\\ \psi_i \end{matrix}\right)\cdot\left(\begin{matrix} U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\ 0 \end{matrix}\right) \end{align*}

and finally expand out the matrix multiplication and dot products, then apply the divergence theorem to obtain a single equation:

\begin{align*} \sum_j u_j \int_\Omega[(1 - kr)\phi_i\phi_j + k\phi_i\psi_j - k\nabla\phi_i\nabla\psi_j + \psi_i\phi_j - \nabla\psi_i\nabla\psi_j - \psi_i\psi_j] &= \int_\Omega\phi_i(U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3) \end{align*}

This last equation represents matrix multiplication of the solution vector by the \(i^{th}\) row of the system matrix, and the left hand side without the summation or \(u_j\) term is what we use to compute the \((i, j)^{th}\) entry of the system matrix.

Boundary Conditions and Choosing a Suitable Domain

This code implements both zero Dirichlet and zero Neumann boundary conditions. Both of these conditions affect the patterns that form. To minimize this effect, we use Neumann boundary conditions and we choose the boundary to be some multiple of the wave number. For example, this code chooses the square mesh to have a side length of \(6\cdot 2\pi\). For all other domains used, we chose a domain size with a similar area to that of the square. For instance, we solve on a torus with outer radius 9 and inner radius 4 because this results in exactly the same area as the square. Note that this is not strictly necessary for the code to function, but does make it easier to compare results between different geometries.

Initial Conditions and Parameters

The code implements two main types of initial conditions: random initial conditions, and creating a small initial hot spot. The SH equation is interesting because it describes pattern formation and self-organization, so choosing random initial conditions allows for this to be observed. Note that the results shown below were all run with the initial seed 314, which was arbitrarily chosen. Setting a fixed seed is useful for comparing pattern formation with different choices of parameters in the SH equation.

The hot spot initial condition is useful for the opposite reason: it is very simple, but it lets us see what happens to a single pattern "wave" as it propagates along our surface. This is particularly useful in distinguishing the effect of curvature and geometry on pattern propagation.

As previously mentioned, we chose \(k_c = 1\) and \(r = 0.3\) for this code. We then let \(g_1\) be the parameter that we change to vary the patterns formed. On the plane, increasing the value of \(g_1\) allows for the formation of hexagonal grids rather than just ripples. Varying \(g_1\) does something similar to patterns on a curved manifold, though with notably different effects in some cases. Increasing \(g_1\) also causes the solution to grow larger in magnitude at certain points.

Checking Convergence

We checked the convergence of this code using 3 tests: we confirmed that a constant initial condition remained constant and converged to a solution that was verified using an ordinary differential equation, we checked that solutions on the square converged across mesh refinements, and we checked that solutions converged over refinements of the timestep on the finest mesh.

Below are the results of several runs of constant initial conditions

image

We also validated that given a fixed random start on a very fine mesh, refining the timestep resulted in the same final solution. The initial condition for each is shown above, While the final solutions are shown in the matrix below. Note that the timestep begins at 1/25 and the denominator increases by 25 across each row, to a max of 1/200 in the bottom right:

image

We validated that solutions converged across mesh refinement by defining psuedorandom functions \(\displaystyle f(x) = \sum_{n=1}^{10} C_n \sin\left(\frac{16x}{3i}\right)\) and \(\displaystyle g(y) = \sum_{n=1}^{10} D_n \sin\left(\frac{16y}{3i}\right)\), where \(C_i\) and \(D_i\) are randomly chosen in the range \((-\sqrt{r}, \sqrt{r})\). The overall pseudorandom function is \(h(x) = f(x)g(y)\). Note that the period of the sine waves was chosen so that the smallest wave could be resolved by a mesh refinement of 7 or higher. The following matrix shows the initial and final solution ranging from a refinement of 0 to a refinement of 7:

image

image

Results

We can see the effects of varying the \(g_1\) parameter and the effects of curvature using the hot spot initial condition. On the plane, an initial hot spot creates one ripple wave, which breaks into discrete pieces as \(g_1\) is increased. In the matrix below, \(g_1\) is increased by 0.2 starting from 0 to a maximum value of 1.4. Note that each final solution is at 100 time units:

image

On the cylinder, the front looks similar to the square, but the back has an overlapping wave pattern:

image

On the sphere, the hot spot generates a single wave. Note that this may be due to the fact that our sphere has a surface area proportional to the period of our pattern wave.

image

On the torus, the pattern propagates similar to the cylinder, with some minor imperfections

image

But on the back side of the torus, we see wave overlapping and spot patterns forming

image

On shapes with stranger curvature, we can see that the pattern wave has a tendency to break apart when crossing lines of curvature. This shape was made by warping the boundary of a cylinder by a cosine wave, and is equivalent to the surface of revolution bounded by \(1 + 0.5\cos(\frac{\pi}{10}x)\)

image

image

Finally, here is a small selection of random initial conditions and the patterns that form. Each image sequence was taken at times 0, 10, 25, 50, and 100:

image

image

image

References

[1] Svetlana Gurevich. Chapter 4: Swift-Hohenberg Equation. url : https://www.uni-muenster.de/Physik.TP/archive/fileadmin/lehre/NumMethoden/SoSe10/Skript/SH.pdf.

Annotated version of Generalized-Swift-Hohenberg-Solver.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2023 by Samuel Scheuerman
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  *
  * Author: Samuel Scheuerman, Colorado State University, 2023
  */
 
 
  #include <deal.II/base/utilities.h>
  #include <deal.II/base/quadrature_lib.h>
  #include <deal.II/base/function.h>
  #include <deal.II/base/logstream.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/full_matrix.h>
  #include <deal.II/lac/dynamic_sparsity_pattern.h>
  #include <deal.II/lac/sparse_matrix.h>
  #include <deal.II/lac/solver_cg.h>
  #include <deal.II/lac/precondition.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/grid/tria.h>
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/grid/grid_refinement.h>
  #include <deal.II/grid/grid_out.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/fe/fe_system.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/numerics/solution_transfer.h>
  #include <deal.II/numerics/matrix_tools.h>
  #include <deal.II/lac/sparse_direct.h>
  #include <deal.II/base/timer.h>
 
  #include <deal.II/grid/manifold_lib.h>
  #include <deal.II/grid/grid_tools.h>
 
  #include <boost/math/special_functions/ellint_1.hpp>
 
  #include <fstream>
  #include <iostream>
  #include <random>
 
  namespace SwiftHohenbergSolver
  {
  using namespace dealii;
 
 
 
 
  enum MeshType {HYPERCUBE, CYLINDER, SPHERE, TORUS, SINUSOID};
 
 
 
  enum InitialConditionType {HOTSPOT, PSUEDORANDOM, RANDOM};
 
 
 
 
 
  template<int spacedim>
  Point<spacedim> transform_function(const Point<spacedim>&p)
  {
Definition point.h:111

Currently this only works for a 3-dimensional embedding space because we are explicitly referencing the x, y, and z coordinates

  Assert(spacedim == 3, ExcNotImplemented());
 
#define Assert(cond, exc)

Returns a point where the x-coordinate is unchanged but the y and z coordinates are adjusted by a cos wave of period 20, amplitude .5, and vertical shift 1

  return Point<spacedim>(p(0), p(1)*(1 + .5*std::cos((3.14159/10)*p(0))), p(2)*(1 + .5*std::cos((3.14159/10)*p(0))));
  }
 
 
 
 
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  class SHEquation
  {
  public:
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)

/

  SHEquation();
 
 
 
  SHEquation(const unsigned int degree
  , double time_step_denominator
  , unsigned int ref_num
  , double r_constant = 0.5
  , double g1_constant = 0.5
  , std::string output_file_name = "solution-"
  , double end_time = 0.5);
  void run();
 
  private:
  void setup_system();
  void solve_time_step();
  void output_results() const;
 
  void make_grid();
 
 
  void make_cylinder();

/ Uses the same process as creating a cylinder, but then also warps the boundary of the cylinder by the function (1 + 0.5*cos(pi*x/10))

  void make_sinusoid();

/ Generates a spherical mesh of radius 6*pi using GridGenerator and refines it refinement_number times.

  void make_sphere();

/ Generates a torus mesh with inner radius 4 and outer radius 9 using GridGenerator and refines it refinement_number times.

  void make_torus();

/ Generates a hypercube mesh with sidelenth 12*pi using GridGenerator and refines it refinement_number times.

  void make_hypercube();
 
 

/ The degree of finite element to be used, default 1

  const unsigned int degree;
 

/ Object holding the mesh

 
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

/ Object which understands which finite elements are at each node

/ Describes the sparsity of the system matrix, allows for more efficient storage

  SparsityPattern sparsity_pattern;
 

/ Object holding the system matrix, stored as a sparse matrix

  SparseMatrix<double> system_matrix;
 
 
  Vector<double> solution;

/ Stores the solution from the previous timestep. Used to compute non-linear terms

  Vector<double> old_solution;
 
  Vector<double> system_rhs;
 

/ Stores the current time, in the units of the problem

  double time;

/ The amount time is increased each iteration/ the denominator of the discretized time derivative

  double time_step;

/ Counts the number of iterations that have elapsed

  unsigned int timestep_number;

/ Used to compute the time_step: time_step = 1/timestep_denominator

  unsigned int timestep_denominator;

/ Determines how much to globally refine each mesh

  unsigned int refinement_number;
 

/ Coefficient of the linear term in the SH equation. This is often taken to be constant and g_1 allowed to vary

  const double r;

/ Coefficient of the quadratic term in the SH equation. Determines whether hexagonal lattices can form

  const double g1;

/ A control parameter for the cubic term. Can be useful for testing, in this code we let k=1 in all cases

  const double k;
 

/ Name used to create output file. Should not include extension

  const std::string output_file_name;
 

/ Determines when the solver terminates, endtime of ~100 are useful to see equilibrium results

  const double end_time;
  };
 
 
 
  template <int spacedim>
  class BoundaryValues : public Function<spacedim>
  {
  public:
  BoundaryValues()
  : Function<spacedim>(2)
  {}
 
  virtual double value(const Point<spacedim> & p,
  const unsigned int component = 0) const override;
  };
 
 
 
 
  template <int spacedim>
  double BoundaryValues<spacedim>::value(const Point<spacedim> & p,
  const unsigned int component) const
  {
  (void)component;
  AssertIndexRange(component, 2);
 
  return 0.;
  }
 
 
  template<int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  class InitialCondition : public Function<spacedim>
  {
  private:
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
#define AssertIndexRange(index, range)

/ The value of the parameter r, used to determine a bound for the magnitude of the initial conditions

  const double r;

/ A center point, used to determine the location of the hot spot for the HotSpot initial condition

  Point<spacedim> center;

/ Radius of the hot spot

  double radius;

/ Stores the randomly generated coefficients for planar sine waves along the x-axis, used for psuedorandom initial conditions

  double x_sin_coefficients[10];

/ Stores the randomly generated coefficients for planar sine waves along the y-axis, used for psuedorandom initial conditions

  double y_sin_coefficients[10];
 
  public:
 
  InitialCondition()
  : Function<spacedim>(2),
  r(0.5),
  radius(.5)
  {
  for(int i = 0; i < 10; ++i){
  x_sin_coefficients[i] = 2*std::sqrt(r)*(std::rand()%1001)/1000 - std::sqrt(r);
  y_sin_coefficients[i] = 2*std::sqrt(r)*(std::rand()%1001)/1000 - std::sqrt(r);
  }
  }
 
 
  InitialCondition(const double r,
  const double radius)
  : Function<spacedim>(2),
  r(r),
  radius(radius)
  {
  for(int i = 0; i < 10; ++i){
  x_sin_coefficients[i] = 2*std::sqrt(r)*(std::rand()%1001)/1000 - std::sqrt(r);
  y_sin_coefficients[i] = 2*std::sqrt(r)*(std::rand()%1001)/1000 - std::sqrt(r);
  }
  }
 
 
  virtual double value(const Point<spacedim> &p, const unsigned int component) const override;
  };
 
 
  template <>
  double InitialCondition<2, HYPERCUBE, HOTSPOT>::value(
  const Point<2> &p,
  const unsigned int component) const
  {
  if(component == 0){
  if(p.square() <= radius){
  return std::sqrt(r);
  }
  else{
  return -std::sqrt(r);
  }
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, CYLINDER, HOTSPOT>::value(
  const Point<3> &p,
  const unsigned int component) const
  {
  if(component == 0){
  const Point<3> center(0, 0, 6);
  const Point<3> compare(p - center);
  if(compare.square() <= radius){
  return std::sqrt(r);
  }
  else{
  return -std::sqrt(r);
  }
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, SPHERE, HOTSPOT>::value(
  const Point<3> &p,
  const unsigned int component) const
  {
  if(component == 0){
  const Point<3> center(18.41988074, 0, 0);
  const Point<3> compare(p - center);
  if(compare.square() <= radius){
  return std::sqrt(r);
  }
  else{
  return -std::sqrt(r);
  }
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, TORUS, HOTSPOT>::value(
  const Point<3> &p,
  const unsigned int component) const
  {
  if(component == 0){
  const Point<3> center(13., 0, 0);
  const Point<3> compare(p - center);
  if(compare.square() <= radius){
  return std::sqrt(r);
  }
  else{
  return -std::sqrt(r);
  }
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, SINUSOID, HOTSPOT>::value(
  const Point<3> &p,
  const unsigned int component) const
  {
  if(component == 0){
  const Point<3> center(0, 0, 9.);
  const Point<3> compare(p - center);
  if(compare.square() <= radius){
  return std::sqrt(r);
  }
  else{
  return -std::sqrt(r);
  }
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<2, HYPERCUBE, PSUEDORANDOM>::value(
  const Point<2> &p,
  const unsigned int component) const
  {
  if(component == 0){
  double x_val = 0;
  double y_val = 0;
  for(int i=0; i < 10; ++i){
  x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
  y_val += y_sin_coefficients[i]*std::sin(2*3.141592653*p(1)/((i+1)*1.178097245));
  }
 
  return x_val*y_val;
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, CYLINDER, PSUEDORANDOM>::value(
  const Point<3> &p,
  const unsigned int component) const
  {
  if(component == 0){
  double x_val = 0;
  double w_val = 0;
  double width = ((std::atan2(p(1),p(2)) - 3.1415926)/3.1415926)*18.84955592;
  for(int i=0; i < 10; ++i){
  x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
  w_val += y_sin_coefficients[i]*std::sin(2*3.141592653*width/((i+1)*1.178097245));
  }
 
  return x_val*w_val;
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, SPHERE, PSUEDORANDOM>::value(
  const Point<3> &p,
  const unsigned int component) const
  {
  if(component == 0){
  double x_val = 0;
  double y_val = 0;
  for(int i=0; i < 10; ++i){
  x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
  y_val += y_sin_coefficients[i]*std::sin(2*3.141592653*p(1)/((i+1)*1.178097245));
  }
 
  return x_val*y_val;
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, TORUS, PSUEDORANDOM>::value(
  const Point<3> &p,
  const unsigned int component) const
  {
  if(component == 0){
  double x_val = 0;
  double z_val = 0;
  for(int i=0; i < 10; ++i){
  x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
  z_val += y_sin_coefficients[i]*std::sin(2*3.141592653*p(2)/((i+1)*1.178097245));
  }
 
  return x_val*z_val;
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, SINUSOID, PSUEDORANDOM>::value(
  const Point<3> &p,
  const unsigned int component) const
  {
  if(component == 0){
  double x_val = 0;
  double w_val = 0;
  double width = ((std::atan2(p(1),p(2)) - 3.1415926)/3.1415926)*18.84955592;
  for(int i=0; i < 10; ++i){
  x_val += x_sin_coefficients[i]*std::sin(2*3.141592653*p(0)/((i+1)*1.178097245));
  w_val += y_sin_coefficients[i]*std::sin(2*3.141592653*width/((i+1)*1.178097245));
  }
 
  return x_val*w_val;
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<2, HYPERCUBE, RANDOM>::value(
  const Point<2> &/*p*/,
  const unsigned int component) const
  {
  if(component == 0){
  return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, CYLINDER, RANDOM>::value(
  const Point<3> &/*p*/,
  const unsigned int component) const
  {
  if(component == 0){
  return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, SPHERE, RANDOM>::value(
  const Point<3> &/*p*/,
  const unsigned int component) const
  {
  if(component == 0){
  return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, TORUS, RANDOM>::value(
  const Point<3> &/*p*/,
  const unsigned int component) const
  {
  if(component == 0){
  return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
  }
  else{
  return 1e18;
  }
  }
 
 
  template <>
  double InitialCondition<3, SINUSOID, RANDOM>::value(
  const Point<3> &/*p*/,
  const unsigned int component) const
  {
  if(component == 0){
  return 2*std::sqrt(r)*(std::rand()%10001)/10000 - std::sqrt(r);
  }
  else{
  return 1e18;
  }
  }
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  SHEquation<dim, spacedim, MESH, ICTYPE>::SHEquation()
  : degree(1)
  , fe(FE_Q<dim, spacedim>(degree), 2)
  , dof_handler(triangulation)
  , time_step(1. / 1500)
  , timestep_denominator(1500)
  , refinement_number(4)
  , r(0.5)
  , g1(0.5)
  , k(1.)
  , output_file_name("solution-")
  , end_time(0.5)
  {}
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  SHEquation<dim, spacedim, MESH, ICTYPE>::SHEquation(const unsigned int degree,
  double time_step_denominator,
  unsigned int ref_num,
  double r_constant,
  double g1_constant,
  std::string output_file_name,
  double end_time)
  : degree(degree)
  , fe(FE_Q<dim, spacedim>(degree), 2)
  , dof_handler(triangulation)
  , time_step(1. / time_step_denominator)
  , timestep_denominator(time_step_denominator)
  , refinement_number(ref_num)
  , r(r_constant)
  , g1(g1_constant)
  , k(1.)
  , output_file_name(output_file_name)
  , end_time(end_time)
  {}
 
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  void SHEquation<dim, spacedim, MESH, ICTYPE>::setup_system()
  {
  dof_handler.distribute_dofs(fe);
 
Definition fe_q.h:554
constexpr numbers::NumberTraits< Number >::real_type square() const
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)

Counts the DoF's for outputting to consolse

  const std::vector<types::global_dof_index> dofs_per_component =
  const unsigned int n_u = dofs_per_component[0],
  n_v = dofs_per_component[1];
 
  std::cout << "Number of active cells: " << triangulation.n_active_cells()
  << std::endl
  << "Total number of cells: " << triangulation.n_cells()
  << std::endl
  << "Number of degrees of freedom: " << dof_handler.n_dofs()
  << " (" << n_u << '+' << n_v << ')' << std::endl;
 
  DynamicSparsityPattern dsp(dof_handler.n_dofs());
 
  dsp);
  sparsity_pattern.copy_from(dsp);
 
  system_matrix.reinit(sparsity_pattern);
 
  solution.reinit(dof_handler.n_dofs());
  old_solution.reinit(dof_handler.n_dofs());
  system_rhs.reinit(dof_handler.n_dofs());
  }
 
 
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  void SHEquation<dim, spacedim, MESH, ICTYPE>::solve_time_step()
  {
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})

std::cout << "Solving linear system" << std::endl; Timer timer;

  SparseDirectUMFPACK direct_solver;
 
  direct_solver.initialize(system_matrix);
 
  direct_solver.vmult(solution, system_rhs);
 
void initialize(const SparsityPattern &sparsity_pattern)

timer.stop(); std::cout << "done (" << timer.cpu_time() << " s)" << std::endl;

  }
 
 
 
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  void SHEquation<dim, spacedim, MESH, ICTYPE>::output_results() const
  {
  std::vector<std::string> solution_names(1, "u");
  solution_names.emplace_back("v");
  std::vector<DataComponentInterpretation::DataComponentInterpretation>
  interpretation(1,
 
  data_out.add_data_vector(dof_handler,
  solution,
  solution_names,
  interpretation /*,
  DataOut<dim, spacedim>::type_dof_data*/);
 
  data_out.build_patches(degree + 1);
 
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})

Takes the output_file_name string and appends timestep_number with up to three leading 0's

  const std::string filename =
  output_file_name + Utilities::int_to_string(timestep_number, 3) + ".vtu";
 
  std::ofstream output(filename);
  data_out.write_vtu(output);
  }
 
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470

Below are all the different template cases for the make_grid() function

  template <>
  void SHEquation<2, 2, HYPERCUBE, HOTSPOT>::make_grid()
  {
  make_hypercube();
  }
 
  template <>
  void SHEquation<2, 3, CYLINDER, HOTSPOT>::make_grid()
  {
  make_cylinder();
  }
 
  template <>
  void SHEquation<2, 3, SPHERE, HOTSPOT>::make_grid()
  {
  make_sphere();
  }
 
  template <>
  void SHEquation<2, 3, TORUS, HOTSPOT>::make_grid()
  {
  make_torus();
  }
 
  template <>
  void SHEquation<2, 3, SINUSOID, HOTSPOT>::make_grid()
  {
  make_sinusoid();
  }
 
  template <>
  void SHEquation<2, 2, HYPERCUBE, PSUEDORANDOM>::make_grid()
  {
  make_hypercube();
  }
 
  template <>
  void SHEquation<2, 3, CYLINDER, PSUEDORANDOM>::make_grid()
  {
  make_cylinder();
  }
 
  template <>
  void SHEquation<2, 3, SPHERE, PSUEDORANDOM>::make_grid()
  {
  make_sphere();
  }
 
  template <>
  void SHEquation<2, 3, TORUS, PSUEDORANDOM>::make_grid()
  {
  make_torus();
  }
 
  template <>
  void SHEquation<2, 3, SINUSOID, PSUEDORANDOM>::make_grid()
  {
  make_sinusoid();
  }
 
  template <>
  void SHEquation<2, 2, HYPERCUBE, RANDOM>::make_grid()
  {
  make_hypercube();
  }
 
  template <>
  void SHEquation<2, 3, CYLINDER, RANDOM>::make_grid()
  {
  make_cylinder();
  }
 
  template <>
  void SHEquation<2, 3, SPHERE, RANDOM>::make_grid()
  {
  make_sphere();
  }
 
  template <>
  void SHEquation<2, 3, TORUS, RANDOM>::make_grid()
  {
  make_torus();
  }
 
  template <>
  void SHEquation<2, 3, SINUSOID, RANDOM>::make_grid()
  {
  make_sinusoid();
  }
 
 
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  void SHEquation<dim, spacedim, MESH, ICTYPE>::run()
  {
  make_grid();
 
  setup_system();
 

Counts total time elapsed

  time = 0.0;

Counts number of iterations

  timestep_number = 0;
 

Sets the random seed so runs are repeatable, remove for varying random initial conditions

  std::srand(314);
 
  InitialCondition<spacedim, MESH, ICTYPE> initial_conditions(r, 0.5);
 

Applies the initial conditions to the old_solution

  VectorTools::interpolate(dof_handler,
  initial_conditions,
  old_solution);
  solution = old_solution;
 
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={}, const unsigned int level=numbers::invalid_unsigned_int)

Outputs initial solution

  output_results();
 

Sets up the quadrature formula and FEValues object

  const QGauss<dim> quadrature_formula(degree + 2);
 
  FEValues<dim, spacedim> fe_values(fe, quadrature_formula,
 
  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 
  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
  Vector<double> cell_rhs(dofs_per_cell);
 
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.

The vector which stores the global indices that each local index connects to

  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 

Extracts the finite elements associated to u and v

Loops over the cells to create the system matrix. We do this only once because the timestep is constant

  for(const auto &cell : dof_handler.active_cell_iterators()){
  cell_matrix = 0;
  cell_rhs = 0;
 
  fe_values.reinit(cell);
 
  cell->get_dof_indices(local_dof_indices);
 
  for(const unsigned int q_index : fe_values.quadrature_point_indices()){
 
  for(const unsigned int i : fe_values.dof_indices()){

These are the ith finite elements associated to u and v

  const double phi_i_u = fe_values[u].value(i, q_index);
  const Tensor<1, spacedim> grad_phi_i_u = fe_values[u].gradient(i, q_index);
  const double phi_i_v = fe_values[v].value(i, q_index);
  const Tensor<1, spacedim> grad_phi_i_v = fe_values[v].gradient(i, q_index);
 
  for(const unsigned int j : fe_values.dof_indices())
  {

These are the jth finite elements associated to u and v

  const double phi_j_u = fe_values[u].value(j, q_index);
  const Tensor<1, spacedim> grad_phi_j_u = fe_values[u].gradient(j, q_index);
  const double phi_j_v = fe_values[v].value(j, q_index);
  const Tensor<1, spacedim> grad_phi_j_v = fe_values[v].gradient(j, q_index);
 

This formula comes from expanding the PDE system

  cell_matrix(i, j) += (phi_i_u*phi_j_u - time_step*r*phi_i_u*phi_j_u
  + time_step*phi_i_u*phi_j_v - time_step*grad_phi_i_u*grad_phi_j_v
  + phi_i_v*phi_j_u - grad_phi_i_v*grad_phi_j_u
  - phi_i_v*phi_j_v)*fe_values.JxW(q_index);
  }
  }
  }
 

Loops over the dof indices to fill the entries of the system_matrix with the local data

  for(unsigned int i : fe_values.dof_indices()){
  for(unsigned int j : fe_values.dof_indices()){
  system_matrix.add(local_dof_indices[i],
  local_dof_indices[j],
  cell_matrix(i, j));
  }
  }
  }
 

Loops over time, incrementing by timestep, to create the RHS, solve the linear system, then output the result

  while (time <= end_time)
  {

Increments time and timestep_number

  time += time_step;
  ++timestep_number;
 

Outputs to console the number of iterations and current time. Currently outputs once every "second"

  if(timestep_number%timestep_denominator == 0){
  std::cout << "Time step " << timestep_number << " at t=" << time
  << std::endl;
  }
 

Resets the system_rhs vector. THIS IS VERY IMPORTANT TO ENSURE THE SYSTEM IS SOLVED CORRECTLY AT EACH TIMESTEP

  system_rhs = 0;
 

Loops over cells, then quadrature points, then dof indices to construct the RHS

  for(const auto &cell : dof_handler.active_cell_iterators()){

Resets the cell_rhs. THIS IS ALSO VERY IMPORTANT TO ENSURE THE SYSTEM IS SOLVED CORRECTLY

  cell_rhs = 0;
 

Resets the FEValues object to only the current cell

  fe_values.reinit(cell);
 
  cell->get_dof_indices(local_dof_indices);
 

Loop over the quadrature points

  for(const unsigned int q_index : fe_values.quadrature_point_indices()){

Stores the value of the previous solution at the quadrature point

  double Un1 = 0;
 

Loops over the dof indices to get the value of Un1

  for(const unsigned int i : fe_values.dof_indices()){
  Un1 += old_solution(local_dof_indices[i])*fe_values[u].value(i, q_index);
  }
 

Loops over the dof indices, using Un1 to construct the RHS for the current timestep. Un1 is used to account for the nonlinear terms in the SH equation

  for(const unsigned int i : fe_values.dof_indices()){
  cell_rhs(i) += (Un1 + time_step*g1*std::pow(Un1, 2) - time_step*k*std::pow(Un1, 3))
  *fe_values[u].value(i, q_index)*fe_values.JxW(q_index);
  }
  }
 
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)

Loops over the dof indices to store the local data in the global RHS vector

  for(unsigned int i : fe_values.dof_indices()){
  system_rhs(local_dof_indices[i]) += cell_rhs(i);
  }
 
 
  }

This is where Dirichlet conditions are applied, or Neumann conditions if the code is commented out

  /* {
  BoundaryValues<spacedim> boundary_values_function;
  boundary_values_function.set_time(time);
 
  std::map<types::global_dof_index, double> boundary_values;
  VectorTools::interpolate_boundary_values(dof_handler,
  0,
  boundary_values_function,
  boundary_values);
 
  MatrixTools::apply_boundary_values(boundary_values,
  system_matrix,
  solution,
  system_rhs);
  } */
 
  solve_time_step();
 

Outputs the solution at regular intervals, currently once every "second" The SH equation evolves slowly in time, so this saves disk space

  if(timestep_number%timestep_denominator == 0){
  output_results();
  }
 
  old_solution = solution;
  }
  }
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  void SHEquation<dim, spacedim, MESH, ICTYPE>::make_cylinder()
  {

Creates a volumetric cylinder

  Triangulation<3> cylinder;
  GridGenerator::cylinder(cylinder, 6, 18.84955592);
 
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)

Extracts the boundary mesh with ID 0, which happens to be the tube part of the cylinder

 
return_type extract_boundary_mesh(const MeshType< dim, spacedim > &volume_mesh, MeshType< dim - 1, spacedim > &surface_mesh, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())

The manifold information is lost upon boundary extraction. This sets the mesh boundary type to be a cylinder again

  triangulation.set_all_manifold_ids(0);
  triangulation.set_manifold(0, boundary);
 
  triangulation.refine_global(refinement_number);
  }
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  void SHEquation<dim, spacedim, MESH, ICTYPE>::make_sinusoid()
  {

Same process as above

  Triangulation<3> cylinder;
  GridGenerator::cylinder(cylinder, 6, 18.84955592);
 
 
  triangulation.set_all_manifold_ids(0);
  triangulation.set_manifold(0, boundary);
 
  triangulation.refine_global(refinement_number);
 

We warp the mesh after refinement to avoid a jagged mesh. We can't tell the code that the boundary should be a perfect sine wave, so we only warp after the mesh is fine enough to resolve this

  GridTools::transform(transform_function<spacedim>, triangulation);
  }
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  void SHEquation<dim, spacedim, MESH, ICTYPE>::make_sphere()
  {
  triangulation.refine_global(refinement_number);
  }
 
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  void SHEquation<dim, spacedim, MESH, ICTYPE>::make_torus()
  {
  triangulation.refine_global(refinement_number);
  }
  template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
  void SHEquation<dim, spacedim, MESH, ICTYPE>::make_hypercube()
  {
  GridGenerator::hyper_cube(triangulation, -18.84955592, 18.84955592);
  triangulation.refine_global(refinement_number);
  }
  } // namespace SwiftHohenbergSolver
 
 
 
  int main()
  {
  using namespace SwiftHohenbergSolver;
 
void hyper_sphere(Triangulation< spacedim - 1, spacedim > &tria, const Point< spacedim > &center=Point< spacedim >(), const double radius=1.)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void torus(Triangulation< dim, spacedim > &tria, const double centerline_radius, const double inner_radius, const unsigned int n_cells_toroidal=6, const double phi=2.0 *numbers::PI)
void transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)

An array of mesh types. We iterate over this to allow for longer runs without having to stop the code

  MeshType mesh_types[5] = {HYPERCUBE, CYLINDER, SPHERE, TORUS, SINUSOID};

An array of initial condition types. We iterate this as well, for the same reason

  InitialConditionType ic_types[3] = {HOTSPOT, PSUEDORANDOM, RANDOM};
 

Controls how long the code runs

  const double end_time = 100.;
 

The number of times we refine the hypercube mesh

  const unsigned int ref_num = 6;
 

The timestep will be 1/timestep_denominator

  const unsigned int timestep_denominator = 25;
 

Loops over mesh types, then initial condition types, then loops over values of g_1

  for(const auto MESH : mesh_types){
  for(const auto ICTYPE: ic_types){
  for(int i = 0; i < 8; ++i){

The value of g_1 passed to the solver object

  const double g_constant = 0.2*i;
 

Used to distinguish the start of each run

  std::cout<< std::endl << std::endl;
 
  try{

Switch statement that determines what template parameters are used by the solver object. Template parameters must be known at compile time, so we cannot pass this as a variable unfortunately. In each case, we create a filename string (named appropriately for the particular case), output to the console what we are running, create the solver object, and call run(). Note that for the cylinder, sphere, and sinusoid we decrease the refinement number by 1. This keeps the number of dofs used in these cases comparable to the number of dofs on the 2D hypercube (otherwise the number of dofs is much larger). For the torus, we decrease the refinement number by 2.

  switch (MESH)
  {
  case HYPERCUBE:
  switch (ICTYPE){
  case HOTSPOT:
  {
  std::string filename = "HYPERCUBE-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 2, HYPERCUBE, HOTSPOT> heat_equation_solver(1, timestep_denominator,
  ref_num, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case PSUEDORANDOM:
  {
  std::string filename = "HYPERCUBE-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 2, HYPERCUBE, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case RANDOM:
  {
  std::string filename = "HYPERCUBE-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 2, HYPERCUBE, RANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
  }
  break;
  case CYLINDER:
  switch (ICTYPE){
  case HOTSPOT:
  {
  std::string filename = "CYLINDER-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, CYLINDER, HOTSPOT> heat_equation_solver(1, timestep_denominator,
  ref_num-1, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case PSUEDORANDOM:
  {
  std::string filename = "CYLINDER-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, CYLINDER, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num-1, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case RANDOM:
  {
  std::string filename = "CYLINDER-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, CYLINDER, RANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num-1, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
  }
  break;
  case SPHERE:
  switch (ICTYPE){
  case HOTSPOT:
  {
  std::string filename = "SPHERE-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, SPHERE, HOTSPOT> heat_equation_solver(1, timestep_denominator,
  ref_num-1, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case PSUEDORANDOM:
  {
  std::string filename = "SPHERE-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, SPHERE, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num-1, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case RANDOM:
  {
  std::string filename = "SPHERE-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, SPHERE, RANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num-1, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
  }
  break;
  case TORUS:
  switch (ICTYPE){
  case HOTSPOT:
  {
  std::string filename = "TORUS-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, TORUS, HOTSPOT> heat_equation_solver(1, timestep_denominator,
  ref_num-2, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case PSUEDORANDOM:
  {
  std::string filename = "TORUS-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, TORUS, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num-2, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case RANDOM:
  {
  std::string filename = "TORUS-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, TORUS, RANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num-2, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
  }
  break;
  case SINUSOID:
  switch (ICTYPE){
  case HOTSPOT:
  {
  std::string filename = "SINUSOID-HOTSPOT-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, SINUSOID, HOTSPOT> heat_equation_solver(1, timestep_denominator,
  ref_num-1, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case PSUEDORANDOM:
  {
  std::string filename = "SINUSOID-PSUEDORANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, SINUSOID, PSUEDORANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num-1, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
 
  case RANDOM:
  {
  std::string filename = "SINUSOID-RANDOM-G1-0.2x" + Utilities::int_to_string(i, 1) + "-";
  std::cout << "Running: " << filename << std::endl << std::endl;
 
  SHEquation<2, 3, SINUSOID, RANDOM> heat_equation_solver(1, timestep_denominator,
  ref_num-1, 0.3, g_constant,
  filename, end_time);
  heat_equation_solver.run();
  }
  break;
  }
  break;
  default:
  break;
  }
  }
  catch (std::exception &exc)
  {
  std::cout << "An error occurred" << std::endl;
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Exception on processing: " << std::endl
  << exc.what() << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
 
  return 1;
  }
  catch (...)
  {
  std::cout << "Error occurred, made it past first catch" << std::endl;
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Unknown exception!" << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  }
  }
  }
  return 0;
  }