v0.14.0
Loading...
Searching...
No Matches
shallow_wave.cpp
Go to the documentation of this file.
1/**
2 * \file shallow_wave.cpp
3 * \example shallow_wave.cpp
4 *
5 * Solving shallow wave equation on manifold
6 *
7 * The inital conditions are set following this paper \cite scott2016test.
8 *
9 */
10
11#include <MoFEM.hpp>
12using namespace MoFEM;
13
14#include <boost/math/quadrature/gauss_kronrod.hpp>
15using namespace boost::math::quadrature;
16
17template <int DIM> struct ElementsAndOps {};
18
19template <> struct ElementsAndOps<2> {
22};
23
24constexpr int FE_DIM = 2;
25
30
33
34// Use forms iterators for Grad-Grad term
39
40// Use forms for Mass term
45
49
51 GAUSS>::OpSource<1, 3>;
53 GAUSS>::OpSource<1, 1>;
54
61
62constexpr double omega = 7.292 * 1e-5;
63constexpr double g = 9.80616;
64constexpr double mu = 1e4;
65constexpr double h0 = 1e4;
66
67constexpr double h_hat = 120;
68constexpr double u_max = 80;
69constexpr double phi_0 = M_PI / 7;
70constexpr double phi_1 = M_PI / 2 - phi_0;
71constexpr double phi_2 = M_PI / 4;
72constexpr double alpha_montain = 1. / 3.;
73constexpr double beta_montain = 1. / 15.;
74
75constexpr double penalty = 1;
76
81
82struct OpURhs : public AssemblyDomainEleOp {
83
84 OpURhs(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
85 boost::shared_ptr<MatrixDouble> u_dot_ptr,
86 boost::shared_ptr<MatrixDouble> grad_u_ptr,
87 boost::shared_ptr<MatrixDouble> grad_h_ptr)
89 uPtr(u_ptr), uDotPtr(u_dot_ptr), uGradPtr(grad_u_ptr),
90 hGradPtr(grad_h_ptr) {}
91
94
95 const double vol = getMeasure();
96 auto t_w = getFTensor0IntegrationWeight();
97 auto t_row_base = row_data.getFTensor0N();
98 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
99 auto t_dot_u = getFTensor1FromMat<3>(*uDotPtr);
100 auto t_u = getFTensor1FromMat<3>(*uPtr);
101 auto t_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
102 auto t_grad_h = getFTensor1FromMat<3>(*hGradPtr);
103 auto t_coords = getFTensor1CoordsAtGaussPts();
104 auto t_normal = getFTensor1NormalsAtGaussPts();
105
106 for (int gg = 0; gg != nbIntegrationPts; gg++) {
107
108 const double alpha = t_w * vol;
109 auto t_nf = getFTensor1FromArray<3, 3>(locF);
110
113
114 const auto a = std::sqrt(t_coords(i) * t_coords(i));
115 const auto sin_fi = t_coords(2) / a;
116 const auto f = 2 * omega * sin_fi;
117
119 t_r(i) = t_normal(i);
120 t_r.normalize();
121
122 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
124 t_P(i, j) = t_r(i) * t_r(j);
125 t_Q(i, j) = t_kd(i, j) - t_P(i, j);
126
128 t_A(i, m) = levi_civita(i, j, m) * t_r(j);
129
130 t_rhs(m) = t_Q(m, i) * (t_dot_u(i) + t_grad_u(i, j) * t_u(j) +
131 f * t_A(i, j) * t_u(j) + g * t_grad_h(i));
132 t_rhs_grad(m, j) = t_Q(m, i) * (mu * t_grad_u(i, j));
133
134 t_rhs(m) += t_P(m, j) * t_u(j);
135
136 int rr = 0;
137 for (; rr != nbRows / 3; ++rr) {
138 t_nf(i) += alpha * t_row_base * t_rhs(i);
139 t_nf(i) += alpha * t_row_diff_base(j) * t_rhs_grad(i, j);
140 ++t_row_base;
141 ++t_row_diff_base;
142 ++t_nf;
143 }
144 for (; rr < nbRowBaseFunctions; ++rr) {
145 ++t_row_base;
146 ++t_row_diff_base;
147 }
148
149 ++t_w;
150 ++t_u;
151 ++t_dot_u;
152 ++t_grad_u;
153 ++t_grad_h;
154 ++t_coords;
155 ++t_normal;
156 }
157
159 }
160
161private:
162 boost::shared_ptr<MatrixDouble> uPtr;
163 boost::shared_ptr<MatrixDouble> uDotPtr;
164 boost::shared_ptr<MatrixDouble> uGradPtr;
165 boost::shared_ptr<MatrixDouble> hGradPtr;
166};
167
169
170 OpULhs_dU(const std::string field_name_row, const std::string field_name_col,
171 boost::shared_ptr<MatrixDouble> u_ptr,
172 boost::shared_ptr<MatrixDouble> grad_u_ptr)
173 : AssemblyDomainEleOp(field_name_row, field_name_col,
174 AssemblyDomainEleOp::OPROWCOL),
175 uPtr(u_ptr), uGradPtr(grad_u_ptr) {
176 this->sYmm = false;
177 }
178
180 EntitiesFieldData::EntData &col_data) {
182
183 const double vol = getMeasure();
184 auto t_w = getFTensor0IntegrationWeight();
185 auto t_row_base = row_data.getFTensor0N();
186 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
187 auto t_coords = getFTensor1CoordsAtGaussPts();
188 auto t_normal = getFTensor1NormalsAtGaussPts();
189
190 auto t_u = getFTensor1FromMat<3>(*uPtr);
191 auto t_grad_u = getFTensor2FromMat<3, 3>(*uGradPtr);
192
193 auto get_t_mat = [&](const int rr) {
195 &locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
196
197 &locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
198
199 &locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
200 };
201
202 const auto ts_a = getFEMethod()->ts_a;
203
204 for (int gg = 0; gg != nbIntegrationPts; gg++) {
205
206 const auto a = std::sqrt(t_coords(i) * t_coords(i));
207 const auto sin_fi = t_coords(2) / a;
208 const auto f = 2 * omega * sin_fi;
209
211 t_r(i) = t_normal(i);
212 t_r.normalize();
213
214 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
216 t_P(i, j) = t_r(i) * t_r(j);
217 t_Q(i, j) = t_kd(i, j) - t_P(i, j);
218
220 t_A(i, m) = levi_civita(i, j, m) * t_r(j);
221
223 t_rhs_du(m, j) =
224 t_Q(m, i) * (ts_a * t_kd(i, j) + t_grad_u(i, j) + f * t_A(i, j)) +
225 t_P(m, j);
226
227 const double alpha = t_w * vol;
228
229 int rr = 0;
230 for (; rr != nbRows / 3; rr++) {
231
232 auto t_col_base = col_data.getFTensor0N(gg, 0);
233 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
234 auto t_mat = get_t_mat(3 * rr);
235
236 for (int cc = 0; cc != nbCols / 3; cc++) {
237 t_mat(i, j) += (alpha * t_row_base * t_col_base) * t_rhs_du(i, j);
238 t_mat(i, j) += (alpha * mu) * t_Q(i, j) *
239 (t_row_diff_base(m) * t_col_diff_base(m));
240 ++t_col_diff_base;
241 ++t_col_base;
242 ++t_mat;
243 }
244 ++t_row_base;
245 ++t_row_diff_base;
246 }
247 for (; rr < nbRowBaseFunctions; ++rr) {
248 ++t_row_base;
249 ++t_row_diff_base;
250 }
251
252 ++t_w;
253 ++t_coords;
254 ++t_normal;
255 ++t_u;
256 ++t_grad_u;
257 }
258
260 }
261
262private:
263 boost::shared_ptr<MatrixDouble> uPtr;
264 boost::shared_ptr<MatrixDouble> uGradPtr;
265};
266
268
269 OpULhs_dH(const std::string field_name_row, const std::string field_name_col)
270 : AssemblyDomainEleOp(field_name_row, field_name_col,
271 AssemblyDomainEleOp::OPROWCOL) {
272 this->sYmm = false;
273 }
274
276 EntitiesFieldData::EntData &col_data) {
278
279 // get element volume
280 const double vol = getMeasure();
281 // get integration weights
282 auto t_w = getFTensor0IntegrationWeight();
283 // get base function gradient on rows
284 auto t_row_base = row_data.getFTensor0N();
285 // normal
286 auto t_normal = getFTensor1NormalsAtGaussPts();
287
288 auto get_t_vec = [&](const int rr) {
290 &locMat(rr + 0, 0),
291
292 &locMat(rr + 1, 0),
293
294 &locMat(rr + 2, 0)};
295 };
296
297 // loop over integration points
298 for (int gg = 0; gg != nbIntegrationPts; gg++) {
299
301 t_r(i) = t_normal(i);
302 t_r.normalize();
303
304 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
306 t_P(i, j) = t_r(i) * t_r(j);
307 t_Q(i, j) = t_kd(i, j) - t_P(i, j);
308
309 const double alpha = t_w * vol;
310
311 int rr = 0;
312 for (; rr != nbRows / 3; rr++) {
313 auto t_vec = get_t_vec(3 * rr);
314 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
315 const double a = alpha * g * t_row_base;
316
317 for (int cc = 0; cc != nbCols; cc++) {
318 t_vec(i) += a * (t_Q(i, m) * t_col_diff_base(m));
319 ++t_vec;
320 ++t_col_diff_base;
321 }
322
323 ++t_row_base;
324 }
325
326 for (; rr < nbRowBaseFunctions; ++rr)
327 ++t_row_base;
328
329 ++t_w;
330 ++t_normal;
331 }
332
334 }
335};
336
337struct Example {
338
339 Example(MoFEM::Interface &m_field) : mField(m_field) {}
340
342
343private:
345
354
355 boost::shared_ptr<FEMethod> domianLhsFEPtr;
356 boost::shared_ptr<FEMethod> domianRhsFEPtr;
357};
358
359//! [Create common data]
362
364}
365//! [Create common data]
366
367//! [Run problem]
379}
380//! [Run problem]
381
382//! [Read mesh]
386 CHKERR simple->getOptions();
387 CHKERR simple->loadFile();
389}
390//! [Read mesh]
391
392//! [Set up problem]
396 // Add field
397 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 3);
398 CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
399 CHKERR simple->addDataField("HO_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3);
400
401 int order = 2;
402 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
403 CHKERR simple->setFieldOrder("U", order);
404 CHKERR simple->setFieldOrder("H", order);
405 CHKERR simple->setFieldOrder("HO_POSITIONS", 3);
406
407 CHKERR simple->setUp();
409}
410//! [Set up problem]
411
412//! [Boundary condition]
415
416 PetscBool is_restart = PETSC_FALSE;
417 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_restart", &is_restart,
418 PETSC_NULL);
419
420 auto restart_vector = [&]() {
423 auto dm = simple->getDM();
424 MOFEM_LOG("SW", Sev::inform)
425 << "reading vector in binary from vector.dat ...";
426 PetscViewer viewer;
427 PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_READ,
428 &viewer);
429 auto T = createDMVector(simple->getDM());
430 VecLoad(T, viewer);
431 CHKERR DMoFEMMeshToLocalVector(dm, T, INSERT_VALUES, SCATTER_REVERSE);
433 };
434
435 if (is_restart) {
436
437 CHKERR restart_vector();
438
439 } else {
440
441 const double e_n = exp(-4 / pow(phi_1 - phi_0, 2));
442 const double u0 = u_max / e_n;
443
444 FTensor::Tensor1<double, 3> t_k{0., 0., 1.};
446 t_A(i, m) = levi_civita(i, j, m) * t_k(j);
447
448 auto get_phi = [&](const double x, const double y, const double z) {
449 FTensor::Tensor1<double, 3> t_r{x, y, 0.};
450 const double r = std::sqrt(t_r(i) * t_r(i));
451 return atan2(z, r);
452 };
453
454 auto init_u_phi = [&](const double phi) {
455 if (phi > phi_0 && phi < phi_1) {
456 return u0 * exp(1. / ((phi - phi_0) * (phi - phi_1)));
457 } else {
458 return 0.;
459 }
460 };
461
462 auto init_u = [&](const double x, const double y, const double z) {
463 FTensor::Tensor1<double, 3> t_u{0., 0., 0.};
464 const double u_phi = init_u_phi(get_phi(x, y, z));
465 if (u_phi > 0) {
466 FTensor::Tensor1<double, 3> t_r{x, y, 0.};
467 t_r.normalize();
468 t_u(i) = ((t_A(i, j) * t_r(j)) * u_phi);
469 }
470 return t_u;
471 };
472
473 auto init_h = [&](const double x, const double y, const double z) {
474 const double a = std::sqrt(x * x + y * y + z * z);
475
476 auto integral = [&](const double fi) {
477 const double u_phi = init_u_phi(fi);
478 const auto f = 2 * omega * sin(fi);
479 return a * u_phi * (f + (tan(fi) / a) * u_phi);
480 };
481
482 auto montain = [&](const double lambda, const double fi) {
483 if (lambda > -M_PI && lambda < M_PI)
484 return h_hat * cos(fi) * exp(-pow(lambda / alpha_montain, 2)) *
485 exp(-pow((phi_2 - fi) / beta_montain, 2));
486 else
487 return 0.;
488 };
489
490 const double fi = get_phi(x, y, z);
491 const double lambda = atan2(x, y);
492
493 double h1 = 0;
494 if (fi > phi_0)
495 h1 = gauss_kronrod<double, 32>::integrate(
496 integral, phi_0, fi, 0, std::numeric_limits<float>::epsilon());
497
498 return h0 + (montain(lambda, fi) - (h1 / g));
499 };
500
501 auto set_domain_general = [&](auto &pipeline) {
502 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
503 pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
504 pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
505 pipeline.push_back(new OpSetHOWeightsOnFace());
506 };
507
508 auto set_domain_rhs = [&](auto &pipeline) {
509 pipeline.push_back(new OpSourceU("U", init_u));
510 pipeline.push_back(new OpSourceH("H", init_h));
511 };
512
513 auto set_domain_lhs = [&](auto &pipeline) {
514 pipeline.push_back(
515 new OpMassUU("U", "U", [](double, double, double) { return 1; }));
516 pipeline.push_back(
517 new OpMassHH("H", "H", [](double, double, double) { return 1; }));
518 };
519
520 auto post_proc = [&]() {
523 auto dm = simple->getDM();
524
525 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
526
527 auto det_ptr = boost::make_shared<VectorDouble>();
528 auto jac_ptr = boost::make_shared<MatrixDouble>();
529 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
530
531 post_proc_fe->getOpPtrVector().push_back(
532 new OpGetHONormalsOnFace("HO_POSITIONS"));
533 post_proc_fe->getOpPtrVector().push_back(
534 new OpCalculateHOCoords("HO_POSITIONS"));
535 post_proc_fe->getOpPtrVector().push_back(
537 post_proc_fe->getOpPtrVector().push_back(
538 new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
539
540 auto u_ptr = boost::make_shared<MatrixDouble>();
541 auto h_ptr = boost::make_shared<VectorDouble>();
542 auto pos_ptr = boost::make_shared<MatrixDouble>();
543
544 post_proc_fe->getOpPtrVector().push_back(
545 new OpCalculateVectorFieldValues<3>("U", u_ptr));
546 post_proc_fe->getOpPtrVector().push_back(
547 new OpCalculateScalarFieldValues("H", h_ptr));
548 post_proc_fe->getOpPtrVector().push_back(
549 new OpCalculateVectorFieldValues<3>("HO_POSITIONS", pos_ptr));
550
552
553 post_proc_fe->getOpPtrVector().push_back(
554
555 new OpPPMap(post_proc_fe->getPostProcMesh(),
556 post_proc_fe->getMapGaussPts(),
557
558 {{"H", h_ptr}},
559
560 {{"U", u_ptr}, {"HO_POSITIONS", pos_ptr}},
561
562 {}, {}
563
564 )
565
566 );
567
568 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
569 CHKERR post_proc_fe->writeFile("out_init.h5m");
570
572 };
573
574 auto solve_init = [&]() {
577 auto pipeline_mng = mField.getInterface<PipelineManager>();
578
579 auto solver = pipeline_mng->createKSP();
580 CHKERR KSPSetFromOptions(solver);
581 PC pc;
582 CHKERR KSPGetPC(solver, &pc);
583 PetscBool is_pcfs = PETSC_FALSE;
584 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
585 if (is_pcfs == PETSC_TRUE) {
586 auto bc_mng = mField.getInterface<BcManager>();
587 auto name_prb = simple->getProblemName();
589 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
590 name_prb, ROW, "U", 0, 3, is_u);
591 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
592 CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_ADDITIVE);
593 }
594
595 CHKERR KSPSetUp(solver);
596
597 auto dm = simple->getDM();
598 auto D = createDMVector(dm);
599 auto F = vectorDuplicate(D);
600
601 CHKERR KSPSolve(solver, F, D);
602 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
603 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
604 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
605
607 };
608
609 auto pipeline_mng = mField.getInterface<PipelineManager>();
610
611 auto integration_rule = [](int, int, int approx_order) {
612 return 2 * approx_order + 4;
613 };
614 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
615 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
616
617 set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
618 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
619 set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
620 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
621
622 CHKERR solve_init();
623 CHKERR post_proc();
624
625 // Clear pipelines
626 pipeline_mng->getOpDomainRhsPipeline().clear();
627 pipeline_mng->getOpDomainLhsPipeline().clear();
628 }
629
631}
632//! [Boundary condition]
633
634//! [Push operators to pipeline]
637
638 // Push element from reference configuration to current configuration in 3d
639 // space
640 auto set_domain_general = [&](auto &pipeline) {
641 auto det_ptr = boost::make_shared<VectorDouble>();
642 auto jac_ptr = boost::make_shared<MatrixDouble>();
643 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
644 pipeline.push_back(new OpGetHONormalsOnFace("HO_POSITIONS"));
645 pipeline.push_back(new OpCalculateHOCoords("HO_POSITIONS"));
646 pipeline.push_back(new OpSetHOWeightsOnFace());
647 pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
648 pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
649 pipeline.push_back(new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
650 };
651
652 auto set_domain_rhs = [&](auto &pipeline) {
653 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
654 auto u_ptr = boost::make_shared<MatrixDouble>();
655 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
656 auto div_u_ptr = boost::make_shared<VectorDouble>();
657 auto dot_h_ptr = boost::make_shared<VectorDouble>();
658 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
659
660 pipeline.push_back(new OpCalculateVectorFieldValuesDot<3>("U", dot_u_ptr));
661 pipeline.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
662
663 pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
664 pipeline.push_back(
665 new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
666 pipeline.push_back(
668 pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
669
670 pipeline.push_back(new OpBaseTimesDotH(
671 "H", dot_h_ptr, [](double, double, double) { return 1.; }));
672 pipeline.push_back(new OpBaseTimesDivU(
673 "H", div_u_ptr, [](double, double, double) { return h0; }));
674 pipeline.push_back(new OpConvectiveH("H", u_ptr, grad_h_ptr));
675 pipeline.push_back(
676 new OpURhs("U", u_ptr, dot_u_ptr, grad_u_ptr, grad_h_ptr));
677 };
678
679 auto set_domain_lhs = [&](auto &pipeline) {
680 auto u_ptr = boost::make_shared<MatrixDouble>();
681 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
682 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
683
684 pipeline.push_back(new OpCalculateVectorFieldValues<3>("U", u_ptr));
685 pipeline.push_back(
686 new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
687 pipeline.push_back(new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
688
689 pipeline.push_back(new OpMassHH("H", "H", [&](double, double, double) {
690 return domianLhsFEPtr->ts_a;
691 }));
692 pipeline.push_back(new OpBaseDivU(
693 "H", "U", [](const double, const double, const double) { return h0; },
694 false, false));
695 pipeline.push_back(
696 new OpConvectiveH_dU("H", "U", grad_h_ptr, []() { return 1; }));
697 pipeline.push_back(
698 new OpConvectiveH_dGradH("H", "H", u_ptr, []() { return 1; }));
699 pipeline.push_back(new OpULhs_dU("U", "U", u_ptr, grad_u_ptr));
700 pipeline.push_back(new OpULhs_dH("U", "H"));
701 };
702
703 auto *pipeline_mng = mField.getInterface<PipelineManager>();
704
705 auto integration_rule = [](int, int, int approx_order) {
706 return 2 * approx_order + 4;
707 };
708 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
709 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
710
711 set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
712 set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
713
714 set_domain_rhs(pipeline_mng->getOpDomainRhsPipeline());
715 set_domain_lhs(pipeline_mng->getOpDomainLhsPipeline());
716
717 domianLhsFEPtr = pipeline_mng->getDomainLhsFE();
718 domianRhsFEPtr = pipeline_mng->getDomainRhsFE();
719
721}
722//! [Push operators to pipeline]
723
724/**
725 * @brief Monitor solution
726 *
727 * This functions is called by TS solver at the end of each step. It is used
728 * to output results to the hard drive.
729 */
730struct Monitor : public FEMethod {
731 Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
732 : dM(dm), postProc(post_proc){};
735 constexpr int save_every_nth_step = 50;
736 if (ts_step % save_every_nth_step == 0) {
738 CHKERR postProc->writeFile(
739 "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
740 MOFEM_LOG("SW", Sev::verbose)
741 << "writing vector in binary to vector.dat ...";
742 PetscViewer viewer;
743 PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_WRITE,
744 &viewer);
745 VecView(ts_u, viewer);
746 PetscViewerDestroy(&viewer);
747 }
749 }
750
751private:
753 boost::shared_ptr<PostProcEle> postProc;
754};
755
756//! [Solve]
759 auto *simple = mField.getInterface<Simple>();
760 auto *pipeline_mng = mField.getInterface<PipelineManager>();
761
762 auto dm = simple->getDM();
763
764 auto set_initial_step = [&](auto ts) {
766 int step = 0;
767 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-step", &step,
768 PETSC_NULL);
769 CHKERR TSSetStepNumber(ts, step);
771 };
772
773 auto set_fieldsplit_preconditioner_ksp = [&](auto ksp) {
775 PC pc;
776 CHKERR KSPGetPC(ksp, &pc);
777 PetscBool is_pcfs = PETSC_FALSE;
778 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
779 if (is_pcfs == PETSC_TRUE) {
780 auto bc_mng = mField.getInterface<BcManager>();
781 auto name_prb = simple->getProblemName();
783 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
784 name_prb, ROW, "U", 0, 3, is_u);
785 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_u);
786 }
788 };
789
790 // Setup postprocessing
791 auto get_fe_post_proc = [&]() {
792 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
793
794 auto det_ptr = boost::make_shared<VectorDouble>();
795 auto jac_ptr = boost::make_shared<MatrixDouble>();
796 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
797
798 post_proc_fe->getOpPtrVector().push_back(
799 new OpGetHONormalsOnFace("HO_POSITIONS"));
800 post_proc_fe->getOpPtrVector().push_back(
801 new OpCalculateHOCoords("HO_POSITIONS"));
802 post_proc_fe->getOpPtrVector().push_back(
804 post_proc_fe->getOpPtrVector().push_back(
805 new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
806 post_proc_fe->getOpPtrVector().push_back(
808
809 auto u_ptr = boost::make_shared<MatrixDouble>();
810 auto h_ptr = boost::make_shared<VectorDouble>();
811 auto pos_ptr = boost::make_shared<MatrixDouble>();
812
813 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
814 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
815
816 post_proc_fe->getOpPtrVector().push_back(
817 new OpCalculateVectorFieldValues<3>("U", u_ptr));
818 post_proc_fe->getOpPtrVector().push_back(
819 new OpCalculateScalarFieldValues("H", h_ptr));
820 post_proc_fe->getOpPtrVector().push_back(
821 new OpCalculateVectorFieldValues<3>("HO_POSITIONS", pos_ptr));
822
823 post_proc_fe->getOpPtrVector().push_back(
824 new OpCalculateVectorFieldGradient<3, 3>("U", grad_u_ptr));
825 post_proc_fe->getOpPtrVector().push_back(
826 new OpCalculateScalarFieldGradient<3>("H", grad_h_ptr));
827
829
830 post_proc_fe->getOpPtrVector().push_back(
831
832 new OpPPMap(
833 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
834
835 {{"H", h_ptr}},
836
837 {{"U", u_ptr}, {"HO_POSITIONS", pos_ptr}, {"GRAD_H", grad_h_ptr}},
838
839 {{"GRAD_U", grad_u_ptr}}, {}
840
841 )
842
843 );
844
845 return post_proc_fe;
846 };
847
848 auto set_fieldsplit_preconditioner_ts = [&](auto solver) {
850 SNES snes;
851 CHKERR TSGetSNES(solver, &snes);
852 KSP ksp;
853 CHKERR SNESGetKSP(snes, &ksp);
854 CHKERR set_fieldsplit_preconditioner_ksp(ksp);
856 };
857
859 ts = pipeline_mng->createTSIM();
860
861 boost::shared_ptr<FEMethod> null_fe;
862 auto monitor_ptr = boost::make_shared<Monitor>(dm, get_fe_post_proc());
863 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
864 null_fe, monitor_ptr);
865
866 // Add monitor to time solver
867 double ftime = 1;
868 // CHKERR TSSetMaxTime(ts, ftime);
869 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
870
871 auto T = createDMVector(simple->getDM());
872 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
873 SCATTER_FORWARD);
874 CHKERR TSSetSolution(ts, T);
875 CHKERR TSSetFromOptions(ts);
876 CHKERR set_fieldsplit_preconditioner_ts(ts);
877 CHKERR TSSetUp(ts);
878 CHKERR set_initial_step(ts);
879 CHKERR TSSolve(ts, NULL);
880 CHKERR TSGetTime(ts, &ftime);
881
883}
884//! [Solve]
885
886//! [Postprocess results]
890}
891//! [Postprocess results]
892
893//! [Check]
897}
898//! [Check]
899
900static char help[] = "...\n\n";
901
902int main(int argc, char *argv[]) {
903
904 // Initialisation of MoFEM/PETSc and MOAB data structures
905 const char param_file[] = "param_file.petsc";
907
908 // Add logging channel for example
909 auto core_log = logging::core::get();
910 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "SW"));
911 LogManager::setLog("SW");
912 MOFEM_LOG_TAG("SW", "example");
913
914 try {
915
916 //! [Register MoFEM discrete manager in PETSc]
917 DMType dm_name = "DMMOFEM";
918 CHKERR DMRegister_MoFEM(dm_name);
919 //! [Register MoFEM discrete manager in PETSc
920
921 //! [Create MoAB]
922 moab::Core mb_instance; ///< mesh database
923 moab::Interface &moab = mb_instance; ///< mesh database interface
924 //! [Create MoAB]
925
926 //! [Create MoFEM]
927 MoFEM::Core core(moab); ///< finite element database
928 MoFEM::Interface &m_field = core; ///< finite element database insterface
929 //! [Create MoFEM]
930
931 //! [Example]
932 Example ex(m_field);
933 CHKERR ex.runProblem();
934 //! [Example]
935 }
937
939}
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
@ ROW
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ 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 ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
const double init_u
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
static double phi
@ F
auto integration_rule
auto init_h
Initialisation function.
constexpr auto t_kd
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1003
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
static double lambda
double D
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::typ levi_civita)(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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 DMMoFEM.cpp:1042
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)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
const double u0
inital vale on blocksets
int r
Definition sdf.py:8
int save_every_nth_step
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< SPACE_DIM > OpBaseDivU
Definition seepage.cpp:79
constexpr int FE_DIM
constexpr double mu
constexpr double phi_0
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseTimesDotH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::TriLinearForm< GAUSS >::OpConvectiveTermLhsDy< 1, 1, 3 > OpConvectiveH_dGradH
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixVectorTimesGrad< 1, 3, 3 > OpBaseGradH
static char help[]
[Check]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::TriLinearForm< GAUSS >::OpConvectiveTermLhsDu< 1, 1, 3 > OpConvectiveH_dU
constexpr double phi_1
constexpr double penalty
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpConvectiveTermRhs< 1, 1, 3 > OpConvectiveH
constexpr double beta_montain
FTensor::Index< 'l', 3 > l
constexpr double h0
OpBaseTimesDotH OpBaseTimesDivU
constexpr double u_max
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMassUU
FTensor::Index< 'j', 3 > j
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixScalarTimesDiv< 3 > OpBaseDivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMassHH
constexpr double omega
FTensor::Index< 'i', 3 > i
constexpr double g
FTensor::Index< 'm', 3 > m
constexpr double h_hat
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpSourceH
constexpr double phi_2
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpSourceU
constexpr double alpha_montain
[Example]
Definition plastic.cpp:228
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition plastic.cpp:235
boost::shared_ptr< FEMethod > domianRhsFEPtr
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
boost::shared_ptr< FEMethod > domianLhsFEPtr
Simple interface for fast problem set-up.
Definition BcManager.hpp:25
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
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
structure for User Loop Methods on finite elements
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
Calculate HO coordinates at gauss points.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Post post-proc data at points from hash maps.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
intrusive_ptr for managing petsc objects
PetscInt ts_step
time step number
Vec & ts_u
state vector
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcFace > postProc
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
OpULhs_dH(const std::string field_name_row, const std::string field_name_col)
boost::shared_ptr< MatrixDouble > uPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > uGradPtr
OpULhs_dU(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > grad_u_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Class dedicated to integrate operator.
boost::shared_ptr< MatrixDouble > uDotPtr
boost::shared_ptr< MatrixDouble > uGradPtr
boost::shared_ptr< MatrixDouble > hGradPtr
boost::shared_ptr< MatrixDouble > uPtr
OpURhs(const std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > u_dot_ptr, boost::shared_ptr< MatrixDouble > grad_u_ptr, boost::shared_ptr< MatrixDouble > grad_h_ptr)