From c36db94fb2f4fd31cb99282f91ab626e0d2a3a79 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Sat, 3 Jun 2023 00:58:07 +0000 Subject: [PATCH] (Hopefully) adjusted the syntax in the README to properly render matrices, fixed the image paths in the README so that images appear inline, changed the comment style in the Generalized-Swift-Hohenberg-Solver.cc file to (hopefully) render multi-line comments better --- .../Generalized-Swift-Hohenberg-Solver.cc | 392 +++++++++--------- Swift-Hohenberg-Solver/README.md | 78 ++-- 2 files changed, 246 insertions(+), 224 deletions(-) diff --git a/Swift-Hohenberg-Solver/Generalized-Swift-Hohenberg-Solver.cc b/Swift-Hohenberg-Solver/Generalized-Swift-Hohenberg-Solver.cc index d88478b..6e9284e 100644 --- a/Swift-Hohenberg-Solver/Generalized-Swift-Hohenberg-Solver.cc +++ b/Swift-Hohenberg-Solver/Generalized-Swift-Hohenberg-Solver.cc @@ -55,36 +55,38 @@ #include #include - 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 Point transform_function(const Point&p) { @@ -98,31 +100,16 @@ namespace SwiftHohenbergSolver } - /// @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 - Point transform_function_2_electric_boogaloo(const Point &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 class SHEquation { @@ -131,14 +118,15 @@ namespace SwiftHohenbergSolver 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 @@ -152,12 +140,14 @@ namespace SwiftHohenbergSolver 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(); @@ -174,8 +164,9 @@ namespace SwiftHohenbergSolver /// @brief Object holding the mesh Triangulation 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 fe; /// @brief Object which understands which finite elements are at each node DoFHandler dof_handler; @@ -186,13 +177,15 @@ namespace SwiftHohenbergSolver /// @brief Object holding the system matrix, stored as a sparse matrix SparseMatrix 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 solution; /// @brief Stores the solution from the previous timestep. Used to compute non-linear terms Vector 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 system_rhs; /// @brief Stores the current time, in the units of the problem @@ -221,10 +214,11 @@ namespace SwiftHohenbergSolver }; - /// @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 class BoundaryValues : public Function { @@ -239,12 +233,13 @@ namespace SwiftHohenbergSolver - /// @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 double BoundaryValues::value(const Point & p, const unsigned int component) const @@ -255,18 +250,19 @@ namespace SwiftHohenbergSolver 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 class InitialCondition : public Function { @@ -283,8 +279,9 @@ namespace SwiftHohenbergSolver 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(2), r(0.5), @@ -296,10 +293,11 @@ namespace SwiftHohenbergSolver } } - /// @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(2), @@ -312,32 +310,34 @@ namespace SwiftHohenbergSolver } } - /// @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 &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, @@ -356,10 +356,11 @@ namespace SwiftHohenbergSolver } } - /// @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, @@ -380,10 +381,11 @@ namespace SwiftHohenbergSolver } } - /// @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, @@ -404,10 +406,11 @@ namespace SwiftHohenbergSolver } } - /// @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, @@ -428,10 +431,11 @@ namespace SwiftHohenbergSolver } } - /// @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, @@ -452,10 +456,11 @@ namespace SwiftHohenbergSolver } } - /// @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, @@ -476,10 +481,11 @@ namespace SwiftHohenbergSolver } } - /// @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, @@ -501,10 +507,11 @@ namespace SwiftHohenbergSolver } } - /// @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, @@ -525,10 +532,11 @@ namespace SwiftHohenbergSolver } } - /// @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, @@ -549,10 +557,11 @@ namespace SwiftHohenbergSolver } } - /// @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, @@ -574,10 +583,11 @@ namespace SwiftHohenbergSolver } } - /// @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*/, @@ -591,10 +601,11 @@ namespace SwiftHohenbergSolver } } - /// @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*/, @@ -608,10 +619,11 @@ namespace SwiftHohenbergSolver } } - /// @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*/, @@ -625,10 +637,11 @@ namespace SwiftHohenbergSolver } } - /// @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*/, @@ -642,10 +655,11 @@ namespace SwiftHohenbergSolver } } - /// @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*/, @@ -695,12 +709,13 @@ namespace SwiftHohenbergSolver , 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 void SHEquation::setup_system() { @@ -733,12 +748,13 @@ namespace SwiftHohenbergSolver } - /// @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 void SHEquation::solve_time_step() { @@ -757,11 +773,12 @@ namespace SwiftHohenbergSolver - /// @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 void SHEquation::output_results() const { @@ -881,12 +898,13 @@ namespace SwiftHohenbergSolver } - /// @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 void SHEquation::run() { diff --git a/Swift-Hohenberg-Solver/README.md b/Swift-Hohenberg-Solver/README.md index a7e0925..87a451b 100755 --- a/Swift-Hohenberg-Solver/README.md +++ b/Swift-Hohenberg-Solver/README.md @@ -2,7 +2,9 @@ This program is used to solve the generalized Swift-Hohenberg equation -$$\frac{\partial u}{\partial t} = ru - (k_c + \Delta)^2 u + g_1 u^2 - u^3$$ +$$\begin{aligned} + \frac{\partial u}{\partial t} = ru - (k_c + \Delta)^2 u + g_1 u^2 - u^3 +\end{aligned}$$ 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. @@ -15,7 +17,9 @@ are interesting behaviors that occur when $g_1$ is smaller or larger than $r$ in magnitude, so this allows us room to vary $g_1$ and explore these behavior. To summarize, this code solves: -$$\frac{\partial u}{\partial t} = 0.3u - (1 + \Delta)^2 u + g_1 u^2 - u^3$$ +$$\begin{aligned} + \frac{\partial u}{\partial t} = 0.3u - (1 + \Delta)^2 u + g_1 u^2 - u^3 +\end{aligned}$$ # Discretization and Solving the Bilaplacian @@ -48,69 +52,69 @@ the previous timestep. We then reframe this system as a vector valued problem $$\begin{aligned} - \begin{pmatrix} + \left(\begin{matrix} 1 - kr & k(1 + \Delta)\\ 1 + \Delta & -1 - \end{pmatrix} - \begin{pmatrix} + \end{matrix}\right) + \left(\begin{matrix} U_n\\ V_n - \end{pmatrix} &= \begin{pmatrix} + \end{matrix}\right) &= \left(\begin{matrix} U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\ 0 - \end{pmatrix} + \end{matrix}\right) \end{aligned}$$ As usual, we multiply each side of the equation by a test function -$$\overrightarrow\varphi_i = \begin{pmatrix} +$$\overrightarrow\varphi_i = \left(\begin{matrix} \phi_i\\ \psi_i -\end{pmatrix}$$ +\end{matrix}\right)$$ and then integrate over the domain $\Omega$ to get the equation $$\begin{aligned} - \int_\Omega \begin{pmatrix} + \int_\Omega \left(\begin{matrix} \phi_i\\ \psi_i - \end{pmatrix}\cdot\begin{pmatrix} + \end{matrix}\right)\cdot\left(\begin{matrix} 1 - kr & k(1 + \Delta)\\ 1 + \Delta & -1 - \end{pmatrix} - \begin{pmatrix} + \end{matrix}\right) + \left(\begin{matrix} U_n\\ V_n - \end{pmatrix} &= \int_\Omega \begin{pmatrix} + \end{matrix}\right) &= \int_\Omega \left(\begin{matrix} \phi_i\\ \psi_i - \end{pmatrix}\cdot\begin{pmatrix} + \end{matrix}\right)\cdot\left(\begin{matrix} U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\ 0 - \end{pmatrix}\\ + \end{matrix}\right)\\ \end{aligned}$$ We can expand our solution vector in this basis $$\begin{aligned} - \int_\Omega \sum_j u_j\begin{pmatrix} + \int_\Omega \sum_j u_j\left(\begin{matrix} \phi_i\\ \psi_i - \end{pmatrix}\cdot\begin{pmatrix} + \end{matrix}\right)\cdot\left(\begin{matrix} 1 - kr & k(1 + \Delta)\\ 1 + \Delta & -1 - \end{pmatrix} - \begin{pmatrix} + \end{matrix}\right) + \left(\begin{matrix} \phi_j\\ \psi_j - \end{pmatrix} &= \int_\Omega\begin{pmatrix} + \end{matrix}\right) &= \int_\Omega\left(\begin{matrix} \phi_i\\ \psi_i - \end{pmatrix}\cdot\begin{pmatrix} + \end{matrix}\right)\cdot\left(\begin{matrix} U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\ 0 - \end{pmatrix} + \end{matrix}\right) \end{aligned}$$ and finally expand out the matrix multiplication @@ -176,7 +180,7 @@ on the finest mesh. Below are the results of several runs of constant initial conditions -![image](doc/images/Figures_1_and_2.png) +![image](./doc/images/Figures_1_and_2.png) 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 @@ -184,7 +188,7 @@ condition for each is shown above, While the final solutions are shown in the ma timestep begins at 1/25 and the denominator increases by 25 across each row, to a max of 1/200 in the bottom right: -![image](doc/images/TC_table.png) +![image](./doc/images/TC_table.png) We validated that solutions converged across mesh refinement by defining psuedorandom functions @@ -198,9 +202,9 @@ so that the smallest wave could be resolved by a mesh refinement of 7 or higher. The following matrix shows the initial and final solution ranging from a refinement of 0 to a refinement of 7: -![image](doc/images/Refinement_Convergence_Table_1.png) +![image](./doc/images/Refinement_Convergence_Table_1.png) -![image](doc/images/Refinement_Convergence_Table_2.png) +![image](./doc/images/Refinement_Convergence_Table_2.png) # Results @@ -211,29 +215,29 @@ pieces as $g_1$ is increased. In the matrix below, $g_1$ is increased by 0.2 starting from 0 to a maximum value of 1.4. Note that each final solution is at 100 time units: -![image](doc/images/Square_Hotspot_Table.png) +![image](./doc/images/Square_Hotspot_Table.png) On the cylinder, the front looks similar to the square, but the back has an overlapping wave pattern: -![image](doc/images/Cylinder_Hotspot_Table.png) +![image](./doc/images/Cylinder_Hotspot_Table.png) On the sphere, the hot spot generates a single wave. Note that this may be due to the fact that our sphere has a surface area proportional to the period of our pattern wave. -![image](doc/images/Sphere_Hotspot_Table.png) +![image](./doc/images/Sphere_Hotspot_Table.png) On the torus, the pattern propagates similar to the cylinder, with some minor imperfections -![image](doc/images/Torus_Hotspot_Front_Table.png) +![image](./doc/images/Torus_Hotspot_Front_Table.png) But on the back side of the torus, we see wave overlapping and spot patterns forming -![image](doc/images/Torus_Hotspot_Back_Table.png) +![image](./doc/images/Torus_Hotspot_Back_Table.png) 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 @@ -241,16 +245,16 @@ was made by warping the boundary of a cylinder by a cosine wave, and is equivalent to the surface of revolution bounded by $1 + 0.5\cos(\frac{\pi}{10}x)$ -![image](doc/images/Sinusoid_Hotspot_Front_Table.png) +![image](./doc/images/Sinusoid_Hotspot_Front_Table.png) -![image](doc/images/Sinusoid_Hotspot_Back_Table.png) +![image](./doc/images/Sinusoid_Hotspot_Back_Table.png) Finally, here is a small selection of random initial conditions and the patterns that form. Each image sequence was taken at times 0, 10, 25, 50, and 100: -![image](doc/images/Square_Random_Table.png) +![image](./doc/images/Square_Random_Table.png) -![image](doc/images/Sphere_Random_Table.png) +![image](./doc/images/Sphere_Random_Table.png) -![image](doc/images/Sinusoid_Random_Table.png) \ No newline at end of file +![image](./doc/images/Sinusoid_Random_Table.png) \ No newline at end of file -- 2.39.5