Reference documentation for deal.II version GIT relicensing-1182-g1782d71672 2024-07-22 00:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
No Matches
The 'Crystal growth problem using phase field modeling' code gallery program

This program was contributed by Umair Hussain <>.
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

Crystal solidification using phase field modeling


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.



[1] Kobayashi, R. (1993). Modeling and numerical simulations of dendritic crystal growth. Physica D: Nonlinear Phenomena, 63(3–4), 410–423.

[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.

Annotated version of InitialValues.cpp

  #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

  #include "PhaseFieldSolver.h"
  : 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.)
  , alpha(0.9)
  , t_eq(1.)
  , a(0.01)
Definition fe_q.h:554
friend class Tensor
Definition tensor.h:882
STL namespace.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

Annotated version of PhaseFieldSolver.h

  #include <deal.II/base/function.h>
  #include <deal.II/base/utilities.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/precondition.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/grid/tria.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/mpi.h>
  #include <deal.II/lac/petsc_vector.h>
  #include <deal.II/grid/grid_tools.h>
  #include <deal.II/dofs/dof_renumbering.h>
  #include <deal.II/distributed/solution_transfer.h>
  using namespace dealii;
  void run();
  void solve();
  void output_results(const unsigned int timestep_number) const;
  double compute_residual();
  void applying_bc();
  MPI_Comm mpi_communicator;
  const unsigned int n_mpi_processes;
  const unsigned int this_mpi_process;
  DoFHandler<2> dof_handler;
  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
  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>
  virtual void vector_value(const Point<2> &p,
  Vector<double> & value) const override;
Definition point.h:111

Annotated version of applying_bc.cpp

  #include "PhaseFieldSolver.h"
  FEValues<2> fe_values(fe,
  ComponentMask p_mask = fe.component_mask(phase_parameter);
  ComponentMask t_mask = fe.component_mask(temperature);
  std::map<types::global_dof_index,double> boundary_values;
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.

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)

void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})

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

Annotated version of assemble_system.cpp

Created by ubuntu on 1/23/21.

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

  FEValues<2> fe_values(fe,
  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);
  for (const auto &cell : dof_handler.active_cell_iterators()){
  if (cell->subdomain_id() == this_mpi_process) {
  cell_matrix = 0;

Copying old solution values

  for (unsigned int q = 0; q < n_q_points; ++q){

Old solution values

  double p_on = old_newton_solution_values_p[q]; //old newton solution
  double p_ot = old_time_solution_values_p[q]; //old time step solution
  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);
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)

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 f2_p = -mpt; // doh f2 by doh p (third Jacobian terms)
  double mt = phi_i*(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 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 mt_t = phi_i*(t_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);
  for (unsigned int i = 0; i < dofs_per_cell; ++i)
  for (unsigned int j = 0; j < dofs_per_cell; ++j)
  system_rhs(local_dof_indices[i]) += cell_rhs(i);
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74

Applying zero BC

  std::map<types::global_dof_index, double> boundary_values;
  system_rhs, false);
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)

Annotated version of get_random_number.cpp

  #include "PhaseFieldSolver.h"
  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

  #include "PhaseFieldSolver.h"

Reading mesh

  std::ifstream f("mesh/Kobayashi_mesh100x400.msh");
  DynamicSparsityPattern dsp(dof_handler.n_dofs());
  const std::vector<IndexSet> locally_owned_dofs_per_proc =
  const IndexSet locally_owned_dofs =
  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);
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void subdomain_wise(DoFHandler< dim, spacedim > &dof_handler)
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
void partition_triangulation(const unsigned int n_partitions, Triangulation< dim, spacedim > &triangulation, const SparsityTools::Partitioner partitioner=SparsityTools::Partitioner::metis)

Annotated version of main.cpp

Annotated version of output_results.cpp

  #include "PhaseFieldSolver.h"
  void PhaseFieldSolver::output_results(const unsigned int timestep_number) const {

using only one process to output the result

  if (this_mpi_process == 0)
  DataOut<2> data_out;
  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";
  vtk_flags.compression_level =
  std::ofstream output(filename);
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)

Annotated version of run.cpp

  #include "PhaseFieldSolver.h"
  #include <time.h>
  void PhaseFieldSolver::run() {
  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,
  pcout << ")" << std::endl;

Initialise the solution

  InitialValues initial_value;
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={})

Applying Boundary Conditions at t=0

Plotting initial solution

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;

Saving parallel vectors as non-parallel ones

Initialise the delta solution as zero

Assemble Jacobian and Residual

Solving to get delta solution


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

output the solution at only specific number of time steps

Annotated version of solve.cpp

  #include "PhaseFieldSolver.h"
  void PhaseFieldSolver::solve(){

Using a direct parallel solver

Updating the solution by adding the delta solution