Crystal_Growth_Phase_Field_Model.h

Go to the documentation of this file.

 
 *   
 *   #include "PhaseFieldSolver.h"
 *   
 *   void InitialValues::vector_value(const Point<2> &p,
 *                                    Vector<double> & values) const
 *   {
 *       values(0)= 0.0; //Initial p value of domain
 *       values(1)= 0.2; //Initial temperature of domain
 *   }
 * @endcode
 
 
<a name="ann-PhaseFieldSolver.cpp"></a>
<h1>Annotated version of PhaseFieldSolver.cpp</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Umair Hussain
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "PhaseFieldSolver.h"
 *   
 *   PhaseFieldSolver::PhaseFieldSolver()
 *   : mpi_communicator(MPI_COMM_WORLD)
 *   , n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
 *   , this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
 *   , pcout(std::cout, (this_mpi_process == 0))
 *   , fe(FE_Q<2>(1), 2)
 *   , dof_handler(triangulation)
 *   , time(0.0)
 *   , final_time(1.)
 *   , time_step(.0002)
 *   , theta(0.5)
 *   , epsilon(0.01)
 *   , tau(0.0003)
 *   , gamma(10.)
 *   , latent_heat(1.4)
 *   , alpha(0.9)
 *   , t_eq(1.)
 *   , a(0.01)
 *   {}
 * @endcode
 
 
<a name="ann-PhaseFieldSolver.h"></a>
<h1>Annotated version of PhaseFieldSolver.h</h1>
 * 
 * 
 * 
 * 
 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *   *
 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *   * Copyright (C) 2024 by Umair Hussain
 *   *
 *   * This file is part of the deal.II code gallery.
 *   *
 *   * -----------------------------------------------------------------------------
 *   */
 *   
 *   #ifndef KOBAYASHI_PARALLEL_PHASEFIELDSOLVER_H
 *   #define KOBAYASHI_PARALLEL_PHASEFIELDSOLVER_H
 *   
 *   #include <deal.II/base/quadrature_lib.h>
 *   #include <deal.II/base/function.h>
 *   #include <deal.II/base/utilities.h>
 *   #include <deal.II/lac/vector.h>
 *   #include <deal.II/lac/full_matrix.h>
 *   #include <deal.II/lac/sparse_matrix.h>
 *   #include <deal.II/lac/sparse_direct.h>
 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
 *   #include <deal.II/lac/solver_cg.h>
 *   #include <deal.II/lac/precondition.h>
 *   #include <deal.II/lac/affine_constraints.h>
 *   #include <deal.II/grid/tria.h>
 *   #include <deal.II/grid/grid_generator.h>
 *   #include <deal.II/dofs/dof_handler.h>
 *   #include <deal.II/dofs/dof_tools.h>
 *   #include <deal.II/fe/fe_q.h>
 *   #include <deal.II/fe/fe_values.h>
 *   #include <deal.II/fe/fe_system.h>
 *   #include <deal.II/numerics/vector_tools.h>
 *   #include <deal.II/numerics/matrix_tools.h>
 *   #include <deal.II/numerics/data_out.h>
 *   #include <deal.II/grid/grid_in.h>
 *   
 * @endcode
 * 
 * For Parallel Computation
 * 
 * @code
 *   #include <deal.II/base/conditional_ostream.h>
 *   #include <deal.II/base/mpi.h>
 *   #include <deal.II/lac/petsc_vector.h>
 *   #include <deal.II/lac/petsc_sparse_matrix.h>
 *   #include <deal.II/lac/petsc_solver.h>
 *   #include <deal.II/lac/petsc_precondition.h>
 *   #include <deal.II/grid/grid_tools.h>
 *   #include <deal.II/dofs/dof_renumbering.h>
 *   
 *   #include <deal.II/distributed/solution_transfer.h>
 *   
 *   #include <fstream>
 *   #include <iostream>
 *   
 *   using namespace dealii;
 *   
 *   class PhaseFieldSolver {
 *   public:
 *       PhaseFieldSolver();
 *       void run();
 *   
 *   private:
 *       void          make_grid_and_dofs();
 *       void          assemble_system();
 *       void          solve();
 *       void          output_results(const unsigned int timestep_number) const;
 *       double        compute_residual();
 *       void          applying_bc();
 *       float         get_random_number();
 *   
 *       MPI_Comm mpi_communicator;
 *       const unsigned int n_mpi_processes;
 *       const unsigned int this_mpi_process;
 *       ConditionalOStream pcout;
 *   
 *       Triangulation<2>     triangulation;
 *       FESystem<2>          fe;
 *       DoFHandler<2>        dof_handler;
 *       GridIn<2>            gridin;
 *   
 *       PETScWrappers::MPI::SparseMatrix jacobian_matrix;
 *   
 *       double       time;
 *       const double final_time, time_step;
 *       const double theta;
 *       const double epsilon, tau, gamma, latent_heat, alpha, t_eq, a; //as given in Ref. [1]
 *   
 *       PETScWrappers::MPI::Vector conv_solution; //solution vector at last newton-raphson iteration
 *       PETScWrappers::MPI::Vector old_solution; //solution vector at last time step
 *       PETScWrappers::MPI::Vector solution_update; //increment in solution or delta solution
 *       PETScWrappers::MPI::Vector system_rhs; //to store residual
 *       Vector<double> conv_solution_np, old_solution_np; //creating non parallel vectors to store data for easy access of old solution values by all processes
 *   
 *   };
 *   
 * @endcode
 * 
 * Initial values class
 * 
 * @code
 *   class InitialValues : public Function<2>
 *   {
 *   public:
 *       InitialValues(): Function<2>(2)
 *       {}
 *       virtual void vector_value(const Point<2> &p,
 *                                 Vector<double> & value) const override;
 *   };
 *   
 *   
 *   #endif //KOBAYASHI_PARALLEL_PHASEFIELDSOLVER_H
 * @endcode
 
 
<a name="ann-applying_bc.cpp"></a>
<h1>Annotated version of applying_bc.cpp</h1>
 * 
 * 
 * 
 * 
 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Umair Hussain
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "PhaseFieldSolver.h"
 *   
 *   void PhaseFieldSolver::applying_bc(){
 *       FEValuesExtractors::Scalar phase_parameter(0);
 *       FEValuesExtractors::Scalar temperature(1);
 *   
 *       QGauss<2> quadrature_formula(fe.degree + 1);
 *       FEValues<2> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients | update_JxW_values);
 *   
 *       ComponentMask p_mask = fe.component_mask(phase_parameter);
 *       ComponentMask t_mask = fe.component_mask(temperature);
 *   
 *       std::map<types::global_dof_index,double> boundary_values;
 *   
 * @endcode
 * 
 * Prescribing p=1 at the left face (this will be maintained in the subsequent iterations when zero BC is applied in the Newton-Raphson iterations)
 * 
 * @code
 *       VectorTools::interpolate_boundary_values (dof_handler,
 *                                                 1,
 *                                                 Functions::ConstantFunction<2>(1., 2),
 *                                                 boundary_values,p_mask);
 *   
 * @endcode
 * 
 * To apply the boundary values only to the solution vector without the Jacobian Matrix and RHS Vector
 * 
 * @code
 *       for (auto &boundary_value : boundary_values)
 *           old_solution(boundary_value.first) = boundary_value.second;
 *   
 *       old_solution.compress(VectorOperation::insert);
 *   
 *   }
 * @endcode
 
 
<a name="ann-assemble_system.cpp"></a>
<h1>Annotated version of assemble_system.cpp</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Umair Hussain
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "PhaseFieldSolver.h"
 *   #include <cmath>
 *   
 *   void PhaseFieldSolver::assemble_system() {
 * @endcode
 * 
 * Separating each variable as a scalar to easily call the respective shape functions
 * 
 * @code
 *       FEValuesExtractors::Scalar phase_parameter(0);
 *       FEValuesExtractors::Scalar temperature(1);
 *   
 *       QGauss<2> quadrature_formula(fe.degree + 1);
 *       FEValues<2> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients | update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *       FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *   
 *       Vector<double>     cell_rhs(dofs_per_cell);
 *   
 * @endcode
 * 
 * To copy values and gradients of solution from previous iteration
 * Old Newton iteration
 * 
 * @code
 *       std::vector<Tensor<1, 2>> old_newton_solution_gradients_p(n_q_points);
 *       std::vector<double> old_newton_solution_values_p(n_q_points);
 *       std::vector<Tensor<1, 2>> old_newton_solution_gradients_t(n_q_points);
 *       std::vector<double> old_newton_solution_values_t(n_q_points);
 * @endcode
 * 
 * Old time step iteration
 * 
 * @code
 *       std::vector<Tensor<1, 2>> old_time_solution_gradients_p(n_q_points);
 *       std::vector<double> old_time_solution_values_p(n_q_points);
 *       std::vector<Tensor<1, 2>> old_time_solution_gradients_t(n_q_points);
 *       std::vector<double> old_time_solution_values_t(n_q_points);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *       jacobian_matrix.operator=(0.0);
 *       system_rhs.operator=(0.0);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators()){
 *           if (cell->subdomain_id() == this_mpi_process) {
 *               cell_matrix = 0;
 *               cell_rhs = 0;
 *   
 *               fe_values.reinit(cell);
 *   
 * @endcode
 * 
 * Copying old solution values
 * 
 * @code
 *               fe_values[phase_parameter].get_function_values(conv_solution_np,old_newton_solution_values_p);
 *               fe_values[phase_parameter].get_function_gradients(conv_solution_np,old_newton_solution_gradients_p);
 *               fe_values[temperature].get_function_values(conv_solution_np,old_newton_solution_values_t);
 *               fe_values[temperature].get_function_gradients(conv_solution_np,old_newton_solution_gradients_t);
 *               fe_values[phase_parameter].get_function_values(old_solution_np,old_time_solution_values_p);
 *               fe_values[phase_parameter].get_function_gradients(old_solution_np,old_time_solution_gradients_p);
 *               fe_values[temperature].get_function_values(old_solution_np,old_time_solution_values_t);
 *               fe_values[temperature].get_function_gradients(old_solution_np,old_time_solution_gradients_t);
 *   
 *               for (unsigned int q = 0; q < n_q_points; ++q){
 *                   double khi = get_random_number();
 * @endcode
 * 
 * Old solution values
 * 
 * @code
 *                   double p_on = old_newton_solution_values_p[q]; //old newton solution
 *                   auto grad_p_on = old_newton_solution_gradients_p[q];
 *                   double p_ot = old_time_solution_values_p[q]; //old time step solution
 *                   auto grad_p_ot = old_time_solution_gradients_p[q];
 *                   double t_on = old_newton_solution_values_t[q];
 *                   auto grad_t_on = old_newton_solution_gradients_t[q];
 *                   double t_ot = old_time_solution_values_t[q];
 *                   auto grad_t_ot = old_time_solution_gradients_t[q];
 *                   for (unsigned int i = 0; i < dofs_per_cell; ++i){
 * @endcode
 * 
 * Shape Functions
 * 
 * @code
 *                       double psi_i = fe_values[phase_parameter].value(i,q);
 *                       auto grad_psi_i = fe_values[phase_parameter].gradient(i,q);
 *                       double phi_i = fe_values[temperature].value(i,q);
 *                       auto grad_phi_i = fe_values[temperature].gradient(i,q);
 *                       for (unsigned int j = 0; j < dofs_per_cell; ++j){
 * @endcode
 * 
 * Shape Functions
 * 
 * @code
 *                           double psi_j = fe_values[phase_parameter].value(j,q);
 *                           auto grad_psi_j = fe_values[phase_parameter].gradient(j,q);
 *                           double phi_j = fe_values[temperature].value(j,q);
 *                           auto grad_phi_j = fe_values[temperature].gradient(j,q);
 *   
 *                           double mp = psi_i*(tau*psi_j);
 *                           double kp = grad_psi_i*(std::pow(epsilon,2)*grad_psi_j);
 *                           double m  = (alpha/M_PI)*std::atan(gamma*(t_eq - t_on));
 *                           double t1 = (1-p_on)*(p_on-0.5+m);
 *                           double t2 = -(p_on)*(p_on-0.5+m);
 *                           double t3 = (p_on)*(1-p_on);
 *                           double nl_p = psi_i*((t1+t2+t3)*psi_j);
 * @endcode
 * 
 * Adding random noise at the interface
 * 
 * @code
 *                           nl_p -= a*khi*psi_i*((1.0 - 2*(p_on))*psi_j);
 *                           double f1_p= mp + time_step*theta*kp - time_step*theta*nl_p; // doh f1 by doh p (first Jacobian terms)
 *   
 *                           double t4 = (p_on)*(1-p_on)*(-(alpha*gamma/(M_PI*(1+std::pow((gamma*(t_eq-t_on)),2)))));
 *                           double nl_t = psi_i*(t4*phi_j);
 *                           double f1_t = -time_step*theta*nl_t; // doh f1 by doh t (second Jacobian terms)
 *   
 *                           double mpt = phi_i*(latent_heat*psi_j);
 *                           double f2_p = -mpt; // doh f2 by doh p (third Jacobian terms)
 *   
 *                           double mt = phi_i*(phi_j);
 *                           double kt = grad_phi_i*(grad_phi_j);
 *                           double f2_t = mt + time_step*theta*kt; // doh f2 by doh t (fourth Jacobian terms)
 *   
 * @endcode
 * 
 * Assembling Jacobian matrix
 * 
 * @code
 *                           cell_matrix(i,j) += (f1_p + f1_t + f2_p + f2_t)*fe_values.JxW(q);
 *   
 *                       }
 * @endcode
 * 
 * Finding f1 and f2 at previous iteration for rhs vector
 * 
 * @code
 *                       double mp_n = psi_i*(tau*p_on);
 *                       double kp_n = grad_psi_i*(std::pow(epsilon,2)*grad_p_on);
 *                       double m_n = (alpha/M_PI)*std::atan(gamma*(t_eq-t_on));
 *                       double nl_n = psi_i*((p_on)*(1-p_on)*(p_on-0.5+m_n));
 *                       double mp_t = psi_i*(tau*p_ot);
 *                       double kp_t = grad_psi_i*(tau*grad_p_ot);
 *                       double m_t = (alpha/M_PI)*std::atan(gamma*(t_eq-t_ot));
 *                       double nl_t = psi_i*(p_ot)*(1-p_ot)*(p_ot-0.5+m_t);
 * @endcode
 * 
 * Adding random noise at the interface
 * 
 * @code
 *                       nl_n -= psi_i*(a*khi*(p_on)*(1-p_on));
 *                       nl_t -= psi_i*(a*khi*(p_ot)*(1-p_ot));
 *   
 *                       double f1n = mp_n + time_step*theta*kp_n - time_step*theta*nl_n - mp_t + time_step*(1-theta)*kp_t - time_step*(1-theta)*nl_t; //f1 at last newton iteration
 *   
 *                       double mt_n = phi_i*(t_on);
 *                       double kt_n = grad_phi_i*(grad_t_on);
 *                       double mpt_n = phi_i*(latent_heat*p_on);
 *                       double mt_t = phi_i*(t_ot);
 *                       double kt_t = grad_phi_i*(grad_t_ot);
 *                       double mpt_t = phi_i*(latent_heat*p_ot);
 *   
 *                       double f2n = mt_n + time_step*theta*kt_n - mpt_n - mt_t + time_step*(1-theta)*kt_t + mpt_t; //f2 at last newton iteration
 *   
 * @endcode
 * 
 * Assembling RHS vector
 * 
 * @code
 *                       cell_rhs(i) -= (f1n + f2n)*fe_values.JxW(q);
 *                   }
 *               }
 *   
 *               cell->get_dof_indices(local_dof_indices);
 *               for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               {
 *                   for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                       jacobian_matrix.add(local_dof_indices[i],
 *                                           local_dof_indices[j],
 *                                           cell_matrix(i, j));
 *                   system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *               }
 *           }
 *       }
 *   
 *       jacobian_matrix.compress(VectorOperation::add);
 *       system_rhs.compress(VectorOperation::add);
 *   
 * @endcode
 * 
 * Applying zero BC
 * 
 * @code
 *       std::map<types::global_dof_index, double> boundary_values;
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                1,
 *                                                Functions::ZeroFunction<2>(2),
 *                                                boundary_values);
 *       MatrixTools::apply_boundary_values(boundary_values,
 *                                          jacobian_matrix,
 *                                          solution_update,
 *                                          system_rhs, false);
 *   
 *       jacobian_matrix.compress(VectorOperation::insert);
 *   
 *       system_rhs.compress(VectorOperation::insert);
 *   }
 * @endcode
 
 
<a name="ann-get_random_number.cpp"></a>
<h1>Annotated version of get_random_number.cpp</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Umair Hussain
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "PhaseFieldSolver.h"
 *   #include <random>
 *   
 *   float PhaseFieldSolver::get_random_number()
 *   {
 *       static std::default_random_engine e;
 *       static std::uniform_real_distribution<> dis(-0.5, 0.5); // returns a random number in the range of -0.5 to 0.5
 *       return dis(e);
 *   }
 * @endcode
 
 
<a name="ann-grid_dof.cpp"></a>
<h1>Annotated version of grid_dof.cpp</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Umair Hussain
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "PhaseFieldSolver.h"
 *   
 *   void PhaseFieldSolver::make_grid_and_dofs() {
 * @endcode
 * 
 * Reading mesh
 * 
 * @code
 *       gridin.attach_triangulation(triangulation);
 *       std::ifstream f("mesh/Kobayashi_mesh100x400.msh");
 *       gridin.read_msh(f);
 *   
 *   
 *       GridTools::partition_triangulation(n_mpi_processes, triangulation);
 *       dof_handler.distribute_dofs(fe);
 *       DoFRenumbering::subdomain_wise(dof_handler);
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *   
 *       const std::vector<IndexSet> locally_owned_dofs_per_proc =
 *               DoFTools::locally_owned_dofs_per_subdomain(dof_handler);
 *       const IndexSet locally_owned_dofs =
 *               locally_owned_dofs_per_proc[this_mpi_process];
 *       jacobian_matrix.reinit(locally_owned_dofs,
 *                            locally_owned_dofs,
 *                            dsp,
 *                            mpi_communicator);
 *       old_solution.reinit(locally_owned_dofs, mpi_communicator);
 *       system_rhs.reinit(locally_owned_dofs, mpi_communicator);
 *       conv_solution.reinit(locally_owned_dofs, mpi_communicator);
 *       solution_update.reinit(locally_owned_dofs, mpi_communicator);
 *   
 *       conv_solution_np.reinit(dof_handler.n_dofs());
 *       old_solution_np.reinit(dof_handler.n_dofs());
 *   }
 * @endcode
 
 
<a name="ann-main.cpp"></a>
<h1>Annotated version of main.cpp</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Umair Hussain
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include <iostream>
 *   #include "PhaseFieldSolver.h"
 *   
 *   int main(int argc, char **argv) {
 *       Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
 *       PhaseFieldSolver phasefieldsolver;
 *       phasefieldsolver.run();
 *       return 0;
 *   }
 * @endcode
 
 
<a name="ann-output_results.cpp"></a>
<h1>Annotated version of output_results.cpp</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Umair Hussain
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "PhaseFieldSolver.h"
 *   
 *   void PhaseFieldSolver::output_results(const unsigned int timestep_number) const {
 *       const Vector<double> localized_solution(old_solution);
 *   
 * @endcode
 * 
 * using only one process to output the result
 * 
 * @code
 *       if (this_mpi_process == 0)
 *       {
 *           DataOut<2> data_out;
 *           data_out.attach_dof_handler(dof_handler);
 *   
 *           std::vector<std::string> solution_names;
 *           solution_names.emplace_back ("p");
 *           solution_names.emplace_back ("T");
 *   
 *           data_out.add_data_vector(localized_solution, solution_names);
 *           const std::string filename =
 *                   "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtk";
 *           DataOutBase::VtkFlags vtk_flags;
 *           vtk_flags.compression_level =
 *                   DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
 *           data_out.set_flags(vtk_flags);
 *           std::ofstream output(filename);
 *   
 *           data_out.build_patches();
 *           data_out.write_vtk(output);
 *       }
 *   }
 * @endcode
 
 
<a name="ann-run.cpp"></a>
<h1>Annotated version of run.cpp</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Umair Hussain
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "PhaseFieldSolver.h"
 *   #include <time.h>
 *   
 *   void PhaseFieldSolver::run() {
 *       make_grid_and_dofs();
 *       pcout << "Processors used: " << n_mpi_processes << std::endl;
 *       pcout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
 *             << " (by partition:";
 *       for (unsigned int p = 0; p < n_mpi_processes; ++p)
 *           pcout << (p == 0 ? ' ' : '+')
 *                 << (DoFTools::count_dofs_with_subdomain_association(dof_handler,
 *                                                                     p));
 *       pcout << ")" << std::endl;
 * @endcode
 * 
 * Initialise the solution
 * 
 * @code
 *       InitialValues initial_value;
 *       VectorTools::interpolate(dof_handler,
 *                                initial_value,
 *                                old_solution);
 *       VectorTools::interpolate(dof_handler,
 *                                initial_value,
 *                                conv_solution);
 * @endcode
 * 
 * Applying Boundary Conditions at t=0
 * 
 * @code
 *       applying_bc();
 * @endcode
 * 
 * Plotting initial solution
 * 
 * @code
 *       output_results(0);
 *   
 * @endcode
 * 
 * Time steps begin here:
 * 
 * @code
 *       unsigned int timestep_number = 1;
 *       for (; time <= final_time; time += time_step, ++timestep_number) {
 *   
 *           pcout << "Time step " << timestep_number << " at t=" << time+time_step
 *                     << std::endl;
 *   
 *           conv_solution.operator=(old_solution); // initialising the newton solution
 *   
 * @endcode
 * 
 * Newton-Raphson iterations begin here:
 * 
 * @code
 *           for (unsigned int it = 1; it <= 100; ++it) {
 *               pcout << "Newton iteration number:" << it << std::endl;
 *   
 *               if (it == 100) {
 *                   pcout << "Convergence Failure!!!!!!!!!!!!!!!" << std::endl;
 *                   std::exit(0);
 *               }
 * @endcode
 * 
 * Saving parallel vectors as non-parallel ones
 * 
 * @code
 *               conv_solution_np = conv_solution;
 *               old_solution_np = old_solution;
 * @endcode
 * 
 * Initialise the delta solution as zero
 * 
 * @code
 *               VectorTools::interpolate(dof_handler,
 *                                        Functions::ZeroFunction<2>(2),
 *                                        solution_update);
 *               solution_update.compress(VectorOperation::insert);
 * @endcode
 * 
 * Assemble Jacobian and Residual
 * 
 * @code
 *               assemble_system();
 * @endcode
 * 
 * Solving to get delta solution
 * 
 * @code
 *               solve();
 * @endcode
 * 
 * Checking for convergence
 * 
 * @code
 *               double residual_norm = system_rhs.l2_norm(); //the norm of residual should converge to zero as the solution converges
 * @endcode
 * 
 * pcout << "Nothing wrong till here!!!!!!" << std::endl;
 * 
 * @code
 *               pcout << "the residual is:" << residual_norm << std::endl;
 *               if (residual_norm <= (1e-4)) {
 *                   pcout << "Solution Converged!" << std::endl;
 *                   break; //Break to next time step if the N-R iterations converge
 *               }
 *           }
 * @endcode
 * 
 * Transfer the converged solution to the old_solution vector to plot output
 * 
 * @code
 *           old_solution.operator=(conv_solution);
 *           old_solution.compress(VectorOperation::insert);
 * @endcode
 * 
 * output the solution at only specific number of time steps
 * 
 * @code
 *           if (timestep_number%10 == 0)
 *               output_results(timestep_number);
 *       }
 *   }
 * @endcode
 
 
<a name="ann-solve.cpp"></a>
<h1>Annotated version of solve.cpp</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Umair Hussain
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "PhaseFieldSolver.h"
 *   
 *   void PhaseFieldSolver::solve(){
 * @endcode
 * 
 * Using a direct parallel solver
 * 
 * @code
 *       SolverControl cn;
 *       PETScWrappers::SparseDirectMUMPS A_direct(cn, mpi_communicator);
 *       A_direct.solve(jacobian_matrix, solution_update, system_rhs);
 * @endcode
 * 
 * Updating the solution by adding the delta solution
 * 
 * @code
 *       conv_solution.add(1, solution_update);
 *       conv_solution.compress(VectorOperation::add);
 *   }
 * @endcode
 
 
*/