 |
deal.II version GIT relicensing-6759-gba22b5843e 2026-09-18 10:10:01+00:00
|
Loading...
Searching...
No Matches
|
deal.II version GIT relicensing-6759-gba22b5843e 2026-09-18 10:10:01+00:00
|
This program was contributed by Samuel Scheuerman <[email protected]>.
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):
 |  |  |  |
 |  |  |  |
 |  |  |  |
 |  |
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*}
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.
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.
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.
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:
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:
[1] Svetlana Gurevich. Chapter 4: Swift-Hohenberg Equation. url : https://www.uni-muenster.de/Physik.TP/archive/fileadmin/lehre/NumMethoden/SoSe10/Skript/SH.pdf.
This enum defines the five mesh types implemented in this program and allows the user to pass which mesh is desired to the solver at runtime. This is useful for looping over different meshes.
This enum defines the three initial conditions used by the program. This allows for the solver class to use a template argument to determine the desired initial condition, which is helpful for setting up loops to solve with a variety of different conditions
This function warps points on a cylindrical mesh by cosine wave along the central axis. We use this function to generate the "sinusoid" mesh, which is the surface of revolution bounded by the cosine wave. spacedim is the dimension of the embedding space, which is where the input point lives. p is the input point to be translated. The return is a translated point in the same dimensional space. This is the new point on the mesh.
Currently this only works for a 3-dimensional embedding space because we are explicitly referencing the x, y, and z coordinates
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
This is the class that holds all the important variables for the solver, as well as the important member functions. This class is based off the HeatEquation class from step-26, so we won't go into full detail on all the features, but we will highlight what has been changed for this problem. dim is the intrinsic dimension of the manifold we are solving on. spacedim is the dimension of the embedding space. MESH determines what manifold we are solving on ICTYPE determines what initial condition we use
Default constructor, initializes all variables and objects with default values
Overloaded constructor, allows user to pass values for important constants. degree is the degree of finite element used, time_step_denominator determines what size timestep we use. The timestep is 1/time_step_denominator. ref_num is the number of times the mesh will be globally refined. r_constant is a constant for the linear component, default 0.5, g1_constant is a constant for the quadratic component, default 0.5. output_file_name is self explanatory, default "solution-" end_time determines when the solver stops, default 0.5, should be ~100 to see equilibrium solutions
This function calls a different grid generation function depending on the template argument MESH. Allows the solver object to generate different mesh types based on the template parameter.
Generates a cylindrical mesh with radius 6 and width 6*pi by first creating a volumetric cylinder, extracting the boundary, and redefining the mesh as a cylinder, then refining the mesh refinement_number times
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 sidelength 12*pi using GridGenerator and refines it refinement_number times.
The degree of finite element to be used, default 1
Object holding the mesh
Object describing the finite element vectors at each node (I believe this gives a basis for the finite elements at each node)
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
Vector of coefficients for the solution in the current timestep. We solve for this in each timestep
Stores the solution from the previous timestep. Used to compute non-linear terms
Stores the coefficients of the right hand side function(in terms of the finite elements). Is the RHS for the linear system
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
Used to compute the time_step: time_step = 1/timestep_denominator
Determines how much to globally refine each mesh
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
Determines when the solver terminates, endtime of ~100 are useful to see equilibrium results
The function which applies zero Dirichlet boundary conditions, and is not being used by the solver currently. Leaving the code in case this is ever needed. spacedim is the dimension of the points which the function takes as input
Returns 0 for all points. This is the output for the boundary spacedim is the dimension of points that are input, p is the input point, component determines whether we are solving for u or v, which determines which part of the system we are solving. Returns 0, which is the boundary value for all points
This class holds the initial condition function we will use for the solver. Note that this class takes both MeshType and InitialConditionType as parameters. This class is capable of producing several different initial conditions without having to change the code each time, which makes it useful for running longer experiments without having to stop the code each time. The downside of this is the code is that the class is rather large, and functions have to be defined multiple times to be compatible with the different configurations of MESH and ICTYPE. Because of this, our implementation is not a good solution if more than a few variations of mesh and initial conditions need to be used. spacedim is the dimension of the input points. MESH is the type of mesh to apply initial conditions to, of type MeshType ICTYPE is the type of initial condition to apply, of type InitialConditionType
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
Stores the randomly generated coefficients for planar sine waves along the y-axis, used for psuedorandom initial conditions
The default constructor for the class. Initializes a function of 2 parameters and sets r and radius to default values. The constructor also loops through the coefficient arrays and stores the random coefficients for the psuedorandom initial condition.
An overloaded constructor, takes r and radius as parameters and uses these for initialization. Also loops through the coefficient arrays and stores the random coefficients for the psuedorandom initial condition. r is the value of the r parameter in the SH equation. radius is the radius of the hot spot
The return value of the initial condition function. This function is highly overloaded to account for a variety of different initial condition and mesh configurations, based on the template parameter given. Note that each initial condition sets the v component to 1e18. The v initial condition should not effect our solutions, and this is a good way to make any bugs causing v's initial condition to affect the solution easy to detect The RANDOM initial condition type does not change from mesh to mesh, it just returns a random number between -sqrt(r) and sqrt(r) The HOTSPOT initial condition changes the center depending on the input mesh type so that the hotspot is on the surface of the mesh The PSEUDORANDOM initial condition generates a function by summing up 10 sine waves in the x and y directions, with periods chosen so that the smallest period wave can still be resolved by a mesh with global refinement 5 or higher. On the plane, the value at each point is the product of the x sine sum and the y sine sum evaluated at the point. On the cylinder and Sinusoid, the x component is still used for the x sine sum, but we use ((arctan(y, z) - pi)/pi)*6*pi for the y sine sum. This wraps the psuedorandom function around the cylinder so that we can compare it to the same initial conditions on the plane. This function will run for the torus and sphere, but it has not been implemented to be comparable to the plane.
Places a small hot spot in the center of the plane on the u solution, and set v to a large number. p is the input point. component determines whether the input is for u or v. The function returns the value of the initial solution at the point
Places the hot spot in the center of the cylinder, on the positive z side. p is the input point. component determines whether the input is for u or v. The functions returns the value of the initial solution at the point.
Places the hot spot on the outside of the sphere, along the positive x axis p is the input point. component determines whether the input is for u or v. The function returns the value of the initial solution at the point.
Places the hot spot on the outside of the torus, along the x axis. p is the input point. component determines whether the input is for u or v. The function returns the value of the initial solution at the point.
Places the hot spot in the center of the sinusoid, on the positive z side. p is the input point. component determines whether the input is for u or v. The function returns the value of the initial solution at the point.
Returns the value of the psuedorandom function at the input point, as described above. p is the input point. component determines whether the input is for u or v. The function returns the value of the initial solution at the point.
Returns the value of the psuedorandom function at the input point, as described above. p is the input point. component determines whether the input is for u or v. The function returns the value of the initial solution at the point.
NOTE: Not particularly useful at the moment. Returns the value of the psuedorandom function at the input point, as described above. p is the input point. component determines whether the input is for u or v. The function returns the value of the initial solution at the point.
NOTE: Not particularly useful at the moment. Returns the value of the psuedorandom function at the input point, as described above. p is the input point. component determines whether the input is for u or v. The function returns the value of the initial solution at the point.
Returns the value of the psuedorandom function at the input point, as described above. p is the input point. component determines whether the input is for u or v. The function returns the value of the initial solution at the point.
Returns a random value between -sqrt(r) and sqrt(r)
Returns a random value between -sqrt(r) and sqrt(r)
Returns a random value between -sqrt(r) and sqrt(r)
Returns a random value between -sqrt(r) and sqrt(r)
Returns a random value between -sqrt(r) and sqrt(r)
Distributes the finite element vectors to each DoF, creates the system matrix, solution, old_solution, and system_rhs vectors, and outputs the number of DoF's to the console. dim is the dimension of the manifold. spacedim is the dimension of the ambient space. MESH is the type of mesh being used, doesn't change how this function works. ICTYPE is the type of initial condition used, doesn't change how this function works.
Counts the DoF's for outputting to consolse
Uses a direct solver to invert the system matrix, then multiplies the RHS vector by the inverted matrix to get the solution. Also includes a timer feature, which is currently commented out, but can be helpful to compute how long a run will take. dim is the dimension of the manifold. spacedim is the dimension of the ambient space. MESH is the type of mesh being used, doesn't change how this function works. ICTYPE is the type of initial condition used, doesn't change how this function works.
std::cout << "Solving linear system" << std::endl; Timer timer;
timer.stop(); std::cout << "done (" << timer.cpu_time() << " s)" << std::endl;
Converts the solution vector into a .vtu file and labels the outputs as u and v. dim is the dimension of the manifold. spacedim is the dimension of the ambient space. MESH is the type of mesh being used, doesn't change how this function works. ICTYPE is the type of initial condition used, doesn't change how this function works.
Takes the output_file_name string and appends timestep_number with up to three leading 0's
Below are all the different template cases for the make_grid() function
Runs the solver. First it creates the mesh and sets up the system, then constructs the system matrix, and finally loops over time to create the RHS vector and solve the system at each step. dim is the dimension of the manifold. spacedim is the dimension of the ambient space. MESH is the type of mesh being used. ICTYPE is the type of initial condition used, doesn't change how this function works.
Counts total time elapsed
Counts number of iterations
Sets the random seed so runs are repeatable, remove for varying random initial conditions
Applies the initial conditions to the old_solution
Outputs initial solution
Sets up the quadrature formula and FEValues object
The vector which stores the global indices that each local index connects to
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
These are the ith finite elements associated to u and v
These are the jth finite elements associated to u and v
This formula comes from expanding the PDE system
Loops over the dof indices to fill the entries of the system_matrix with the local data
Loops over time, incrementing by timestep, to create the RHS, solve the linear system, then output the result
Increments time and timestep_number
Outputs to console the number of iterations and current time. Currently outputs once every "second"
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
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
Loop over the quadrature points
Stores the value of the previous solution at the quadrature point
Loops over the dof indices to get the value of Un1
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
Loops over the dof indices to store the local data in the global RHS vector
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
Creates a volumetric cylinder
Extracts the boundary mesh with ID 0, which happens to be the tube part of the cylinder
The manifold information is lost upon boundary extraction. This sets the mesh boundary type to be a cylinder again
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
An array of mesh types. We iterate over this to allow for longer runs without having to stop the code
An array of initial condition types. We iterate this as well, for the same reason
Controls how long the code runs
The number of times we refine the hypercube mesh
The timestep will be 1/timestep_denominator
Loops over mesh types, then initial condition types, then loops over values of g_1
The value of g_1 passed to the solver object
Used to distinguish the start of each run
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.
1.9.8