Reference documentation for deal.II version 9.5.0
Searching...
No Matches
The 'Distributed moving laser heating' code gallery program

This program was contributed by Hongfeng Ma <hongfeng.mark.ma@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):

# Introduction

In this tutorial, the studied problem is to simulate temperature distributions of a sample under a moving laser. Light penetrates the substrate without loss. The top-covered thin-film is, however, a light absorbing material. For simplicity, the thin-film is assumed to be TiO $$_2$$ mixed with silver nanoparticles, which supports light heating of the material by absorbing light energy. For this tutorial, we only consider the isotropic absorption. Figure fgr:s1} illustrates the sketch of the problem. The absorption coefficient is assumed to be $$10^4 m^{-1}$$. The substrate is glass. The thickness of the thin-film is assumed to be 400 nm. The spot size at the top of thin film is $$20 \mu m$$ at $$e^{-2}$$. The writing speed is assumed to be 10 mm/s. The laser power is 0.4 W. The time step is set as 10 $$\mu s$$. The initial position of laser center is $$-50 \mu m$$ with 50 $$\mu m$$ distance away from the left boundary to avoid boundary effects.

# Discretization of the non-uniform isotropic heat equation

In general, the non-uniform isotropic heat equation is as following

\begin{align*} \rho C_m \frac{\partial T}{\partial t} -\nabla \cdot (k\nabla T) = f(\textbf{x}) \end{align*}

Now, we discretize the equation in time with the theta-scheme as

\begin{align*} \rho C_m \frac{T^{n} - T^{n-1}}{dt} - [(1-\theta) \nabla \cdot (k \nabla T^{n-1}) + \theta \nabla \cdot (k \nabla T^n)] = (1-\theta) f^{n-1}(\textbf{x}) + \theta f^n(\textbf{x}) \end{align*}

where $$\theta$$ is a parameter; if $$\theta = 0 (\text{or} = 1)$$, it becomes forward (backward) explicit Euler method; the Crank-Nicolson method is when $$\theta = 0.5$$. Integrating by the test function $$T^{'}$$ and do integration by parts

\begin{align*} \int T^{'} \nabla \cdot (k \nabla T) dV = T^{'} k \nabla T |^a_b - \int k \nabla T \cdot \nabla T^{'} dV \end{align*}

since the test function is time invariant (because the grid is not changed), we have $$T^{'n}$$ = $$T^{'n-1}$$.

\begin{align*} T &= \sum{_i} {u_i} T^{'}_{i} \\ \int T T^{'}_{i} dV &= u_{i} \end{align*}

\noindent let

\begin{align*} M &= \int \rho C_m T^{'_i} T^{'_j} dV \\ A & = \int k \nabla T^{'_i} \cdot \nabla T^{'_j} dV \\ F^n & = \int T' f^{n}(\textbf{x}) \end{align*}

we have the following term,

\begin{align*} \int T^{'} \rho C_m [T^{n} - T^{n-1}] - dt \int T^{'} [(1-\theta) \nabla \cdot (k \nabla T^{n-1}) + \theta \nabla \cdot (k \nabla T^n)] \\ = dt \int T^{'} (1-\theta) f^{n-1}(\textbf{x}) + dt \int T^{'} \theta f^n(\textbf{x}) \end{align*}

\begin{align*} M U^n - M U^{n-1} - dt \int T^{'} [(1-\theta) \nabla \cdot (k \nabla T^{n-1}) + \theta \nabla \cdot (k \nabla T^n)] \\ = dt \int T^{'} (1-\theta) f^{n-1}(\textbf{x}) + dt \int T^{'} \theta f^n(\textbf{x}) \end{align*}

\begin{align*} M U^n - M U^{n-1} + dt [(1-\theta) A U^{n-1} + \theta A U^n] \\ = dt (1-\theta) F^{n-1} + dt \theta F^{n} \end{align*}

the final equation becomes

\begin{align*} (M + dt \theta A) U^n = [M - dt (1-\theta) A] U^{n-1} + dt (1-\theta) F^{n-1} + dt \theta F^{n} \end{align*}

# Initial temperature

The initial temperature can be interpolated over each vertex as follows,

\begin{align*} M_0 &= \int T^{'_i} T^{'_j} dV \\ T_0 &= \sum_i u_i T^{'i} = g_0(x) \\ M_0 U &= \int g_0(\textbf{x}) T^{'i} dV \end{align*}

which is robust for general use. In fact, Deal.II provides a function (VectorTools::interpolate) doing the same thing, which is, however, may not necessary work for parallel version.

## Results

To simplify the question, the heat equation is solved in two-dimensions (x-y) by assuming that the z-axis is homogeneous. Following is part of the running results in 4-threads:



Solving problem in 2 space dimensions.
Number of active cells: 6567
Total number of cells: 6567
Number of degrees of freedom: 11185
9 CG iterations needed to obtain convergence.
initial convergence value = nan
final convergence value = 2.31623e-20

Time step 1 at t=1e-05 time_step = 1e-05
80 CG iterations needed to obtain convergence.
initial convergence value = nan
final convergence value = 1.66925e-13

+------------------------------------------+------------+------+
| Total wallclock time elapsed since start  |     2.98s |      |
|                               |           |           |      |
| Section                       | no. calls |  wall time| %    |
+------------------------------+-----------+------------+------+
| assemble_rhs_T()              |         1 |    0.107s | 3.6% |
| assemble_system_matrix_init() |         1 |    0.245s | 8.2% |
| make_grid()                   |         1 |     1.11s |  37% |
| refine_mesh_at_beginning      |         1 |    0.652s |  22% |
| setup_system()                |         1 |    0.276s | 9.3% |
| solve_T()                     |         2 |    0.426s |  14% |
+-------------------------------+-----------+-----------+------+

Time step 2 at t=2e-05 time_step = 1e-05
79 CG iterations needed to obtain convergence.
initial convergence value = nan
final convergence value = 2.06942e-13

+------------------------------------------+------------+------+
| Total wallclock time elapsed since start |     0.293s |      |
|                                          |            |      |
| Section                      | no. calls |  wall time | %    |
+---------------------------------+--------+------------+------+
| assemble_rhs_T()             |         1 |    0.0969s |  33% |
| solve_T()                    |         1 |     0.161s |  55% |
+------------------------------+-----------+------------+------+

Time step 3 at t=3e-05 time_step = 1e-05
80 CG iterations needed to obtain convergence.
initial convergence value = nan
final convergence value = 1.71207e-13




## Temperature distribution

The colors stand for different cores.

# References

@article{ma2021numerical, title={Numerical study of laser micro-and nano-processing of nanocomposite porous materials}, author={Ma, Hongfeng}, journal={arXiv preprint arXiv:2103.07334}, year={2021} }

@article{ma2019laser, title={Laser-generated ag nanoparticles in mesoporous tio2 films: Formation processes and modeling-based size prediction}, author={Ma, Hongfeng and Bakhti, Said and Rudenko, Anton and Vocanson, Francis and Slaughter, Daniel S and Destouches, Nathalie and Itina, Tatiana E}, journal={The Journal of Physical Chemistry C}, volume={123}, number={42}, pages={25898–25907}, year={2019}, publisher={ACS Publications} }

@article{dealII93, title = {The \texttt{deal.II} Library, Version 9.3}, author = {Daniel Arndt and Wolfgang Bangerth and Bruno Blais and Marc Fehling and Rene Gassm{"o}ller and Timo Heister and Luca Heltai and Uwe K{"o}cher and Martin Kronbichler and Matthias Maier and Peter Munch and Jean-Paul Pelteret and Sebastian Proell and Konrad Simon and Bruno Turcksin and David Wells and Jiaqi Zhang}, journal = {Journal of Numerical Mathematics}, year = {2021, accepted for publication}, url = {https://dealii.org/deal93-preprint.pdf} }

@inproceedings{crank1947practical, title={A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type}, author={Crank, John and Nicolson, Phyllis}, booktitle={Mathematical Proceedings of the Cambridge Philosophical Society}, volume={43}, number={1}, pages={50–67}, year={1947}, organization={Cambridge University Press} }

# Annotated version of Distributed_Moving_Laser_Heating.cc

/* ---------------------------------------------------------------------
*
* Copyright (C) 2021 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Authors: Hongfeng Ma, 2021.
*/

### Include files

All the include files have already been discussed in previous tutorials.

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/base/function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>

#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/trilinos_precondition.h>

#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/generic_linear_algebra.h>

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/distributed/tria.h>

#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>

#include <deal.II/base/logstream.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/base/timer.h>

const ::Triangulation< dim, spacedim > & tria

grid refine

#include <deal.II/numerics/error_estimator.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>

remember to use the namespace dealii before defining a class that is inheritaged from Function<dim>

using namespace dealii;

In general, it is more clear to separate boundary and initial condictions, as well as the right hand side function from a file holding all the things. To do so, in this work, the globalPara.h file defines physical constants, laser parameters, and heat characteristics of materials involved. Boundary and initial conditions are defined in boundaryInit.h. The rightHandSide.h defines the heat source, which in this work is a moving Gaussian beam.

#ifndef GLOBAL_PARA
#define GLOBAL_PARA
#include "./globalPara.h"
#include "./boundaryInit.h"
#include "./rightHandSide.h"
#endif

Now the main class is defined as following

template <int dim>
class LaserHeating
{
public:
LaserHeating ();
~LaserHeating ();
void run ();

private:
void make_grid ();
void setup_system();

void assemble_system_matrix_init (double time_step);
void dynamic_assemble_rhs_T (double time, double time_step);

void solve_T ();

void refine_mesh();
void output_results (int output_num) const;

MPI_Comm mpi_communicator;

FE_Q<dim> fe;
DoFHandler<dim> dof_handler;

AffineConstraints<double> constraints_T;

Definition fe_q.h:551
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

system_matrix

TrilinosWrappers::SparseMatrix system_matrix_T;

for storing left matrix

TrilinosWrappers::SparseMatrix left_system_matrix_T;

for storing right matrix

TrilinosWrappers::SparseMatrix right_system_matrix_T;

System_rhs, only locally owned cells

Solutions Old Solutions with ghost cells, for output

TrilinosWrappers::MPI::Vector old_solution_T;

Old Solutions only with locally owned cells

TrilinosWrappers::MPI::Vector old_solution_T_cal;

New Solutions only with locally owned cells

TrilinosWrappers::MPI::Vector new_solution_T;

Dynamic assembling of the righthandside terms

IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;

TimerOutput computing_timer;

double theta;

};

the constructor

template <int dim>
LaserHeating<dim>::LaserHeating ()
:
mpi_communicator (MPI_COMM_WORLD),
triangulation (mpi_communicator),
fe (1),
dof_handler (triangulation),
pcout (std::cout,Utilities::MPI::this_mpi_process(mpi_communicator) == 0),
computing_timer (mpi_communicator, pcout, TimerOutput::summary, TimerOutput::wall_times),
theta(0.5)
{}

STL namespace.

the destructor

template <int dim>
LaserHeating<dim>::~LaserHeating ()
{
dof_handler.clear();
}

#### LaserHeating::make_grid

make grid by importing msh file, and rescale

template <int dim>
void LaserHeating<dim>::make_grid ()
{
TimerOutput::Scope t(computing_timer,"make_grid()");

GridIn<dim> grid_in;
std::ifstream input_file ("geometry.msh");

pcout << " Number of active cells: "
<< triangulation.n_global_active_cells()
<< std::endl
<< " Total number of cells: "
<< triangulation.n_cells()
<< std::endl;
}

void attach_triangulation(Triangulation< dim, spacedim > &tria)
Definition grid_in.cc:156
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)

#### LaserHeating::setup_system

initialization

template <int dim>
void LaserHeating<dim>::setup_system ()
{

TimerOutput::Scope t(computing_timer,"setup_system()");

dof_handler.distribute_dofs (fe);

pcout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;

locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_dofs);

IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)

we want to output solution, so here should have ghost cells

old_solution_T.reinit (locally_owned_dofs,locally_relevant_dofs,mpi_communicator);

locally owned cells

old_solution_T_cal.reinit (locally_owned_dofs,mpi_communicator);
new_solution_T.reinit (locally_owned_dofs,mpi_communicator);
dynamic_rhs_T.reinit (locally_owned_dofs,mpi_communicator);
system_rhs_T.reinit (locally_owned_dofs,mpi_communicator);

constraints_T.clear();
constraints_T.reinit (locally_relevant_dofs);
DoFTools::make_hanging_node_constraints (dof_handler, constraints_T);
constraints_T.close();

const unsigned int myid = Utilities::MPI::this_mpi_process (mpi_communicator);
DynamicSparsityPattern dsp_T(locally_relevant_dofs);
DoFTools::make_sparsity_pattern (dof_handler, dsp_T, constraints_T, false, myid);

locally_owned_dofs,
mpi_communicator,
locally_relevant_dofs);

left_system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
right_system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);
system_matrix_T.reinit (locally_owned_dofs,locally_owned_dofs,dsp_T,mpi_communicator);

}

void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:161

#### LaserHeating::assemble_system_matrix

template <int dim>
void LaserHeating<dim>::assemble_system_matrix_init (double time_step)
{

TimerOutput::Scope t(computing_timer,"assemble_system_matrix_init()");

const InitialValues<dim> initial_value_func_T;

const RhoC<dim> rho_C_fun_T;
const K_T<dim> k_fun_T;

const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();

FullMatrix<double> local_init_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> local_init_T_rhs (dofs_per_cell);

FullMatrix<double> local_rho_c_T_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> local_k_T_matrix (dofs_per_cell, dofs_per_cell);

std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

system_matrix_T = 0;
left_system_matrix_T = 0;
right_system_matrix_T = 0;

system_rhs_T = 0;

for (const auto &cell : dof_handler.active_cell_iterators())
if(cell->is_locally_owned())
{

fe_values.reinit (cell);
local_init_matrix = 0;

local_rho_c_T_matrix = 0;
local_k_T_matrix = 0;

local_init_T_rhs = 0;

for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim> div_phi_i_u = fe_values.shape_grad (i,q);
const double phi_i_u = fe_values.shape_value (i,q);

for (unsigned int j=0; j<dofs_per_cell; ++j)
{

const double phi_j_u = fe_values.shape_value (j,q);

local_init_matrix(i,j) += (phi_i_u *
phi_j_u *
fe_values.JxW (q));

phi_i_u *
phi_j_u *
fe_values.JxW (q))
+
time_step * (theta) *
div_phi_i_u *
div_phi_j_u *
fe_values.JxW (q));

phi_i_u *
phi_j_u *
fe_values.JxW (q))
-
time_step * (1.0-theta) *
div_phi_i_u *
div_phi_j_u *
fe_values.JxW (q));

}

local_init_T_rhs(i) += (phi_i_u *
fe_values.JxW (q));

}

cell->get_dof_indices (local_dof_indices);

@ update_values
Shape function values.
@ update_JxW_values

copy to system_matrix_T and system_rhs_T for projecting initial values

constraints_T.distribute_local_to_global(local_init_matrix,
local_init_T_rhs,
local_dof_indices,
system_matrix_T,
system_rhs_T);

store M + dt*theta*A as the left_system_matrix

constraints_T.distribute_local_to_global(local_rho_c_T_matrix,
local_dof_indices,
left_system_matrix_T);

store M - dt*(1-theta)*A as the right_system_matrix

constraints_T.distribute_local_to_global(local_k_T_matrix,
local_dof_indices,
right_system_matrix_T);

}

}

#### LaserHeating::dynamic_assemble_rhs_T

The right hand side is assembled each time during running, which is necessary as the laser source is moving. To separate the heat source and the right hand side assembling, the right hand side function is defined as RightHandside<dim>.

template <int dim>
void LaserHeating<dim>::dynamic_assemble_rhs_T (double time, double time_step)
{

TimerOutput::Scope t(computing_timer,"assemble_rhs_T()");

RightHandside<dim> rhs_func_T_1;
rhs_func_T_1.set_time(time);

RightHandside<dim> rhs_func_T_2;
rhs_func_T_2.set_time(time-time_step);

const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();

Vector<double> local_rhs_vector_T (dofs_per_cell);

std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

std::vector<double> Rnp_cal_assemble (n_q_points);

dynamic_rhs_T = 0 ;

cell = dof_handler.begin_active(),
endc = dof_handler.end();

for (; cell!=endc; ++cell)
if(cell->is_locally_owned())
{
fe_values.reinit (cell);
local_rhs_vector_T = 0;

for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double phi_i_u = fe_values.shape_value (i,q);

local_rhs_vector_T(i) += time_step * theta *
(phi_i_u *
fe_values.JxW (q))
+
time_step * (1.0 - theta) *
(phi_i_u *
fe_values.JxW (q));
}

cell->get_dof_indices (local_dof_indices);

constraints_T.distribute_local_to_global(local_rhs_vector_T,
local_dof_indices,
dynamic_rhs_T);

}

}

typename ActiveSelector::active_cell_iterator active_cell_iterator

#### LaserHeating::solve

Solving the equation is direct. Recall that we have defined several matrices and vectors, to avoid ambiguous, here, we only use system_matrix_T as the system matrix, system_rhs_T as the system right hand side. The vector completely_distributed_solution is used to store the obtained solution.

template <int dim>
void LaserHeating<dim>::solve_T ()
{
TimerOutput::Scope t(computing_timer,"solve_T()");

TrilinosWrappers::MPI::Vector completely_distributed_solution (locally_owned_dofs,mpi_communicator);
SolverControl solver_control (1*system_rhs_T.size(),1e-12*system_rhs_T.l2_norm(),true);

TrilinosWrappers::SolverCG solver (solver_control);

preconditioner.initialize(system_matrix_T,data);

solver.solve (system_matrix_T,completely_distributed_solution,system_rhs_T,preconditioner);

Print the number of iterations by hand.

pcout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence." << std::endl
<< "\t initial convergence value = " << solver_control.initial_value() << std::endl
<< "\t final convergence value = " << solver_control.last_value() << std::endl
<< std::endl;

constraints_T.distribute (completely_distributed_solution);
new_solution_T = completely_distributed_solution;

}

#### LaserHeating::refine_mesh

template <int dim>
void LaserHeating<dim>::refine_mesh()
{

TimerOutput::Scope t(computing_timer,"refine_mesh_at_beginning");

only refine mesh 5um above the TiO2 and glass interface

cell = triangulation.begin_active();
cell != triangulation.end(); ++cell)
if(cell->is_locally_owned())
{
fe_values.reinit(cell);
{
cell->set_refine_flag();
}
else
{
cell->clear_refine_flag();
}
}
triangulation.execute_coarsening_and_refinement();

}

::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)

#### LaserHeatingh::output_results

template <int dim>
void LaserHeating<dim>::output_results (int output_num) const
{

DataOut<dim> data_out;

data_out.attach_dof_handler (dof_handler);

void attach_dof_handler(const DoFHandler< dim, spacedim > &)

the length of output numbering

int step_N = 7;

data_out.build_patches ();

const std::string filename = ("solution-" +
Utilities::int_to_string (output_num, step_N) +
"." +
Utilities::int_to_string (triangulation.locally_owned_subdomain(),4) +
".vtu");
std::ofstream output (filename.c_str());
data_out.write_vtu (output);

std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:471

output the overall solution

if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
{
std::vector<std::string> filenames;
for (unsigned int i=0;i<Utilities::MPI::n_mpi_processes(mpi_communicator);++i)
filenames.push_back ("solution-" +
Utilities::int_to_string (output_num, step_N) +
"." +
".vtu");
std::ofstream master_output (("solution-" +
Utilities::int_to_string (output_num,step_N)+
".pvtu").c_str());
data_out.write_pvtu_record (master_output,filenames);

}

}

unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:150

#### LaserHeating::run

This is the function which has the top-level control over everything. Apart from one line of additional output, it is the same as for the previous example.

template <int dim>
void LaserHeating<dim>::run ()
{
pcout << "Solving problem in " << dim << " space dimensions." << std::endl;

make_grid();
refine_mesh();
setup_system ();
assemble_system_matrix_init (global_simulation_time_step);

projection of initial conditions by solving. solution stored in new_solution_T;

solve_T ();

old_solution_T = new_solution_T;
old_solution_T_cal = new_solution_T;

reinitialization

system_matrix_T = 0;
system_rhs_T = 0;

double time_step = global_simulation_time_step;
double time = 0;
int timestep_number = 0;

output initial values; need ghost cells

output_results (0);

while(time < global_simulation_end_time)
{

time += time_step;
timestep_number ++;

pcout << "Time step " << timestep_number
<< " at t=" << time
<< " time_step = " << time_step
<< std::endl;

the dynamic solving part

{

right_system_matrix_T.vmult(system_rhs_T,old_solution_T_cal);

dynamic_assemble_rhs_T (time,time_step);

system_matrix_T.copy_from (left_system_matrix_T);

{
BoundaryValues<dim> boundary_values_function;
std::map<types::global_dof_index,double> boundary_values;

BOUNDARY_NUM,
boundary_values_function,
boundary_values);

system_matrix_T,
new_solution_T,
system_rhs_T,
false);
}

solve_T ();

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)
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=ComponentMask())

old_solution_T is used for output, holding ghost cells old_solution_T_cal is used for calculation, holding only locally owned cells.

old_solution_T = new_solution_T;
old_solution_T_cal = new_solution_T;

if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 96 && (timestep_number % 50 == 0 ))
{
TimerOutput::Scope t(computing_timer,"output");
output_results (timestep_number);
}

computing_timer.print_summary ();
computing_timer.reset();

pcout << std::endl;

}
}

}

### The main function

int main (int argc, char *argv[])
{
try
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);

LaserHeating<2> laserHeating_2d;
laserHeating_2d.run ();

}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;

return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}

return 0;
}

# Annotated version of boundaryInit.h

#-------------------------------------------------------— #

# This file defines the boundary and initial conditions

# #-------------------------------------------------------—

#ifndef GLOBAL_PARA
#define GLOBAL_PARA
#include "./globalPara.h"
#endif

#-------------------------------------------------------—

# Declaration

template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues () : Function<dim>() {}

virtual double value (const Point<dim> &p,
const unsigned int component = 0) const override;
};

template <int dim>
class InitialValues : public Function<dim>
{
public:
InitialValues () : Function<dim>() {}

virtual double value (const Point<dim> &p,
const unsigned int component = 0) const override;
};

virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Definition point.h:112

#-------------------------------------------------------—

# Implementation

template <int dim>
double BoundaryValues<dim>::value (const Point<dim> &/*p*/,
const unsigned int /*component*/) const
{
return 293;
}

template <int dim>
double InitialValues<dim>::value (const Point<dim> &/*p*/,
const unsigned int /*component*/) const
{
return 293;
}

#-------------------------------------------------------—

# mass density and heat capacity

template <int dim>
class RhoC : public Function<dim>
{
public:
RhoC () : Function<dim>() {}

virtual double value (const Point<dim> &p,
const unsigned int component = 0) const override;
};

template <int dim>
double RhoC<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{

# for 2D problems in xy, we assume the non-uniform is in y-axis.

if ( p[1] >= -global_film_thickness )
{
return global_rho_Tio2 * global_C_Tio2;
}
else
return global_rho_glass * global_C_glass;
}

#-------------------------------------------------------—

# thermal conductivity

template <int dim>
class K_T : public Function<dim>
{
public:
K_T () : Function<dim>() {}

virtual double value (const Point<dim> &p,
const unsigned int component = 0) const override;
};

template <int dim>
double K_T<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{

# for 2D problems in xy, we assume the non-uniform is in y-axis.

if ( p[1] >= -global_film_thickness)
{
return global_k_Tio2;
}
else
return global_k_glass;
}

# physics constants

double global_PI = 3.1415927;

# Laser

double global_Pow_laser = 0.4; // power [W]
double global_spotsize_at_e_2 = 20e-6; // laser spot size at e^(-2) [m]
double global_c_laser = global_spotsize_at_e_2 / 4.0; // C parameter in Gaussian func, [m]
double global_c_hwhm = global_c_laser * 2.35482 / 2; // HWHM, [m]
double global_V_scan_x = 10e-3; // scan speed, [m/s]

double global_init_position_x0 = -50e-6; // initial spot center position

# material

thin film

double global_rho_Tio2 = 4200; // mass density, [kg/m^3]
double global_C_Tio2 = 690; // heat capacity, [J/kg/K]
double global_k_Tio2 = 4.8; // thermal conductivity, [W/m/K]

substrate

double global_rho_glass = 2200;
double global_C_glass = 700;
double global_k_glass = 1.8;

double global_film_thickness = 400e-9; // film thickness, [m]

# simulation time

double global_simulation_time_step = 1e-5; // 10 [us]
double global_simulation_end_time = 100e-6 / global_V_scan_x; // 100 [um] / scan speed

#define BOUNDARY_NUM 11

# Annotated version of rightHandSide.h

#-------------------------------------------------------— #

# This file defines the boundary and initial conditions

# #-------------------------------------------------------—

#ifndef GLOBAL_PARA
#define GLOBAL_PARA
#include "./globalPara.h"
#endif

#-------------------------------------------------------—

# Declaration

template <int dim>
class RightHandside : public Function<dim>
{
public:
RightHandside () : Function<dim>() {}

double value_v2 (const Point<dim> &p);
};

#-------------------------------------------------------—

# Implementation

template <int dim>
double RightHandside<dim>::value_v2 (const Point<dim> &p)
{

double alpha_abs = 1e4;

if(p[1] >= -global_film_thickness)
{

double P00 = global_Pow_laser / global_PI / global_c_laser / global_c_laser / 2.0;

double I00 = P00 * std::exp(-
(
(p[0] - global_V_scan_x * this->get_time()-global_init_position_x0) *
(p[0] - global_V_scan_x * this->get_time()-global_init_position_x0)
) /
(2.0 * global_c_laser * global_c_laser) );

return alpha_abs * I00 * std::exp(-alpha_abs*(0 - p[1]));

}

else
{
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)

return 0;
}

}