v0.13.1
initial_diffusion.cpp
Go to the documentation of this file.
1/**
2 * \file inital_diffusion.cpp
3 * \example inital_diffusion.cpp
4 *
5 **/
6
7#include <stdlib.h>
8#include <cmath>
9#include <BasicFiniteElements.hpp>
10
11#define BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND
12
13using namespace MoFEM;
14
15static char help[] = "...\n\n";
16
17template <int DIM> struct ElementsAndOps {};
18
19//! [Define dimension]
20constexpr int SPACE_DIM = 3; //< Space dimension of problem, mesh
21//! [Define dimension]
22
26
37
38const double n = 1.44; ///< refractive index of diffusive medium
39const double c = 30.; ///< speed of light (cm/ns)
40const double v = c / n; ///< phase velocity of light in medium (cm/ns)
41
42double mu_a; ///< absorption coefficient (cm^-1)
43double mu_sp; ///< scattering coefficient (cm^-1)
44double D;
45
47double beam_radius; //< spot radius
50double flux_magnitude = 1e3; ///< impulse magnitude
52
53char out_file_name[255] = "init_file.dat";;
55
56PetscBool output_volume = PETSC_FALSE;
57
58#include <boost/math/quadrature/gauss_kronrod.hpp>
59using namespace boost::math::quadrature;
60
61struct OpError : public DomainEleOp {
62 OpError(boost::shared_ptr<VectorDouble> u_at_pts_ptr, double &l2_error)
63 : DomainEleOp("PHOTON_FLUENCE_RATE", OPROW), uAtPtsPtr(u_at_pts_ptr),
64 l2Error(l2_error) {
65 // Only will be executed once on vertices
66 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
67 }
68
69private:
70 boost::shared_ptr<VectorDouble> uAtPtsPtr;
71 double &l2Error;
72};
73
75public:
77
78 // Declaration of the main function to run analysis
80
81 /**
82 * @brief Pulse is infinitely short.
83 *
84 * \note It is good approximation of pulse in femtosecond scale. To make it
85 * longer one can apply third integral over time.
86 *
87 * \note Note analysis in photon_diffusion is shifted in time by initial_time.
88 *
89 * @param x
90 * @param y
91 * @param z
92 * @return double
93 */
94 static double sourceFunction(const double x, const double y, const double z) {
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 };
121
122private:
123 // Declaration of other main functions called in runProgram()
133
134 // Main interfaces
136};
137
139
142
143 auto *simple = mField.getInterface<Simple>();
145 CHKERR simple->getOptions();
146 CHKERR simple->loadFile();
147
149}
150
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}
235
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}
251
255}
256
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}
283
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}
317
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}
348
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}
382
385
393 if (output_volume)
395
397}
398
399int main(int argc, char *argv[]) {
400
401 // Initialisation of MoFEM/PETSc and MOAB data structures
402 const char param_file[] = "param_file.petsc";
404
405 // Add logging channel for example
406 auto core_log = logging::core::get();
407 core_log->add_sink(
408 LogManager::createSink(LogManager::getStrmWorld(), "INITIAL"));
409 LogManager::setLog("INITIAL");
410 MOFEM_LOG_TAG("INITIAL", "initial_diffusion")
411
412 // Error handling
413 try {
414 // Register MoFEM discrete manager in PETSc
415 DMType dm_name = "DMMOFEM";
416 CHKERR DMRegister_MoFEM(dm_name);
417
418 // Create MOAB instance
419 moab::Core mb_instance; // mesh database
420 moab::Interface &moab = mb_instance; // mesh database interface
421
422 // Create MoFEM instance
423 MoFEM::Core core(moab); // finite element database
424 MoFEM::Interface &m_field = core; // finite element interface
425
426 // Run the main analysis
427 PhotonDiffusion heat_problem(m_field);
428 CHKERR heat_problem.runProgram();
429 }
431
432 // Finish work: cleaning memory, getting statistics, etc.
434
435 return 0;
436}
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
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.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
auto integration_rule
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
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
#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.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:27
int main(int argc, char *argv[])
double mu_sp
scattering coefficient (cm^-1)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
static char help[]
double flux_magnitude
impulse magnitude
char out_file_name[255]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
double beam_centre_y
int numHoLevels
constexpr int SPACE_DIM
[Define dimension]
double beam_centre_x
const double c
speed of light (cm/ns)
double slab_thickness
double beam_radius
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
double initial_time
double D
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
PetscBool output_volume
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
const double T
auto f
Definition: HenckyOps.hpp:5
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
OpCalculateHOJacForVolume OpCalculateHOJacVolume
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
double A
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static constexpr int approx_order
constexpr double g
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:165
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
@ OPROW
operator doWork function is executed on FE rows
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< VectorDouble > uAtPtsPtr
double & l2Error
OpError(boost::shared_ptr< VectorDouble > u_at_pts_ptr, double &l2_error)
MoFEMErrorCode assembleSystem()
MoFEMErrorCode solveSystem()
MoFEM::Interface & mField
MoFEMErrorCode readMesh()
MoFEMErrorCode outputResults()
MoFEMErrorCode initialCondition()
PhotonDiffusion(MoFEM::Interface &m_field)
MoFEMErrorCode checkResults()
MoFEMErrorCode runProgram()
MoFEMErrorCode boundaryCondition()
static double sourceFunction(const double x, const double y, const double z)
Pulse is infinitely short.
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode setupProblem()