v0.13.1
Loading...
Searching...
No Matches
Classes | Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | List of all members
PhotonDiffusion Struct Reference
Collaboration diagram for PhotonDiffusion:
[legend]

Classes

struct  CommonData
 
struct  Monitor
 
struct  OpCameraInteg
 
struct  OpGetScalarFieldGradientValuesOnSkin
 

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 createCommonData ()
 
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< FEMethoddomainLhsFEPtr
 
boost::shared_ptr< FEMethodboundaryLhsFEPtr
 
boost::shared_ptr< FEMethodboundaryRhsFEPtr
 
boost::shared_ptr< CommonDatacommonDataPtr
 

Detailed Description

Examples
photon_diffusion.cpp.

Definition at line 74 of file initial_diffusion.cpp.

Constructor & Destructor Documentation

◆ PhotonDiffusion() [1/2]

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

Definition at line 138 of file initial_diffusion.cpp.

138: 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 284 of file initial_diffusion.cpp.

284 {
286
287 auto bc_mng = mField.getInterface<BcManager>();
288 auto &bc_map = bc_mng->getBcMapByBlockName();
289
290 auto add_domain_base_ops = [&](auto &pipeline) {
291 auto jac_ptr = boost::make_shared<MatrixDouble>();
292 auto det_ptr = boost::make_shared<VectorDouble>();
293 pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
294 pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, nullptr));
295 pipeline.push_back(new OpSetHOWeights(det_ptr));
296 };
297
298 auto add_domain_lhs_ops = [&](auto &pipeline) {
299 pipeline.push_back(new OpDomainMass(
300 "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
301 [](const double, const double, const double) { return 1; }));
302 };
303
304 auto add_domain_rhs_ops = [&](auto &pipeline) {
305 pipeline.push_back(
306 new OpDomainSource("PHOTON_FLUENCE_RATE", sourceFunction));
307 };
308
309 auto pipeline_mng = mField.getInterface<PipelineManager>();
310 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
311 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
312 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
313 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
314
316}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
OpCalculateHOJacForVolume OpCalculateHOJacVolume
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:23
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:165
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 257 of file initial_diffusion.cpp.

257 {
259
260 auto *simple = mField.getInterface<Simple>();
261
262 // Get boundary edges marked in block named "BOUNDARY_CONDITION"
263 Range boundary_ents;
265 std::string entity_name = it->getName();
266 if (entity_name.compare(0, 3, "INT") == 0) {
267 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
268 boundary_ents, true);
269 }
270 }
271 // Add vertices to boundary entities
272 Range boundary_verts;
273 CHKERR mField.get_moab().get_connectivity(boundary_ents, boundary_verts,
274 true);
275 boundary_ents.merge(boundary_verts);
276
277 // Remove DOFs as homogeneous boundary condition is used
278 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
279 simple->getProblemName(), "PHOTON_FLUENCE_RATE", boundary_ents);
280
282}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
@ BLOCKSET
Definition: definitions.h:148
#define CHKERR
Inline error check.
Definition: definitions.h:535
#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
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26

◆ boundaryCondition() [2/2]

MoFEMErrorCode PhotonDiffusion::boundaryCondition ( )
private

◆ checkResults()

MoFEMErrorCode PhotonDiffusion::checkResults ( )
private

◆ createCommonData()

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

Definition at line 223 of file photon_diffusion.cpp.

223 {
225 commonDataPtr = boost::make_shared<CommonData>();
226 PetscInt ghosts[1] = {0};
227 if (!mField.get_comm_rank())
228 commonDataPtr->petscVec =
229 createSmartGhostVector(mField.get_comm(), 1, 1, 0, ghosts);
230 else
231 commonDataPtr->petscVec =
232 createSmartGhostVector(mField.get_comm(), 0, 1, 1, ghosts);
233 commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
235}
auto createSmartGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
boost::shared_ptr< CommonData > commonDataPtr

◆ initialCondition() [1/2]

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

Definition at line 252 of file initial_diffusion.cpp.

252 {
255}

◆ initialCondition() [2/2]

MoFEMErrorCode PhotonDiffusion::initialCondition ( )
private

◆ outputResults() [1/2]

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

Definition at line 349 of file initial_diffusion.cpp.

349 {
351 auto *pipeline_mng = mField.getInterface<PipelineManager>();
352 pipeline_mng->getDomainLhsFE().reset();
353 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
354
355 auto u_ptr = boost::make_shared<VectorDouble>();
356 post_proc_fe->getOpPtrVector().push_back(
357 new OpCalculateScalarFieldValues("U", u_ptr));
358
360
361 post_proc_fe->getOpPtrVector().push_back(
362
363 new OpPPMap(
364
365 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
366
367 {{"PHOTON_FLUENCE_RATE", u_ptr}},
368
369 {},
370
371 {},
372
373 {})
374
375 );
376
377 pipeline_mng->getDomainRhsFE() = post_proc_fe;
378 CHKERR pipeline_mng->loopFiniteElements();
379 CHKERR post_proc_fe->writeFile("out_initial.h5m");
381}
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
boost::shared_ptr< FEMethod > & getDomainLhsFE()

◆ outputResults() [2/2]

MoFEMErrorCode PhotonDiffusion::outputResults ( )
private

◆ readMesh() [1/2]

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

Definition at line 140 of file initial_diffusion.cpp.

140 {
142
143 auto *simple = mField.getInterface<Simple>();
144 CHKERR mField.getInterface(simple);
145 CHKERR simple->getOptions();
146 CHKERR simple->loadFile();
147
149}

◆ readMesh() [2/2]

MoFEMErrorCode PhotonDiffusion::readMesh ( )
private

◆ runProgram() [1/2]

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

Definition at line 383 of file initial_diffusion.cpp.

383 {
385
393 if (output_volume)
395
397}
PetscBool output_volume
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 236 of file initial_diffusion.cpp.

236 {
238
239 auto integration_rule = [](int o_row, int o_col, int approx_order) {
240 return 2 * approx_order;
241 };
242
243 auto *pipeline_mng = mField.getInterface<PipelineManager>();
244 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
245 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
246 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
247 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
248
250}
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 151 of file initial_diffusion.cpp.

151 {
153
154 auto *simple = mField.getInterface<Simple>();
155 CHKERR simple->addDomainField("PHOTON_FLUENCE_RATE", H1,
157 CHKERR simple->addBoundaryField("PHOTON_FLUENCE_RATE", H1,
159
160 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-flux_magnitude",
161 &flux_magnitude, PETSC_NULL);
162 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-slab_thickness",
163 &slab_thickness, PETSC_NULL);
164 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-beam_radius", &beam_radius,
165 PETSC_NULL);
166 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-beam_centre_x", &beam_centre_x,
167 PETSC_NULL);
168 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-beam_centre_y", &beam_centre_y,
169 PETSC_NULL);
170 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-mu_a", &mu_a, PETSC_NULL);
171 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-mu_sp", &mu_sp, PETSC_NULL);
172 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-initial_time", &initial_time,
173 PETSC_NULL);
174
175 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-output_file", out_file_name,
176 255, PETSC_NULL);
177 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-output_volume", &output_volume,
178 PETSC_NULL);
179
180 D = 1. / (3. * (mu_a + mu_sp));
181
182 MOFEM_LOG("INITIAL", Sev::inform) << "Refractive index: " << n;
183 MOFEM_LOG("INITIAL", Sev::inform) << "Speed of light (cm/ns): " << c;
184 MOFEM_LOG("INITIAL", Sev::inform)
185 << "Phase velocity in medium (cm/ns): " << v;
186 MOFEM_LOG("INITIAL", Sev::inform)
187 << "Absorption coefficient (cm^-1): " << mu_a;
188 MOFEM_LOG("INITIAL", Sev::inform)
189 << "Scattering coefficient (cm^-1): " << mu_sp;
190 MOFEM_LOG("INITIAL", Sev::inform) << "Diffusion coefficient D : " << D;
191 MOFEM_LOG("INITIAL", Sev::inform) << "Impulse magnitude: " << flux_magnitude;
192 MOFEM_LOG("INITIAL", Sev::inform) << "Compute time (ns): " << initial_time;
193 MOFEM_LOG("INITIAL", Sev::inform) << "Slab thickness: " << slab_thickness;
194
195 int order = 2;
196 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
197
198 MOFEM_LOG("INITIAL", Sev::inform) << "Approximation order: " << order;
199
200 CHKERR simple->setFieldOrder("PHOTON_FLUENCE_RATE", order);
201
202 // if (numHoLevels > 0) {
203
204 // Range ho_ents;
205 // for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
206 // if (it->getName().compare(0, 3, "CAM") == 0) {
207 // CHKERR mField.get_moab().get_entities_by_dimension(it->getMeshset(), 2,
208 // ho_ents, true);
209 // }
210 // }
211
212 // EntityHandle meshset;
213 // CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
214 // CHKERR mField.get_moab().add_entities(meshset, ho_ents);
215 // std::string field_name;
216 // field_name = "out_test_" +
217 // boost::lexical_cast<std::string>(mField.get_comm_rank()) +
218 // ".vtk";
219 // CHKERR mField.get_moab().write_file(field_name.c_str(), "VTK", "", &meshset,
220 // 1);
221 // CHKERR mField.get_moab().delete_entities(&meshset, 1);
222
223 // CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
224
225 // CHKERR simple->setFieldOrder("PHOTON_FLUENCE_RATE", order + 1, &ho_ents);
226
227 // CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
228 // "PHOTON_FLUENCE_RATE");
229 // }
230
231 CHKERR simple->setUp();
232
234}
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
double flux_magnitude
impulse magnitude
char out_file_name[255]
double beam_centre_y
double beam_centre_x
double slab_thickness
double beam_radius
double initial_time
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)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
double mu_sp
scattering coefficient (cm^-1)
const double c
speed of light (cm/ns)
int order
double D
double mu_a
absorption coefficient (cm^-1)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium

◆ setupProblem() [2/2]

MoFEMErrorCode PhotonDiffusion::setupProblem ( )
private

◆ solveSystem() [1/2]

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

Definition at line 318 of file initial_diffusion.cpp.

318 {
320 auto *simple = mField.getInterface<Simple>();
321 auto *pipeline_mng = mField.getInterface<PipelineManager>();
322 auto solver = pipeline_mng->createKSP();
323
324 CHKERR KSPSetFromOptions(solver);
325 // CHKERR KSPSetUp(solver);
326
327 auto dm = simple->getDM();
328 auto X = smartCreateDMVector(dm);
329 auto F = smartVectorDuplicate(X);
330
331 MOFEM_LOG("INITIAL", Sev::inform) << "Solver start";
332 CHKERR KSPSolve(solver, F, X);
333 CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
334 CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
335 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
336
337 MOFEM_LOG("INITIAL", Sev::inform)
338 << "writing vector in binary to " << out_file_name << " ...";
339 PetscViewer viewer;
340 PetscViewerBinaryOpen(PETSC_COMM_WORLD, out_file_name, FILE_MODE_WRITE,
341 &viewer);
342 VecView(X, viewer);
343 PetscViewerDestroy(&viewer);
344
345 MOFEM_LOG("INITIAL", Sev::inform) << "Solver done";
347}
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:470
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
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 
)
inlinestatic

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 94 of file initial_diffusion.cpp.

94 {
95 const double A = 4. * D * v * initial_time;
96 const double T =
97 (v / pow(M_PI * A, 3. / 2.)) * exp(-mu_a * v * initial_time);
98
99 auto phi_pulse = [&](const double r_s, const double phi_s) {
100 const double xs = r_s * cos(phi_s);
101 const double ys = r_s * sin(phi_s);
102 const double xp = x - xs - beam_centre_x;
103 const double yp = y - ys - beam_centre_y;
104 const double zp1 = z + slab_thickness / 2. - 1. / mu_sp;
105 const double zp2 = z + slab_thickness / 2. + 1. / mu_sp;
106 const double P1 = xp * xp + yp * yp + zp1 * zp1;
107 const double P2 = xp * xp + yp * yp + zp2 * zp2;
108 return r_s * (exp(-P1 / A) - exp(-P2 / A));
109 };
110
111 auto f = [&](const double r_s) {
112 auto g = [&](const double phi_s) { return phi_pulse(r_s, phi_s); };
113 return gauss_kronrod<double, 15>::integrate(
114 g, 0, 2 * M_PI, 0, std::numeric_limits<float>::epsilon());
115 };
116
117 return T * flux_magnitude *
118 gauss_kronrod<double, 15>::integrate(
119 f, 0, beam_radius, 0, std::numeric_limits<float>::epsilon());
120 };
const double T
auto f
Definition: HenckyOps.hpp:5
double A
constexpr double g

Member Data Documentation

◆ boundaryLhsFEPtr

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

Definition at line 99 of file photon_diffusion.cpp.

◆ boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > PhotonDiffusion::boundaryMarker
private
Examples
photon_diffusion.cpp.

Definition at line 96 of file photon_diffusion.cpp.

◆ boundaryRhsFEPtr

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

Definition at line 100 of file photon_diffusion.cpp.

◆ commonDataPtr

boost::shared_ptr<CommonData> PhotonDiffusion::commonDataPtr
private
Examples
photon_diffusion.cpp.

Definition at line 112 of file photon_diffusion.cpp.

◆ domainLhsFEPtr

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

Definition at line 98 of file photon_diffusion.cpp.

◆ mField

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

Definition at line 135 of file initial_diffusion.cpp.


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