This program was contributed by Umair Hussain <husain.umair2010@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Pictures from this code gallery program

Annotated version of README.md

Crystal solidification using phase field modeling

Overview

This code solves the solidification problem based in the famous work by Ryo Kobayashi (1993) [1]. The model is based on the Allen-Cahn [2] phase field equation coupled with the transient heat equation Though we have covered only the isotropic directional solidification from the paper in the results, the same code can be modified and used for other types of solidification problems. Let us quickly go through the governing equations and the boundary conditions solved in this problem.

Problem Definition

\begin{align*} \tau \frac{\partial p}{\partial t} = \nabla \cdot( \left(\epsilon^2\right)\nabla p) + p(1-p)(p-\frac{1}{2}+m) \label{heateq} \\ \frac{\partial T}{\partial t}=\nabla^2T+K\frac{\partial p}{\partial t} \end{align*}

where \(m(T) = \frac{a}{\pi}\tan^{-1}(\gamma(T_e-T))\)

The problem is subjected to the boundary conditions:

\begin{align*} p(0,y,t)= 1 \\ T(0,y,t)= T_\gamma -\Delta T \end{align*}

and the initial conditions:

\begin{align*} p(x,y,0)= 0 \\ T(x,y,0)= T_\gamma -\Delta T \end{align*}

Here, \(\Delta T\) is the degree of undercooling.

Dendritic Growth

Using this code, we have reproduced one of the study from Kobayashi's work regarding the dendritic behaviour during directional solidification. The latent heat parameter 'K' in the equation determines the amount of heat released as the phase solidifies. If this value is high enough, we would observe an unstable interface between the solid and liquid phase, which would lead to the formation of dendrites as shown in these images. To assist this growth we need to add a random perturbation term ' \(a \chi p (1-p)\)' to the dynamic term of the phase field equation.

Screenshot

References

[1] Kobayashi, R. (1993). Modeling and numerical simulations of dendritic crystal growth. Physica D: Nonlinear Phenomena, 63(3–4), 410–423. https://doi.org/10.1016/0167-2789(93)90120-P

[2] Allen, S. M., & Cahn, J. W. (1979). A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metallurgica, 27(6), 1085–1095. https://doi.org/10.1016/0001-6160(79)90196-2

Annotated version of InitialValues.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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 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
  }

Annotated version of PhaseFieldSolver.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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)
  {}

Annotated version of PhaseFieldSolver.h

  /* -----------------------------------------------------------------------------
  *
  * 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>
  

For Parallel Computation

  #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
  
  };
  

Initial values class

  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

Annotated version of applying_bc.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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;
  

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)

      VectorTools::interpolate_boundary_values (dof_handler,
                                                1,
                                                Functions::ConstantFunction<2>(1., 2),
                                                boundary_values,p_mask);
  

To apply the boundary values only to the solution vector without the Jacobian Matrix and RHS Vector

      for (auto &boundary_value : boundary_values)
          old_solution(boundary_value.first) = boundary_value.second;
  
      old_solution.compress(VectorOperation::insert);
  
  }

Annotated version of assemble_system.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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() {

Separating each variable as a scalar to easily call the respective shape functions

      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);
  

To copy values and gradients of solution from previous iteration Old Newton iteration

      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);

Old time step iteration

      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);
  

Copying old solution values

              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();

Old solution values

                  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){

Shape Functions

                      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){

Shape Functions

                          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);

Adding random noise at the interface

                          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)
  

Assembling Jacobian matrix

                          cell_matrix(i,j) += (f1_p + f1_t + f2_p + f2_t)*fe_values.JxW(q);
  
                      }

Finding f1 and f2 at previous iteration for rhs vector

                      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);

Adding random noise at the interface

                      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
  

Assembling RHS vector

                      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);
  

Applying zero BC

      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);
  }

Annotated version of get_random_number.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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);
  }

Annotated version of grid_dof.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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() {

Reading mesh

      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());
  }

Annotated version of main.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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;
  }

Annotated version of output_results.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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);
  

using only one process to output the result

      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::CompressionLevel::best_speed;
          data_out.set_flags(vtk_flags);
          std::ofstream output(filename);
  
          data_out.build_patches();
          data_out.write_vtk(output);
      }
  }

Annotated version of run.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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;

Initialise the solution

      InitialValues initial_value;
      VectorTools::interpolate(dof_handler,
                               initial_value,
                               old_solution);
      VectorTools::interpolate(dof_handler,
                               initial_value,
                               conv_solution);

Applying Boundary Conditions at t=0

applying_bc();

Plotting initial solution

output_results(0);

Time steps begin here:

      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
  

Newton-Raphson iterations begin here:

          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);
              }

Saving parallel vectors as non-parallel ones

conv_solution_np = conv_solution;

old_solution_np = old_solution;

Initialise the delta solution as zero

              VectorTools::interpolate(dof_handler,
                                       Functions::ZeroFunction<2>(2),
                                       solution_update);
              solution_update.compress(VectorOperation::insert);

Assemble Jacobian and Residual

assemble_system();

Solving to get delta solution

solve();

Checking for convergence

double residual_norm = system_rhs.l2_norm(); //the norm of residual should converge to zero as the solution converges

pcout << "Nothing wrong till here!!!!!!" << std::endl;

              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
              }
          }

Transfer the converged solution to the old_solution vector to plot output

old_solution.operator=(conv_solution);

old_solution.compress(VectorOperation::insert);

output the solution at only specific number of time steps

          if (timestep_number%10 == 0)
              output_results(timestep_number);
      }
  }

Annotated version of solve.cpp

  /* -----------------------------------------------------------------------------
   *
   * 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(){

Using a direct parallel solver

      SolverControl cn;
      PETScWrappers::SparseDirectMUMPS A_direct(cn, mpi_communicator);
      A_direct.solve(jacobian_matrix, solution_update, system_rhs);

Updating the solution by adding the delta solution

      conv_solution.add(1, solution_update);
      conv_solution.compress(VectorOperation::add);
  }