#include <iostream>
#include <random>
-
namespace SwiftHohenbergSolver
{
using namespace dealii;
- /// @brief 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.
+ /** @brief 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.
+ */
enum MeshType {HYPERCUBE, CYLINDER, SPHERE, TORUS, SINUSOID};
- /// @brief 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
+ /** @brief 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
+ */
enum InitialConditionType {HOTSPOT, PSUEDORANDOM, RANDOM};
- /// @brief This function warps points on a cyclindrical 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.
- /// @tparam spacedim This is the dimension of the embedding space, which is where the input point lives
- /// @param p This is thel input point to be translated.
- /// @return The return as a tranlated point in the same dimensional space. This is the new point on the mesh.
+ /** @brief This function warps points on a cyclindrical 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.
+ * @tparam spacedim This is the dimension of the embedding space, which is where the input point lives
+ * @param p This is the input point to be translated.
+ * @return The return as a tranlated point in the same dimensional space. This is the new point on the mesh.
+ */
template<int spacedim>
Point<spacedim> transform_function(const Point<spacedim>&p)
{
}
- /// @brief Not currently implemented, but will function the same as above only with and undulary boundary curve rather
- /// than a cosine boundary curve.
- /// @tparam spacedim See above
- /// @param p See above
- /// @return See above
- template<int spacedim>
- Point<spacedim> transform_function_2_electric_boogaloo(const Point<spacedim> &p)
- {
- Assert(spacedim == 3, ExcNotImplemented());
- return 0;
- }
-
-
-
-
-
- /// @brief 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.
- /// @tparam dim This is the intrinsic dimension of the manifold we are solving on.
- /// @tparam spacedim This is the dimension of the embedding space.
- /// @tparam MESH This determines what manifold we are solving on
- /// @tparam ICTYPE This determines what initial condition we use
+ /** @brief 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.
+ * @tparam dim This is the intrinsic dimension of the manifold we are solving on.
+ * @tparam spacedim This is the dimension of the embedding space.
+ * @tparam MESH This determines what manifold we are solving on
+ * @tparam ICTYPE This determines what initial condition we use
+ */
template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
class SHEquation
{
SHEquation();
- /// @brief Overloaded constructor, allows user to pass values for important constants
- /// @param degree This is the degree of finite element used
- /// @param time_step_denominator This determines what size timestep we use. The timestep is 1/time_step_denominator
- /// @param ref_num The number of times the mesh will be globally refined.
- /// @param r_constant Constant for linear component, default 0.5
- /// @param g1_constant Constant for quadratic component, default 0.5
- /// @param output_file_name Self explanatory, default "solution-"
- /// @param end_time Determines when the solver stops, default 0.5, should be ~100 to see equilibrium solutions
+ /** @brief Overloaded constructor, allows user to pass values for important constants
+ * @param degree This is the degree of finite element used
+ * @param time_step_denominator This determines what size timestep we use. The timestep is 1/time_step_denominator
+ * @param ref_num The number of times the mesh will be globally refined.
+ * @param r_constant Constant for linear component, default 0.5
+ * @param g1_constant Constant for quadratic component, default 0.5
+ * @param output_file_name Self explanatory, default "solution-"
+ * @param end_time Determines when the solver stops, default 0.5, should be ~100 to see equilibrium solutions
+ */
SHEquation(const unsigned int degree
, double time_step_denominator
, unsigned int ref_num
void setup_system();
void solve_time_step();
void output_results() const;
- /// @brief 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.
+ /** @brief 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.
+ */
void make_grid();
- /// @brief 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
+ /** @brief 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
+ */
void make_cylinder();
/// @brief 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();
/// @brief Object holding the mesh
Triangulation<dim, spacedim> triangulation;
- /// @brief Object describing the finite element vectors at each node
- /// (I believe this gives a basis for the finite elements at each node)
+ /** @brief Object describing the finite element vectors at each node
+ * (I believe this gives a basis for the finite elements at each node)
+ */
FESystem<dim, spacedim> fe;
/// @brief Object which understands which finite elements are at each node
DoFHandler<dim, spacedim> dof_handler;
/// @brief Object holding the system matrix, stored as a sparse matrix
SparseMatrix<double> system_matrix;
- /// @brief Vector of coefficients for the solution in the current timestep
- /// We solve for this in each timestep
+ /** @brief Vector of coefficients for the solution in the current timestep
+ * We solve for this in each timestep
+ */
Vector<double> solution;
/// @brief Stores the solution from the previous timestep. Used to compute non-linear terms
Vector<double> old_solution;
- /// @brief Stores the coefficients of the right hand side function(in terms of the finite elements)
- /// Is the RHS for the linear system
+ /** @brief Stores the coefficients of the right hand side function(in terms of the finite elements)
+ * Is the RHS for the linear system
+ */
Vector<double> system_rhs;
/// @brief Stores the current time, in the units of the problem
};
- /// @brief 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.
- /// @tparam spacedim The dimension of the points which the function takes as input
+ /** @brief 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.
+ * @tparam spacedim The dimension of the points which the function takes as input
+ */
template <int spacedim>
class BoundaryValues : public Function<spacedim>
{
- /// @brief Returns 0 for all points. This is the output for the boundary
- /// @tparam spacedim The dimension of points that are input
- /// @param p The input point
- /// @param component Determines whether we are solving for u or v.
- /// This determines which part of the system we are solving
- /// @return 0; This is the boundary value for all points
+ /** @brief Returns 0 for all points. This is the output for the boundary
+ * @tparam spacedim The dimension of points that are input
+ * @param p The input point
+ * @param component Determines whether we are solving for u or v.
+ * This determines which part of the system we are solving
+ * @return 0; This is the boundary value for all points
+ */
template <int spacedim>
double BoundaryValues<spacedim>::value(const Point<spacedim> & p,
const unsigned int component) const
return 0.;
}
- /// @brief 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.
- /// @tparam spacedim The dimension of the input points
- /// @tparam MESH The type of mesh to apply initial conditions to, of type MeshType
- /// @tparam ICTYPE The type of initial condition to apply, of type InitialConditionType
+ /** @brief 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.
+ * @tparam spacedim The dimension of the input points
+ * @tparam MESH The type of mesh to apply initial conditions to, of type MeshType
+ * @tparam ICTYPE The type of initial condition to apply, of type InitialConditionType
+ */
template<int spacedim, MeshType MESH, InitialConditionType ICTYPE>
class InitialCondition : public Function<spacedim>
{
double y_sin_coefficients[10];
public:
- /// @brief 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.
+ /** @brief 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.
+ */
InitialCondition()
: Function<spacedim>(2),
r(0.5),
}
}
- /// @brief 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.
- /// @param r The value of the r parameter in the SH equation
- /// @param radius The radius of the hot spot
+ /** @brief 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.
+ * @param r The value of the r parameter in the SH equation
+ * @param radius The radius of the hot spot
+ */
InitialCondition(const double r,
const double radius)
: Function<spacedim>(2),
}
}
- /// @brief 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.
- /// @param p
- /// @param component
- /// @return
+ /** @brief 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.
+ * @param p
+ * @param component
+ * @return
+ */
virtual double value(const Point<spacedim> &p, const unsigned int component) const override;
};
- /// @brief Places a small hot spot in the center of the plane on the u solution, and set v to a large number
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Places a small hot spot in the center of the plane on the u solution, and set v to a large number
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<2, HYPERCUBE, HOTSPOT>::value(
const Point<2> &p,
}
}
- /// @brief Places the hot spot in the center of the cylinder, on the positive z side
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Places the hot spot in the center of the cylinder, on the positive z side
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, CYLINDER, HOTSPOT>::value(
const Point<3> &p,
}
}
- /// @brief Places the hot spot on the outside of the sphere, along the positive x axis
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Places the hot spot on the outside of the sphere, along the positive x axis
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, SPHERE, HOTSPOT>::value(
const Point<3> &p,
}
}
- /// @brief Places the hot spot on the outside of the torus, along the x axis
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Places the hot spot on the outside of the torus, along the x axis
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, TORUS, HOTSPOT>::value(
const Point<3> &p,
}
}
- /// @brief Places the hot spot in the center of the sinusoid, on the positive z side
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Places the hot spot in the center of the sinusoid, on the positive z side
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, SINUSOID, HOTSPOT>::value(
const Point<3> &p,
}
}
- /// @brief Returns the value of the psuedorandom function at the input point, as described above
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Returns the value of the psuedorandom function at the input point, as described above
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<2, HYPERCUBE, PSUEDORANDOM>::value(
const Point<2> &p,
}
}
- /// @brief Returns the value of the psuedorandom function at the input point, as described above
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Returns the value of the psuedorandom function at the input point, as described above
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, CYLINDER, PSUEDORANDOM>::value(
const Point<3> &p,
}
}
- /// @brief NOTE: Not particularly useful at the moment. Returns the value of the psuedorandom function at the input point, as described above
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief NOTE: Not particularly useful at the moment. Returns the value of the psuedorandom function at the input point, as described above
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, SPHERE, PSUEDORANDOM>::value(
const Point<3> &p,
}
}
- /// @brief NOTE: Not particularly useful at the moment. Returns the value of the psuedorandom function at the input point, as described above
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief NOTE: Not particularly useful at the moment. Returns the value of the psuedorandom function at the input point, as described above
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, TORUS, PSUEDORANDOM>::value(
const Point<3> &p,
}
}
- /// @brief Returns the value of the psuedorandom function at the input point, as described above
- /// @param p The input point
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Returns the value of the psuedorandom function at the input point, as described above
+ * @param p The input point
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, SINUSOID, PSUEDORANDOM>::value(
const Point<3> &p,
}
}
- /// @brief Returns a random value between -sqrt(r) and sqrt(r)
- /// @param p The input point, not used in this function
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Returns a random value between -sqrt(r) and sqrt(r)
+ * @param p The input point, not used in this function
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<2, HYPERCUBE, RANDOM>::value(
const Point<2> &/*p*/,
}
}
- /// @brief Returns a random value between -sqrt(r) and sqrt(r)
- /// @param p The input point, not used in this function
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Returns a random value between -sqrt(r) and sqrt(r)
+ * @param p The input point, not used in this function
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, CYLINDER, RANDOM>::value(
const Point<3> &/*p*/,
}
}
- /// @brief Returns a random value between -sqrt(r) and sqrt(r)
- /// @param p The input point, not used in this function
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Returns a random value between -sqrt(r) and sqrt(r)
+ * @param p The input point, not used in this function
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, SPHERE, RANDOM>::value(
const Point<3> &/*p*/,
}
}
- /// @brief Returns a random value between -sqrt(r) and sqrt(r)
- /// @param p The input point, not used in this function
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Returns a random value between -sqrt(r) and sqrt(r)
+ * @param p The input point, not used in this function
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, TORUS, RANDOM>::value(
const Point<3> &/*p*/,
}
}
- /// @brief Returns a random value between -sqrt(r) and sqrt(r)
- /// @param p The input point, not used in this function
- /// @param component Determines whether the input is for u or v
- /// @return The value of the initial solution at the point
+ /** @brief Returns a random value between -sqrt(r) and sqrt(r)
+ * @param p The input point, not used in this function
+ * @param component Determines whether the input is for u or v
+ * @return The value of the initial solution at the point
+ */
template <>
double InitialCondition<3, SINUSOID, RANDOM>::value(
const Point<3> &/*p*/,
, end_time(end_time)
{}
- /// @brief Distrubutes 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.
- /// @tparam dim The dimension of the manifold
- /// @tparam spacedim The dimension of the ambient space
- /// @tparam MESH The type of mesh being used, doesn't change how this function works
- /// @tparam ICTYPE The type of initial condition used, doesn't change how this function works
+ /** @brief Distrubutes 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.
+ * @tparam dim The dimension of the manifold
+ * @tparam spacedim The dimension of the ambient space
+ * @tparam MESH The type of mesh being used, doesn't change how this function works
+ * @tparam ICTYPE The type of initial condition used, doesn't change how this function works
+ */
template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
void SHEquation<dim, spacedim, MESH, ICTYPE>::setup_system()
{
}
- /// @brief 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
- /// @tparam dim The dimension of the manifold
- /// @tparam spacedim The dimension of the ambient space
- /// @tparam MESH The type of mesh being used, doesn't change how this function works
- /// @tparam ICTYPE The type of initial condition used, doesn't change how this function works
+ /** @brief 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
+ * @tparam dim The dimension of the manifold
+ * @tparam spacedim The dimension of the ambient space
+ * @tparam MESH The type of mesh being used, doesn't change how this function works
+ * @tparam ICTYPE The type of initial condition used, doesn't change how this function works
+ */
template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
void SHEquation<dim, spacedim, MESH, ICTYPE>::solve_time_step()
{
- /// @brief Converts the solution vector into a .vtu file and labels the outputs as u and v
- /// @tparam dim The dimension of the manifold
- /// @tparam spacedim The dimension of the ambient space
- /// @tparam MESH The type of mesh being used, doesn't change how this function works
- /// @tparam ICTYPE The type of initial condition used, doesn't change how this function works
+ /** @brief Converts the solution vector into a .vtu file and labels the outputs as u and v
+ * @tparam dim The dimension of the manifold
+ * @tparam spacedim The dimension of the ambient space
+ * @tparam MESH The type of mesh being used, doesn't change how this function works
+ * @tparam ICTYPE The type of initial condition used, doesn't change how this function works
+ */
template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
void SHEquation<dim, spacedim, MESH, ICTYPE>::output_results() const
{
}
- /// @brief 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
- /// @tparam dim The dimension of the manifold
- /// @tparam spacedim The dimension of the ambient space
- /// @tparam MESH The type of mesh being used
- /// @tparam ICTYPE The type of initial condition used, doesn't change how this function works
+ /** @brief 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
+ * @tparam dim The dimension of the manifold
+ * @tparam spacedim The dimension of the ambient space
+ * @tparam MESH The type of mesh being used
+ * @tparam ICTYPE The type of initial condition used, doesn't change how this function works
+ */
template <int dim, int spacedim, MeshType MESH, InitialConditionType ICTYPE>
void SHEquation<dim, spacedim, MESH, ICTYPE>::run()
{
\r
This program is used to solve the generalized Swift-Hohenberg equation\r
\r
-$$\frac{\partial u}{\partial t} = ru - (k_c + \Delta)^2 u + g_1 u^2 - u^3$$\r
+$$\begin{aligned}\r
+ \frac{\partial u}{\partial t} = ru - (k_c + \Delta)^2 u + g_1 u^2 - u^3\r
+\end{aligned}$$\r
\r
where $k_c$ is the wave number, $r$ is some fixed constant, and\r
$g_1$ is a parameter which determines the behavior of the solutions.\r
than $r$ in magnitude, so this allows us room to vary $g_1$ and\r
explore these behavior. To summarize, this code solves:\r
\r
-$$\frac{\partial u}{\partial t} = 0.3u - (1 + \Delta)^2 u + g_1 u^2 - u^3$$\r
+$$\begin{aligned}\r
+ \frac{\partial u}{\partial t} = 0.3u - (1 + \Delta)^2 u + g_1 u^2 - u^3\r
+\end{aligned}$$\r
\r
# Discretization and Solving the Bilaplacian\r
\r
problem\r
\r
$$\begin{aligned}\r
- \begin{pmatrix}\r
+ \left(\begin{matrix}\r
1 - kr & k(1 + \Delta)\\\r
1 + \Delta & -1\r
- \end{pmatrix}\r
- \begin{pmatrix}\r
+ \end{matrix}\right)\r
+ \left(\begin{matrix}\r
U_n\\\r
V_n\r
- \end{pmatrix} &= \begin{pmatrix}\r
+ \end{matrix}\right) &= \left(\begin{matrix}\r
U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\\r
0\r
- \end{pmatrix}\r
+ \end{matrix}\right)\r
\end{aligned}$$\r
\r
As usual, we multiply each side of the equation by a\r
test function \r
\r
-$$\overrightarrow\varphi_i = \begin{pmatrix}\r
+$$\overrightarrow\varphi_i = \left(\begin{matrix}\r
\phi_i\\\r
\psi_i\r
-\end{pmatrix}$$\r
+\end{matrix}\right)$$\r
\r
and then integrate over the domain $\Omega$ to get the equation\r
\r
$$\begin{aligned}\r
- \int_\Omega \begin{pmatrix}\r
+ \int_\Omega \left(\begin{matrix}\r
\phi_i\\ \r
\psi_i\r
- \end{pmatrix}\cdot\begin{pmatrix}\r
+ \end{matrix}\right)\cdot\left(\begin{matrix}\r
1 - kr & k(1 + \Delta)\\\r
1 + \Delta & -1\r
- \end{pmatrix}\r
- \begin{pmatrix}\r
+ \end{matrix}\right)\r
+ \left(\begin{matrix}\r
U_n\\\r
V_n\r
- \end{pmatrix} &= \int_\Omega \begin{pmatrix}\r
+ \end{matrix}\right) &= \int_\Omega \left(\begin{matrix}\r
\phi_i\\\r
\psi_i\r
- \end{pmatrix}\cdot\begin{pmatrix}\r
+ \end{matrix}\right)\cdot\left(\begin{matrix}\r
U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\\r
0\r
- \end{pmatrix}\\\r
+ \end{matrix}\right)\\\r
\end{aligned}$$\r
\r
We can expand our solution vector in this basis\r
\r
$$\begin{aligned}\r
- \int_\Omega \sum_j u_j\begin{pmatrix}\r
+ \int_\Omega \sum_j u_j\left(\begin{matrix}\r
\phi_i\\\r
\psi_i\r
- \end{pmatrix}\cdot\begin{pmatrix}\r
+ \end{matrix}\right)\cdot\left(\begin{matrix}\r
1 - kr & k(1 + \Delta)\\\r
1 + \Delta & -1\r
- \end{pmatrix}\r
- \begin{pmatrix}\r
+ \end{matrix}\right)\r
+ \left(\begin{matrix}\r
\phi_j\\\r
\psi_j\r
- \end{pmatrix} &= \int_\Omega\begin{pmatrix}\r
+ \end{matrix}\right) &= \int_\Omega\left(\begin{matrix}\r
\phi_i\\\r
\psi_i\r
- \end{pmatrix}\cdot\begin{pmatrix}\r
+ \end{matrix}\right)\cdot\left(\begin{matrix}\r
U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\\r
0\r
- \end{pmatrix}\r
+ \end{matrix}\right)\r
\end{aligned}$$\r
\r
and finally expand out the matrix multiplication\r
\r
Below are the results of several runs of constant initial conditions\r
\r
-\r
+\r
\r
We also validated that given a fixed random start on a very fine mesh,\r
refining the timestep resulted in the same final solution. The initial\r
timestep begins at 1/25 and the denominator increases by 25 across each\r
row, to a max of 1/200 in the bottom right:\r
\r
-\r
+\r
\r
We validated that solutions converged across mesh refinement by defining\r
psuedorandom functions\r
higher. The following matrix shows the initial and final solution\r
ranging from a refinement of 0 to a refinement of 7:\r
\r
-\r
+\r
\r
-\r
+\r
\r
# Results\r
\r
increased by 0.2 starting from 0 to a maximum value of 1.4. Note that\r
each final solution is at 100 time units:\r
\r
-\r
+\r
\r
On the cylinder, the front looks similar to the square, but the back has\r
an overlapping wave pattern:\r
\r
-\r
+\r
\r
\r
On the sphere, the hot spot generates a single wave. Note that this may\r
be due to the fact that our sphere has a surface area proportional to\r
the period of our pattern wave.\r
\r
-\r
+\r
\r
On the torus, the pattern propagates similar to the cylinder, with some\r
minor imperfections\r
\r
-\r
+\r
\r
But on the back side of the torus, we see wave overlapping and spot\r
patterns forming\r
\r
-\r
+\r
\r
On shapes with stranger curvature, we can see that the pattern wave has\r
a tendency to break apart when crossing lines of curvature. This shape\r
equivalent to the surface of revolution bounded by\r
$1 + 0.5\cos(\frac{\pi}{10}x)$\r
\r
-\r
+\r
\r
-\r
+\r
\r
Finally, here is a small selection of random initial conditions and the\r
patterns that form. Each image sequence was taken at times 0, 10, 25,\r
50, and 100:\r
\r
-\r
+\r
\r
-\r
+\r
\r
-
\ No newline at end of file
+
\ No newline at end of file
projects / code-gallery.git / commitdiff