v0.13.0
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 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #include <stdlib.h>
22 #include <cmath>
23 #include <BasicFiniteElements.hpp>
24 
25 using namespace MoFEM;
26 
27 static char help[] = "...\n\n";
28 
29 template <int DIM> struct ElementsAndOps {};
30 
31 //! [Define dimension]
32 constexpr int SPACE_DIM = 3; //< Space dimension of problem, mesh
33 //! [Define dimension]
34 
38 
49 
50 double n = 1.44; ///< refractive index of diffusive medium
51 double c = 30.; ///< speed of light (cm/ns)
52 double v = c / n; ///< phase velocity of light in medium (cm/ns)
53 double mu_a = 0.09; ///< absorption coefficient (cm^-1)
54 double mu_sp = 16.5; ///< scattering coefficient (cm^-1)
55 double flux_magnitude = 1e3; ///< impulse magnitude
56 
57 double slab_thickness = 5.0;
58 double spot_radius = 2.5; //< spot spot_radius
59 double initial_time = 0.1;
60 
61 double D = 1. / (3. * (mu_a + mu_sp));
62 
63 #include <boost/math/quadrature/gauss_kronrod.hpp>
64 using namespace boost::math::quadrature;
65 
66 struct OpError : public DomainEleOp {
67  OpError(boost::shared_ptr<VectorDouble> u_at_pts_ptr, double &l2_error)
68  : DomainEleOp("U", OPROW), uAtPtsPtr(u_at_pts_ptr), l2Error(l2_error) {
69  // Only will be executed once on vertices
70  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
71  }
72 
73  // MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
74  // MoFEMFunctionBegin;
75  // SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
76  // MogEMFunctionReturn(0);
77  // }
78 
79 private:
80  boost::shared_ptr<VectorDouble> uAtPtsPtr;
81  double &l2Error;
82 };
83 
85 public:
87 
88  // Declaration of the main function to run analysis
89  MoFEMErrorCode runProgram();
90 
91  /**
92  * @brief Pulse is infinitely short.
93  *
94  * \note It is good approximation of pulse in femtosecond scale. To make it
95  * longer one can apply third integral over time.
96  *
97  * \note Note analysis in photon_diffusion is shifted in time by initial_time.
98  *
99  * @param x
100  * @param y
101  * @param z
102  * @return double
103  */
104  static double sourceFunction(const double x, const double y, const double z) {
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  };
131 
132 private:
133  // Declaration of other main functions called in runProgram()
134  MoFEMErrorCode readMesh();
135  MoFEMErrorCode setupProblem();
136  MoFEMErrorCode setIntegrationRules();
137  MoFEMErrorCode initialCondition();
138  MoFEMErrorCode boundaryCondition();
139  MoFEMErrorCode assembleSystem();
140  MoFEMErrorCode solveSystem();
141  MoFEMErrorCode checkResults();
142  MoFEMErrorCode outputResults();
143 
144  // Main interfaces
146 
147  // Object to mark boundary entities for the assembling of domain elements
148  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
149 };
150 
152 
155 
156  auto *simple = mField.getInterface<Simple>();
158  CHKERR simple->getOptions();
159  CHKERR simple->loadFile();
160 
162 }
163 
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 }
202 
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 }
218 
222 }
223 
227 }
228 
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 }
260 
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 }
291 
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 }
308 
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 }
321 
324 
325  CHKERR readMesh();
332  // CHKERR checkResults();
334 
336 }
337 
338 int main(int argc, char *argv[]) {
339 
340  // Initialisation of MoFEM/PETSc and MOAB data structures
341  const char param_file[] = "param_file.petsc";
342  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
343 
344  // Add logging channel for example
345  auto core_log = logging::core::get();
346  core_log->add_sink(
347  LogManager::createSink(LogManager::getStrmWorld(), "INITIAL"));
348  LogManager::setLog("INITIAL");
349  MOFEM_LOG_TAG("INITIAL", "initial_diffusion")
350 
351  // Error handling
352  try {
353  // Register MoFEM discrete manager in PETSc
354  DMType dm_name = "DMMOFEM";
355  CHKERR DMRegister_MoFEM(dm_name);
356 
357  // Create MOAB instance
358  moab::Core mb_instance; // mesh database
359  moab::Interface &moab = mb_instance; // mesh database interface
360 
361  // Create MoFEM instance
362  MoFEM::Core core(moab); // finite element database
363  MoFEM::Interface &m_field = core; // finite element interface
364 
365  // Run the main analysis
366  PhotonDiffusion heat_problem(m_field);
367  CHKERR heat_problem.runProgram();
368  }
369  CATCH_ERRORS;
370 
371  // Finish work: cleaning memory, getting statistics, etc.
373 
374  return 0;
375 }
std::string param_file
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:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#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
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto integration_rule
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
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
Definition: helmholtz.cpp:37
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpDomainTimesScalarField
double c
speed of light (cm/ns)
double v
phase velocity of light in medium (cm/ns)
constexpr int SPACE_DIM
[Define dimension]
double slab_thickness
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
double spot_radius
double n
refractive index of diffusive medium
double initial_time
double D
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
double mu_a
absorption coefficient (cm^-1)
const double T
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
double A
static constexpr int approx_order
constexpr double g
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:131
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Calculate jacobian on Hex or other volume which is not simplex.
Get value at integration points for scalar field.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
boost::shared_ptr< VectorDouble > uAtPtsPtr
double & l2Error
OpError(boost::shared_ptr< VectorDouble > u_at_pts_ptr, double &l2_error)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
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()
Post processing.