v0.13.0
Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | List of all members
PhotonDiffusion Struct Reference
Collaboration diagram for PhotonDiffusion:
[legend]

Public Member Functions

 PhotonDiffusion (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProgram ()
 
 PhotonDiffusion (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProgram ()
 

Static Public Member Functions

static double sourceFunction (const double x, const double y, const double z)
 Pulse is infinitely short. More...
 

Private Member Functions

MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 
MoFEMErrorCode initialCondition ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 
MoFEMErrorCode initialCondition ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 

Private Attributes

MoFEM::InterfacemField
 
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
 
boost::shared_ptr< FEMethoddomianLhsFEPtr
 
boost::shared_ptr< FEMethodboundaryLhsFEPtr
 
boost::shared_ptr< FEMethodboundaryRhsFEPtr
 

Detailed Description

Examples
photon_diffusion.cpp.

Definition at line 84 of file initial_diffusion.cpp.

Constructor & Destructor Documentation

◆ PhotonDiffusion() [1/2]

PhotonDiffusion::PhotonDiffusion ( MoFEM::Interface m_field)
Examples
photon_diffusion.cpp.

Definition at line 151 of file initial_diffusion.cpp.

151 : mField(m_field) {}
MoFEM::Interface & mField

◆ PhotonDiffusion() [2/2]

PhotonDiffusion::PhotonDiffusion ( MoFEM::Interface m_field)

Member Function Documentation

◆ assembleSystem() [1/2]

MoFEMErrorCode PhotonDiffusion::assembleSystem ( )
private
Examples
photon_diffusion.cpp.

Definition at line 229 of file initial_diffusion.cpp.

229  {
231 
232  auto bc_mng = mField.getInterface<BcManager>();
233  auto &bc_map = bc_mng->getBcMapByBlockName();
234 
235  auto add_domain_base_ops = [&](auto &pipeline) {
236  auto jac_ptr = boost::make_shared<MatrixDouble>();
237  auto det_ptr = boost::make_shared<VectorDouble>();
238  pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
239  pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, nullptr));
240  pipeline.push_back(new OpSetHOWeights(det_ptr));
241  };
242 
243  auto add_domain_lhs_ops = [&](auto &pipeline) {
244  pipeline.push_back(new OpDomainMass(
245  "U", "U", [](const double, const double, const double) { return 1; }));
246  };
247 
248  auto add_domain_rhs_ops = [&](auto &pipeline) {
249  pipeline.push_back(new OpDomainSource("U", sourceFunction));
250  };
251 
252  auto pipeline_mng = mField.getInterface<PipelineManager>();
253  add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
254  add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
255  add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
256  add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
257 
259 }
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:131
Calculate jacobian on Hex or other volume which is not simplex.
Set inverse jacobian to base functions.
PipelineManager interface.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
static double sourceFunction(const double x, const double y, const double z)
Pulse is infinitely short.

◆ assembleSystem() [2/2]

MoFEMErrorCode PhotonDiffusion::assembleSystem ( )
private

◆ boundaryCondition() [1/2]

MoFEMErrorCode PhotonDiffusion::boundaryCondition ( )
private
Examples
photon_diffusion.cpp.

Definition at line 224 of file initial_diffusion.cpp.

224  {
227 }

◆ boundaryCondition() [2/2]

MoFEMErrorCode PhotonDiffusion::boundaryCondition ( )
private

◆ checkResults()

MoFEMErrorCode PhotonDiffusion::checkResults ( )
private

Definition at line 292 of file initial_diffusion.cpp.

292  {
295  pipeline_mng->getDomainLhsFE().reset();
296  pipeline_mng->getDomainRhsFE().reset();
297  pipeline_mng->getOpDomainRhsPipeline().clear();
298  auto u_vals_at_gauss_pts = boost::make_shared<VectorDouble>();
299  double l2_error = 0;
300  pipeline_mng->getOpDomainRhsPipeline().push_back(
301  new OpCalculateScalarFieldValues("U", u_vals_at_gauss_pts));
302  pipeline_mng->getOpDomainRhsPipeline().push_back(
303  new OpError(u_vals_at_gauss_pts, l2_error));
304  CHKERR pipeline_mng->loopFiniteElements();
305  MOFEM_LOG("INITIAL", Sev::inform) << "L2 error " << l2_error;
307 }
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
Get value at integration points for scalar field.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()

◆ initialCondition() [1/2]

MoFEMErrorCode PhotonDiffusion::initialCondition ( )
private
Examples
photon_diffusion.cpp.

Definition at line 219 of file initial_diffusion.cpp.

219  {
222 }

◆ initialCondition() [2/2]

MoFEMErrorCode PhotonDiffusion::initialCondition ( )
private

◆ outputResults() [1/2]

MoFEMErrorCode PhotonDiffusion::outputResults ( )
private
Examples
photon_diffusion.cpp.

Definition at line 309 of file initial_diffusion.cpp.

309  {
311  auto *pipeline_mng = mField.getInterface<PipelineManager>();
312  pipeline_mng->getDomainLhsFE().reset();
313  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
314  post_proc_fe->generateReferenceElementMesh();
315  post_proc_fe->addFieldValuesPostProc("U");
316  pipeline_mng->getDomainRhsFE() = post_proc_fe;
317  CHKERR pipeline_mng->loopFiniteElements();
318  CHKERR post_proc_fe->writeFile("out_initial.h5m");
320 }

◆ outputResults() [2/2]

MoFEMErrorCode PhotonDiffusion::outputResults ( )
private

◆ readMesh() [1/2]

MoFEMErrorCode PhotonDiffusion::readMesh ( )
private
Examples
photon_diffusion.cpp.

Definition at line 153 of file initial_diffusion.cpp.

153  {
155 
156  auto *simple = mField.getInterface<Simple>();
158  CHKERR simple->getOptions();
159  CHKERR simple->loadFile();
160 
162 }
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Simple interface for fast problem set-up.
Definition: Simple.hpp:33

◆ readMesh() [2/2]

MoFEMErrorCode PhotonDiffusion::readMesh ( )
private

◆ runProgram() [1/2]

MoFEMErrorCode PhotonDiffusion::runProgram ( )
Examples
photon_diffusion.cpp.

Definition at line 322 of file initial_diffusion.cpp.

322  {
324 
325  CHKERR readMesh();
332  // CHKERR checkResults();
334 
336 }
MoFEMErrorCode assembleSystem()
MoFEMErrorCode solveSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode outputResults()
MoFEMErrorCode initialCondition()
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode setupProblem()

◆ runProgram() [2/2]

MoFEMErrorCode PhotonDiffusion::runProgram ( )

◆ setIntegrationRules() [1/2]

MoFEMErrorCode PhotonDiffusion::setIntegrationRules ( )
private
Examples
photon_diffusion.cpp.

Definition at line 203 of file initial_diffusion.cpp.

203  {
205 
206  auto integration_rule = [](int o_row, int o_col, int approx_order) {
207  return 2 * approx_order;
208  };
209 
210  auto *pipeline_mng = mField.getInterface<PipelineManager>();
211  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
212  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
213  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
214  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
215 
217 }
auto integration_rule
static constexpr int approx_order

◆ setIntegrationRules() [2/2]

MoFEMErrorCode PhotonDiffusion::setIntegrationRules ( )
private

◆ setupProblem() [1/2]

MoFEMErrorCode PhotonDiffusion::setupProblem ( )
private
Examples
photon_diffusion.cpp.

Definition at line 164 of file initial_diffusion.cpp.

164  {
166 
167  auto *simple = mField.getInterface<Simple>();
168  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
169  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
170 
171  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-flux_magnitude",
172  &flux_magnitude, PETSC_NULL);
173  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-slab_thickness",
174  &slab_thickness, PETSC_NULL);
175  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-spot_radius", &spot_radius,
176  PETSC_NULL);
177  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-initial_time", &initial_time,
178  PETSC_NULL);
179 
180  MOFEM_LOG("INITIAL", Sev::inform) << "Refractive index: " << n;
181  MOFEM_LOG("INITIAL", Sev::inform) << "Speed of light (cm/ns): " << c;
182  MOFEM_LOG("INITIAL", Sev::inform)
183  << "Phase velocity in medium (cm/ns): " << v;
184  MOFEM_LOG("INITIAL", Sev::inform)
185  << "Absorption coefficient (cm^-1): " << mu_a;
186  MOFEM_LOG("INITIAL", Sev::inform)
187  << "Scattering coefficient (cm^-1): " << mu_sp;
188  MOFEM_LOG("INITIAL", Sev::inform) << "Impulse magnitude: " << flux_magnitude;
189  MOFEM_LOG("INITIAL", Sev::inform) << "Initial time (ns): " << initial_time;
190 
191  int order = 3;
192  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
193 
194  MOFEM_LOG("INITIAL", Sev::inform) << "Approximation order: " << order;
195 
196  CHKERR simple->setFieldOrder("U", order);
197 
198  CHKERR simple->setUp();
199 
201 }
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
double mu_sp
scattering coefficient (cm^-1)
double flux_magnitude
impulse magnitude
double c
speed of light (cm/ns)
double v
phase velocity of light in medium (cm/ns)
double slab_thickness
double spot_radius
double n
refractive index of diffusive medium
double initial_time
double mu_a
absorption coefficient (cm^-1)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)

◆ setupProblem() [2/2]

MoFEMErrorCode PhotonDiffusion::setupProblem ( )
private

◆ solveSystem() [1/2]

MoFEMErrorCode PhotonDiffusion::solveSystem ( )
private
Examples
photon_diffusion.cpp.

Definition at line 261 of file initial_diffusion.cpp.

261  {
263  auto *simple = mField.getInterface<Simple>();
264  auto *pipeline_mng = mField.getInterface<PipelineManager>();
265  auto solver = pipeline_mng->createKSP();
266 
267  CHKERR KSPSetFromOptions(solver);
268  CHKERR KSPSetUp(solver);
269 
270  auto dm = simple->getDM();
271  auto D = smartCreateDMVector(dm);
272  auto F = smartVectorDuplicate(D);
273 
274  MOFEM_LOG("INITIAL", Sev::inform) << "Solver start";
275  CHKERR KSPSolve(solver, F, D);
276  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
277  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
278  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
279 
280  MOFEM_LOG("INITIAL", Sev::inform)
281  << "writing vector in binary to vector.dat ...";
282  PetscViewer viewer;
283  PetscViewerBinaryOpen(PETSC_COMM_WORLD, "initial_vector.dat", FILE_MODE_WRITE,
284  &viewer);
285  VecView(D, viewer);
286  PetscViewerDestroy(&viewer);
287 
288  MOFEM_LOG("INITIAL", Sev::inform) << "Solver done";
290 }
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
double D
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.

◆ solveSystem() [2/2]

MoFEMErrorCode PhotonDiffusion::solveSystem ( )
private

◆ sourceFunction()

static double PhotonDiffusion::sourceFunction ( const double  x,
const double  y,
const double  z 
)
static

Pulse is infinitely short.

Note
It is good approximation of pulse in femtosecond scale. To make it longer one can apply third integral over time.
Note analysis in photon_diffusion is shifted in time by initial_time.
Parameters
x
y
z
Returns
double

Definition at line 104 of file initial_diffusion.cpp.

104  {
105  const double A = 4 * D * v * initial_time;
106  const double T =
107  (v / pow(M_PI * A, 3. / 2.)) * exp(-mu_a * v * initial_time);
108 
109  auto phi_pulse = [&](const double r_s, const double phi_s) {
110  const double xs = r_s * cos(phi_s);
111  const double ys = r_s * sin(phi_s);
112  const double xp = x - xs;
113  const double yp = y - ys;
114  const double zp1 = z + slab_thickness / 2. - 1. / mu_sp;
115  const double zp2 = z + slab_thickness / 2. + 1. / mu_sp;
116  const double P1 = xp * xp + yp * yp + zp1 * zp1;
117  const double P2 = xp * xp + yp * yp + zp2 * zp2;
118  return r_s * (exp(-P1 / A) - exp(-P2 / A));
119  };
120 
121  auto f = [&](const double r_s) {
122  auto g = [&](const double phi_s) { return phi_pulse(r_s, phi_s); };
123  return gauss_kronrod<double, 42>::integrate(
124  g, 0, 2 * M_PI, 0, std::numeric_limits<float>::epsilon());
125  };
126 
127  return T * flux_magnitude *
128  gauss_kronrod<double, 32>::integrate(
129  f, 0, spot_radius, 0, std::numeric_limits<float>::epsilon());
130  };
const double T
double A
constexpr double g

Member Data Documentation

◆ boundaryLhsFEPtr

boost::shared_ptr<FEMethod> PhotonDiffusion::boundaryLhsFEPtr
private
Examples
photon_diffusion.cpp.

Definition at line 142 of file photon_diffusion.cpp.

◆ boundaryMarker

boost::shared_ptr< std::vector< unsigned char > > PhotonDiffusion::boundaryMarker
private

Definition at line 148 of file initial_diffusion.cpp.

◆ boundaryRhsFEPtr

boost::shared_ptr<FEMethod> PhotonDiffusion::boundaryRhsFEPtr
private
Examples
photon_diffusion.cpp.

Definition at line 143 of file photon_diffusion.cpp.

◆ domianLhsFEPtr

boost::shared_ptr<FEMethod> PhotonDiffusion::domianLhsFEPtr
private
Examples
photon_diffusion.cpp.

Definition at line 141 of file photon_diffusion.cpp.

◆ mField

MoFEM::Interface & PhotonDiffusion::mField
private
Examples
photon_diffusion.cpp.

Definition at line 145 of file initial_diffusion.cpp.


The documentation for this struct was generated from the following files: