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
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:
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
#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
{
enum MeshType {HYPERCUBE, CYLINDER, SPHERE, TORUS, SINUSOID};
enum InitialConditionType {HOTSPOT, PSUEDORANDOM, RANDOM};
template<int spacedim>
{
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
}
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))
/ Generates a spherical mesh of radius 6*pi using GridGenerator and refines it refinement_number times.
/ Generates a torus mesh with inner radius 4 and outer radius 9 using GridGenerator and refines it refinement_number times.
/ Generates a hypercube mesh with sidelenth 12*pi using GridGenerator and refines it refinement_number times.
/ 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
/ Object holding the system matrix, stored as a sparse matrix
/ Stores the solution from the previous timestep. Used to compute non-linear terms
/ Stores the current time, in the units of the problem
/ The amount time is increased each iteration/ the denominator of the discretized time derivative
/ 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
/ Coefficient of the quadratic term in the SH equation. Determines whether hexagonal lattices can form
/ A control parameter for the cubic term. Can be useful for testing, in this code we let k=1 in all cases
/ 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()
{}
const unsigned int component = 0) const override;
};
template <int spacedim>
const unsigned int component) const
{
(void)component;
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
/ A center point, used to determine the location of the hot spot for the HotSpot initial condition
/ 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[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()
r(0.5),
radius(.5)
{
for(int i = 0; i < 10; ++i){
}
}
InitialCondition(const double r,
const double radius)
r(r),
radius(radius)
{
for(
int i = 0; i < 10; ++i){
}
}
};
template <>
double InitialCondition<2, HYPERCUBE, HOTSPOT>::value(
const unsigned int component) const
{
if(component == 0){
}
else{
}
}
else{
return 1e18;
}
}
template <>
double InitialCondition<3, CYLINDER, HOTSPOT>::value(
const unsigned int component) const
{
if(component == 0){
if(compare.square() <= radius){
}
else{
}
}
else{
return 1e18;
}
}
template <>
double InitialCondition<3, SPHERE, HOTSPOT>::value(
const unsigned int component) const
{
if(component == 0){
const Point<3> center(18.41988074, 0, 0);
if(compare.square() <= radius){
}
else{
}
}
else{
return 1e18;
}
}
template <>
double InitialCondition<3, TORUS, HOTSPOT>::value(
const unsigned int component) const
{
if(component == 0){
if(compare.square() <= radius){
}
else{
}
}
else{
return 1e18;
}
}
template <>
double InitialCondition<3, SINUSOID, HOTSPOT>::value(
const unsigned int component) const
{
if(component == 0){
if(compare.square() <= radius){
}
else{
}
}
else{
return 1e18;
}
}
template <>
double InitialCondition<2, HYPERCUBE, PSUEDORANDOM>::value(
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 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 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 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 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 unsigned int component) const
{
if(component == 0){
}
else{
return 1e18;
}
}
template <>
double InitialCondition<3, CYLINDER, RANDOM>::value(
const unsigned int component) const
{
if(component == 0){
}
else{
return 1e18;
}
}
template <>
double InitialCondition<3, SPHERE, RANDOM>::value(
const unsigned int component) const
{
if(component == 0){
}
else{
return 1e18;
}
}
template <>
double InitialCondition<3, TORUS, RANDOM>::value(
const unsigned int component) const
{
if(component == 0){
}
else{
return 1e18;
}
}
template <>
double InitialCondition<3, SINUSOID, RANDOM>::value(
const unsigned int component) const
{
if(component == 0){
}
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)
, 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)
, 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);
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
<< std::endl
<< "Number of degrees of freedom: " << dof_handler.n_dofs()
<< " (" << n_u << '+' << n_v << ')' << std::endl;
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::cout << "Solving linear system" << std::endl; Timer timer;
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
);
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 =
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)
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
Counts number of iterations
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
initial_conditions,
old_solution);
solution = old_solution;
Outputs initial solution
Sets up the quadrature formula and FEValues object
const unsigned int dofs_per_cell = fe.n_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 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
+ 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
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
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
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
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
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
}
template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
void SHEquation<dim, spacedim, MESH, ICTYPE>::make_sinusoid()
{
Same process as above
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
}
template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
void SHEquation<dim, spacedim, MESH, ICTYPE>::make_sphere()
{
}
template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
void SHEquation<dim, spacedim, MESH, ICTYPE>::make_torus()
{
}
template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
void SHEquation<dim, spacedim, MESH, ICTYPE>::make_hypercube()
{
}
}
int main()
{
using namespace SwiftHohenbergSolver;
void hyper_sphere(Triangulation< spacedim - 1, spacedim > &tria, const Point< spacedim > ¢er=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)
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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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;
}
1.9.8