- Note
- Prerequisite of this tutorial is SCL-6: Heat equation
- Note
- Intended learning outcome:
- first time-dependent problem in MoFEM
- TS solver with implicit scheme
Introduction
This tutorial covers the numerical modelling of propagation of photons in a slab with a hidden object in its centre. An initial, analytical pulse is generated on the mesh and its diffusion is then simulated using the PETSc SNES time solver. The code calculates the fluence rate of photons over an area at each timestep (the photon time-of-flight data) at the opposite end of the incident beam, see Fig. 1. These photons can then be detected by a camera and the hidden object reconstructed via Bayesian machine learning methods. The simulation complements a physical experiment where a camera is used to count photons by providing a wide range of simulations that can be used as a training set for the machine learning algorithms.
Figure 1: Experimental setup for imaging through turbid media lyons_computational_2019
Light propagation in turbid media is associated with two processes: scattering (interaction of photons with the media causing secondary radiation in every direction) and absorption (transformation of the electromagnetic (EM) wave into other forms of energy) [82] . The average distance a photon travels before its trajectory is randomised by interactions with the medium is known as the transport mean free path (TFMP). TFMP can be defined as:
\[
\begin{equation}
\label{Equation:TFMP}
l^* = \frac{1}{\mu_a + \mu'_s},
\end{equation}
\]
where \(\mu_a\) and \(\mu'_s\) are constants known as the absorption and reduced scattering coefficient. We denote the thickness of the medium through which light propagates as \(L\). For this experiment, we will typically use \(L = 5 \, \textrm{cm}\), \(\mu_a = 0.09 \, \textrm{cm}^{-1}\), and \(\mu'_s = 16.5 \, \textrm{cm}^{-1}\), which give \(l^* = 600 \, \mu\textrm{m}\). For \(L \gg l^*\), light transport can be modelled using the photon diffusion equation [82], [81] :
\[
\begin{equation}
\label{Equation:pde}
\frac{1}{c} \frac{\partial \Phi(\mathbf{r}, t)}{\partial t} - D \nabla^2\Phi(\mathbf{r},t)+\mu_a\Phi(\mathbf{r},t) = S(\mathbf{r}, t),
\end{equation}
\]
where \(\Phi(\mathbf{r}, t)\) is the photon fluence rate, i.e. the number of photons crossing a sphere with unit area in unit time, \(\mathbf{r} = [x,y,z]^\intercal\) is the spatial coordinate vector, \(t\) is time, \(c\) is the speed of light in the media, \(D = l^*/3\) is the diffusion coefficient, and \(S(\mathbf{r}, t)\) is the photon source.
For our experiment, we illuminate the slab with a \(120 \ \text{fs} \) ( \(10^{-15} \ \text{s}\)) pulse and record for at least \(5 \ \textrm{ns}\) . Considering an initial incident beam as a boundary condition would require extremely small time steps and a sufficiently fine mesh near the surface illuminated by the source. Instead, we use the analytical solution of the PDE in Eq. \eqref{Equation:pde} in a semi-infinite domain to provide an initial distribution of the photon fluence rate. This then allows the numerical solution to be reliably computed on a relatively coarse mesh with time steps of around \(0.1 \text{ns}\). For an isotropic point source, where we can represent \(S(\mathbf{r}, t)\) as a \(\delta\)-function, the solution of Eq. \eqref{Equation:pde} is [71] :
\[
\begin{equation}
\label{Equation:phi_star}
\Phi^*(\mathbf{r}, t) = c (4\pi D c t)^{-3/2}\exp{\left(-\frac{|\mathbf{r}|^2}{4Dct} - \mu_a ct\right)}.
\end{equation}
\]
Figure 2: (a) Schematic of the image method in a semi-infinite domain. The black circle represents the positive source and the white circle represents the 'negative' one. (b) Plot of the fluence rate generated by the image method at different times along the z=0 line.
To obtain the analytical solution for a point source located in the boundary of the medium, we use the method of images, see Fig. 2. This method sums contributions from a 'positive' source inside the body at depth \(z_0 = 1/\mu'_s\) and a 'negative' source at a distance \(z_0\) from the surface [71] . We denote this fluence rate as \(\Phi^{**}(\mathbf{r}, t)\):
\[
\begin{equation}
\label{eq:phi_double_star}
\Phi^{**}(\mathbf{r}, t) = \Phi^{*}(\mathbf{r} - z_0\mathbf{e}_z, t) - \Phi^{*}(\mathbf{r} + z_0\mathbf{e}_z, t),
\end{equation}
\]
Convolving \(\Phi^{**}(\mathbf{r}, t)\) with an initial condition
\[
\begin{equation}
\label{eq:phi_0}
\Phi_0(\textbf{r}) = \begin{cases}
I, \quad z = 0 \; \text{and}\; \sqrt{x^2+y^2} \leq R \\
0, \quad \text{otherwise},
\end{cases}
\end{equation}
\]
we arrive at the strong form of the problem given by:
\[
{\label{eq:ibvp}}
\begin{cases}
\frac{1}{c} \frac{\partial \Phi(\mathbf{r}, t)}{\partial t} - D \nabla^2\Phi(\mathbf{r},t)+\mu_a\Phi(\mathbf{r},t) = 0, \quad & \mathbf{r} \in \Omega, \; t \in (t_0, T]\\
\Phi(\mathbf{r}, t_0) = \int_{\Omega} \Phi^{**}(\mathbf{r}-\mathbf{q}, t_0)\, \Phi_0(\mathbf{q})\, \textrm{d}\textbf{q}, \quad & \mathbf{r} \in \Omega \\
\Phi(\mathbf{r}, t) = 0, \quad & \mathbf{r} \in \Gamma^\text{int}, \; t \in (t_0, T] \\
\Phi(\mathbf{r}, t) + 2DA\mathbf{n}\cdot\nabla\Phi(\mathbf{r},t) = 0, \quad & \mathbf{r} \in \Gamma^\text{ext}, \; t \in (t_0, T]
\end{cases}
\]
where the third equation is the Dirichlet boundary condition representing the internal perfectly absorbing object,the fourth equation is the Robin boundary condition imposed on the external surface of the domain, \(\mathbf{n}\) is the outward normal, \(t_0\) is the time until which the analytical solution has been calculated, and \(A\approx3\) is a coefficient accounting for the difference of the refractive index between the inside of the turbid domain and non-scattering surrounding medium (air), see [81] .
Using test functions \(\delta\Phi \equiv \delta\Phi(\mathbf{r}, t)\) satisfying \(\delta\Phi = 0\) at \(\Gamma^{\text{int}}\), we arrive at the following weak function after integrating the first line of Eq. \eqref{eq:ibvp} by parts and applying the Robin Boundary condition:
\[
\begin{equation}
\label{Equation:weakform}
\frac{1}{c} \int_\Omega \delta\Phi\frac{\partial \Phi}{\partial t}\, \textrm{d}\Omega + D\int_\Omega \nabla\delta\Phi \cdot \nabla\Phi\, \textrm{d}\Omega + \frac{1}{2A}\int_{\Gamma^\text{ext}} \delta\Phi \Phi\, \textrm{d}\Gamma + \mu_a \int_\Omega \delta\Phi\Phi\, \textrm{d}\Omega = 0 \quad \forall\, \delta\Phi\in H^1_0,
\end{equation}
\]
where we have omitted the spatial and temporal dependencies of \(\Phi(\mathbf{r},t)\) by writing \(\Phi \equiv \Phi(\mathbf{r}, t)\). To solve the initial step, we will require a solution to the least squares problem of the second equation in Eq. \eqref{eq:ibvp}, which we write as follows:
\[
\begin{equation}
\label{Equation:weakform_init}
\int_\Omega \delta\Phi\Phi(t_0)\, \textrm{d}\Omega - \int_\Omega \delta\Phi s\, \textrm{d}\Omega = 0 \quad \forall\, \delta\Phi\in H^1_0.
\end{equation}
\]
Here, we have substituted the integral on the right hand side by \(s\) for notational simplicity. We require \(\Phi\) and \(\delta\Phi\) to be members of a subset of the Sobolev space \(H^1(\Omega)\). This is the space of all square integrable functions with square integrable derivatives. The spaces for \(\Phi\) and \(\delta\Phi\) are defined in Table 1.
Table 1: Functional spaces for the photon fluence rate and test functions
| Field | Space |
| \(\Phi\) | \( \Phi \in \ H^1_0(\Omega) = \{u \in H^1(\Omega) | u = 0 \ \text{on} \ \Gamma_\textrm{int} \} \) |
| \(\delta\Phi\) | \( \delta\Phi \in \ H^1_0(\Omega) \) |
Linearisation
To solve the time-dependent problem, we need to express the time derivative in terms of a discrete change in photon fluence rate, \(\Delta\Phi \in H^1_0(\Omega)\), over an increase in time, \(\Delta t\). To this end, we linearise and solve the equation using the Newton-Raphson method [17]. For this particular problem, the Newton-Raphson method is not mathematically needed, but the linearised form is required for the time solver. We write a residual:
\[
\begin{equation}
\label{Equation:residual}
R(\Phi,\, t,\, \delta\Phi) = 0 \quad \forall\, \delta\Phi\in H^1_0
\end{equation}
\]
and equate the left-hand side of Eq. \eqref{Equation:weakform} to the residual. A first order Taylor series expansion of the residual around \(\Phi_0\) gives:
\[
\begin{equation}
\label{Equation:linearexp}
R(\Phi_0 + \Delta\Phi,\, t,\, \delta\Phi) \approx R(\Phi_0,\, t,\, \delta\Phi) + DR(\Phi_0,\, t, \,\delta\Phi) \left[ \Delta\Phi \right] = 0.
\end{equation}
\]
In the expression above, the second term on the right hand side corresponds to the Gateaux, or directional, derivative of the residual along the direction of \(\Delta\Phi\). It can be expressed as:
\[
\begin{equation}
\label{Equation:gateauxder}
\begin{split}
DR(\Phi_0,\, t, \,\delta\Phi) \left[ \Delta\Phi \right] & = \lim_{\varepsilon\to 0} \frac{R(\Phi_0 + \varepsilon\Delta\Phi,\, t,\, \delta\Phi)-R(\Phi_0, \, t,\, \delta\Phi)}{\varepsilon} = \frac{\textrm{d}}{\textrm{d}\varepsilon} R(\Phi_0 + \varepsilon\Delta\Phi,\, t,\, \delta\Phi) \big|_{\varepsilon = 0} \\
& = \frac{1}{c} \int_\Omega \delta\Phi\frac{\Delta \Phi}{\Delta t}\, \textrm{d}\Omega + D\int_\Omega \nabla\delta\Phi \cdot \nabla\Delta\Phi\, \textrm{d}\Omega + \frac{1}{2A}\int_{\Gamma^\text{ext}} \delta\Phi \Delta\Phi\, \textrm{d}\Gamma + \mu_a \int_\Omega \delta\Phi\Delta\Phi\, \textrm{d}\Omega,
\end{split}
\end{equation}
\]
where we have expressed the time derivative as \(\Delta\Phi/\Delta t\), where \(\Delta t\) can be obtained using PETSc. Rearranging Eq. \eqref{Equation:linearexp}, we can express the problem as:
\[
\begin{equation}
\label{Equation:NRproblem}
DR(\Phi_0,\, t, \,\delta\Phi) \left[ \Delta\Phi \right] = - R(\Phi_0,\, t,\, \delta\Phi) \quad \forall\, \delta\Phi\in H^1_0.
\end{equation}
\]
The Newton-Raphson algorithm then consists of finding the increment \(\Delta\Phi\) which solves Eq. \eqref{Equation:residual}.
Discretisation
Before moving onto the MoFEM implementation of the problem, we first discretise the fields \(\delta\Phi\) and \(\Delta\Phi\) as:
\[
\begin{equation}
\label{Equation:shapefunc}
\delta\Phi^h = \delta\Phi_\beta N_\beta \quad \textrm{and} \quad \Delta\Phi^h = \Delta\Phi_\gamma N_\gamma,
\end{equation}
\]
where \(h\) denotes the element on which the field is being approximated, \(N\) are hierarchical shape functions [94] , and we sum over repeated indices.
We now recast Eq. \eqref{Equation:NRproblem} into matrix form:
\[
\begin{equation}
\label{Equation:NRproblemdisc}
\left[ K_{\beta\gamma} \right] \left[ \Delta\Phi_\gamma \right] = - \left[ R_\beta \right],
\end{equation}
\]
where \(\left[ \Delta\Phi_\gamma\right]\) is the vector of unkown coefficients, \(\left[K_{\beta\gamma}\right]\) is the tangent matrix [17] constructed from Eq. \eqref{Equation:gateauxder} and \(\left[ R_\beta \right]\) is the residual vector. The expressions for the tangent matrix and residual vector are:
\[
\begin{equation}
\label{Equation:tangentmat}
\left[ K_{\beta\gamma} \right] = \frac{1}{c} \int_\Omega N_\beta N_\gamma \Delta t^{-1}\, \textrm{d}\Omega + D\int_\Omega \nabla N_\beta \cdot \nabla N_\gamma\, \textrm{d}\Omega + \frac{1}{2A}\int_{\Gamma^\text{ext}} N_\beta N_\gamma \, \textrm{d}\Gamma + \mu_a \int_\Omega N_\beta N_\gamma \, \textrm{d}\Omega,
\end{equation}
\]
\[
\begin{equation}
\label{Equation:resvec}
\left[ R_\beta \right] = \frac{1}{c} \int_\Omega N_\beta \dot{\Phi}_0 \, \textrm{d}\Omega + D\int_\Omega \nabla N_\beta \cdot \nabla\Phi_0 \, \textrm{d}\Omega + \frac{1}{2A}\int_{\Gamma^\text{ext}} N_\beta \Phi_0 \, \textrm{d}\Gamma + \mu_a \int_\Omega N_\beta \Phi_0 \, \textrm{d}\Omega = 0,
\end{equation}
\]
where \(\dot{\Phi}_0\) is the time derivative of \(\Phi_0\). \(\dot{\Phi}_0\) is calculated by PETSc through a finite difference scheme. For the initial time step, we will also need to discretise the initial condition Eq. \eqref{Equation:weakform_init}. Approximating the initial configuration of the fluence rate as \(\Phi^h(t_0)= \Phi_\alpha (t_0) N_\alpha\), the least squares problem of the initial condition can be written as:
\[
\begin{equation}
\label{eq:init_cond_disc}
\left[ \int_\Omega N_\beta N_\alpha \, \textrm{d}\Omega \right] [\Phi_\alpha{}(t_0)] = \left[ \int_\Omega N_\beta s\, \textrm{d}\Omega \right], \\
\end{equation}
\]
where we input the integrand of the second equation in Eq. \eqref{eq:ibvp}, \(s\), analytically as in Eq. \eqref{Equation:weakform_init}. Solving this equation will give us the initial value of \(\Phi_0\) that can be input to calculate Eqs. \eqref{Equation:tangentmat} and \eqref{Equation:resvec}.
Implementation
Initial Diffusion
We first set up an initial analytical pulse on the mesh using the method of images described above. At this stage, the code generates an initial pulse at on the z face of the mesh using the second equation in Eq. \eqref{eq:ibvp}. The code uses Gauss-Kronrod quadrature from the boost library to perform the integration for the analytical source function. This is done as follows:
double sourceFunctionEval(
const double x,
const double y,
const double z,
const double beam_radius,
const double beam_centre_x,
const double T =
auto phi_pulse = [&](const double r_s, const double phi_s) {
const double xs = r_s * cos(phi_s);
const double ys = r_s * sin(phi_s);
const double P1 = xp * xp + yp * yp + zp1 * zp1;
const double P2 = xp * xp + yp * yp + zp2 * zp2;
return r_s * (exp(-P1 /
A) - exp(-P2 /
A));
};
auto f = [&](const double r_s) {
auto g = [&](
const double phi_s) {
return phi_pulse(r_s, phi_s); };
return gauss_kronrod<double, 15>::integrate(
g, 0, 2 * M_PI, 0, std::numeric_limits<float>::epsilon());
};
gauss_kronrod<double, 15>::integrate(
f, 0,
beam_radius, 0, std::numeric_limits<float>::epsilon());
};
double mu_sp
scattering coefficient (cm^-1)
double flux_magnitude
impulse magnitude
double mu_a
absorption coefficient (cm^-1)
const double v
phase velocity of light in medium (cm/ns)
The operators for solving Eq. \eqref{eq:init_cond_disc} are then assembled and pushed to the pipeline:
};
auto set_domain = [&]() {
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainMass(
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](const double, const double, const double) { return 1; }));
pipeline_mng->getOpDomainRhsPipeline().push_back(
};
auto set_boundary = [&]() {
};
}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
static constexpr int approx_order
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
MoFEMErrorCode assembleSystem()
[interior_boundary]
MoFEM::Interface & mField
static double sourceFunction(const double x, const double y, const double z)
Table 2 shows how the operators match the terms in Eq. \eqref{eq:init_cond_disc}
Table 2: Operators used to solve the initial timestep.
| Term | Operator |
| \( \int_\Omega N_\beta N_\alpha \, \textrm{d}\Omega \) | OpDomainMass("PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", ... {return 1;}) |
| \( \int_\Omega N_\beta s\, \textrm{d}\Omega \) | OpDomainSource("PHOTON_FLUENCE_RATE", sourcefunction) |
Finally, Eq. \eqref{eq:init_cond_disc} is solved using a Krylov Subspace Solver (KSP) [80] and the output data is written into a binary file which will be taken as the initial configuration for the numerical solver later on:
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Solver start";
CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD,
out_file_name, FILE_MODE_WRITE,
&viewer);
VecView(X, viewer);
PetscViewerDestroy(&viewer);
MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Solver done";
}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
#define MOFEM_LOG(channel, severity)
Log.
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
MoFEMErrorCode solveSystem()
[assembleSystem]
The function that runs the program is defined as:
}
}
MoFEMErrorCode readMesh()
[sourceFunction]
MoFEMErrorCode outputResults()
[solveSystem]
MoFEMErrorCode initialCondition()
MoFEMErrorCode checkResults()
[runProgram]
MoFEMErrorCode runProgram()
[runProgram]
MoFEMErrorCode createCommonData()
MoFEMErrorCode boundaryCondition()
[interior_boundary]
MoFEMErrorCode setupProblem()
readMesh() loads the mesh and the functions initialCondition() and boundaryCondition() set the initial and boundary condition. In this case, initialCondition() does not do anything as we start from an empty mesh. outputVolume() adds the output data into a mesh for visualization.
Photon diffusion
The operators used to implement Eqs. \eqref{Equation:tangentmat} and \eqref{Equation:resvec} are the following:
};
auto &bc_map = bc_mng->getBcMapByBlockName();
auto set_domain = [&]() {
pipeline_mng->getOpDomainLhsPipeline(), {H1});
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](
double,
double,
double) ->
double {
return D; }));
};
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainMass(
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", get_mass_coefficient));
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE",
grad_u_at_gauss_pts));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
u_at_gauss_pts));
pipeline_mng->getOpDomainRhsPipeline().push_back(
dot_u_at_gauss_pts));
"PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts,
[](
double,
double,
double) ->
double {
return D; }));
pipeline_mng->getOpDomainRhsPipeline().push_back(
"PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts,
[](
const double,
const double,
const double) {
return inv_v; }));
pipeline_mng->getOpDomainRhsPipeline().push_back(
"PHOTON_FLUENCE_RATE", u_at_gauss_pts,
[](
const double,
const double,
const double) {
return mu_a; }));
};
auto set_boundary = [&]() {
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpSetHOWeightsOnFace());
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpSetHOWeightsOnFace());
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
u_at_gauss_pts));
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpBoundaryMass(
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](
const double,
const double,
const double) {
return h; },
b.second->getBcEntsPtr()));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
"PHOTON_FLUENCE_RATE", u_at_gauss_pts,
[](
const double,
const double,
const double) {
return h; },
b.second->getBcEntsPtr()));
}
}
};
}
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
double mu_a
absorption coefficient (cm^-1)
boost::shared_ptr< FEMethod > boundaryRhsFEPtr
boost::shared_ptr< FEMethod > domainLhsFEPtr
boost::shared_ptr< FEMethod > boundaryLhsFEPtr
where these lines implement the exterior boundary:
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
new OpBoundaryMass(
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](
const double,
const double,
const double) {
return h; },
b.second->getBcEntsPtr()));
"PHOTON_FLUENCE_RATE", u_at_gauss_pts,
[](
const double,
const double,
const double) {
return h; },
b.second->getBcEntsPtr()));
}
}
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
and the degrees of freedom on the interior boundary, which represents the hidden object, are removed with the following code:
std::string entity_name = it->getName();
if (entity_name.compare(0, 3, "INT") == 0) {
boundary_faces, true);
}
}
if (boundary_faces.empty()) {
std::string entity_name = it->getName();
boundary_faces, true);
}
}
boundary_faces, 1, false, boundary_ents, moab::Interface::UNION);
boundary_faces, 0, false, boundary_verts, moab::Interface::UNION);
boundary_faces.merge(boundary_verts);
boundary_faces.merge(boundary_ents);
boundary_faces);
EntityHandle meshset;
simple->getProblemName(),
"PHOTON_FLUENCE_RATE", boundary_faces);
}
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual moab::Interface & get_moab()=0
Table 3 shows how the operators match the terms in Eqs. \eqref{Equation:tangentmat} and \eqref{Equation:resvec}:
Table 3: Operators used to solve the diffusion equation at each timestep.
| Term | Operator |
| \( \frac{1}{c} \int_\Omega N_\beta N_\gamma \Delta t^{-1}\, \textrm{d}\Omega + \mu_a \int_\Omega N_\beta N_\gamma \, \textrm{d}\Omega \) | OpDomainMass("PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", ... {return get_mass_coefficient;}) |
| \( D\int_\Omega \nabla N_\beta \cdot \nabla N_\gamma\, \textrm{d}\Omega \) | OpDomainGradGrad("PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", ... {return D;}) |
| \( \frac{1}{2A}\int_{\Gamma^\text{ext}} N_\beta N_\gamma \, \textrm{d}\Gamma \) | OpBoundaryMass("PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", ... {return h;}) |
| \( \frac{1}{c} \int_\Omega N_\beta \dot{\Phi}_0 \, \textrm{d}\Omega \) | OpDomainTimesScalarField("PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts, ... {return inv_v;}) |
| \( D\int_\Omega \nabla N_\beta \cdot \nabla\Phi_0, \textrm{d}\Omega \) | OpDomainGradTimesVec("PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts, ... {return D;}) |
| \( \mu_a \int_\Omega N_\beta \Phi_0 \, \textrm{d}\Omega \) | OpDomainTimesScalarField("PHOTON_FLUENCE_RATE", u_at_gauss_pts, ... {return mu_a;}) |
| \( \frac{1}{2A}\int_{\Gamma^\text{ext}} N_\beta \Phi_0 \, \textrm{d}\Gamma \) | OpDomainBoundaryTimesScalarField("PHOTON_FLUENCE_RATE", u_at_gauss_pts, ... {return h;}) |
The following code in the solveSystem function creates the post-processing element, which will be used to calculate the integral of the photon fluence rate over the area of the camera at each timestep.
auto create_post_process_element = [&]() {
auto post_froc_fe = boost::make_shared<PostProcEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
post_froc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
post_froc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE",
grad_ptr));
post_froc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_froc_fe->getPostProcMesh(), post_froc_fe->getMapGaussPts(),
{{"PHOTON_FLUENCE_RATE", u_ptr}},
{{"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
return post_froc_fe;
};
auto create_post_process_camera_element = [&]() {
if (mField.check_finite_element("CAMERA_FE")) {
auto post_proc_skin = boost::make_shared<PostProcFaceEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto op_loop_side = new OpLoopSide<FaceElementForcesAndSourcesCoreOnSide>(
op_loop_side->getOpPtrVector(), {H1});
op_loop_side->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
op_loop_side->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<SPACE_DIM>("PHOTON_FLUENCE_RATE",
grad_ptr));
post_proc_skin->getOpPtrVector().push_back(op_loop_side);
post_proc_skin->getOpPtrVector().push_back(
new OpPPMap(
post_proc_skin->getPostProcMesh(), post_proc_skin->getMapGaussPts(),
{{"PHOTON_FLUENCE_RATE", u_ptr}},
{{"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
return post_proc_skin;
} else {
return boost::shared_ptr<PostProcFaceEle>();
}
};
auto create_post_process_integ_camera_element = [&]() {
if (mField.check_finite_element("CAMERA_FE")) {
auto post_proc_integ_skin = boost::make_shared<BoundaryEle>(mField);
<< "Creating post process integ camera element";
post_proc_integ_skin->getOpPtrVector().push_back(
new OpSetHOWeightsOnFace());
post_proc_integ_skin->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
commonDataPtr->approxVals));
post_proc_integ_skin->getOpPtrVector().push_back(
new OpCameraInteg(commonDataPtr));
} else {
<< "Creating testing camera error element";
post_proc_integ_skin->getOpPtrVector().push_back(
new OpSetHOWeightsOnFace());
};
post_proc_integ_skin->getOpPtrVector(), {NOSPACE});
post_proc_integ_skin->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
commonDataPtr->uAtPtsPtr));
post_proc_integ_skin->getOpPtrVector().push_back(
}
return post_proc_integ_skin;
} else {
return boost::shared_ptr<BoundaryEle>();
}
};
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
The Time Solver is called at the end of solveSystem:
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(
dm, create_post_process_element(), create_post_process_camera_element(),
create_post_process_integ_camera_element(), commonDataPtr, mField));
boost::shared_ptr<ForcesAndSourcesCore> null;
CHKERR DMMoFEMTSSetMonitor(dm, solver,
simple->getDomainFEName(),
monitor_ptr, null, null);
};
auto X = createDMVector(dm);
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"reading vector in binary from file "
PetscViewer viewer;
&viewer);
VecLoad(X, viewer);
CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
}
auto solver = pipeline_mng->createTSIM();
CHKERR TSSetSolution(solver, X);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, X);
CHKERR TSSetFromOptions(solver);
CHKERR TSSetIJacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
}
auto createDMMatrix(DM dm)
Get smart matrix from DM.
char init_data_file_name[255]
[Push operators to pipeline]
In particular, createTS() creates a PETSc implicit time solver.
Finally, the photon fluence integral at a given timestep is calculated by summing the fluence (the solution to Eq. \eqref{Equation:weakform}) at each gauss point at said timestep over the camera surface:
MoFEMErrorCode
EntitiesFieldData::EntData &data) {
const int nb_integration_pts = getGaussPts().size2();
const double area = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
double values_integ = 0;
#ifdef ENABLE_PYTHON_BINDING
sens_vals_vec =
}
#endif
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double sens = 1;
#ifdef ENABLE_PYTHON_BINDING
sens = sens_vals_vec(gg);
}
#endif
const double alpha = t_w * area * sens;
values_integ += alpha * t_val;
++t_w;
++t_val;
}
std::array<double, 1> values;
values[0] = values_integ;
ADD_VALUES);
}
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
static auto getFTensor0FromVec(V &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
VectorDouble interp_function(const std::string sens_image, MatrixDouble &m_ref_coords, int nb_gauss_pts, double cam_len_x, double cam_len_y, const std::string block_name)
char interp_image_name[255]
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
[Integral_calc]
boost::shared_ptr< CommonData > commonDataPtr
The code loops through all entities in the camera and adds the contribution of each one into a commondataptr() which manages the global calculation:
constexpr std::array<int, 1> indices = {CommonData::VALUES_INTEG};
std::array<double, 1> values;
values[0] = values_integ;
CHKERR VecSetValues(commonDataPtr->petscVec, 1, indices.data(), values.data(),
ADD_VALUES);
Visualizing the results
We can use a python script to run the simulation for a given amount of time. The accuracy of the result will depend on the timestep, order, and number of elements. For a relatively coarse mesh, a timestep of \( 0.1 \text{s} \) and setting the order of the field to \(3\) yields a smooth approximation of the solution. Capturing the integral of the photon rate over the area of the camera for each timestep, we can plot the evolution of the photon fluence at each timestep. This can be seen on Fig. 3 alongside an example of the behaviour of the simulation around the hidden object on Fig. 4.
Figure 3: Example result of the evolution over time of the photon fluence rate at the camera surface.
Figure 4: Clip of the mesh used to produce the results in Fig. 3 showing the hidden object at the centre.
References
Plain Code
#include <stdlib.h>
#include <cmath>
#define BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND
static char help[] =
"...\n\n";
15;
;
#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace boost::math::quadrature;
public:
const double z) {
}
private:
};
};
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
if (
const size_t nb_dofs = data.
getIndices().size()) {
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
nf.clear();
const double volume = getMeasure();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
t_coords(0), t_coords(1), t_coords(2));
error += alpha * pow(diff, 2);
for (
size_t r = 0;
r != nb_dofs; ++
r) {
nf[
r] += alpha * t_row_base * diff;
++t_row_base;
}
++t_w;
++t_val;
++t_coords;
}
const int index = 0;
}
}
};
}
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
255, PETSC_NULLPTR);
PETSC_NULLPTR);
MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Refractive index: " <<
n;
MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Speed of light (cm/ns): " <<
c;
<<
"Phase velocity in medium (cm/ns): " <<
v;
<<
"Absorption coefficient (cm^-1): " <<
mu_a;
<<
"Scattering coefficient (cm^-1): " <<
mu_sp;
MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Diffusion coefficient D : " <<
D;
MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Approximation order: " <<
order;
}
commonDataPtr = boost::make_shared<PhotonDiffusion::CommonData>();
}
}
std::string entity_name = it->getName();
if (entity_name.compare(0, 3, "INT") == 0) {
boundary_faces, true);
}
}
if (boundary_faces.empty()) {
std::string entity_name = it->getName();
boundary_faces, true);
}
}
boundary_faces, 1, false, boundary_ents, moab::Interface::UNION);
boundary_faces, 0, false, boundary_verts, moab::Interface::UNION);
boundary_faces.merge(boundary_verts);
boundary_faces.merge(boundary_ents);
boundary_faces);
EntityHandle meshset;
simple->getProblemName(),
"PHOTON_FLUENCE_RATE", boundary_faces);
}
};
auto set_domain = [&]() {
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](const double, const double, const double) { return 1; }));
};
auto set_boundary = [&]() {
};
}
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Solver start";
CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
PetscViewer viewer;
PetscViewerBinaryOpen(PETSC_COMM_WORLD,
out_file_name, FILE_MODE_WRITE,
&viewer);
VecView(X, viewer);
PetscViewerDestroy(&viewer);
MOFEM_LOG(
"INITIAL", Sev::inform) <<
"Solver done";
}
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"PHOTON_FLUENCE_RATE", u_ptr}},
{},
{},
{})
);
pipeline_mng->getDomainPostProcFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElementsPostProc();
CHKERR post_proc_fe->writeFile(
"out_initial.h5m");
}
}
}
double nrm2;
const double *array;
PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n",
std::sqrt(array[0]), nrm2);
constexpr double eps = 1e-8;
"Not converged solution");
}
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "INITIAL"));
LogManager::setLog("INITIAL");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
CHKERR heat_problem.runProgram();
}
return 0;
}
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ MOFEM_DATA_INCONSISTENCY
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
const int kronrod_points
This has been tested and gives the same result for any number of points. Increasing the number of poi...
const double c
speed of light (cm/ns)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
double sourceFunctionEval(const double x, const double y, const double z, const double beam_radius, const double beam_centre_x, const double beam_centre_y, const double slab_thickness, const double mu_a, const double mu_sp, const double flux_magnitude, double initial_time, const double v, const double D)
Pulse is infinitely short.
double mu_sp
scattering coefficient (cm^-1)
double flux_magnitude
impulse magnitude
const int kronrod_points
This has been tested and gives the same result for any number of points. Increasing the number of poi...
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Specialization for double precision scalar field values calculation.
Post post-proc data at points from hash maps.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Get domain right-hand side finite element.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Get domain left-hand side finite element.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Set integration rule for domain right-hand side finite element.
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
Set integration rule for boundary left-hand side finite element.
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
Set integration rule for boundary right-hand side finite element.
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Set integration rule for domain left-hand side finite element.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
intrusive_ptr for managing petsc objects
Volume finite element base.
boost::shared_ptr< VectorDouble > uAtPtsPtr
SmartPetscObj< Vec > L2Vec
SmartPetscObj< Vec > resVec
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
PhotonDiffusion(MoFEM::Interface &m_field)
boost::shared_ptr< CommonData > commonDataPtr
#include <stdlib.h>
#include <cmath>
#include <boost/math/constants/constants.hpp>
#ifdef ENABLE_PYTHON_BINDING
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
static char help[] =
"...\n\n";
15;
#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace boost::math::quadrature;
public:
struct InterpPython;
private:
#ifdef ENABLE_PYTHON_BINDING
boost::shared_ptr<InterpPython> interpPythonPtr;
#endif
};
std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
};
struct OpGetScalarFieldGradientValuesOnSkin :
public BoundaryEleOp {
DataForcesAndSourcesCore::EntData &data) {
}
};
boost::shared_ptr<PostProcFaceEle> skin_post_proc,
boost::shared_ptr<BoundaryEle> skin_post_proc_integ,
boost::shared_ptr<CommonData> common_data_ptr,
}
}
auto vector_update = [&](auto vec) {
CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(vec, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(vec, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
};
const double *array;
<< "Fluence rate integral: " << array[0];
} else {
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Testing";
double nrm2;
const double *array2;
<< "Error " << array2[0] << " N of Entities " << nrm2;
constexpr double eps = 0.02;
if (array2[0] / nrm2 >
eps)
"Not converged solution");
}
boost::lexical_cast<std::string>(
ts_step) +
".h5m");
}
"out_camera_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
}
}
}
private:
boost::shared_ptr<PostProcEle>
postProc;
};
};
const double z, double time) {
};
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
double nf = 1;
const double volume = getMeasure();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
double analytical =
double diff = t_val - analytical;
error += alpha * pow(diff, 2) / std::abs(t_val);
++t_w;
++t_val;
++t_coords;
}
const int index = 0;
ADD_VALUES);
}
};
#ifdef ENABLE_PYTHON_BINDING
struct PhotonDiffusion::InterpPython {
InterpPython() = default;
virtual ~InterpPython() = default;
np::ndarray gauss_coords_x,
np::ndarray gauss_coords_y,
double cam_len_x,
template <typename T>
inline std::vector<T>
py_list_to_std_vector(const boost::python::object &iterable) {
return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
boost::python::stl_input_iterator<T>());
}
private:
bp::object mainNamespace;
bp::object InterpFun;
};
static boost::weak_ptr<PhotonDiffusion::InterpPython> interpPythonWeakPtr;
#endif
#ifdef ENABLE_PYTHON_BINDING
PhotonDiffusion::InterpPython::InterpInit(const std::string py_file) {
try {
bp::object main_module = bp::import("__main__");
mainNamespace = main_module.attr("__dict__");
bp::object ignored = bp::exec_file(py_file.c_str(), mainNamespace);
InterpFun = mainNamespace
["py_Interpolate"];
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
const std::string sens_image, np::ndarray gauss_coords_x,
np::ndarray &sens_vals) {
try {
sens_vals = bp::extract<np::ndarray>(InterpFun(
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
inline np::ndarray convert_to_numpy(
VectorDouble &data,
int nb_gauss_pts,
int id) {
auto dtype = np::dtype::get_builtin<double>();
auto size = bp::make_tuple(nb_gauss_pts);
auto stride = bp::make_tuple(3 * sizeof(double));
return (np::from_data(&data[id], dtype, size, stride, bp::object()));
}
#endif
const std::string block_name) {
#ifdef ENABLE_PYTHON_BINDING
if (auto interp_ptr = interpPythonWeakPtr.lock()) {
bp::list python_coords;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(convert_to_numpy(v_ref_coords, nb_gauss_pts, idx));
}
np::ndarray np_interp = np::empty(bp::make_tuple(nb_gauss_pts, 3),
np::dtype::get_builtin<double>());
auto interp_block_name = "(.*)INTERPOLATION(.*)";
std::regex reg_interp_name(interp_block_name);
if (std::regex_match(block_name, reg_interp_name)) {
sens_image, bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
cam_len_x,
"Failed py_Interp() python call");
} else {
}
if (np_interp.get_shape()[0] != nb_gauss_pts ||
np_interp.get_shape()[1] != 1) {
"Wrong shape of analytical expression returned from "
"python, expected: (" +
std::to_string(nb_gauss_pts) + ", 1), got: (" +
std::to_string(np_interp.get_shape()[0]) + ", " +
std::to_string(np_interp.get_shape()[1]) + ")");
}
double *interp_val_ptr = reinterpret_cast<double *>(np_interp.get_data());
v_interp.resize(nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
v_interp(gg) = *(interp_val_ptr + gg);
}
return v_interp;
} else {
"InterpPython pointer is expired");
}
#endif
}
}
PetscInt ghosts[1] = {0};
} else {
}
}
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
#ifdef ENABLE_PYTHON_BINDING
auto file_exists = [](std::string myfile) {
std::ifstream
file(myfile.c_str());
if (file) {
return true;
}
return false;
};
#ifdef ENABLE_PYTHON_BINDING
interpPythonPtr = boost::make_shared<PhotonDiffusion::InterpPython>();
interpPythonWeakPtr = interpPythonPtr;
} else {
}
}
#endif
#endif
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Refractive index: " <<
n;
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Speed of light (cm/ns): " <<
c;
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Phase velocity in medium (cm/ns): " <<
v;
<<
"Absorption coefficient (cm^-1): " <<
mu_a;
<<
"Scattering coefficient (cm^-1): " <<
mu_sp;
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Diffusion coefficient D : " <<
D;
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Coefficient A : " <<
A;
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Coefficient h : " <<
h;
auto set_camera_skin_fe = [&]() {
const std::string block_name = "CAM";
bool add_fe = false;
if (
bit->getName().compare(0, block_name.size(), block_name) == 0) {
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Found CAM block";
bit->getMeshset(), 2, camera_surface,
true);
add_fe = true;
}
}
MOFEM_LOG(
"PHOTON", Sev::noisy) <<
"CAM block entities:\n"
<< camera_surface;
if (add_fe) {
"PHOTON_FLUENCE_RATE");
"CAMERA_FE");
}
};
auto my_simple_set_up = [&]() {
}
};
}
}
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"EXT",
"PHOTON_FLUENCE_RATE", 0, 0, false);
std::string entity_name = it->getName();
if (entity_name.compare(0, 3, "INT") == 0) {
boundary_faces, true);
}
}
if (boundary_faces.empty()) {
std::string entity_name = it->getName();
boundary_faces, true);
}
}
boundary_faces, 1, false, boundary_ents, moab::Interface::UNION);
boundary_faces, 0, false, boundary_verts, moab::Interface::UNION);
boundary_faces.merge(boundary_verts);
boundary_faces.merge(boundary_ents);
boundary_faces);
simple->getProblemName(),
"PHOTON_FLUENCE_RATE", boundary_faces);
}
};
auto set_domain = [&]() {
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](
double,
double,
double) ->
double {
return D; }));
};
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE", get_mass_coefficient));
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
grad_u_at_gauss_pts));
u_at_gauss_pts));
dot_u_at_gauss_pts));
"PHOTON_FLUENCE_RATE", grad_u_at_gauss_pts,
[](
double,
double,
double) ->
double {
return D; }));
"PHOTON_FLUENCE_RATE", dot_u_at_gauss_pts,
[](
const double,
const double,
const double) {
return inv_v; }));
"PHOTON_FLUENCE_RATE", u_at_gauss_pts,
[](
const double,
const double,
const double) {
return mu_a; }));
};
auto set_boundary = [&]() {
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
u_at_gauss_pts));
for (auto b : bc_map) {
if (std::regex_match(b.first, std::regex("(.*)EXT(.*)"))) {
"PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
[](
const double,
const double,
const double) {
return h; },
b.second->getBcEntsPtr()));
"PHOTON_FLUENCE_RATE", u_at_gauss_pts,
[](
const double,
const double,
const double) {
return h; },
b.second->getBcEntsPtr()));
}
}
};
}
auto create_post_process_element = [&]() {
auto post_froc_fe = boost::make_shared<PostProcEle>(
mField);
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
post_froc_fe->getOpPtrVector().push_back(
post_froc_fe->getOpPtrVector().push_back(
grad_ptr));
post_froc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_froc_fe->getPostProcMesh(), post_froc_fe->getMapGaussPts(),
{{"PHOTON_FLUENCE_RATE", u_ptr}},
{{"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
return post_froc_fe;
};
auto create_post_process_camera_element = [&]() {
auto post_proc_skin = boost::make_shared<PostProcFaceEle>(mField);
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
op_loop_side->getOpPtrVector(), {H1});
op_loop_side->getOpPtrVector().push_back(
op_loop_side->getOpPtrVector().push_back(
grad_ptr));
post_proc_skin->getOpPtrVector().push_back(op_loop_side);
post_proc_skin->getOpPtrVector().push_back(
new OpPPMap(
post_proc_skin->getPostProcMesh(), post_proc_skin->getMapGaussPts(),
{{"PHOTON_FLUENCE_RATE", u_ptr}},
{{"GRAD_PHOTON_FLUENCE_RATE", grad_ptr}}, {}, {}));
return post_proc_skin;
} else {
return boost::shared_ptr<PostProcFaceEle>();
}
};
auto create_post_process_integ_camera_element = [&]() {
auto post_proc_integ_skin = boost::make_shared<BoundaryEle>(mField);
<< "Creating post process integ camera element";
post_proc_integ_skin->getOpPtrVector().push_back(
post_proc_integ_skin->getOpPtrVector().push_back(
commonDataPtr->approxVals));
post_proc_integ_skin->getOpPtrVector().push_back(
new OpCameraInteg(commonDataPtr));
} else {
<< "Creating testing camera error element";
post_proc_integ_skin->getOpPtrVector().push_back(
};
post_proc_integ_skin->getOpPtrVector(), {NOSPACE});
post_proc_integ_skin->getOpPtrVector().push_back(
commonDataPtr->uAtPtsPtr));
post_proc_integ_skin->getOpPtrVector().push_back(
}
return post_proc_integ_skin;
} else {
return boost::shared_ptr<BoundaryEle>();
}
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(
dm, create_post_process_element(), create_post_process_camera_element(),
create_post_process_integ_camera_element(), commonDataPtr, mField));
boost::shared_ptr<ForcesAndSourcesCore> null;
monitor_ptr, null, null);
};
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"reading vector in binary from file "
PetscViewer viewer;
&viewer);
VecLoad(X, viewer);
}
auto solver = pipeline_mng->createTSIM();
CHKERR TSSetSolution(solver, X);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetSolution(solver, X);
CHKERR TSSetFromOptions(solver);
CHKERR TSSetIJacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
}
}
}
const int nb_integration_pts = getGaussPts().size2();
const double area = getMeasure();
auto t_w = getFTensor0IntegrationWeight();
double values_integ = 0;
#ifdef ENABLE_PYTHON_BINDING
sens_vals_vec =
}
#endif
for (int gg = 0; gg != nb_integration_pts; ++gg) {
double sens = 1;
#ifdef ENABLE_PYTHON_BINDING
sens = sens_vals_vec(gg);
}
#endif
const double alpha = t_w * area * sens;
values_integ += alpha * t_val;
++t_w;
++t_val;
}
std::array<double, 1> values;
values[0] = values_integ;
ADD_VALUES);
}
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "PHOTON"));
LogManager::setLog("PHOTON");
#ifdef ENABLE_PYTHON_BINDING
PETSC_NULLPTR);
Py_Initialize();
np::initialize();
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Python initialised";
} else {
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Python NOT initialised";
}
#endif
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
CHKERR heat_problem.runProgram();
}
#ifdef ENABLE_PYTHON_BINDING
MOFEM_LOG(
"PHOTON", Sev::inform) <<
"Finalizing Python";
if (Py_FinalizeEx() < 0) {
exit(120);
}
}
#endif
return 0;
}
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore > PostProcFaceEle
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
BcMapByBlockName & getBcMapByBlockName()
Get the boundary condition map.
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
char interp_file_name[255]
Boundary condition manager for finite element problem setup.
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Element used to execute operators on side of the element.
Modify integration weights on face to take into account higher-order geometry.
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
Get boundary left-hand side finite element.
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Get boundary right-hand side finite element.
PetscInt ts_step
Current time step number.
Base volume element used to integrate on skeleton.
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > petscVec
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
boost::shared_ptr< PostProcFaceEle > skinPostProc
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
MoFEMErrorCode operator()()
Main operator function executed for each loop iteration.
MoFEM::Interface & mField
boost::shared_ptr< BoundaryEle > skinPostProcInteg
boost::shared_ptr< PostProcEle > postProc
OpCameraInteg(boost::shared_ptr< CommonData > common_data_ptr)
static double sourceFunction(const double x, const double y, const double z, double time)
boost::shared_ptr< VolSideFe > sideOpFe
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetScalarFieldGradientValuesOnSkin(boost::shared_ptr< VolSideFe > side_fe)
MoFEMErrorCode setIntegrationRules()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
