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):
Pictures from this code gallery program
Annotated version of README.md
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.
Illustration of the problem
Numerical results
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.
Mesh
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
8-threads
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
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/quadrature_lib.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>
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>
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;
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
system_matrix
for storing left matrix
for storing right matrix
System_rhs, only locally owned cells
Solutions Old Solutions with ghost cells, for output
Old Solutions only with locally owned cells
New Solutions only with locally owned cells
Dynamic assembling of the righthandside terms
the constructor
template <int dim>
LaserHeating<dim>::LaserHeating ()
:
mpi_communicator (MPI_COMM_WORLD),
fe (1),
theta(0.5)
{}
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 ()
{
std::ifstream input_file ("geometry.msh");
grid_in.read_msh (input_file);
pcout << " Number of active cells: "
<< std::endl
<< " Total number of cells: "
<< std::endl;
}
void attach_triangulation(Triangulation< dim, spacedim > &tria)
LaserHeating::setup_system
initialization
template <int dim>
void LaserHeating<dim>::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();
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);
constraints_T.close();
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={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
LaserHeating::assemble_system_matrix
template <int dim>
void LaserHeating<dim>::assemble_system_matrix_init (double time_step)
{
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();
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 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));
local_rho_c_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
phi_i_u *
phi_j_u *
fe_values.JxW (q))
+
time_step * (theta) *
(k_fun_T.value(fe_values.quadrature_point(q)) *
div_phi_i_u *
div_phi_j_u *
fe_values.JxW (q));
local_k_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
phi_i_u *
phi_j_u *
fe_values.JxW (q))
-
time_step * (1.0-theta) *
(k_fun_T.value(fe_values.quadrature_point(q)) *
div_phi_i_u *
div_phi_j_u *
fe_values.JxW (q));
}
local_init_T_rhs(i) += (phi_i_u *
initial_value_func_T.value (fe_values.quadrature_point (q)) *
fe_values.JxW (q));
}
cell->get_dof_indices (local_dof_indices);
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
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)
{
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();
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 *
rhs_func_T_1.value_v2 (fe_values.quadrature_point (q)) *
fe_values.JxW (q))
+
time_step * (1.0 - theta) *
(phi_i_u *
rhs_func_T_2.value_v2 (fe_values.quadrature_point (q)) *
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 ()
{
SolverControl solver_control (1*system_rhs_T.size(),1e-12*system_rhs_T.l2_norm(),
true);
preconditioner.initialize(system_matrix_T,
data);
solver.solve (system_matrix_T,completely_distributed_solution,system_rhs_T,preconditioner);
std::vector< index_type > data
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()
{
only refine mesh 5um above the TiO2 and glass interface
if(cell->is_locally_owned())
{
fe_values.reinit(cell);
if(
std::abs(fe_values.quadrature_point(0)[1]) <= global_film_thickness+5e-6)
{
cell->set_refine_flag();
}
else
{
cell->clear_refine_flag();
}
}
}
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
LaserHeatingh::output_results
template <int dim>
void LaserHeating<dim>::output_results (int output_num) const
{
data_out.add_data_vector (old_solution_T, "T");
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-" +
"." +
".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)
output the overall solution
{
std::vector<std::string> filenames;
filenames.push_back ("solution-" +
"." +
".vtu");
std::ofstream master_output (("solution-" +
".pvtu").c_str());
data_out.write_pvtu_record (master_output,filenames);
}
}
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
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_rhs_T.add(1,dynamic_rhs_T);
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 ();
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;
{
output_results (timestep_number);
}
computing_timer.print_summary ();
computing_timer.reset();
pcout << std::endl;
}
}
}
The main function
int main (int argc, char *argv[])
{
try
{
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:
const unsigned int component = 0) const override;
};
template <int dim>
class InitialValues :
public Function<dim>
{
public:
const unsigned int component = 0) const override;
};
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
#-------------------------------------------------------—
Implementation
template <int dim>
double BoundaryValues<dim>::value (
const Point<dim> &,
const unsigned int ) const
{
return 293;
}
template <int dim>
double InitialValues<dim>::value (
const Point<dim> &,
const unsigned int ) const
{
return 293;
}
#-------------------------------------------------------—
Declaration and Implementation
mass density and heat capacity
template <int dim>
{
public:
const unsigned int component = 0) const override;
};
template <int dim>
const unsigned int ) const
{
p stores the xyz coordinates at each vertex
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;
}
#-------------------------------------------------------—
Declaration and Implementation
thermal conductivity
template <int dim>
{
public:
const unsigned int component = 0) const override;
};
template <int dim>
const unsigned int ) const
{
p stores the xyz coordinates at each vertex
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;
}
Annotated version of globalPara.h
physics constants
double global_PI = 3.1415927;
Laser
double global_Pow_laser = 0.4;
double global_spotsize_at_e_2 = 20e-6;
double global_c_laser = global_spotsize_at_e_2 / 4.0;
double global_c_hwhm = global_c_laser * 2.35482 / 2;
double global_V_scan_x = 10e-3;
double global_init_position_x0 = -50e-6;
material
thin film
double global_rho_Tio2 = 4200;
double global_C_Tio2 = 690;
double global_k_Tio2 = 4.8;
substrate
double global_rho_glass = 2200;
double global_C_glass = 700;
double global_k_glass = 1.8;
double global_film_thickness = 400e-9;
simulation time
double global_simulation_time_step = 1e-5;
double global_simulation_end_time = 100e-6 / global_V_scan_x;
about the MESH
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:
};
#-------------------------------------------------------—
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;
(
(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 > &)
no absorption
1.9.8