 Reference documentation for deal.II version GIT edc7d5c3ce 2023-09-25 07:10:03+00:00
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              # 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. 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 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: 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:  # 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: On the cylinder, the front looks similar to the square, but the back has an overlapping wave pattern: 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. On the torus, the pattern propagates similar to the cylinder, with some minor imperfections But on the back side of the torus, we see wave overlapping and spot patterns forming 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)$$  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:   # Annotated version of Generalized-Swift-Hohenberg-Solver.cc

/* ---------------------------------------------------------------------
*   *
*   * Copyright (C) 2023 by the deal.II authors
*   *
*   * This file is part of the deal.II library.
*   *
*   * The deal.II library is free software; you can use it, redistribute
*   * it, and/or modify it under the terms of the GNU Lesser General
*   * version 2.1 of the License, or (at your option) any later version.
*   * The full text of the license can be found in the file LICENSE.md at
*   * the top level directory of deal.II.
*   *
*   * ---------------------------------------------------------------------
*
*   *
*   * Author: Samuel Scheuerman, Colorado State University, 2023
*   */

#include <deal.II/base/utilities.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)
{
const ::Triangulation< dim, spacedim > & tria

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)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()

Retruns 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:

/

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();
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition: work_stream.h:472

/ 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

/ 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 ellapsed

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)
Definition: exceptions.h:1857

/ 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< 3 > center

/ Radius of the hot spot

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

double x_sin_coefficients;

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

double y_sin_coefficients;

public:

InitialCondition()
: Function<spacedim>(2),
r(0.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,
: Function<spacedim>(2),
r(r),
{
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){
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);
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);
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);
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);
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:551
numbers::NumberTraits< Number >::real_type square() const
Expression atan2(const Expression &y, const Expression &x)
::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,
n_v = dofs_per_component;

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={})
Definition: dof_tools.cc:1925

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,

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:471

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>
{
make_grid();

setup_system();

Counts total time ellapsed

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={})

Outputs initial solution

output_results();

Sets up the quadrature formula and FEValues object

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
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double >> &velocity, const double factor=1.)

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 becase 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 double phi_i_v = fe_values[v].value(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 double phi_j_v = fe_values[v].value(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
- 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()){
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);

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

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

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 torus(Triangulation< dim, spacedim > &tria, const double R, const double r, const unsigned int n_cells_toroidal=6, const double phi=2.0 *numbers::PI)
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 transform(const Transformation &transformation, Triangulation< dim, spacedim > &triangulation)

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

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

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

InitialConditionType ic_types = {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 varible 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 occured" << 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 occured, 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;
}