v0.9.1
mixed_reac_diff.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
3 #include <RDOperators.hpp>
4 
5 using namespace MoFEM;
6 using namespace ReactionDiffusion;
7 
8 static char help[] = "...\n\n";
9 
10 // #define M_PI 3.14159265358979323846 /* pi */
11 
12 struct RDProblem {
13 public:
14  RDProblem(MoFEM::Core &core, const int order)
15  : m_field(core)
16  , order(order)
17  , cOmm(m_field.get_comm())
18  , rAnk(m_field.get_comm_rank()) {
19  vol_ele_slow_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
20  natural_bdry_ele_slow_rhs =
21  boost::shared_ptr<BoundaryEle>(new BoundaryEle(m_field));
22  vol_ele_stiff_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
23  vol_ele_stiff_lhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
24  post_proc = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
25  new PostProcFaceOnRefinedMesh(m_field));
26 
27  data1 = boost::shared_ptr<PreviousData>(new PreviousData());
28  data2 = boost::shared_ptr<PreviousData>(new PreviousData());
29  data3 = boost::shared_ptr<PreviousData>(new PreviousData());
30 
31  flux_values_ptr1 =
32  boost::shared_ptr<MatrixDouble>(data1, &data1->flux_values);
33 
34  flux_divs_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->flux_divs);
35 
36  mass_values_ptr1 =
37  boost::shared_ptr<VectorDouble>(data1, &data1->mass_values);
38 
39  mass_dots_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->mass_dots);
40 
41  flux_values_ptr2 =
42  boost::shared_ptr<MatrixDouble>(data2, &data2->flux_values);
43  flux_divs_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->flux_divs);
44 
45  mass_values_ptr2 =
46  boost::shared_ptr<VectorDouble>(data2, &data2->mass_values);
47  mass_dots_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->mass_dots);
48 
49  flux_values_ptr3 =
50  boost::shared_ptr<MatrixDouble>(data3, &data3->flux_values);
51  flux_divs_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->flux_divs);
52 
53  mass_values_ptr3 =
54  boost::shared_ptr<VectorDouble>(data3, &data3->mass_values);
55 
56  mass_dots_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->mass_dots);
57  }
58 
59  // RDProblem(const int order) : order(order){}
60  MoFEMErrorCode run_analysis(int nb_sp);
61 
62  double global_error;
63 
64 private:
65  MoFEMErrorCode setup_system();
66  MoFEMErrorCode add_fe(std::string mass_field, std::string flux_field);
67  MoFEMErrorCode set_blockData(std::map<int, BlockData> &block_data_map);
68  MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL, std::string internal);
69  MoFEMErrorCode extract_initial_ents(int block_id, Range &surface);
70  MoFEMErrorCode update_slow_rhs(std::string mass_fiedl,
71  boost::shared_ptr<VectorDouble> &mass_ptr);
72  MoFEMErrorCode push_slow_rhs(std::string mass_field, std::string flux_field,
73  boost::shared_ptr<PreviousData> &data);
74  MoFEMErrorCode update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
75  boost::shared_ptr<PreviousData> &data);
77  update_stiff_rhs(std::string mass_field, std::string flux_field,
78  boost::shared_ptr<VectorDouble> &mass_ptr,
79  boost::shared_ptr<MatrixDouble> &flux_ptr,
80  boost::shared_ptr<VectorDouble> &mass_dot_ptr,
81  boost::shared_ptr<VectorDouble> &flux_div_ptr);
82  MoFEMErrorCode push_stiff_rhs(std::string mass_field, std::string flux_field,
83  boost::shared_ptr<PreviousData> &data,
84  std::map<int, BlockData> &block_map);
85  MoFEMErrorCode update_stiff_lhs(std::string mass_fiedl,
86  std::string flux_field,
87  boost::shared_ptr<VectorDouble> &mass_ptr,
88  boost::shared_ptr<MatrixDouble> &flux_ptr);
89  MoFEMErrorCode push_stiff_lhs(std::string mass_field, std::string flux_field,
90  boost::shared_ptr<PreviousData> &data,
91  std::map<int, BlockData> &block_map);
92 
93  MoFEMErrorCode set_integration_rule();
94  MoFEMErrorCode apply_IC(std::string mass_field, Range &surface,
95  boost::shared_ptr<FaceEle> &initial_ele);
96  MoFEMErrorCode apply_BC(std::string flux_field);
97  MoFEMErrorCode loop_fe();
98  MoFEMErrorCode post_proc_fields(std::string mass_field,
99  std::string flux_field);
100  MoFEMErrorCode output_result();
101  MoFEMErrorCode solve();
102 
103  MoFEM::Interface &m_field;
104  Simple *simple_interface;
107 
108  Range essential_bdry_ents;
109  Range natural_bdry_ents;
110 
112 
113  Range inner_surface1; // nb_species times
114  Range inner_surface2;
115  Range inner_surface3;
116 
117  MPI_Comm cOmm;
118  const int rAnk;
119 
120  int order;
121  int nb_species;
122 
123  std::map<int, BlockData> material_blocks;
124 
125  boost::shared_ptr<FaceEle> vol_ele_slow_rhs;
126  boost::shared_ptr<FaceEle> vol_ele_stiff_rhs;
127  boost::shared_ptr<FaceEle> vol_ele_stiff_lhs;
128  boost::shared_ptr<BoundaryEle> natural_bdry_ele_slow_rhs;
129 
130  boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc;
131  boost::shared_ptr<Monitor> monitor_ptr;
132 
133  boost::shared_ptr<PreviousData> data1; // nb_species times
134  boost::shared_ptr<PreviousData> data2;
135  boost::shared_ptr<PreviousData> data3;
136 
137  boost::shared_ptr<MatrixDouble> flux_values_ptr1; // nb_species times
138  boost::shared_ptr<MatrixDouble> flux_values_ptr2;
139  boost::shared_ptr<MatrixDouble> flux_values_ptr3;
140 
141  boost::shared_ptr<VectorDouble> flux_divs_ptr1; // nb_species times
142  boost::shared_ptr<VectorDouble> flux_divs_ptr2;
143  boost::shared_ptr<VectorDouble> flux_divs_ptr3;
144 
145  boost::shared_ptr<VectorDouble> mass_values_ptr1; // nb_species times
146  boost::shared_ptr<VectorDouble> mass_values_ptr2;
147  boost::shared_ptr<VectorDouble> mass_values_ptr3;
148 
149  boost::shared_ptr<VectorDouble> mass_dots_ptr1; // nb_species times
150  boost::shared_ptr<VectorDouble> mass_dots_ptr2;
151  boost::shared_ptr<VectorDouble> mass_dots_ptr3;
152 
153  boost::shared_ptr<ForcesAndSourcesCore> null;
154 };
155 
156 const double ramp_t = 1.0;
157 const double sml = 0.0;
158 const double T = M_PI / 2.0;
159 
160 struct KinkFunction {
161  double operator() (const double x) const {
162  return 1 - abs(x);
163  }
164 };
165 
167  double operator()(const double x) const {
168  if(x > 0)
169  {return -1;}
170  else
171  {return 1;}
172 }
173 };
174 
175 struct ExactFunction {
176  double operator()(const double x, const double y, const double t) const {
177  double g = cos(T * x) * cos(T * y);
178  if (t <= ramp_t) {
179  return g * t;
180  } else {
181  return g * ramp_t;
182  }
183  }
184 };
185 
186 struct ExactFunctionGrad {
187  FTensor::Tensor1<double, 3> operator()(const double x, const double y,
188  const double t) const {
190  double mx = - T * sin(T * x) * cos(T * y);
191  double my = - T * cos(T * x) * sin(T * y);
192 
193 
194  if (t <= ramp_t) {
195  grad(0) = mx * t;
196  grad(1) = my * t;
197  } else {
198  grad(0) = mx * ramp_t;
199  grad(1) = my * ramp_t;
200  }
201  grad(2) = 0.0;
202  return grad;
203  }
204 };
205 
207  double operator()(const double x, const double y, const double t) const {
208  double glap = -2.0 * pow(T, 2) * cos(T * x) * cos(T * y);
209  if (t <= ramp_t) {
210  return glap * t;
211  } else {
212  return glap * ramp_t;
213  }
214  }
215 };
216 
218  double operator()(const double x, const double y, const double t) const {
219  double gdot = cos(T * x) * cos(T * y);
220  if (t <= ramp_t) {
221  return gdot;
222  } else {
223  return 0;
224  }
225  }
226 };
227 
230  CHKERR m_field.getInterface(simple_interface);
231  CHKERR simple_interface->getOptions();
232  CHKERR simple_interface->loadFile();
234 }
235 
236 MoFEMErrorCode RDProblem::add_fe(std::string mass_field,
237  std::string flux_field) {
239  CHKERR simple_interface->addDomainField(mass_field, L2,
241 
242  CHKERR simple_interface->addDomainField(flux_field, HCURL,
244 
245  CHKERR simple_interface->addBoundaryField(flux_field, HCURL,
247  CHKERR simple_interface->addDataField("ERROR", L2, AINSWORTH_LEGENDRE_BASE,
248  1);
249 
250  CHKERR simple_interface->setFieldOrder(mass_field, order - 1);
251  CHKERR simple_interface->setFieldOrder(flux_field, order);
252  CHKERR simple_interface->setFieldOrder("ERROR", 0); // approximation order for error
253 
255 }
256 MoFEMErrorCode RDProblem::set_blockData(std::map<int, BlockData> &block_map) {
259  string name = it->getName();
260  const int id = it->getMeshsetId();
261  if (name.compare(0, 14, "REGION1") == 0) {
262  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
263  id, BLOCKSET, 2, block_map[id].block_ents, true);
264  block_map[id].B0 = 1e-4;
265  block_map[id].block_id = id;
266  } else if (name.compare(0, 14, "REGION2") == 0) {
267  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
268  id, BLOCKSET, 2, block_map[id].block_ents, true);
269  block_map[id].B0 = 1e-4;
270  block_map[id].block_id = id;
271  } else if (name.compare(0, 14, "REGION3") == 0) {
272  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
273  id, BLOCKSET, 2, block_map[id].block_ents, true);
274  block_map[id].B0 = 1e-4;
275  block_map[id].block_id = id;
276  } else if (name.compare(0, 14, "REGION4") == 0) {
277  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
278  id, BLOCKSET, 2, block_map[id].block_ents, true);
279  block_map[id].B0 = 1e-4;
280  block_map[id].block_id = id;
281  } else if (name.compare(0, 14, "REGION5") == 0) {
282  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
283  id, BLOCKSET, 2, block_map[id].block_ents, true);
284  block_map[id].B0 = 1e-4;
285  block_map[id].block_id = id;
286  }
287  }
289 }
290 
292  std::string natural,
293  std::string internal) {
296  string name = it->getName();
297  if (name.compare(0, 14, natural) == 0) {
298 
299  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
300  natural_bdry_ents, true);
301  } else if (name.compare(0, 14, essential) == 0) {
302  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
303  essential_bdry_ents, true);
304  } else if (name.compare(0, 14, internal) == 0) {
305  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
306  internal_edge_ents, true);
307  }
308  }
310 }
311 
312 MoFEMErrorCode RDProblem::extract_initial_ents(int block_id, Range &surface) {
314  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(block_id,
315  BLOCKSET)) {
316  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
317  block_id, BLOCKSET, 2, surface, true);
318  }
320 }
322 RDProblem::update_slow_rhs(std::string mass_field,
323  boost::shared_ptr<VectorDouble> &mass_ptr) {
325  vol_ele_slow_rhs->getOpPtrVector().push_back(
326  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
328 }
329 
331  std::string flux_field,
332  boost::shared_ptr<PreviousData> &data) {
334 
335  vol_ele_slow_rhs->getOpPtrVector().push_back(
336  new OpAssembleSlowRhsV(mass_field, data, ExactFunction(), ExactFunctionDot(), ExactFunctionLap()));
337 
338  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
339  new OpSkeletonSource(mass_field,
340  KinkFunction(),
341  ExactFunction(),
342  internal_edge_ents));
343 
344  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
345  new OpAssembleNaturalBCRhsTau(flux_field, natural_bdry_ents));
346 
347  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
348  new OpEssentialBC(flux_field, essential_bdry_ents));
350 }
351 
352 MoFEMErrorCode RDProblem::update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
353  boost::shared_ptr<PreviousData> &data) {
355  vol_ele->getOpPtrVector().push_back(new OpCalculateJacForFace(data->jac));
356  vol_ele->getOpPtrVector().push_back(
357  new OpCalculateInvJacForFace(data->inv_jac));
358  vol_ele->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
359 
360  vol_ele->getOpPtrVector().push_back(
361  new OpSetContravariantPiolaTransformFace(data->jac));
362 
363  vol_ele->getOpPtrVector().push_back(new OpSetInvJacHcurlFace(data->inv_jac));
365 }
367 RDProblem::update_stiff_rhs(std::string mass_field, std::string flux_field,
368  boost::shared_ptr<VectorDouble> &mass_ptr,
369  boost::shared_ptr<MatrixDouble> &flux_ptr,
370  boost::shared_ptr<VectorDouble> &mass_dot_ptr,
371  boost::shared_ptr<VectorDouble> &flux_div_ptr) {
372 
374 
375  vol_ele_stiff_rhs->getOpPtrVector().push_back(
376  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
377 
378  vol_ele_stiff_rhs->getOpPtrVector().push_back(
379  new OpCalculateHdivVectorField<3>(flux_field, flux_ptr));
380 
381  vol_ele_stiff_rhs->getOpPtrVector().push_back(
382  new OpCalculateScalarValuesDot(mass_field, mass_dot_ptr));
383 
384  vol_ele_stiff_rhs->getOpPtrVector().push_back(
385  new OpCalculateHdivVectorDivergence<3, 2>(flux_field, flux_div_ptr));
387 }
388 
390  std::string flux_field,
391  boost::shared_ptr<PreviousData> &data,
392  std::map<int, BlockData> &block_map) {
394  vol_ele_stiff_rhs->getOpPtrVector().push_back(
395  new OpAssembleStiffRhsTau<3>(flux_field, data, block_map));
396 
397  vol_ele_stiff_rhs->getOpPtrVector().push_back(
398  new OpAssembleStiffRhsV<3>(mass_field, data, block_map, ExactFunction(),
401 }
402 
404 RDProblem::update_stiff_lhs(std::string mass_field, std::string flux_field,
405  boost::shared_ptr<VectorDouble> &mass_ptr,
406  boost::shared_ptr<MatrixDouble> &flux_ptr) {
408  vol_ele_stiff_lhs->getOpPtrVector().push_back(
409  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
410 
411  vol_ele_stiff_lhs->getOpPtrVector().push_back(
412  new OpCalculateHdivVectorField<3>(flux_field, flux_ptr));
414 }
415 
417  std::string flux_field,
418  boost::shared_ptr<PreviousData> &data,
419  std::map<int, BlockData> &block_map) {
421  vol_ele_stiff_lhs->getOpPtrVector().push_back(
422  new OpAssembleLhsTauTau<3>(flux_field, data, block_map));
423 
424  vol_ele_stiff_lhs->getOpPtrVector().push_back(
425  new OpAssembleLhsVV(mass_field));
426 
427  vol_ele_stiff_lhs->getOpPtrVector().push_back(
428  new OpAssembleLhsTauV<3>(flux_field, mass_field, data, block_map));
429 
430  vol_ele_stiff_lhs->getOpPtrVector().push_back(
431  new OpAssembleLhsVTau(mass_field, flux_field));
433 }
434 
437  auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
438  vol_ele_slow_rhs->getRuleHook = vol_rule;
439  natural_bdry_ele_slow_rhs->getRuleHook = vol_rule;
440 
441  vol_ele_stiff_rhs->getRuleHook = vol_rule;
442 
443  vol_ele_stiff_lhs->getRuleHook = vol_rule;
445 }
446 
447 MoFEMErrorCode RDProblem::apply_IC(std::string mass_field, Range &surface,
448  boost::shared_ptr<FaceEle> &initial_ele) {
450  initial_ele->getOpPtrVector().push_back(
451  new OpInitialMass(mass_field, surface));
453 }
454 
455 MoFEMErrorCode RDProblem::apply_BC(std::string flux_field) {
457  CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
458  "SimpleProblem", flux_field, essential_bdry_ents);
459 
461 }
464 
465  // Vec R, U;
466  // CHKERR DMCreateGlobalVector(dm, &R);
467 
468  // CHKERR VecDuplicate(R, &U);
469 
470  // DM dm_sub_ff, dm_sub_fm;
471  // ublas::matrix<Mat> nested_matrices(2, 2);
472  // ublas::vector<IS> nested_is(2);
473 
474  // CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_ff);
475  // CHKERR DMSetType(dm_sub_ff, "DMMOFEM");
476  // CHKERR DMMoFEMCreateSubDM(dm_sub_hh, dm, "SUB_FF");
477 
478  // CHKERR DMMoFEMSetSquareProblem(dm_sub_ff, PETSC_TRUE);
479  // CHKERR DMMoFEMAddSubFieldRow(dm_sub_ff, "FLUX1");
480  // CHKERR DMMoFEMAddSubFieldCol(dm_sub_ff, "FLUX1");
481  // CHKERR DMMoFEMAddElement(dm_sub_ff,
482  // simple_interface->getDomainFEName().c_str());
483  // CHKERR DMSetUp(dm_sub_ff);
484 
485  // CHKERR DMMoFEMGetSubRowIS(dm_sub_ff, &nested_is[0]);
486  // CHKERR DMCreateMatrix(dm_sub_ff, &nested_matrices(0, 0));
487 
488  // vol_ele_stiff_lhs->getFEMethod()->ts_B = nested_matrices(0, 0);
489  // CHKERR TSSetType(ts, TSARKIMEX);
490  // CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
491 
492  // CHKERR DMMoFEMTSSetIJacobian(dm_sub_ff, simple_interface->getDomainFEName(),
493  // vol_ele_stiff_lhs, null, null);
494 
495  // CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
496  // CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
497 
498  CHKERR TSSetType(ts, TSARKIMEX);
499  CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
500 
501  CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
502  vol_ele_stiff_lhs, null, null);
503 
504 
505  CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
506  vol_ele_stiff_rhs, null, null);
507 
508  CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
509  vol_ele_slow_rhs, null, null);
510  CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getBoundaryFEName(),
511  natural_bdry_ele_slow_rhs, null, null);
512 
513 
515 }
516 
518  std::string flux_field) {
520  post_proc->addFieldValuesPostProc(mass_field);
521  post_proc->addFieldValuesPostProc(flux_field);
523 }
524 
527  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
528  monitor_ptr, null, null);
530 }
533  // Create solution vector
536  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
537  // Solve problem
538  double ftime = 1;
539  CHKERR TSSetDM(ts, dm);
540  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
541  CHKERR TSSetSolution(ts, X);
542  CHKERR TSSetFromOptions(ts);
543 
544  if (1) {
545  SNES snes;
546  CHKERR TSGetSNES(ts, &snes);
547  KSP ksp;
548  CHKERR SNESGetKSP(snes, &ksp);
549  PC pc;
550  CHKERR KSPGetPC(ksp, &pc);
551  PetscBool is_pcfs = PETSC_FALSE;
552  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
553  // Set up FIELDSPLIT
554  // Only is user set -pc_type fieldsplit
555  if (is_pcfs == PETSC_TRUE) {
556  IS is_mass, is_flux;
557  const MoFEM::Problem *problem_ptr;
558  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
559  CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
560  problem_ptr->getName(), ROW, "MASS1", 0, 1, &is_mass);
561  CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
562  problem_ptr->getName(), ROW, "FLUX1", 0, 1, &is_flux);
563  // CHKERR ISView(is_flux, PETSC_VIEWER_STDOUT_SELF);
564  // CHKERR ISView(is_mass, PETSC_VIEWER_STDOUT_SELF);
565 
566  CHKERR PCFieldSplitSetIS(pc, NULL, is_flux);
567  CHKERR PCFieldSplitSetIS(pc, NULL, is_mass);
568 
569 
570  CHKERR ISDestroy(&is_flux);
571  CHKERR ISDestroy(&is_mass);
572  }
573  }
574 
575  CHKERR TSSolve(ts, X);
577 }
578 
581  global_error = 0;
582  // set nb_species
583  CHKERR setup_system(); // only once
584  nb_species = nb_sp;
585  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
586  CHKERR add_fe("MASS1", "FLUX1"); // nb_species times
587  if (nb_species == 2 || nb_species == 3) {
588  CHKERR add_fe("MASS2", "FLUX2");
589  if (nb_species == 3) {
590  CHKERR add_fe("MASS3", "FLUX3");
591  }
592  }
593  }
594 
595  CHKERR simple_interface->setUp();
596 
597  CHKERR set_blockData(material_blocks);
598 
599  CHKERR extract_bd_ents("ESSENTIAL", "NATURAL", "INTERNAL"); // nb_species times
600 
601  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
602  CHKERR extract_initial_ents(2, inner_surface1);
603  CHKERR update_slow_rhs("MASS1", mass_values_ptr1);
604  if (nb_species == 1) {
605  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
606  "MASS1", data1, data1, data1, material_blocks));
607  } else if (nb_species == 2 || nb_species == 3) {
608  CHKERR extract_initial_ents(3, inner_surface2);
609  CHKERR update_slow_rhs("MASS2", mass_values_ptr2);
610  if (nb_species == 2) {
611  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
612  "MASS1", data1, data2, data2, material_blocks));
613  } else if (nb_species == 3) {
614  CHKERR extract_initial_ents(4, inner_surface3);
615  CHKERR update_slow_rhs("MASS3", mass_values_ptr3);
616  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
617  "MASS1", data1, data2, data3, material_blocks));
618  }
619  }
620  }
621  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
623 
624  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
625  CHKERR push_slow_rhs("MASS1", "FLUX1", data1); // nb_species times
626  if (nb_species == 2 || nb_species == 3) {
627  CHKERR push_slow_rhs("MASS2", "FLUX2", data2);
628  if (nb_species == 3) {
629  CHKERR push_slow_rhs("MASS3", "FLUX3", data3);
630  }
631  }
632  }
633 
634  CHKERR update_vol_fe(vol_ele_stiff_rhs, data1);
635 
636  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
637  CHKERR update_stiff_rhs("MASS1", "FLUX1", mass_values_ptr1,
638  flux_values_ptr1, mass_dots_ptr1, flux_divs_ptr1);
639  CHKERR push_stiff_rhs("MASS1", "FLUX1", data1,
640  material_blocks); // nb_species times
641  if (nb_species == 2 || nb_species == 3) {
642  CHKERR update_stiff_rhs("MASS2", "FLUX2", mass_values_ptr2,
643  flux_values_ptr2, mass_dots_ptr2, flux_divs_ptr2);
644  CHKERR push_stiff_rhs("MASS2", "FLUX2", data2, material_blocks);
645  if (nb_species == 3) {
646  CHKERR update_stiff_rhs("MASS3", "FLUX3", mass_values_ptr3,
647  flux_values_ptr3, mass_dots_ptr3,
648  flux_divs_ptr3);
649  CHKERR push_stiff_rhs("MASS3", "FLUX3", data3, material_blocks);
650  }
651  }
652  }
653 
654  CHKERR update_vol_fe(vol_ele_stiff_lhs, data1);
655 
656  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
657  CHKERR update_stiff_lhs("MASS1", "FLUX1", mass_values_ptr1,
658  flux_values_ptr1);
659  CHKERR push_stiff_lhs("MASS1", "FLUX1", data1,
660  material_blocks); // nb_species times
661  if (nb_species == 2 || nb_species == 3) {
662  CHKERR update_stiff_lhs("MASS2", "FLUX2", mass_values_ptr2,
663  flux_values_ptr2);
664  CHKERR push_stiff_lhs("MASS2", "FLUX2", data2, material_blocks);
665  if (nb_species == 3) {
666  CHKERR update_stiff_lhs("MASS3", "FLUX3", mass_values_ptr3,
667  flux_values_ptr3);
668  CHKERR push_stiff_lhs("MASS3", "FLUX3", data3, material_blocks);
669  }
670  }
671  }
672 
673 
674  // vol_ele_stiff_lhs->getOpPtrVector().push_back(
675  // new OpError(ExactFunction(), ExactFunctionLap(), ExactFunctionGrad(), data1, material_blocks, global_error));
676  CHKERR set_integration_rule();
677  dm = simple_interface->getDM();
678  ts = createTS(m_field.get_comm());
679  boost::shared_ptr<FaceEle> initial_mass_ele(new FaceEle(m_field));
680 
681  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
682  CHKERR apply_IC("MASS1", inner_surface1,
683  initial_mass_ele); // nb_species times
684  if (nb_species == 2 || nb_species == 3) {
685  CHKERR apply_IC("MASS2", inner_surface2, initial_mass_ele);
686  if (nb_species == 3) {
687  CHKERR apply_IC("MASS3", inner_surface3, initial_mass_ele);
688  }
689  }
690  }
691  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
692  initial_mass_ele);
693 
694  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
695  CHKERR apply_BC("FLUX1"); // nb_species times
696  if (nb_species == 2 || nb_species == 3) {
697  CHKERR apply_BC("FLUX2");
698  if (nb_species == 3) {
699  CHKERR apply_BC("FLUX3");
700  }
701  }
702  }
703 
704  CHKERR loop_fe(); // only once
705  post_proc->generateReferenceElementMesh(); // only once
706 
707 
708 
709  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
710  CHKERR post_proc_fields("MASS1", "FLUX1");
711  post_proc->addFieldValuesPostProc("ERROR");
712  if (nb_species == 2 || nb_species == 3) {
713  CHKERR post_proc_fields("MASS2", "FLUX2");
714  if (nb_species == 3) {
715  CHKERR post_proc_fields("MASS3", "FLUX3");
716  }
717  }
718  }
719 
720  //the double global_vector is a global_vector which sums errors in each element
721 
722  // Vec gVec;
723  // CHKERR VecCreateMPI(PETSC_COMM_WORLD, 1, 1, &gVec);
724  // double err_per_step = 0;
725 
726  // auto compute_global_error = [&gVec](double &g_err, double &err_out) {
727  // MoFEMFunctionBegin;
728 
729  // CHKERR VecZeroEntries(gVec);
730  // CHKERR VecAssemblyBegin(gVec);
731  // CHKERR VecAssemblyEnd(gVec);
732 
733  // int ind[1] = {0};
734  // CHKERR VecSetValues(gVec, 1, ind, &err_out, ADD_VALUES);
735  // CHKERR VecAssemblyBegin(gVec);
736  // CHKERR VecAssemblyEnd(gVec);
737 
738  // CHKERR VecGetValues(gVec, 1, ind, &err_out);
739 
740  // MoFEMFunctionReturn(0);
741  // };
742 
743  // Vec error_per_proc;
744  // CHKERR VecCreateMPI(cOmm, 1, PETSC_DECIDE, &error_per_proc);
745  // auto get_global_error = [&]() {
746  // MoFEMFunctionBegin;
747  // CHKERR VecSetValue(error_per_proc, m_field.get_comm_rank(), global_error, INSERT_VALUES);
748  // MoFEMFunctionReturn(0);
749  // };
750 
751  // auto reinitialize_global = [&]() {
752  // MoFEMFunctionBegin;
753  // CHKERR get_global_error();
754  // CHKERR VecAssemblyBegin(error_per_proc);
755  // CHKERR VecAssemblyEnd(error_per_proc);
756  // double error_sum;
757  // CHKERR VecSum(error_per_proc, &error_sum);
758  // CHKERR PetscPrintf(PETSC_COMM_WORLD, "Error : %3.4e \n",
759  // error_sum);
760  // global_error = 0;
761  // MoFEMFunctionReturn(0);
762  // };
763 
764  // vol_ele_stiff_lhs->postProcessHook = reinitialize_global;
765 
766 
767 
768  monitor_ptr = boost::shared_ptr<Monitor>(
769  new Monitor(cOmm, rAnk, dm, post_proc, global_error)); // nb_species times
770  CHKERR output_result(); // only once
771  CHKERR solve(); // only once
773 }
774 
775 int main(int argc, char *argv[]) {
776  const char param_file[] = "param_file.petsc";
777  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
778  try {
779  moab::Core mb_instance;
780  moab::Interface &moab = mb_instance;
781  MoFEM::Core core(moab);
782  DMType dm_name = "DMMOFEM";
783  CHKERR DMRegister_MoFEM(dm_name);
784 
785  int order = 1;
786  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
787  int nb_species = 1;
788  RDProblem reac_diff_problem(core, order+1);
789  CHKERR reac_diff_problem.run_analysis(nb_species);
790  }
791  CATCH_ERRORS;
793  return 0;
794 }
Calculate divergence of vector field.
const double sml
boost::shared_ptr< VectorDouble > flux_divs_ptr3
Deprecated interface functions.
Problem manager is used to build and partition problems.
MoFEMErrorCode add_fe()
const double ramp_t
MoFEMErrorCode update_stiff_lhs()
MoFEMErrorCode setup_system()
auto createTS
Definition: AuxPETSc.hpp:275
Section manager is used to create indexes and sectionsFIXME: ISManager is not properly testsed by ato...
Definition: ISManager.hpp:36
boost::shared_ptr< MatrixDouble > flux_values_ptr3
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:445
double operator()(const double x, const double y, const double t) const
MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
int main(int argc, char *argv[])
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode set_blockData(std::map< int, BlockData > &block_data_map)
Core (interface) class.
Definition: Core.hpp:50
MoFEMErrorCode update_slow_rhs(std::string mass_fiedl, boost::shared_ptr< VectorDouble > &mass_ptr)
keeps basic data about problemThis is low level structure with information about problem,...
EdgeElementForcesAndSourcesCoreBase BoundaryEle
MoFEMErrorCode push_slow_rhs()
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, 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 TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:772
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
MoFEMErrorCode apply_BC(std::string flux_field)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVectorsets the routine to create a global vector associated with the shell DM.
Definition: DMMMoFEM.cpp:1026
MoFEMErrorCode push_stiff_lhs()
MoFEMErrorCode update_vol_fe(boost::shared_ptr< FaceEle > &vol_ele, boost::shared_ptr< PreviousData > &data)
double operator()(const double x, const double y, const double t) const
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
MoFEMErrorCode post_proc_fields()
PetscErrorCode DMMoFEMTSSetRHSFunction(DM dm, 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 TS the right hand side function
Definition: DMMMoFEM.cpp:801
double operator()(const double x, const double y, const double t) const
MoFEMErrorCode push_stiff_rhs()
RDProblem(MoFEM::Core &core, const int order)
MoFEMErrorCode extract_initial_ents(int block_id, Range &surface)
Calculate inverse of jacobian for face element.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
MoFEMErrorCode set_integration_rule()
MoFEMErrorCode loop_fe()
field with continuous tangents
Definition: definitions.h:178
Get vector field for H-div approximation.
boost::shared_ptr< VectorDouble > mass_values_ptr3
MoFEMErrorCode run_analysis()
boost::shared_ptr< PreviousData > data3
#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.
#define CHKERR
Inline error check.
Definition: definitions.h:602
MoFEMErrorCode update_stiff_rhs()
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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.
Definition: DMMMoFEM.cpp:915
Range internal_edge_ents
constexpr int order
MoFEMErrorCode output_result()
Postprocess on face.
MoFEMErrorCode solve()
brief Transform local reference derivatives of shape function to global derivatives for face
MoFEMErrorCode apply_IC(std::string mass_field, Range &surface, boost::shared_ptr< FaceEle > &initial_ele, double &init_val)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Apply contravariant (Piola) transfer to Hdiv space on face.
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:719
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
Calculate jacobian for face element.
std::string getName() const
Get value at integration points for scalar field.
FTensor::Tensor1< double, 3 > operator()(const double x, const double y, const double t) const
OpCalculateScalarFieldValuesDot OpCalculateScalarValuesDot
double operator()(const double x) const
boost::shared_ptr< VectorDouble > mass_dots_ptr3
static char help[]
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:67
Make Hdiv space from Hcurl space in 2d.
field with C-1 continuity
Definition: definitions.h:180