Merge pull request #185 from umairhussaincmm/my-code-gallery-project

Request to add my 'Crystal_Growth_Phase_Field_Model' code to the code gallery.

Crystal_Growth_Phase_Field_Model/CMakeLists.txt

@@ -0,0 +1,33 @@
+cmake_minimum_required(VERSION 3.16)
+# Set the name of the project and target:
+SET(TARGET "Kobayashi_Parallel")
+
+#set(CMAKE_CXX_STANDARD 14)
+
+SET(TARGET_SRC
+        main.cpp
+        PhaseFieldSolver.cpp
+        applying_bc.cpp
+        assemble_system.cpp
+        grid_dof.cpp
+        output_results.cpp
+        run.cpp
+        solve.cpp
+        InitialValues.cpp
+        get_random_number.cpp
+        )
+
+FIND_PACKAGE(deal.II 9.5.0 QUIET
+        HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+        )
+IF(NOT ${deal.II_FOUND})
+    MESSAGE(FATAL_ERROR "\n"
+            "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+            "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+            "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+            )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()

Crystal_Growth_Phase_Field_Model/InitialValues.cpp

@@ -0,0 +1,8 @@
+#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
+}

Crystal_Growth_Phase_Field_Model/PhaseFieldSolver.cpp

@@ -0,0 +1,21 @@
+#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)
+{}

Crystal_Growth_Phase_Field_Model/PhaseFieldSolver.h

@@ -0,0 +1,94 @@
+#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

Crystal_Growth_Phase_Field_Model/README.md

@@ -0,0 +1,42 @@
+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
+@f{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}
+@f}
+where $m(t) = \frac{a}{\pi}\tan^{-1}(\gamma(T_e-T))$
+
+The problem is subjected to the boundary conditions:
+@f{align*}{
+p(0,y,t)= 1 \\
+T(0,y,t)= T_\gamma -\Delta T
+@f}
+and the initial conditions:
+@f{align*}{
+p(x,y,0)= 0 \\
+T(x,y,0)= T_\gamma -\Delta T
+@f}
+Here, $\Delta T$ is the degree of undercooling.
+
+### Dendritic Growth
+Using this code, we have reproduced one of the study from Koabayashi'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 perterbation term 
+'$a \chi p (1-p)$' to the dynamic term of the phase field equation.
+
+![Screenshot](./doc/images/K=1v1.2v1.4.gif)
+
+### 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

Crystal_Growth_Phase_Field_Model/applying_bc.cpp

@@ -0,0 +1,29 @@
+#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);
+
+}

Crystal_Growth_Phase_Field_Model/assemble_system.cpp

@@ -0,0 +1,165 @@
+//
+// Created by ubuntu on 1/23/21.
+//
+
+#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);
+}