v0.15.0
Loading...
Searching...
No Matches
free_surface.cpp
Go to the documentation of this file.
1/**
2 * \file free_surface.cpp
3 * \example mofem/tutorials/vec-5/free_surface.cpp
4 *
5 * Using PipelineManager interface calculate the divergence of base functions,
6 * and integral of flux on the boundary. Since the h-div space is used, volume
7 * integral and boundary integral should give the same result.
8 *
9 * Implementation based on \cite lovric2019low
10 */
11
12#include <MoFEM.hpp>
13#include <petsc/private/petscimpl.h>
14
15using namespace MoFEM;
16
17static char help[] = "...\n\n";
18
19#undef PYTHON_INIT_SURFACE
20
21#ifdef PYTHON_INIT_SURFACE
22#include <boost/python.hpp>
23#include <boost/python/def.hpp>
24namespace bp = boost::python;
25
26struct SurfacePython {
27 SurfacePython() = default;
28 virtual ~SurfacePython() = default;
29
30 /**
31 * @brief Initialize Python surface evaluation from file
32 *
33 * @param py_file Path to Python file containing surface function
34 * Must contain a function named "surface" that takes
35 * (x, y, z, eta) coordinates and returns surface value
36 * @return MoFEMErrorCode Success/failure code
37 *
38 * @note The Python file must define a function called "surface" that
39 * evaluates the initial surface configuration at given coordinates
40 */
41 MoFEMErrorCode surfaceInit(const std::string py_file) {
43 try {
44
45 // create main module
46 auto main_module = bp::import("__main__");
47 mainNamespace = main_module.attr("__dict__");
48 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
49 // create a reference to python function
50 surfaceFun = mainNamespace["surface"];
51
52 } catch (bp::error_already_set const &) {
53 // print all other errors to stderr
54 PyErr_Print();
56 }
58 };
59
60 /**
61 * @brief Evaluate surface function at given coordinates
62 *
63 * @param x X-coordinate
64 * @param y Y-coordinate
65 * @param z Z-coordinate
66 * @param eta Interface thickness parameter
67 * @param s Output surface value (phase field value)
68 * @return MoFEMErrorCode Success/failure code
69 *
70 * @note This function calls the Python "surface" function with the given
71 * coordinates and returns the evaluated surface value
72 */
73 MoFEMErrorCode evalSurface(
74
75 double x, double y, double z, double eta, double &s
76
77 ) {
79 try {
80
81 // call python function
82 s = bp::extract<double>(surfaceFun(x, y, z, eta));
83
84 } catch (bp::error_already_set const &) {
85 // print all other errors to stderr
86 PyErr_Print();
88 }
90 }
91
92private:
93 bp::object mainNamespace;
94 bp::object surfaceFun;
95};
96
97static boost::weak_ptr<SurfacePython> surfacePythonWeakPtr;
98
99#endif
100
101int coord_type = EXECUTABLE_COORD_TYPE;
102
103constexpr int BASE_DIM = 1;
104constexpr int SPACE_DIM = 2;
105constexpr int U_FIELD_DIM = SPACE_DIM;
106
107constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
109 IntegrationType::GAUSS; //< selected integration type
110
111template <int DIM>
113
116using DomainEleOp = DomainEle::UserDataOperator;
119using BoundaryEleOp = BoundaryEle::UserDataOperator;
121using SideOp = SideEle::UserDataOperator;
122
126
128
132
141
148
149template <CoordinateTypes COORD_TYPE>
152
153// Flux is applied by Lagrange Multiplie
157
162
164
165// mesh refinement
166int order = 3; ///< approximation order
167int nb_levels = 4; //< number of refinement levels
168int refine_overlap = 4; //< mesh overlap while refine
169
170constexpr bool debug = true;
171
172auto get_start_bit = []() {
173 return nb_levels + 1;
174}; //< first refinement level for computational mesh
175auto get_current_bit = []() {
176 return 2 * get_start_bit() + 1;
177}; ///< dofs bit used to do calculations
178auto get_skin_parent_bit = []() { return 2 * get_start_bit() + 2; };
179auto get_skin_child_bit = []() { return 2 * get_start_bit() + 3; };
180auto get_projection_bit = []() { return 2 * get_start_bit() + 4; };
181auto get_skin_projection_bit = []() { return 2 * get_start_bit() + 5; };
182
183// FIXME: Set parameters from command line
184
185// Physical parameters
186double a0 = 980;
187double rho_m = 0.998;
188double mu_m = 0.010101 * 1e2;
189double rho_p = 0.0012;
190double mu_p = 0.000182 * 1e2;
191double lambda = 73; ///< surface tension
192double W = 0.25;
193
194// Model parameters
195double h = 0.03; // mesh size
196double eta = h;
197double eta2 = eta * eta;
198
199// Numerical parameters
200double md = 1e-2; // mobility
201double eps = 1e-12;
202double tol = std::numeric_limits<float>::epsilon();
203
204double rho_ave = (rho_p + rho_m) / 2;
205double rho_diff = (rho_p - rho_m) / 2;
206double mu_ave = (mu_p + mu_m) / 2;
207double mu_diff = (mu_p - mu_m) / 2;
208
209double kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
210
211auto integration_rule = [](int, int, int) { return 2 * order + 1; };
212
213/**
214 * @brief Coordinate system scaling factor
215 * @param r Radial coordinate
216 * @return Scaling factor for integration
217 *
218 * Returns 2πr for cylindrical coordinates, 1 for Cartesian.
219 * Used to scale integration measures and force contributions.
220 */
221auto cylindrical = [](const double r) {
222 if (coord_type == CYLINDRICAL)
223 return 2 * M_PI * r;
224 else
225 return 1.;
226};
227
228/**
229 * @brief Wetting angle sub-stepping for gradual application
230 * @param ts_step Current time step number
231 * @return Scaling factor [0,1] for gradual wetting angle application
232 *
233 * Gradually applies wetting angle boundary condition over first 16 steps
234 * to avoid sudden changes that could destabilize the solution.
235 */
236auto wetting_angle_sub_stepping = [](auto ts_step) {
237 constexpr int sub_stepping = 16;
238 return std::min(1., static_cast<double>(ts_step) / sub_stepping);
239};
240
241// Phase field cutoff and smoothing functions
242
243/**
244 * @brief Maximum function with smooth transition
245 * @param x Input value
246 * @return max(x, -1) with smooth transition
247 */
248auto my_max = [](const double x) { return (x - 1 + std::abs(x + 1)) / 2; };
249
250/**
251 * @brief Minimum function with smooth transition
252 * @param x Input value
253 * @return min(x, 1) with smooth transition
254 */
255auto my_min = [](const double x) { return (x + 1 - std::abs(x - 1)) / 2; };
256
257/**
258 * @brief Phase field cutoff function
259 * @param h Phase field value
260 * @return Constrained value in [-1, 1]
261 *
262 * Constrains phase field values to physical range [-1,1] with smooth cutoff
263 */
264auto cut_off = [](const double h) { return my_max(my_min(h)); };
265
266/**
267 * @brief Derivative of cutoff function
268 * @param h Phase field value
269 * @return Derivative of cutoff function
270 */
271auto d_cut_off = [](const double h) {
272 if (h >= -1 && h < 1)
273 return 1.;
274 else
275 return 0.;
276};
277
278/**
279 * @brief Phase-dependent material property interpolation
280 * @param h Phase field value (-1 to 1)
281 * @param diff Difference between phase values (phase_plus - phase_minus)
282 * @param ave Average of phase values (phase_plus + phase_minus)/2
283 * @return Interpolated material property
284 *
285 * Linear interpolation: property = diff * h + ave
286 * Used for density and viscosity interpolation between phases
287 */
288auto phase_function = [](const double h, const double diff, const double ave) {
289 return diff * h + ave;
290};
291
292/**
293 * @brief Derivative of phase function with respect to h
294 * @param h Phase field value (unused in linear case)
295 * @param diff Difference between phase values
296 * @return Derivative value
297 */
298auto d_phase_function_h = [](const double h, const double diff) {
299 return diff;
300};
301
302// Free energy and mobility functions for phase field model
303
304/**
305 * @brief Double-well potential function
306 * @param h Phase field value
307 * @return Free energy density f(h) = 4W*h*(h²-1)
308 *
309 * Double-well potential with minima at h = ±1 (pure phases)
310 * and maximum at h = 0 (interface). Controls interface structure.
311 */
312auto get_f = [](const double h) { return 4 * W * h * (h * h - 1); };
313
314/**
315 * @brief Derivative of double-well potential
316 * @param h Phase field value
317 * @return f'(h) = 4W*(3h²-1)
318 */
319auto get_f_dh = [](const double h) { return 4 * W * (3 * h * h - 1); };
320
321/**
322 * @brief Constant mobility function M₀
323 * @param h Phase field value (unused)
324 * @return Constant mobility value
325 */
326auto get_M0 = [](auto h) { return md; };
327
328/**
329 * @brief Derivative of constant mobility
330 * @param h Phase field value (unused)
331 * @return Zero (constant mobility)
332 */
333auto get_M0_dh = [](auto h) { return 0; };
334
335/**
336 * @brief Degenerate mobility function M₂
337 * @param h Phase field value
338 * @return M₂(h) = md*(1-h²)
339 *
340 * Mobility that vanishes at pure phases (h = ±1)
341 */
342auto get_M2 = [](auto h) {
343 return md * (1 - h * h);
344};
345
346/**
347 * @brief Derivative of degenerate mobility M₂
348 * @param h Phase field value
349 * @return M₂'(h) = -2*md*h
350 */
351auto get_M2_dh = [](auto h) { return -md * 2 * h; };
352
353/**
354 * @brief Non-linear mobility function M₃
355 * @param h Phase field value
356 * @return Piecewise cubic mobility function
357 *
358 * Smooth mobility function with different behavior for h >= 0 and h < 0
359 */
360auto get_M3 = [](auto h) {
361 const double h2 = h * h;
362 const double h3 = h2 * h;
363 if (h >= 0)
364 return md * (2 * h3 - 3 * h2 + 1);
365 else
366 return md * (-2 * h3 - 3 * h2 + 1);
367};
368
369/**
370 * @brief Derivative of non-linear mobility M₃
371 * @param h Phase field value
372 * @return M₃'(h) with sign-dependent cubic behavior
373 */
374auto get_M3_dh = [](auto h) {
375 if (h >= 0)
376 return md * (6 * h * (h - 1));
377 else
378 return md * (-6 * h * (h + 1));
379};
380
381// Select mobility model (currently using constant M₀)
382auto get_M = [](auto h) { return get_M0(h); };
383auto get_M_dh = [](auto h) { return get_M0_dh(h); };
384
385/**
386 * @brief Create deviatoric stress tensor
387 * @param A Coefficient for tensor scaling
388 * @return Fourth-order deviatoric tensor D_ijkl
389 *
390 * Constructs deviatoric part of fourth-order identity tensor:
391 * D_ijkl = A * (δ_ik δ_jl + δ_il δ_jk)/2
392 * Used in viscous stress calculations for fluid flow
393 */
394auto get_D = [](const double A) {
396 t_D(i, j, k, l) = A * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
397 return t_D;
398};
399
400// Initial condition functions for phase field
401
402/**
403 * @brief Oscillating interface initialization
404 * @param r Radial coordinate
405 * @param y Vertical coordinate
406 * @param unused Z-coordinate (unused in 2D)
407 * @return Phase field value (-1 to 1)
408 *
409 * Creates circular interface with sinusoidal perturbations:
410 * - Base radius R with amplitude A
411 * - n-fold symmetry oscillations
412 * - Uses tanh profile for smooth interface
413 */
414auto kernel_oscillation = [](double r, double y, double) {
415 constexpr int n = 3;
416 constexpr double R = 0.0125;
417 constexpr double A = R * 0.2;
418 const double theta = atan2(r, y);
419 const double w = R + A * cos(n * theta);
420 const double d = std::sqrt(r * r + y * y);
421 return tanh((w - d) / (eta * std::sqrt(2)));
422};
423
424/**
425 * @brief Eye-shaped interface initialization
426 * @param r Radial coordinate
427 * @param y Vertical coordinate
428 * @param unused Z-coordinate (unused in 2D)
429 * @return Phase field value (-1 to 1)
430 *
431 * Creates circular droplet centered at (0, y0) with radius R
432 */
433auto kernel_eye = [](double r, double y, double) {
434 constexpr double y0 = 0.5;
435 constexpr double R = 0.4;
436 y -= y0;
437 const double d = std::sqrt(r * r + y * y);
438 return tanh((R - d) / (eta * std::sqrt(2)));
439};
440
441/**
442 * @brief Capillary tube initialization
443 * @param x X-coordinate (unused)
444 * @param y Y-coordinate
445 * @param z Z-coordinate (unused)
446 * @return Phase field value (-1 to 1)
447 *
448 * Creates horizontal interface at specified water height
449 */
450auto capillary_tube = [](double x, double y, double z) {
451 constexpr double water_height = 0.;
452 return tanh((water_height - y) / (eta * std::sqrt(2)));
453 ;
454};
455
456/**
457 * @brief Bubble device initialization
458 * @param x X-coordinate
459 * @param y Y-coordinate (unused)
460 * @param z Z-coordinate (unused)
461 * @return Phase field value (-1 to 1)
462 *
463 * Creates vertical interface for bubble formation device
464 */
465auto bubble_device = [](double x, double y, double z) {
466 return -tanh((-0.039 - x) / (eta * std::sqrt(2)));
467};
468
469/**
470 * @brief Initialisation function
471 *
472 * @note If UMs are compiled with Python to initialise phase field "H"
473 * surface.py function is used, which has to be present in execution folder.
474 *
475 */
476auto init_h = [](double r, double y, double theta) {
477#ifdef PYTHON_INIT_SURFACE
478 double s = 1;
479 if (auto ptr = surfacePythonWeakPtr.lock()) {
480 CHK_THROW_MESSAGE(ptr->evalSurface(r, y, theta, eta, s),
481 "error eval python");
482 }
483 return s;
484#else
485 // return bubble_device(r, y, theta);
486 return capillary_tube(r, y, theta);
487 // return kernel_eye(r, y, theta);
488#endif
489};
490
491/**
492 * @brief Wetting angle function (placeholder)
493 * @param water_level Current water level
494 * @return Wetting angle value
495 *
496 * Currently returns input value directly. Can be modified to implement
497 * complex wetting angle dependencies on interface position or time.
498 */
499auto wetting_angle = [](double water_level) { return water_level; };
500
501/**
502 * @brief Create bit reference level
503 * @param b Bit number to set
504 * @return BitRefLevel with specified bit set
505 */
506auto bit = [](auto b) { return BitRefLevel().set(b); };
507
508/**
509 * @brief Create marker bit reference level
510 * @param b Bit number from end
511 * @return BitRefLevel with bit set from end
512 */
513auto marker = [](auto b) { return BitRefLevel().set(BITREFLEVEL_SIZE - b); };
514
515/**
516 * @brief Get bit reference level from finite element
517 * @param fe_ptr Pointer to finite element method
518 * @return Current bit reference level of the element
519 */
520auto get_fe_bit = [](FEMethod *fe_ptr) {
521 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel();
522};
523
524/**
525 * @brief Get global size across all processors
526 * @param l_size Local size on current processor
527 * @return Global sum of sizes across all processors
528 */
529auto get_global_size = [](int l_size) {
530 int g_size;
531 MPI_Allreduce(&l_size, &g_size, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
532 return g_size;
533};
534
535/**
536 * @brief Save range of entities to file
537 * @param moab MOAB interface for mesh operations
538 * @param name Output filename
539 * @param r Range of entities to save
540 * @return MoFEMErrorCode Success/failure code
541 *
542 * Saves entities to HDF5 file if range is non-empty globally
543 */
544auto save_range = [](moab::Interface &moab, const std::string name,
545 const Range r) {
547 if (get_global_size(r.size())) {
548 auto out_meshset = get_temp_meshset_ptr(moab);
549 CHKERR moab.add_entities(*out_meshset, r);
550 CHKERR moab.write_file(name.c_str(), "MOAB", "PARALLEL=WRITE_PART",
551 out_meshset->get_ptr(), 1);
552 }
554};
555
556/**
557 * @brief Get entities of DOFs by field name - used for debugging
558 * @param dm PETSc DM object containing the problem
559 * @param field_name Name of field to extract entities for
560 * @return Range of entities containing DOFs for specified field
561 *
562 * Extracts all mesh entities that have degrees of freedom for the
563 * specified field. Useful for debugging and visualization of field
564 * distributions across the mesh.
565 */
566auto get_dofs_ents_by_field_name = [](auto dm, auto field_name) {
567 auto prb_ptr = getProblemPtr(dm);
568 std::vector<EntityHandle> ents_vec;
569
570 MoFEM::Interface *m_field_ptr;
571 CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
572
573 auto bit_number = m_field_ptr->get_field_bit_number(field_name);
574 auto dofs = prb_ptr->numeredRowDofsPtr;
575 auto lo_it = dofs->lower_bound(FieldEntity::getLoBitNumberUId(bit_number));
576 auto hi_it = dofs->upper_bound(FieldEntity::getHiBitNumberUId(bit_number));
577 ents_vec.reserve(std::distance(lo_it, hi_it));
578
579 for (; lo_it != hi_it; ++lo_it) {
580 ents_vec.push_back((*lo_it)->getEnt());
581 }
582
583 std::sort(ents_vec.begin(), ents_vec.end());
584 auto it = std::unique(ents_vec.begin(), ents_vec.end());
585 Range r;
586 r.insert_list(ents_vec.begin(), it);
587 return r;
588};
589
590/**
591 * @brief Get all entities with DOFs in the problem - used for debugging
592 * @param dm PETSc DM object containing the problem
593 * @return Range of all entities containing any DOFs
594 *
595 * Extracts all mesh entities that have degrees of freedom for any field.
596 * Useful for debugging mesh partitioning and DOF distribution.
597 */
598auto get_dofs_ents_all = [](auto dm) {
599 auto prb_ptr = getProblemPtr(dm);
600 std::vector<EntityHandle> ents_vec;
601
602 MoFEM::Interface *m_field_ptr;
603 CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
604
605 auto dofs = prb_ptr->numeredRowDofsPtr;
606 ents_vec.reserve(dofs->size());
607
608 for (auto d : *dofs) {
609 ents_vec.push_back(d->getEnt());
610 }
611
612 std::sort(ents_vec.begin(), ents_vec.end());
613 auto it = std::unique(ents_vec.begin(), ents_vec.end());
614 Range r;
615 r.insert_list(ents_vec.begin(), it);
616 return r;
617};
618
619#include <FreeSurfaceOps.hpp>
620using namespace FreeSurfaceOps;
621
622struct FreeSurface;
623
624enum FR { F, R }; // F - forward, and reverse
625
626/**
627 * @brief Set of functions called by PETSc solver used to refine and update
628 * mesh.
629 *
630 * @note Currently theta method is only handled by this code.
631 *
632 */
633struct TSPrePostProc {
634
635 TSPrePostProc() = default;
636 virtual ~TSPrePostProc() = default;
637
638 /**
639 * @brief Used to setup TS solver
640 *
641 * @param ts PETSc time stepping object to configure
642 * @return MoFEMErrorCode Success/failure code
643 *
644 * Sets up time stepping solver with custom preconditioner, monitors,
645 * and callback functions for mesh refinement and solution projection
646 */
648
649 /**
650 * @brief Get scatter context for vector operations
651 *
652 * @param x Local sub-vector
653 * @param y Global vector
654 * @param fr Direction flag (F=forward, R=reverse)
655 * @return SmartPetscObj<VecScatter> Scatter context for vector operations
656 *
657 * Creates scatter context to transfer data between global and sub-problem vectors
658 */
659 SmartPetscObj<VecScatter> getScatter(Vec x, Vec y, enum FR fr);
660
661 /**
662 * @brief Create sub-problem vector
663 *
664 * @return SmartPetscObj<Vec> Vector compatible with sub-problem DM
665 *
666 * Creates a vector with proper size and ghost structure for sub-problem
667 */
669
673
674private:
675 /**
676 * @brief Pre process time step
677 *
678 * Refine mesh and update fields before each time step
679 *
680 * @param ts PETSc time stepping object
681 * @return MoFEMErrorCode Success/failure code
682 *
683 * This function is called before each time step to:
684 * - Apply phase field cutoff constraints
685 * - Refine mesh based on interface location
686 * - Project solution data to new mesh
687 * - Update solver operators
688 */
689 static MoFEMErrorCode tsPreProc(TS ts);
690
691 /**
692 * @brief Post process time step
693 *
694 * Currently this function does not make anything major
695 *
696 * @param ts PETSc time stepping object
697 * @return MoFEMErrorCode Success/failure code
698 *
699 * Called after each time step completion for cleanup operations
700 */
701 static MoFEMErrorCode tsPostProc(TS ts);
702
703 /**
704 * @brief Pre-stage processing for time stepping
705 *
706 * @param ts PETSc time stepping object
707 * @return MoFEMErrorCode Success/failure code
708 */
710
711 /**
712 * @brief Set implicit function for time stepping
713 *
714 * @param ts PETSc time stepping object
715 * @param t Current time
716 * @param u Solution vector at current time
717 * @param u_t Time derivative of solution
718 * @param f Output residual vector F(t,u,u_t)
719 * @param ctx User context (unused)
720 * @return MoFEMErrorCode Success/failure code
721 *
722 * Wrapper that scatters global vectors to sub-problem and evaluates residual
723 */
724 static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
725 Vec f,
726 void *ctx); //< Wrapper for SNES Rhs
727
728 /**
729 * @brief Set implicit Jacobian for time stepping
730 *
731 * @param ts PETSc time stepping object
732 * @param t Current time
733 * @param u Solution vector at current time
734 * @param u_t Time derivative of solution
735 * @param a Shift parameter for implicit methods
736 * @param A Jacobian matrix (input)
737 * @param B Jacobian matrix (output, often same as A)
738 * @param ctx User context (unused)
739 * @return MoFEMErrorCode Success/failure code
740 *
741 * Wrapper that assembles Jacobian matrix for sub-problem
742 */
743 static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
744 PetscReal a, Mat A, Mat B,
745 void *ctx); ///< Wrapper for SNES Lhs
746
747 /**
748 * @brief Monitor solution during time stepping
749 *
750 * @param ts PETSc time stepping object
751 * @param step Current time step number
752 * @param t Current time value
753 * @param u Current solution vector
754 * @param ctx User context (pointer to FreeSurface)
755 * @return MoFEMErrorCode Success/failure code
756 *
757 * Called after each time step to monitor solution and save output
758 */
759 static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u,
760 void *ctx); ///< Wrapper for TS monitor
761
762 /**
763 * @brief Setup preconditioner
764 *
765 * @param pc PETSc preconditioner object
766 * @return MoFEMErrorCode Success/failure code
767 *
768 * Initializes KSP solver for shell preconditioner
769 */
770 static MoFEMErrorCode pcSetup(PC pc);
771
772 /**
773 * @brief Apply preconditioner
774 *
775 * @param pc PETSc preconditioner object
776 * @param pc_f Input vector (right-hand side)
777 * @param pc_x Output vector (preconditioned solution)
778 * @return MoFEMErrorCode Success/failure code
779 *
780 * Applies preconditioner by solving sub-problem with KSP
781 */
782 static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x);
783
784 SmartPetscObj<Vec> globRes; //< global residual
785 SmartPetscObj<Mat> subB; //< sub problem tangent matrix
786 SmartPetscObj<KSP> subKSP; //< sub problem KSP solver
787
788 boost::shared_ptr<SnesCtx>
789 snesCtxPtr; //< internal data (context) for MoFEM SNES functions
790 boost::shared_ptr<TsCtx>
791 tsCtxPtr; //< internal data (context) for MoFEM TS functions.
792};
793
794static boost::weak_ptr<TSPrePostProc> tsPrePostProc;
795
797
798 /**
799 * @brief Constructor
800 *
801 * @param m_field Reference to MoFEM interface for finite element operations
802 */
803 FreeSurface(MoFEM::Interface &m_field) : mField(m_field) {}
804
805 /**
806 * @brief Main function to run the complete free surface simulation
807 *
808 * @return MoFEMErrorCode Success/failure code
809 *
810 * Executes the complete simulation workflow:
811 * 1. Read mesh from file
812 * 2. Setup problem (fields, approximation orders, parameters)
813 * 3. Apply boundary conditions and initialize fields
814 * 4. Assemble system matrices and operators
815 * 5. Solve time-dependent problem
816 */
818
819 /**
820 * @brief Create refined problem for mesh adaptation
821 *
822 * @return MoFEMErrorCode Success/failure code
823 *
824 * Creates mesh refinement data structures needed for adaptive meshing
825 */
827
829
830private:
831 /**
832 * @brief Read mesh from input file
833 *
834 * @return MoFEMErrorCode Success/failure code
835 *
836 * Loads mesh using Simple interface and sets up parent adjacencies
837 * for hierarchical mesh refinement
838 */
840
841 /**
842 * @brief Setup problem fields and parameters
843 *
844 * @return MoFEMErrorCode Success/failure code
845 *
846 * Creates finite element fields:
847 * - U: Velocity field (H1, vector)
848 * - P: Pressure field (H1, scalar)
849 * - H: Phase/order field (H1, scalar)
850 * - G: Chemical potential (H1, scalar)
851 * - L: Lagrange multiplier for boundary conditions (H1, scalar)
852 *
853 * Sets approximation orders and reads physical parameters from command line
854 */
856
857 /**
858 * @brief Apply boundary conditions and initialize fields
859 *
860 * @return MoFEMErrorCode Success/failure code
861 *
862 * - Initializes phase field using analytical or Python functions
863 * - Solves initialization problem for consistent initial conditions
864 * - Performs mesh refinement based on interface location
865 * - Removes DOFs for boundary conditions (symmetry, fixed, etc.)
866 */
868
869 /**
870 * @brief Project solution data between mesh levels
871 *
872 * @return MoFEMErrorCode Success/failure code
873 *
874 * Projects field data from coarse to fine mesh levels during refinement.
875 * Handles both time stepping vectors (for theta method) and regular fields.
876 * Uses L2 projection with mass matrix assembly.
877 */
879
880 /**
881 * @brief Assemble system operators and matrices
882 *
883 * @return MoFEMErrorCode Success/failure code
884 *
885 * Sets up finite element operators for:
886 * - Domain integration (momentum, phase field, mass conservation)
887 * - Boundary integration (surface tension, wetting angle, Lagrange multipliers)
888 * - Parent-child mesh hierarchies for adaptive refinement
889 */
891
892 /**
893 * @brief Solve the time-dependent free surface problem
894 *
895 * @return MoFEMErrorCode Success/failure code
896 *
897 * Creates and configures time stepping solver with:
898 * - Sub-problem DM for refined mesh
899 * - Post-processing for output visualization
900 * - Custom preconditioner and monitors
901 * - TSTheta implicit time integration
902 */
904
905 /// @brief Find entities on refinement levels
906 /// @param overlap Level of overlap around phase interface
907 /// @return Vector of entity ranges for each refinement level
908 ///
909 /// Identifies mesh entities that cross the phase interface by analyzing
910 /// phase field values at vertices. Returns entities that need refinement
911 /// at each hierarchical level with specified overlap.
912 std::vector<Range> findEntitiesCrossedByPhaseInterface(size_t overlap);
913
914 /// @brief Find parent entities that need refinement
915 /// @param ents Child entities requiring refinement
916 /// @param level Bit level for entity marking
917 /// @param mask Bit mask for filtering
918 /// @return Range of parent entities to refine
919 ///
920 /// Traverses mesh hierarchy to find parent entities that should be
921 /// refined to accommodate interface tracking requirements.
923
924 /// @brief Perform adaptive mesh refinement
925 /// @param overlap Number of element layers around interface to refine
926 /// @return MoFEMErrorCode Success/failure code
927 ///
928 /// Performs hierarchical mesh refinement around the phase interface:
929 /// - Identifies entities crossing interface
930 /// - Creates refinement hierarchy
931 /// - Sets up skin elements between levels
932 /// - Updates bit level markings for computation
933 MoFEMErrorCode refineMesh(size_t overlap);
934
935 /// @brief Create hierarchy of elements run on parent levels
936 /// @param fe_top Pipeline element to which hierarchy is attached
937 /// @param field_name Name of field for DOF extraction ("" for all fields)
938 /// @param op Type of operator OPSPACE, OPROW, OPCOL or OPROWCOL
939 /// @param child_ent_bit Bit level marking child entities
940 /// @param get_elem Lambda function to create element on hierarchy
941 /// @param verbosity Verbosity level for debugging output
942 /// @param sev Severity level for logging
943 /// @return MoFEMErrorCode Success/failure code
944 ///
945 /// Sets up parent-child relationships for hierarchical mesh refinement.
946 /// Allows access to field data from parent mesh levels during computation
947 /// on refined child meshes. Essential for projection and interpolation.
949 boost::shared_ptr<FEMethod> fe_top, std::string field_name,
950 ForcesAndSourcesCore::UserDataOperator::OpType op,
951 BitRefLevel child_ent_bit,
952 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
953 int verbosity, LogManager::SeverityLevel sev);
954
955 friend struct TSPrePostProc;
956};
957
958//! [Run programme]
968//! [Run programme]
969
970//! [Read mesh]
973 MOFEM_LOG("FS", Sev::inform) << "Read mesh for problem";
975
977 simple->getBitRefLevel() = BitRefLevel();
978
979 CHKERR simple->getOptions();
980 CHKERR simple->loadFile();
981
983}
984//! [Read mesh]
985
986//! [Set up problem]
989
990 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
991 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-nb_levels", &nb_levels,
992 PETSC_NULLPTR);
993 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-refine_overlap", &refine_overlap,
994 PETSC_NULLPTR);
995
996 const char *coord_type_names[] = {"cartesian", "polar", "cylindrical",
997 "spherical"};
998 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-coords", coord_type_names,
999 LAST_COORDINATE_SYSTEM, &coord_type, PETSC_NULLPTR);
1000
1001 MOFEM_LOG("FS", Sev::inform) << "Approximation order = " << order;
1002 MOFEM_LOG("FS", Sev::inform)
1003 << "Number of refinement levels nb_levels = " << nb_levels;
1004 nb_levels += 1;
1005
1006 auto simple = mField.getInterface<Simple>();
1007
1008 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "Acceleration", "-a0", &a0, PETSC_NULLPTR);
1009 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "\"Minus\" phase density rho_m",
1010 "-rho_m", &rho_m, PETSC_NULLPTR);
1011 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "\"Minus\" phase viscosity", "-mu_m",
1012 &mu_m, PETSC_NULLPTR);
1013 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "\"Plus\" phase density", "-rho_p",
1014 &rho_p, PETSC_NULLPTR);
1015 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "\"Plus\" phase viscosity", "-mu_p",
1016 &mu_p, PETSC_NULLPTR);
1017 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "Surface tension", "-lambda",
1018 &lambda, PETSC_NULLPTR);
1019 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR,
1020 "Height of the well in energy functional", "-W",
1021 &W, PETSC_NULLPTR);
1022 rho_ave = (rho_p + rho_m) / 2;
1023 rho_diff = (rho_p - rho_m) / 2;
1024 mu_ave = (mu_p + mu_m) / 2;
1025 mu_diff = (mu_p - mu_m) / 2;
1026
1027 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-h", &h, PETSC_NULLPTR);
1028 eta = h;
1029 eta2 = eta * eta;
1030 kappa = (3. / (4. * std::sqrt(2. * W))) * (lambda / eta);
1031
1032 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-md", &eta, PETSC_NULLPTR);
1033
1034 MOFEM_LOG("FS", Sev::inform) << "Acceleration a0 = " << a0;
1035 MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase density rho_m = " << rho_m;
1036 MOFEM_LOG("FS", Sev::inform) << "\"Minus\" phase viscosity mu_m = " << mu_m;
1037 MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase density rho_p = " << rho_p;
1038 MOFEM_LOG("FS", Sev::inform) << "\"Plus\" phase viscosity mu_p = " << mu_p;
1039 MOFEM_LOG("FS", Sev::inform) << "Surface tension lambda = " << lambda;
1040 MOFEM_LOG("FS", Sev::inform)
1041 << "Height of the well in energy functional W = " << W;
1042 MOFEM_LOG("FS", Sev::inform) << "Characteristic mesh size h = " << h;
1043 MOFEM_LOG("FS", Sev::inform) << "Mobility md = " << md;
1044
1045 MOFEM_LOG("FS", Sev::inform) << "Average density rho_ave = " << rho_ave;
1046 MOFEM_LOG("FS", Sev::inform) << "Difference density rho_diff = " << rho_diff;
1047 MOFEM_LOG("FS", Sev::inform) << "Average viscosity mu_ave = " << mu_ave;
1048 MOFEM_LOG("FS", Sev::inform) << "Difference viscosity mu_diff = " << mu_diff;
1049 MOFEM_LOG("FS", Sev::inform) << "kappa = " << kappa;
1050
1051
1052 // Fields on domain
1053
1054 // Velocity field
1055 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, U_FIELD_DIM);
1056 // Pressure field
1057 CHKERR simple->addDomainField("P", H1, AINSWORTH_LEGENDRE_BASE, 1);
1058 // Order/phase fild
1059 CHKERR simple->addDomainField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
1060 // Chemical potential
1061 CHKERR simple->addDomainField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
1062
1063 // Field on boundary
1064 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE,
1065 U_FIELD_DIM);
1066 CHKERR simple->addBoundaryField("H", H1, AINSWORTH_LEGENDRE_BASE, 1);
1067 CHKERR simple->addBoundaryField("G", H1, AINSWORTH_LEGENDRE_BASE, 1);
1068 // Lagrange multiplier which constrains slip conditions
1069 CHKERR simple->addBoundaryField("L", H1, AINSWORTH_LEGENDRE_BASE, 1);
1070
1071 CHKERR simple->setFieldOrder("U", order);
1072 CHKERR simple->setFieldOrder("P", order - 1);
1073 CHKERR simple->setFieldOrder("H", order);
1074 CHKERR simple->setFieldOrder("G", order);
1075 CHKERR simple->setFieldOrder("L", order);
1076
1077 // Initialise bit ref levels
1078 auto set_problem_bit = [&]() {
1080 // Set bits to build adjacencies between parents and children. That is
1081 // used by simple interface.
1082 simple->getBitAdjEnt() = BitRefLevel().set();
1083 simple->getBitAdjParent() = BitRefLevel().set();
1084 simple->getBitRefLevel() = BitRefLevel().set();
1085 simple->getBitRefLevelMask() = BitRefLevel().set();
1087 };
1088
1089 CHKERR set_problem_bit();
1090
1091 CHKERR simple->setUp();
1092
1094}
1095//! [Set up problem]
1096
1097//! [Boundary condition]
1100
1101#ifdef PYTHON_INIT_SURFACE
1102 auto get_py_surface_init = []() {
1103 auto py_surf_init = boost::make_shared<SurfacePython>();
1104 CHKERR py_surf_init->surfaceInit("surface.py");
1105 surfacePythonWeakPtr = py_surf_init;
1106 return py_surf_init;
1107 };
1108 auto py_surf_init = get_py_surface_init();
1109#endif
1110
1111 auto simple = mField.getInterface<Simple>();
1112 auto pip_mng = mField.getInterface<PipelineManager>();
1113 auto bc_mng = mField.getInterface<BcManager>();
1114 auto bit_mng = mField.getInterface<BitRefManager>();
1115 auto dm = simple->getDM();
1116
1118
1119 auto reset_bits = [&]() {
1121 BitRefLevel start_mask;
1122 for (auto s = 0; s != get_start_bit(); ++s)
1123 start_mask[s] = true;
1124 // reset bit ref levels
1125 CHKERR bit_mng->lambdaBitRefLevel(
1126 [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
1127 Range level0;
1128 CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
1129 CHKERR bit_mng->setNthBitRefLevel(level0, get_current_bit(), true);
1130 CHKERR bit_mng->setNthBitRefLevel(level0, get_projection_bit(), true);
1132 };
1133
1134 auto add_parent_field = [&](auto fe, auto op, auto field) {
1135 return setParentDofs(
1136 fe, field, op, bit(get_skin_parent_bit()),
1137
1138 [&]() {
1139 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1140 new DomainParentEle(mField));
1141 return fe_parent;
1142 },
1143
1144 QUIET, Sev::noisy);
1145 };
1146
1147 auto h_ptr = boost::make_shared<VectorDouble>();
1148 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
1149 auto g_ptr = boost::make_shared<VectorDouble>();
1150 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
1151
1152 auto set_generic = [&](auto fe) {
1154 auto &pip = fe->getOpPtrVector();
1155
1157
1159 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1160
1161 [&]() {
1162 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1163 new DomainParentEle(mField));
1165 fe_parent->getOpPtrVector(), {H1});
1166 return fe_parent;
1167 },
1168
1169 QUIET, Sev::noisy);
1170
1171 CHKERR add_parent_field(fe, UDO::OPROW, "H");
1172 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1173 pip.push_back(
1174 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
1175
1176 CHKERR add_parent_field(fe, UDO::OPROW, "G");
1177 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1178 pip.push_back(
1179 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
1180
1182 };
1183
1184 auto post_proc = [&](auto exe_test) {
1186 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
1187 post_proc_fe->exeTestHook = exe_test;
1188
1189 CHKERR set_generic(post_proc_fe);
1190
1192
1193 post_proc_fe->getOpPtrVector().push_back(
1194
1195 new OpPPMap(post_proc_fe->getPostProcMesh(),
1196 post_proc_fe->getMapGaussPts(),
1197
1198 {{"H", h_ptr}, {"G", g_ptr}},
1199
1200 {{"GRAD_H", grad_h_ptr}, {"GRAD_G", grad_g_ptr}},
1201
1202 {},
1203
1204 {}
1205
1206 )
1207
1208 );
1209
1210 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
1211 CHKERR post_proc_fe->writeFile("out_init.h5m");
1212
1214 };
1215
1216 auto solve_init = [&](auto exe_test) {
1218
1219 pip_mng->getOpDomainRhsPipeline().clear();
1220 pip_mng->getOpDomainLhsPipeline().clear();
1221
1222 auto set_domain_rhs = [&](auto fe) {
1224 CHKERR set_generic(fe);
1225 auto &pip = fe->getOpPtrVector();
1226
1227 CHKERR add_parent_field(fe, UDO::OPROW, "H");
1228 pip.push_back(new OpRhsH<true>("H", nullptr, nullptr, h_ptr, grad_h_ptr,
1229 grad_g_ptr));
1230 CHKERR add_parent_field(fe, UDO::OPROW, "G");
1231 pip.push_back(new OpRhsG<true>("G", h_ptr, grad_h_ptr, g_ptr));
1233 };
1234
1235 auto set_domain_lhs = [&](auto fe) {
1237
1238 CHKERR set_generic(fe);
1239 auto &pip = fe->getOpPtrVector();
1240
1241 CHKERR add_parent_field(fe, UDO::OPROW, "H");
1242 CHKERR add_parent_field(fe, UDO::OPCOL, "H");
1243 pip.push_back(new OpLhsH_dH<true>("H", nullptr, h_ptr, grad_g_ptr));
1244
1245 CHKERR add_parent_field(fe, UDO::OPCOL, "G");
1246 pip.push_back(new OpLhsH_dG<true>("H", "G", h_ptr));
1247
1248 CHKERR add_parent_field(fe, UDO::OPROW, "G");
1249 pip.push_back(new OpLhsG_dG("G"));
1250
1251 CHKERR add_parent_field(fe, UDO::OPCOL, "H");
1252 pip.push_back(new OpLhsG_dH<true>("G", "H", h_ptr));
1254 };
1255
1256 auto create_subdm = [&]() {
1257 auto level_ents_ptr = boost::make_shared<Range>();
1258 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
1259 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
1260
1261 DM subdm;
1262 CHKERR DMCreate(mField.get_comm(), &subdm);
1263 CHKERR DMSetType(subdm, "DMMOFEM");
1264 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_INIT");
1265 CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
1266 CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
1267 CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_TRUE);
1268
1269 for (auto f : {"H", "G"}) {
1270 CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
1271 CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
1272 }
1273 CHKERR DMSetUp(subdm);
1274
1275 if constexpr (debug) {
1276 if (mField.get_comm_size() == 1) {
1277 auto dm_ents = get_dofs_ents_all(SmartPetscObj<DM>(subdm, true));
1278 CHKERR save_range(mField.get_moab(), "sub_dm_init_ents_verts.h5m",
1279 dm_ents.subset_by_type(MBVERTEX));
1280 dm_ents = subtract(dm_ents, dm_ents.subset_by_type(MBVERTEX));
1281 CHKERR save_range(mField.get_moab(), "sub_dm_init_ents.h5m", dm_ents);
1282 }
1283 }
1284
1285 return SmartPetscObj<DM>(subdm);
1286 };
1287
1288 auto subdm = create_subdm();
1289
1290 pip_mng->getDomainRhsFE().reset();
1291 pip_mng->getDomainLhsFE().reset();
1292 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
1293 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
1294 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1295 pip_mng->getDomainRhsFE()->exeTestHook = exe_test;
1296
1297 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1298 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1299
1300 auto D = createDMVector(subdm);
1301 auto snes = pip_mng->createSNES(subdm);
1302 auto snes_ctx_ptr = getDMSnesCtx(subdm);
1303
1304 auto set_section_monitor = [&](auto solver) {
1306 CHKERR SNESMonitorSet(solver,
1307 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1308 void *))MoFEMSNESMonitorFields,
1309 (void *)(snes_ctx_ptr.get()), nullptr);
1310
1312 };
1313
1314 auto print_section_field = [&]() {
1316
1317 auto section =
1318 mField.getInterface<ISManager>()->sectionCreate("SUB_INIT");
1319 PetscInt num_fields;
1320 CHKERR PetscSectionGetNumFields(section, &num_fields);
1321 for (int f = 0; f < num_fields; ++f) {
1322 const char *field_name;
1323 CHKERR PetscSectionGetFieldName(section, f, &field_name);
1324 MOFEM_LOG("FS", Sev::inform)
1325 << "Field " << f << " " << std::string(field_name);
1326 }
1327
1329 };
1330
1331 CHKERR set_section_monitor(snes);
1332 CHKERR print_section_field();
1333
1334 for (auto f : {"U", "P", "H", "G", "L"}) {
1335 CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
1336 }
1337
1338 CHKERR SNESSetFromOptions(snes);
1339 CHKERR SNESSolve(snes, PETSC_NULLPTR, D);
1340
1341 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1342 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1343 CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
1344
1346 };
1347
1348 CHKERR reset_bits();
1349 CHKERR solve_init(
1350 [](FEMethod *fe_ptr) { return get_fe_bit(fe_ptr).test(0); });
1351 CHKERR refineMesh(refine_overlap);
1352
1353 for (auto f : {"U", "P", "H", "G", "L"}) {
1354 CHKERR mField.getInterface<FieldBlas>()->setField(0, f);
1355 }
1356 CHKERR solve_init([](FEMethod *fe_ptr) {
1357 return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
1358 });
1359
1360 CHKERR post_proc([](FEMethod *fe_ptr) {
1361 return get_fe_bit(fe_ptr).test(get_start_bit() + nb_levels - 1);
1362 });
1363
1364 pip_mng->getOpDomainRhsPipeline().clear();
1365 pip_mng->getOpDomainLhsPipeline().clear();
1366
1367 // Remove DOFs where boundary conditions are set
1368 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
1369 "U", 0, 0);
1370 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_X",
1371 "L", 0, 0);
1372 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
1373 "U", 1, 1);
1374 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "SYM_Y",
1375 "L", 1, 1);
1376 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "U",
1377 0, SPACE_DIM);
1378 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX", "L",
1379 0, 0);
1380 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "ZERO",
1381 "L", 0, 0);
1382
1384}
1385//! [Boundary condition]
1386
1387//! [Data projection]
1390
1391 // FIXME: Functionality in this method should be improved (projection should
1392 // be done field by field), generalized, and move to become core
1393 // functionality.
1394
1395 auto simple = mField.getInterface<Simple>();
1396 auto pip_mng = mField.getInterface<PipelineManager>();
1397 auto bit_mng = mField.getInterface<BitRefManager>();
1398 auto field_blas = mField.getInterface<FieldBlas>();
1399
1400 // Store all existing elements pipelines, replace them by data projection
1401 // pipelines, to put them back when projection is done.
1402 auto fe_domain_lhs = pip_mng->getDomainLhsFE();
1403 auto fe_domain_rhs = pip_mng->getDomainRhsFE();
1404 auto fe_bdy_lhs = pip_mng->getBoundaryLhsFE();
1405 auto fe_bdy_rhs = pip_mng->getBoundaryRhsFE();
1406
1407 pip_mng->getDomainLhsFE().reset();
1408 pip_mng->getDomainRhsFE().reset();
1409 pip_mng->getBoundaryLhsFE().reset();
1410 pip_mng->getBoundaryRhsFE().reset();
1411
1413
1414 // extract field data for domain parent element
1415 auto add_parent_field_domain = [&](auto fe, auto op, auto field, auto bit) {
1416 return setParentDofs(
1417 fe, field, op, bit,
1418
1419 [&]() {
1420 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1421 new DomainParentEle(mField));
1422 return fe_parent;
1423 },
1424
1425 QUIET, Sev::noisy);
1426 };
1427
1428 // extract field data for boundary parent element
1429 auto add_parent_field_bdy = [&](auto fe, auto op, auto field, auto bit) {
1430 return setParentDofs(
1431 fe, field, op, bit,
1432
1433 [&]() {
1434 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1436 return fe_parent;
1437 },
1438
1439 QUIET, Sev::noisy);
1440 };
1441
1442 // run this element on element with given bit, otherwise run other nested
1443 // element
1444 auto create_run_parent_op = [&](auto parent_fe_ptr, auto this_fe_ptr,
1445 auto fe_bit) {
1446 auto parent_mask = fe_bit;
1447 parent_mask.flip();
1448 return new OpRunParent(parent_fe_ptr, BitRefLevel().set(), parent_mask,
1449 this_fe_ptr, fe_bit, BitRefLevel().set(), QUIET,
1450 Sev::inform);
1451 };
1452
1453 // create hierarchy of nested elements
1454 auto get_parents_vel_fe_ptr = [&](auto this_fe_ptr, auto fe_bit) {
1455 std::vector<boost::shared_ptr<DomainParentEle>> parents_elems_ptr_vec;
1456 for (int l = 0; l < nb_levels; ++l)
1457 parents_elems_ptr_vec.emplace_back(
1458 boost::make_shared<DomainParentEle>(mField));
1459 for (auto l = 1; l < nb_levels; ++l) {
1460 parents_elems_ptr_vec[l - 1]->getOpPtrVector().push_back(
1461 create_run_parent_op(parents_elems_ptr_vec[l], this_fe_ptr, fe_bit));
1462 }
1463 return parents_elems_ptr_vec[0];
1464 };
1465
1466 // solve projection problem
1467 auto solve_projection = [&](auto exe_test) {
1469
1470 auto set_domain_rhs = [&](auto fe) {
1472 auto &pip = fe->getOpPtrVector();
1473
1474 auto u_ptr = boost::make_shared<MatrixDouble>();
1475 auto p_ptr = boost::make_shared<VectorDouble>();
1476 auto h_ptr = boost::make_shared<VectorDouble>();
1477 auto g_ptr = boost::make_shared<VectorDouble>();
1478
1479 auto eval_fe_ptr = boost::make_shared<DomainParentEle>(mField);
1480 {
1481 auto &pip = eval_fe_ptr->getOpPtrVector();
1484 eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
1485
1486 [&]() {
1487 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1488 new DomainParentEle(mField));
1489 return fe_parent;
1490 },
1491
1492 QUIET, Sev::noisy);
1493 // That can be done much smarter, by block, field by field. For
1494 // simplicity is like that.
1495 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "U",
1497 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1498 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "P",
1500 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1501 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "H",
1503 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1504 CHKERR add_parent_field_domain(eval_fe_ptr, UDO::OPROW, "G",
1506 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1507 }
1508 auto parent_eval_fe_ptr =
1509 get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
1510 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1512
1513 auto assemble_fe_ptr = boost::make_shared<DomainParentEle>(mField);
1514 {
1515 auto &pip = assemble_fe_ptr->getOpPtrVector();
1518 assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1519
1520 [&]() {
1521 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1522 new DomainParentEle(mField));
1523 return fe_parent;
1524 },
1525
1526 QUIET, Sev::noisy);
1527 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "U",
1529 pip.push_back(new OpDomainAssembleVector("U", u_ptr));
1530 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "P",
1532 pip.push_back(new OpDomainAssembleScalar("P", p_ptr));
1533 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "H",
1535 pip.push_back(new OpDomainAssembleScalar("H", h_ptr));
1536 CHKERR add_parent_field_domain(assemble_fe_ptr, UDO::OPROW, "G",
1538 pip.push_back(new OpDomainAssembleScalar("G", g_ptr));
1539 }
1540 auto parent_assemble_fe_ptr =
1541 get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
1542 pip.push_back(create_run_parent_op(
1543 parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
1544
1546 };
1547
1548 auto set_domain_lhs = [&](auto fe) {
1550
1551 auto &pip = fe->getOpPtrVector();
1552
1554
1556 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1557
1558 [&]() {
1559 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1560 new DomainParentEle(mField));
1561 return fe_parent;
1562 },
1563
1564 QUIET, Sev::noisy);
1565
1566 // That can be done much smarter, by block, field by field. For simplicity
1567 // is like that.
1568 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U",
1570 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U",
1572 pip.push_back(new OpDomainMassU("U", "U"));
1573 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P",
1575 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P",
1577 pip.push_back(new OpDomainMassP("P", "P"));
1578 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H",
1580 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H",
1582 pip.push_back(new OpDomainMassH("H", "H"));
1583 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G",
1585 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G",
1587 pip.push_back(new OpDomainMassG("G", "G"));
1588
1590 };
1591
1592 auto set_bdy_rhs = [&](auto fe) {
1594 auto &pip = fe->getOpPtrVector();
1595
1596 auto l_ptr = boost::make_shared<VectorDouble>();
1597
1598 auto eval_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
1599 {
1600 auto &pip = eval_fe_ptr->getOpPtrVector();
1602 eval_fe_ptr, "", UDO::OPSPACE, bit(get_skin_projection_bit()),
1603
1604 [&]() {
1605 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1607 return fe_parent;
1608 },
1609
1610 QUIET, Sev::noisy);
1611 // That can be done much smarter, by block, field by field. For
1612 // simplicity is like that.
1613 CHKERR add_parent_field_bdy(eval_fe_ptr, UDO::OPROW, "L",
1615 pip.push_back(new OpCalculateScalarFieldValues("L", l_ptr));
1616 }
1617 auto parent_eval_fe_ptr =
1618 get_parents_vel_fe_ptr(eval_fe_ptr, bit(get_projection_bit()));
1619 pip.push_back(create_run_parent_op(parent_eval_fe_ptr, eval_fe_ptr,
1621
1622 auto assemble_fe_ptr = boost::make_shared<BoundaryParentEle>(mField);
1623 {
1624 auto &pip = assemble_fe_ptr->getOpPtrVector();
1626 assemble_fe_ptr, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1627
1628 [&]() {
1629 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1631 return fe_parent;
1632 },
1633
1634 QUIET, Sev::noisy);
1635
1636 struct OpLSize : public BoundaryEleOp {
1637 OpLSize(boost::shared_ptr<VectorDouble> l_ptr)
1638 : BoundaryEleOp(NOSPACE, DomainEleOp::OPSPACE), lPtr(l_ptr) {}
1639 MoFEMErrorCode doWork(int, EntityType, EntData &) {
1641 if (lPtr->size() != getGaussPts().size2()) {
1642 lPtr->resize(getGaussPts().size2());
1643 lPtr->clear();
1644 }
1646 }
1647
1648 private:
1649 boost::shared_ptr<VectorDouble> lPtr;
1650 };
1651
1652 pip.push_back(new OpLSize(l_ptr));
1653
1654 CHKERR add_parent_field_bdy(assemble_fe_ptr, UDO::OPROW, "L",
1656 pip.push_back(new OpBoundaryAssembleScalar("L", l_ptr));
1657 }
1658 auto parent_assemble_fe_ptr =
1659 get_parents_vel_fe_ptr(assemble_fe_ptr, bit(get_current_bit()));
1660 pip.push_back(create_run_parent_op(
1661 parent_assemble_fe_ptr, assemble_fe_ptr, bit(get_current_bit())));
1662
1664 };
1665
1666 auto set_bdy_lhs = [&](auto fe) {
1668
1669 auto &pip = fe->getOpPtrVector();
1670
1672 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1673
1674 [&]() {
1675 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1677 return fe_parent;
1678 },
1679
1680 QUIET, Sev::noisy);
1681
1682 // That can be done much smarter, by block, field by field. For simplicity
1683 // is like that.
1684 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L",
1686 CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "L",
1688 pip.push_back(new OpBoundaryMassL("L", "L"));
1689
1691 };
1692
1693 auto create_subdm = [&]() -> SmartPetscObj<DM> {
1694 auto level_ents_ptr = boost::make_shared<Range>();
1695 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
1696 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
1697
1698 auto get_prj_ents = [&]() {
1699 Range prj_mesh;
1700 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_projection_bit()),
1701 BitRefLevel().set(),
1702 SPACE_DIM, prj_mesh);
1703 auto common_ents = intersect(prj_mesh, *level_ents_ptr);
1704 prj_mesh = subtract(unite(*level_ents_ptr, prj_mesh), common_ents)
1705 .subset_by_dimension(SPACE_DIM);
1706
1707 return prj_mesh;
1708 };
1709
1710 auto prj_ents = get_prj_ents();
1711
1712 if (get_global_size(prj_ents.size())) {
1713
1714 auto rebuild = [&]() {
1715 auto prb_mng = mField.getInterface<ProblemsManager>();
1717
1718 std::vector<std::string> fields{"U", "P", "H", "G", "L"};
1719 std::map<std::string, boost::shared_ptr<Range>> range_maps{
1720
1721 {"U", level_ents_ptr},
1722 {"P", level_ents_ptr},
1723 {"H", level_ents_ptr},
1724 {"G", level_ents_ptr},
1725 {"L", level_ents_ptr}
1726
1727 };
1728
1729 CHKERR prb_mng->buildSubProblem("SUB_SOLVER", fields, fields,
1730 simple->getProblemName(), PETSC_TRUE,
1731 &range_maps, &range_maps);
1732
1733 // partition problem
1734 CHKERR prb_mng->partitionFiniteElements("SUB_SOLVER", true, 0,
1736 // set ghost nodes
1737 CHKERR prb_mng->partitionGhostDofsOnDistributedMesh("SUB_SOLVER");
1738
1740 };
1741
1742 MOFEM_LOG("FS", Sev::verbose) << "Create projection problem";
1743
1744 CHKERR rebuild();
1745
1746 auto dm = simple->getDM();
1747 DM subdm;
1748 CHKERR DMCreate(mField.get_comm(), &subdm);
1749 CHKERR DMSetType(subdm, "DMMOFEM");
1750 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
1751 return SmartPetscObj<DM>(subdm);
1752 }
1753
1754 MOFEM_LOG("FS", Sev::inform) << "Nothing to project";
1755
1756 return SmartPetscObj<DM>();
1757 };
1758
1759 auto create_dummy_dm = [&]() {
1760 auto dummy_dm = createDM(mField.get_comm(), "DMMOFEM");
1762 simple->getProblemName().c_str(),
1764 "create dummy dm");
1765 return dummy_dm;
1766 };
1767
1768 auto subdm = create_subdm();
1769 if (subdm) {
1770
1771 pip_mng->getDomainRhsFE().reset();
1772 pip_mng->getDomainLhsFE().reset();
1773 pip_mng->getBoundaryRhsFE().reset();
1774 pip_mng->getBoundaryLhsFE().reset();
1775 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
1776 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
1777 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
1778 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
1779 pip_mng->getDomainLhsFE()->exeTestHook = exe_test;
1780 pip_mng->getDomainRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
1781 return get_fe_bit(fe_ptr).test(nb_levels - 1);
1782 };
1783 pip_mng->getBoundaryLhsFE()->exeTestHook = exe_test;
1784 pip_mng->getBoundaryRhsFE()->exeTestHook = [](FEMethod *fe_ptr) {
1785 return get_fe_bit(fe_ptr).test(nb_levels - 1);
1786 };
1787
1788 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
1789 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
1790 CHKERR set_bdy_rhs(pip_mng->getCastBoundaryRhsFE());
1791 CHKERR set_bdy_lhs(pip_mng->getCastBoundaryLhsFE());
1792
1793 auto D = createDMVector(subdm);
1794 auto F = vectorDuplicate(D);
1795
1796 auto ksp = pip_mng->createKSP(subdm);
1797 CHKERR KSPSetFromOptions(ksp);
1798 CHKERR KSPSetUp(ksp);
1799
1800 // get vector norm
1801 auto get_norm = [&](auto x) {
1802 double nrm;
1803 CHKERR VecNorm(x, NORM_2, &nrm);
1804 return nrm;
1805 };
1806
1807 /**
1808 * @brief Zero DOFs, used by FieldBlas
1809 */
1810 auto zero_dofs = [](boost::shared_ptr<FieldEntity> ent_ptr) {
1812 for (auto &v : ent_ptr->getEntFieldData()) {
1813 v = 0;
1814 }
1816 };
1817
1818 auto solve = [&](auto S) {
1820 CHKERR VecZeroEntries(S);
1821 CHKERR VecZeroEntries(F);
1822 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
1823 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
1824 CHKERR KSPSolve(ksp, F, S);
1825 CHKERR VecGhostUpdateBegin(S, INSERT_VALUES, SCATTER_FORWARD);
1826 CHKERR VecGhostUpdateEnd(S, INSERT_VALUES, SCATTER_FORWARD);
1827
1828
1829
1831 };
1832
1833 MOFEM_LOG("FS", Sev::inform) << "Solve projection";
1834 CHKERR solve(D);
1835
1836 auto glob_x = createDMVector(simple->getDM());
1837 auto sub_x = createDMVector(subdm);
1838 auto dummy_dm = create_dummy_dm();
1839
1840 /**
1841 * @brief get TSTheta data operators
1842 */
1843 auto apply_restrict = [&]() {
1844 auto get_is = [](auto v) {
1845 IS iy;
1846 auto create = [&]() {
1848 int n, ystart;
1849 CHKERR VecGetLocalSize(v, &n);
1850 CHKERR VecGetOwnershipRange(v, &ystart, NULL);
1851 CHKERR ISCreateStride(PETSC_COMM_SELF, n, ystart, 1, &iy);
1853 };
1854 CHK_THROW_MESSAGE(create(), "create is");
1855 return SmartPetscObj<IS>(iy);
1856 };
1857
1858 auto iy = get_is(glob_x);
1859 auto s = createVecScatter(glob_x, PETSC_NULLPTR, glob_x, iy);
1860
1862 DMSubDomainRestrict(simple->getDM(), s, PETSC_NULLPTR, dummy_dm),
1863 "restrict");
1864 Vec X0, Xdot;
1865 CHK_THROW_MESSAGE(DMGetNamedGlobalVector(dummy_dm, "TSTheta_X0", &X0),
1866 "get X0");
1868 DMGetNamedGlobalVector(dummy_dm, "TSTheta_Xdot", &Xdot),
1869 "get Xdot");
1870
1871 auto forward_ghost = [](auto g) {
1873 CHKERR VecGhostUpdateBegin(g, INSERT_VALUES, SCATTER_FORWARD);
1874 CHKERR VecGhostUpdateEnd(g, INSERT_VALUES, SCATTER_FORWARD);
1876 };
1877
1878 CHK_THROW_MESSAGE(forward_ghost(X0), "");
1879 CHK_THROW_MESSAGE(forward_ghost(Xdot), "");
1880
1881 if constexpr (debug) {
1882 MOFEM_LOG("FS", Sev::inform)
1883 << "Reverse restrict: X0 " << get_norm(X0) << " Xdot "
1884 << get_norm(Xdot);
1885 }
1886
1887 return std::vector<Vec>{X0, Xdot};
1888 };
1889
1890 auto ts_solver_vecs = apply_restrict();
1891
1892 if (ts_solver_vecs.size()) {
1893
1894 for (auto v : ts_solver_vecs) {
1895 MOFEM_LOG("FS", Sev::inform) << "Solve projection vector";
1896
1897 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1898 SCATTER_REVERSE);
1899 CHKERR solve(sub_x);
1900
1901 for (auto f : {"U", "P", "H", "G", "L"}) {
1902 MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1903 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1904 }
1905
1906 CHKERR DMoFEMMeshToLocalVector(subdm, sub_x, INSERT_VALUES,
1907 SCATTER_REVERSE);
1908 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), v, INSERT_VALUES,
1909 SCATTER_FORWARD);
1910
1911 MOFEM_LOG("FS", Sev::inform) << "Norm V " << get_norm(v);
1912 }
1913
1914 CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_X0",
1915 &ts_solver_vecs[0]);
1916 CHKERR DMRestoreNamedGlobalVector(dummy_dm, "TSTheta_Xdot",
1917 &ts_solver_vecs[1]);
1918 }
1919
1920 for (auto f : {"U", "P", "H", "G", "L"}) {
1921 MOFEM_LOG("WORLD", Sev::verbose) << "Zero field " << f;
1922 CHKERR field_blas->fieldLambdaOnEntities(zero_dofs, f);
1923 }
1924 CHKERR DMoFEMMeshToLocalVector(subdm, D, INSERT_VALUES, SCATTER_REVERSE);
1925 }
1926
1928 };
1929
1930 // postprocessing (only for debugging purposes)
1931 auto post_proc = [&](auto exe_test) {
1933 auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
1934 auto &pip = post_proc_fe->getOpPtrVector();
1935 post_proc_fe->exeTestHook = exe_test;
1936
1938
1940 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
1941
1942 [&]() {
1943 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
1944 new DomainParentEle(mField));
1946 fe_parent->getOpPtrVector(), {H1});
1947 return fe_parent;
1948 },
1949
1950 QUIET, Sev::noisy);
1951
1953
1954 auto u_ptr = boost::make_shared<MatrixDouble>();
1955 auto p_ptr = boost::make_shared<VectorDouble>();
1956 auto h_ptr = boost::make_shared<VectorDouble>();
1957 auto g_ptr = boost::make_shared<VectorDouble>();
1958
1959 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U",
1961 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1962 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P",
1964 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
1965 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H",
1967 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
1968 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G",
1970 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
1971
1972 post_proc_fe->getOpPtrVector().push_back(
1973
1974 new OpPPMap(post_proc_fe->getPostProcMesh(),
1975 post_proc_fe->getMapGaussPts(),
1976
1977 {{"P", p_ptr}, {"H", h_ptr}, {"G", g_ptr}},
1978
1979 {{"U", u_ptr}},
1980
1981 {},
1982
1983 {}
1984
1985 )
1986
1987 );
1988
1989 auto dm = simple->getDM();
1990 CHKERR DMoFEMLoopFiniteElements(dm, "dFE", post_proc_fe);
1991 CHKERR post_proc_fe->writeFile("out_projection.h5m");
1992
1994 };
1995
1996 CHKERR solve_projection([](FEMethod *fe_ptr) {
1997 return get_fe_bit(fe_ptr).test(get_current_bit());
1998 });
1999
2000 if constexpr (debug) {
2001 CHKERR post_proc([](FEMethod *fe_ptr) {
2002 return get_fe_bit(fe_ptr).test(get_current_bit());
2003 });
2004 }
2005
2006 fe_domain_lhs.swap(pip_mng->getDomainLhsFE());
2007 fe_domain_rhs.swap(pip_mng->getDomainRhsFE());
2008 fe_bdy_lhs.swap(pip_mng->getBoundaryLhsFE());
2009 fe_bdy_rhs.swap(pip_mng->getBoundaryRhsFE());
2010
2012}
2013//! [Data projection]
2014
2015//! [Push operators to pip]
2018 auto simple = mField.getInterface<Simple>();
2019
2021
2022 auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
2023 return setParentDofs(
2024 fe, field, op, bit(get_skin_parent_bit()),
2025
2026 [&]() {
2027 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2028 new DomainParentEle(mField));
2029 return fe_parent;
2030 },
2031
2032 QUIET, Sev::noisy);
2033 };
2034
2035 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2036 return setParentDofs(
2037 fe, field, op, bit(get_skin_parent_bit()),
2038
2039 [&]() {
2040 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2042 return fe_parent;
2043 },
2044
2045 QUIET, Sev::noisy);
2046 };
2047
2048 auto test_bit_child = [](FEMethod *fe_ptr) {
2049 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2050 get_start_bit() + nb_levels - 1);
2051 };
2052
2053 auto dot_u_ptr = boost::make_shared<MatrixDouble>();
2054 auto u_ptr = boost::make_shared<MatrixDouble>();
2055 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2056 auto dot_h_ptr = boost::make_shared<VectorDouble>();
2057 auto h_ptr = boost::make_shared<VectorDouble>();
2058 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2059 auto g_ptr = boost::make_shared<VectorDouble>();
2060 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2061 auto lambda_ptr = boost::make_shared<VectorDouble>();
2062 auto p_ptr = boost::make_shared<VectorDouble>();
2063 auto div_u_ptr = boost::make_shared<VectorDouble>();
2064
2065 // Push element from reference configuration to current configuration in 3d
2066 // space
2067 auto set_domain_general = [&](auto fe) {
2069 auto &pip = fe->getOpPtrVector();
2070
2072
2074 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2075
2076 [&]() {
2077 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2078 new DomainParentEle(mField));
2080 fe_parent->getOpPtrVector(), {H1});
2081 return fe_parent;
2082 },
2083
2084 QUIET, Sev::noisy);
2085
2086 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
2087 pip.push_back(
2089 pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
2091 "U", grad_u_ptr));
2092
2093 switch (coord_type) {
2094 case CARTESIAN:
2095 pip.push_back(
2097 "U", div_u_ptr));
2098 break;
2099 case CYLINDRICAL:
2100 pip.push_back(
2102 "U", div_u_ptr));
2103 break;
2104 default:
2105 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
2106 "Coordinate system not implemented");
2107 }
2108
2109 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
2110 pip.push_back(new OpCalculateScalarFieldValuesDot("H", dot_h_ptr));
2111 pip.push_back(new OpCalculateScalarFieldValues("H", h_ptr));
2112 pip.push_back(
2113 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
2114
2115 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
2116 pip.push_back(new OpCalculateScalarFieldValues("G", g_ptr));
2117 pip.push_back(
2118 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
2119
2120 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
2121 pip.push_back(new OpCalculateScalarFieldValues("P", p_ptr));
2123 };
2124
2125 auto set_domain_rhs = [&](auto fe) {
2127 auto &pip = fe->getOpPtrVector();
2128
2129 CHKERR set_domain_general(fe);
2130
2131 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
2132 pip.push_back(new OpRhsU("U", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr,
2133 grad_h_ptr, g_ptr, p_ptr));
2134
2135 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
2136 pip.push_back(new OpRhsH<false>("H", u_ptr, dot_h_ptr, h_ptr, grad_h_ptr,
2137 grad_g_ptr));
2138
2139 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
2140 pip.push_back(new OpRhsG<false>("G", h_ptr, grad_h_ptr, g_ptr));
2141
2142 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
2143 pip.push_back(new OpDomainAssembleScalar(
2144 "P", div_u_ptr, [](const double r, const double, const double) {
2145 return cylindrical(r);
2146 }));
2147 pip.push_back(new OpDomainAssembleScalar(
2148 "P", p_ptr, [](const double r, const double, const double) {
2149 return eps * cylindrical(r);
2150 }));
2152 };
2153
2154 auto set_domain_lhs = [&](auto fe) {
2156 auto &pip = fe->getOpPtrVector();
2157
2158 CHKERR set_domain_general(fe);
2159
2160 CHKERR add_parent_field_domain(fe, UDO::OPROW, "U");
2161 {
2162 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
2163 pip.push_back(new OpLhsU_dU("U", u_ptr, grad_u_ptr, h_ptr));
2164 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
2165 pip.push_back(
2166 new OpLhsU_dH("U", "H", dot_u_ptr, u_ptr, grad_u_ptr, h_ptr, g_ptr));
2167
2168 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
2169 pip.push_back(new OpLhsU_dG("U", "G", grad_h_ptr));
2170 }
2171
2172 CHKERR add_parent_field_domain(fe, UDO::OPROW, "H");
2173 {
2174 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
2175 pip.push_back(new OpLhsH_dU("H", "U", grad_h_ptr));
2176 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
2177 pip.push_back(new OpLhsH_dH<false>("H", u_ptr, h_ptr, grad_g_ptr));
2178 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
2179 pip.push_back(new OpLhsH_dG<false>("H", "G", h_ptr));
2180 }
2181
2182 CHKERR add_parent_field_domain(fe, UDO::OPROW, "G");
2183 {
2184 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "H");
2185 pip.push_back(new OpLhsG_dH<false>("G", "H", h_ptr));
2186 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "G");
2187 pip.push_back(new OpLhsG_dG("G"));
2188 }
2189
2190 CHKERR add_parent_field_domain(fe, UDO::OPROW, "P");
2191 {
2192 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "U");
2193
2194 switch (coord_type) {
2195 case CARTESIAN:
2196 pip.push_back(new OpMixScalarTimesDiv<CARTESIAN>(
2197 "P", "U",
2198 [](const double r, const double, const double) {
2199 return cylindrical(r);
2200 },
2201 true, false));
2202 break;
2203 case CYLINDRICAL:
2204 pip.push_back(new OpMixScalarTimesDiv<CYLINDRICAL>(
2205 "P", "U",
2206 [](const double r, const double, const double) {
2207 return cylindrical(r);
2208 },
2209 true, false));
2210 break;
2211 default:
2212 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
2213 "Coordinate system not implemented");
2214 }
2215
2216 CHKERR add_parent_field_domain(fe, UDO::OPCOL, "P");
2217 pip.push_back(new OpDomainMassP("P", "P", [](double r, double, double) {
2218 return eps * cylindrical(r);
2219 }));
2220 }
2221
2223 };
2224
2225 auto get_block_name = [](auto name) {
2226 return boost::format("%s(.*)") % "WETTING_ANGLE";
2227 };
2228
2229 auto get_blocks = [&](auto &&name) {
2230 return mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
2231 std::regex(name.str()));
2232 };
2233
2234 auto set_boundary_rhs = [&](auto fe) {
2236 auto &pip = fe->getOpPtrVector();
2237
2239 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2240
2241 [&]() {
2242 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2244 return fe_parent;
2245 },
2246
2247 QUIET, Sev::noisy);
2248
2249 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
2250 pip.push_back(new OpCalculateVectorFieldValues<U_FIELD_DIM>("U", u_ptr));
2251
2252 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
2253 pip.push_back(new OpCalculateScalarFieldValues("L", lambda_ptr));
2254 pip.push_back(new OpNormalConstrainRhs("L", u_ptr));
2255
2257 pip, mField, "L", {}, "INFLUX",
2258 [](double r, double, double) { return cylindrical(r); }, Sev::inform);
2259
2260 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "U");
2261 pip.push_back(new OpNormalForceRhs("U", lambda_ptr));
2262
2263 auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
2264 if (wetting_block.size()) {
2265 // push operators to the side element which is called from op_bdy_side
2266 auto op_bdy_side =
2267 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
2268 op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
2269
2271 op_bdy_side->getOpPtrVector(), {H1});
2272
2274 op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
2276
2277 [&]() {
2278 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2279 new DomainParentEle(mField));
2281 fe_parent->getOpPtrVector(), {H1});
2282 return fe_parent;
2283 },
2284
2285 QUIET, Sev::noisy);
2286
2287 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
2288 "H");
2289 op_bdy_side->getOpPtrVector().push_back(
2290 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
2291 // push bdy side op
2292 pip.push_back(op_bdy_side);
2293
2294 // push operators for rhs wetting angle
2295 for (auto &b : wetting_block) {
2296 Range force_edges;
2297 std::vector<double> attr_vec;
2298 CHKERR b->getMeshsetIdEntitiesByDimension(
2299 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
2300 b->getAttributes(attr_vec);
2301 if (attr_vec.size() != 1)
2302 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
2303 "Should be one attribute");
2304 MOFEM_LOG("FS", Sev::inform) << "Wetting angle: " << attr_vec.front();
2305 // need to find the attributes and pass to operator
2306 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
2307 pip.push_back(new OpWettingAngleRhs(
2308 "G", grad_h_ptr, boost::make_shared<Range>(force_edges),
2309 attr_vec.front()));
2310 }
2311 }
2312
2314 };
2315
2316 auto set_boundary_lhs = [&](auto fe) {
2318 auto &pip = fe->getOpPtrVector();
2319
2321 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2322
2323 [&]() {
2324 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2326 return fe_parent;
2327 },
2328
2329 QUIET, Sev::noisy);
2330
2331 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "L");
2332 CHKERR add_parent_field_bdy(fe, UDO::OPCOL, "U");
2333 pip.push_back(new OpNormalConstrainLhs("L", "U"));
2334
2335 auto wetting_block = get_blocks(get_block_name("WETTING_ANGLE"));
2336 if (wetting_block.size()) {
2337 auto col_ind_ptr = boost::make_shared<std::vector<VectorInt>>();
2338 auto col_diff_base_ptr = boost::make_shared<std::vector<MatrixDouble>>();
2339
2340 // push operators to the side element which is called from op_bdy_side
2341 auto op_bdy_side =
2342 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
2343 op_bdy_side->getSideFEPtr()->exeTestHook = test_bit_child;
2344
2346 op_bdy_side->getOpPtrVector(), {H1});
2347
2349 op_bdy_side->getSideFEPtr(), "", UDO::OPSPACE,
2351
2352 [&]() {
2353 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2354 new DomainParentEle(mField));
2356 fe_parent->getOpPtrVector(), {H1});
2357 return fe_parent;
2358 },
2359
2360 QUIET, Sev::noisy);
2361
2362 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPROW,
2363 "H");
2364 op_bdy_side->getOpPtrVector().push_back(
2365 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
2366 CHKERR add_parent_field_domain(op_bdy_side->getSideFEPtr(), UDO::OPCOL,
2367 "H");
2368 op_bdy_side->getOpPtrVector().push_back(
2369 new OpLoopSideGetDataForSideEle("H", col_ind_ptr, col_diff_base_ptr));
2370
2371 // push bdy side op
2372 pip.push_back(op_bdy_side);
2373
2374 // push operators for lhs wetting angle
2375 for (auto &b : wetting_block) {
2376 Range force_edges;
2377 std::vector<double> attr_vec;
2378 CHKERR b->getMeshsetIdEntitiesByDimension(
2379 mField.get_moab(), SPACE_DIM - 1, force_edges, true);
2380 b->getAttributes(attr_vec);
2381 if (attr_vec.size() != 1)
2382 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
2383 "Should be one attribute");
2384 MOFEM_LOG("FS", Sev::inform)
2385 << "wetting angle edges size " << force_edges.size();
2386
2387 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "G");
2388 pip.push_back(new OpWettingAngleLhs(
2389 "G", grad_h_ptr, col_ind_ptr, col_diff_base_ptr,
2390 boost::make_shared<Range>(force_edges), attr_vec.front()));
2391 }
2392 }
2393
2395 };
2396
2397 auto *pip_mng = mField.getInterface<PipelineManager>();
2398
2399 CHKERR set_domain_rhs(pip_mng->getCastDomainRhsFE());
2400 CHKERR set_domain_lhs(pip_mng->getCastDomainLhsFE());
2401 CHKERR set_boundary_rhs(pip_mng->getCastBoundaryRhsFE());
2402 CHKERR set_boundary_lhs(pip_mng->getCastBoundaryLhsFE());
2403
2404 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule);
2405 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule);
2406 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule);
2407 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule);
2408
2409 pip_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
2410 pip_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
2411 pip_mng->getBoundaryLhsFE()->exeTestHook = test_bit_child;
2412 pip_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
2413
2415}
2416//! [Push operators to pip]
2417
2418/**
2419 * @brief Monitor solution
2420 *
2421 * This function is called by TS solver at the end of each step. It is used
2422 */
2423struct Monitor : public FEMethod {
2425 SmartPetscObj<DM> dm, boost::shared_ptr<moab::Core> post_proc_mesh,
2426 boost::shared_ptr<PostProcEleDomainCont> post_proc,
2427 boost::shared_ptr<PostProcEleBdyCont> post_proc_edge,
2428 std::pair<boost::shared_ptr<BoundaryEle>, boost::shared_ptr<VectorDouble>>
2429 p)
2430 : dM(dm), postProcMesh(post_proc_mesh), postProc(post_proc),
2431 postProcEdge(post_proc_edge), liftFE(p.first), liftVec(p.second) {}
2434
2435 MOFEM_LOG("FS", Sev::verbose) << "Monitor";
2436 constexpr int save_every_nth_step = 1;
2437 if (ts_step % save_every_nth_step == 0) {
2438 MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc";
2439 MoFEM::Interface *m_field_ptr;
2440 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
2441 auto post_proc_begin =
2442 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
2443 postProcMesh);
2444 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
2445 *m_field_ptr, postProcMesh);
2446 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
2448 this->getCacheWeakPtr());
2450 this->getCacheWeakPtr());
2451 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
2452 CHKERR post_proc_end->writeFile(
2453 "out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
2454 MOFEM_LOG("FS", Sev::verbose) << "Mesh pre proc done";
2455 }
2456
2457 liftVec->resize(SPACE_DIM, false);
2458 liftVec->clear();
2460 MPI_Allreduce(MPI_IN_PLACE, &(*liftVec)[0], SPACE_DIM, MPI_DOUBLE, MPI_SUM,
2461 MPI_COMM_WORLD);
2462 MOFEM_LOG("FS", Sev::inform)
2463 << "Step " << ts_step << " time " << ts_t
2464 << " lift vec x: " << (*liftVec)[0] << " y: " << (*liftVec)[1];
2465
2467 }
2468
2469private:
2471 boost::shared_ptr<moab::Core> postProcMesh;
2472 boost::shared_ptr<PostProcEleDomainCont> postProc;
2473 boost::shared_ptr<PostProcEleBdyCont> postProcEdge;
2474 boost::shared_ptr<BoundaryEle> liftFE;
2475 boost::shared_ptr<VectorDouble> liftVec;
2476};
2477
2478//! [Solve]
2481
2483
2484 auto simple = mField.getInterface<Simple>();
2485 auto pip_mng = mField.getInterface<PipelineManager>();
2486
2487 auto create_solver_dm = [&](auto dm) -> SmartPetscObj<DM> {
2488 DM subdm;
2489
2490 auto setup_subdm = [&](auto dm) {
2492 auto simple = mField.getInterface<Simple>();
2493 auto bit_mng = mField.getInterface<BitRefManager>();
2494 auto dm = simple->getDM();
2495 auto level_ents_ptr = boost::make_shared<Range>();
2496 CHKERR bit_mng->getEntitiesByRefLevel(
2497 bit(get_current_bit()), BitRefLevel().set(), *level_ents_ptr);
2498 CHKERR DMCreate(mField.get_comm(), &subdm);
2499 CHKERR DMSetType(subdm, "DMMOFEM");
2500 CHKERR DMMoFEMCreateSubDM(subdm, dm, "SUB_SOLVER");
2501 CHKERR DMMoFEMAddElement(subdm, simple->getDomainFEName());
2502 CHKERR DMMoFEMAddElement(subdm, simple->getBoundaryFEName());
2503 CHKERR DMMoFEMSetSquareProblem(subdm, PETSC_TRUE);
2504 CHKERR DMMoFEMSetDestroyProblem(subdm, PETSC_FALSE);
2505 for (auto f : {"U", "P", "H", "G", "L"}) {
2506 CHKERR DMMoFEMAddSubFieldRow(subdm, f, level_ents_ptr);
2507 CHKERR DMMoFEMAddSubFieldCol(subdm, f, level_ents_ptr);
2508 }
2509 CHKERR DMSetUp(subdm);
2511 };
2512
2513 CHK_THROW_MESSAGE(setup_subdm(dm), "create subdm");
2514
2515 return SmartPetscObj<DM>(subdm);
2516 };
2517
2518 auto get_fe_post_proc = [&](auto post_proc_mesh) {
2519 auto add_parent_field_domain = [&](auto fe, auto op, auto field) {
2520 return setParentDofs(
2521 fe, field, op, bit(get_skin_parent_bit()),
2522
2523 [&]() {
2524 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2525 new DomainParentEle(mField));
2526 return fe_parent;
2527 },
2528
2529 QUIET, Sev::noisy);
2530 };
2531
2532 auto post_proc_fe =
2533 boost::make_shared<PostProcEleDomainCont>(mField, post_proc_mesh);
2534 post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2535 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2536 get_start_bit() + nb_levels - 1);
2537 };
2538
2539 auto u_ptr = boost::make_shared<MatrixDouble>();
2540 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
2541 auto h_ptr = boost::make_shared<VectorDouble>();
2542 auto grad_h_ptr = boost::make_shared<MatrixDouble>();
2543 auto p_ptr = boost::make_shared<VectorDouble>();
2544 auto g_ptr = boost::make_shared<VectorDouble>();
2545 auto grad_g_ptr = boost::make_shared<MatrixDouble>();
2546
2548 post_proc_fe->getOpPtrVector(), {H1});
2549
2551 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2552
2553 [&]() {
2554 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2555 new DomainParentEle(mField));
2557 fe_parent->getOpPtrVector(), {H1});
2558 return fe_parent;
2559 },
2560
2561 QUIET, Sev::noisy);
2562
2563 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "U");
2564 post_proc_fe->getOpPtrVector().push_back(
2566 post_proc_fe->getOpPtrVector().push_back(
2568 grad_u_ptr));
2569
2570 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "H");
2571 post_proc_fe->getOpPtrVector().push_back(
2572 new OpCalculateScalarFieldValues("H", h_ptr));
2573 post_proc_fe->getOpPtrVector().push_back(
2574 new OpCalculateScalarFieldGradient<SPACE_DIM>("H", grad_h_ptr));
2575
2576 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "P");
2577 post_proc_fe->getOpPtrVector().push_back(
2578 new OpCalculateScalarFieldValues("P", p_ptr));
2579
2580 CHKERR add_parent_field_domain(post_proc_fe, UDO::OPROW, "G");
2581 post_proc_fe->getOpPtrVector().push_back(
2582 new OpCalculateScalarFieldValues("G", g_ptr));
2583 post_proc_fe->getOpPtrVector().push_back(
2584 new OpCalculateScalarFieldGradient<SPACE_DIM>("G", grad_g_ptr));
2585
2587
2588 post_proc_fe->getOpPtrVector().push_back(
2589
2590 new OpPPMap(
2591 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
2592
2593 {{"H", h_ptr}, {"P", p_ptr}, {"G", g_ptr}},
2594
2595 {{"U", u_ptr}, {"H_GRAD", grad_h_ptr}, {"G_GRAD", grad_g_ptr}},
2596
2597 {{"GRAD_U", grad_u_ptr}},
2598
2599 {}
2600
2601 )
2602
2603 );
2604
2605 return post_proc_fe;
2606 };
2607
2608 auto get_bdy_post_proc_fe = [&](auto post_proc_mesh) {
2609 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2610 return setParentDofs(
2611 fe, field, op, bit(get_skin_parent_bit()),
2612
2613 [&]() {
2614 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2615 new BoundaryParentEle(mField));
2616 return fe_parent;
2617 },
2618
2619 QUIET, Sev::noisy);
2620 };
2621
2622 auto post_proc_fe =
2623 boost::make_shared<PostProcEleBdyCont>(mField, post_proc_mesh);
2624 post_proc_fe->exeTestHook = [](FEMethod *fe_ptr) {
2625 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2626 get_start_bit() + nb_levels - 1);
2627 };
2628
2629 CHKERR setParentDofs(
2630 post_proc_fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2631
2632 [&]() {
2633 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2634 new BoundaryParentEle(mField));
2635 return fe_parent;
2636 },
2637
2638 QUIET, Sev::noisy);
2639
2640 struct OpGetNormal : public BoundaryEleOp {
2641
2642 OpGetNormal(boost::shared_ptr<VectorDouble> l_ptr,
2643 boost::shared_ptr<MatrixDouble> n_ptr)
2644 : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), ptrL(l_ptr),
2645 ptrNormal(n_ptr) {}
2646
2647 MoFEMErrorCode doWork(int side, EntityType type,
2650 auto t_l = getFTensor0FromVec(*ptrL);
2651 auto t_n_fe = getFTensor1NormalsAtGaussPts();
2652 ptrNormal->resize(SPACE_DIM, getGaussPts().size2(), false);
2653 auto t_n = getFTensor1FromMat<SPACE_DIM>(*ptrNormal);
2654 for (auto gg = 0; gg != getGaussPts().size2(); ++gg) {
2655 t_n(i) = t_n_fe(i) * t_l / std::sqrt(t_n_fe(i) * t_n_fe(i));
2656 ++t_n_fe;
2657 ++t_l;
2658 ++t_n;
2659 }
2661 };
2662
2663 protected:
2664 boost::shared_ptr<VectorDouble> ptrL;
2665 boost::shared_ptr<MatrixDouble> ptrNormal;
2666 };
2667
2668 auto u_ptr = boost::make_shared<MatrixDouble>();
2669 auto p_ptr = boost::make_shared<VectorDouble>();
2670 auto lambda_ptr = boost::make_shared<VectorDouble>();
2671 auto normal_l_ptr = boost::make_shared<MatrixDouble>();
2672
2673 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "U");
2674 post_proc_fe->getOpPtrVector().push_back(
2676
2677 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "L");
2678 post_proc_fe->getOpPtrVector().push_back(
2679 new OpCalculateScalarFieldValues("L", lambda_ptr));
2680 post_proc_fe->getOpPtrVector().push_back(
2681 new OpGetNormal(lambda_ptr, normal_l_ptr));
2682
2683 CHKERR add_parent_field_bdy(post_proc_fe, UDO::OPROW, "P");
2684 post_proc_fe->getOpPtrVector().push_back(
2685 new OpCalculateScalarFieldValues("P", p_ptr));
2686
2687 auto op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
2688 op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
2689 EntityType type,
2693 };
2694
2696
2697 post_proc_fe->getOpPtrVector().push_back(
2698
2699 new OpPPMap(post_proc_fe->getPostProcMesh(),
2700 post_proc_fe->getMapGaussPts(),
2701
2702 OpPPMap::DataMapVec{{"P", p_ptr}},
2703
2704 OpPPMap::DataMapMat{{"U", u_ptr}, {"L", normal_l_ptr}},
2705
2707
2709
2710 )
2711
2712 );
2713
2714 return post_proc_fe;
2715 };
2716
2717 auto get_lift_fe = [&]() {
2718 auto add_parent_field_bdy = [&](auto fe, auto op, auto field) {
2719 return setParentDofs(
2720 fe, field, op, bit(get_skin_parent_bit()),
2721
2722 [&]() {
2723 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2724 new BoundaryParentEle(mField));
2725 return fe_parent;
2726 },
2727
2728 QUIET, Sev::noisy);
2729 };
2730
2731 auto fe = boost::make_shared<BoundaryEle>(mField);
2732 fe->exeTestHook = [](FEMethod *fe_ptr) {
2733 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
2734 get_start_bit() + nb_levels - 1);
2735 };
2736
2737 auto lift_ptr = boost::make_shared<VectorDouble>();
2738 auto p_ptr = boost::make_shared<VectorDouble>();
2739 auto ents_ptr = boost::make_shared<Range>();
2740
2741 CHKERR setParentDofs(
2742 fe, "", UDO::OPSPACE, bit(get_skin_parent_bit()),
2743
2744 [&]() {
2745 boost::shared_ptr<ForcesAndSourcesCore> fe_parent(
2746 new BoundaryParentEle(mField));
2747 return fe_parent;
2748 },
2749
2750 QUIET, Sev::noisy);
2751
2752 std::vector<const CubitMeshSets *> vec_ptr;
2753 CHKERR mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
2754 std::regex("LIFT"), vec_ptr);
2755 for (auto m_ptr : vec_ptr) {
2756 auto meshset = m_ptr->getMeshset();
2757 Range ents;
2758 CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
2759 ents, true);
2760 ents_ptr->merge(ents);
2761 }
2762
2763 MOFEM_LOG("FS", Sev::noisy) << "Lift ents " << (*ents_ptr);
2764
2765 CHKERR add_parent_field_bdy(fe, UDO::OPROW, "P");
2766 fe->getOpPtrVector().push_back(
2767 new OpCalculateScalarFieldValues("L", p_ptr));
2768 fe->getOpPtrVector().push_back(
2769 new OpCalculateLift("L", p_ptr, lift_ptr, ents_ptr));
2770
2771 return std::make_pair(fe, lift_ptr);
2772 };
2773
2774 auto set_post_proc_monitor = [&](auto ts) {
2776 DM dm;
2777 CHKERR TSGetDM(ts, &dm);
2778 boost::shared_ptr<FEMethod> null_fe;
2779 auto post_proc_mesh = boost::make_shared<moab::Core>();
2780 auto monitor_ptr = boost::make_shared<Monitor>(
2781 SmartPetscObj<DM>(dm, true), post_proc_mesh,
2782 get_fe_post_proc(post_proc_mesh), get_bdy_post_proc_fe(post_proc_mesh),
2783 get_lift_fe());
2784 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
2785 null_fe, monitor_ptr);
2787 };
2788
2789 auto dm = simple->getDM();
2790 auto ts = createTS(mField.get_comm());
2791 CHKERR TSSetDM(ts, dm);
2792
2793 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
2794 tsPrePostProc = ts_pre_post_proc;
2795
2796 if (auto ptr = tsPrePostProc.lock()) {
2797
2798 ptr->fsRawPtr = this;
2799 ptr->solverSubDM = create_solver_dm(simple->getDM());
2800 ptr->globSol = createDMVector(dm);
2801 CHKERR DMoFEMMeshToLocalVector(dm, ptr->globSol, INSERT_VALUES,
2802 SCATTER_FORWARD);
2803 CHKERR VecAssemblyBegin(ptr->globSol);
2804 CHKERR VecAssemblyEnd(ptr->globSol);
2805
2806 auto sub_ts = pip_mng->createTSIM(ptr->solverSubDM);
2807
2808 CHKERR set_post_proc_monitor(sub_ts);
2809
2810 // Add monitor to time solver
2811 CHKERR TSSetFromOptions(ts);
2812 CHKERR ptr->tsSetUp(ts);
2813 CHKERR TSSetUp(ts);
2814
2815 auto print_fields_in_section = [&]() {
2817 auto section = mField.getInterface<ISManager>()->sectionCreate(
2818 simple->getProblemName());
2819 PetscInt num_fields;
2820 CHKERR PetscSectionGetNumFields(section, &num_fields);
2821 for (int f = 0; f < num_fields; ++f) {
2822 const char *field_name;
2823 CHKERR PetscSectionGetFieldName(section, f, &field_name);
2824 MOFEM_LOG("FS", Sev::inform)
2825 << "Field " << f << " " << std::string(field_name);
2826 }
2828 };
2829
2830 CHKERR print_fields_in_section();
2831
2832 CHKERR TSSolve(ts, ptr->globSol);
2833 }
2834
2836}
2837
2838/**
2839 * @brief Main function for free surface simulation
2840 * @param argc Number of command line arguments
2841 * @param argv Array of command line argument strings
2842 * @return Exit code (0 for success)
2843 *
2844 * Main driver function that:
2845 * 1. Initializes MoFEM/PETSc and MOAB data structures
2846 * 2. Sets up logging channels for debugging output
2847 * 3. Registers MoFEM discrete manager with PETSc
2848 * 4. Creates mesh database (MOAB) and finite element database (MoFEM)
2849 * 5. Runs the complete free surface simulation
2850 * 6. Handles cleanup and finalization
2851 *
2852 * Command line options are read from param_file.petsc and command line.
2853 * Python initialization is optional (controlled by PYTHON_INIT_SURFACE).
2854 */
2855int main(int argc, char *argv[]) {
2856
2857#ifdef PYTHON_INIT_SURFACE
2858 Py_Initialize();
2859#endif
2860
2861 // Initialisation of MoFEM/PETSc and MOAB data structures
2862 const char param_file[] = "param_file.petsc";
2863 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
2864
2865 // Add logging channel for example
2866 auto core_log = logging::core::get();
2867 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FS"));
2868 LogManager::setLog("FS");
2869 MOFEM_LOG_TAG("FS", "free surface");
2870
2871 try {
2872
2873 //! [Register MoFEM discrete manager in PETSc]
2874 DMType dm_name = "DMMOFEM";
2875 CHKERR DMRegister_MoFEM(dm_name);
2876 //! [Register MoFEM discrete manager in PETSc
2877
2878 //! [Create MoAB]
2879 moab::Core mb_instance; ///< mesh database
2880 moab::Interface &moab = mb_instance; ///< mesh database interface
2881 //! [Create MoAB]
2882
2883 //! [Create MoFEM]
2884 MoFEM::Core core(moab); ///< finite element database
2885 MoFEM::Interface &m_field = core; ///< finite element database interface
2886 //! [Create MoFEM]
2887
2888 //! [FreeSurface]
2889 FreeSurface ex(m_field);
2890 CHKERR ex.runProblem();
2891 //! [FreeSurface]
2892 }
2894
2896
2897#ifdef PYTHON_INIT_SURFACE
2898 if (Py_FinalizeEx() < 0) {
2899 exit(120);
2900 }
2901#endif
2902}
2903
2904std::vector<Range>
2906
2907 auto &moab = mField.get_moab();
2908 auto bit_mng = mField.getInterface<BitRefManager>();
2909 auto comm_mng = mField.getInterface<CommInterface>();
2910
2911 Range vertices;
2912 CHK_THROW_MESSAGE(bit_mng->getEntitiesByTypeAndRefLevel(
2913 bit(0), BitRefLevel().set(), MBVERTEX, vertices),
2914 "can not get vertices on bit 0");
2915
2916 auto &dofs_mi = mField.get_dofs()->get<Unique_mi_tag>();
2917 auto field_bit_number = mField.get_field_bit_number("H");
2918
2919 Range plus_range, minus_range;
2920 std::vector<EntityHandle> plus, minus;
2921
2922 // get vertices on level 0 on plus and minus side
2923 for (auto p = vertices.pair_begin(); p != vertices.pair_end(); ++p) {
2924
2925 const auto f = p->first;
2926 const auto s = p->second;
2927
2928 // Lowest Dof UId for given field (field bit number) on entity f
2929 const auto lo_uid = DofEntity::getLoFieldEntityUId(field_bit_number, f);
2930 const auto hi_uid = DofEntity::getHiFieldEntityUId(field_bit_number, s);
2931 auto it = dofs_mi.lower_bound(lo_uid);
2932 const auto hi_it = dofs_mi.upper_bound(hi_uid);
2933
2934 plus.clear();
2935 minus.clear();
2936 plus.reserve(std::distance(it, hi_it));
2937 minus.reserve(std::distance(it, hi_it));
2938
2939 for (; it != hi_it; ++it) {
2940 const auto v = (*it)->getFieldData();
2941 if (v > 0)
2942 plus.push_back((*it)->getEnt());
2943 else
2944 minus.push_back((*it)->getEnt());
2945 }
2946
2947 plus_range.insert_list(plus.begin(), plus.end());
2948 minus_range.insert_list(minus.begin(), minus.end());
2949 }
2950
2951 MOFEM_LOG_CHANNEL("SYNC");
2952 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2953 << "Plus range " << plus_range << endl;
2954 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "FS")
2955 << "Minus range " << minus_range << endl;
2957
2958 auto get_elems = [&](auto &ents, auto bit, auto mask) {
2959 Range adj;
2960 CHK_MOAB_THROW(moab.get_adjacencies(ents, SPACE_DIM, false, adj,
2961 moab::Interface::UNION),
2962 "can not get adjacencies");
2963 CHK_THROW_MESSAGE(bit_mng->filterEntitiesByRefLevel(bit, mask, adj),
2964 "can not filter elements with bit 0");
2965 return adj;
2966 };
2967
2968 CHKERR comm_mng->synchroniseEntities(plus_range);
2969 CHKERR comm_mng->synchroniseEntities(minus_range);
2970
2971 std::vector<Range> ele_plus(nb_levels), ele_minus(nb_levels);
2972 ele_plus[0] = get_elems(plus_range, bit(0), BitRefLevel().set());
2973 ele_minus[0] = get_elems(minus_range, bit(0), BitRefLevel().set());
2974 auto common = intersect(ele_plus[0], ele_minus[0]);
2975 ele_plus[0] = subtract(ele_plus[0], common);
2976 ele_minus[0] = subtract(ele_minus[0], common);
2977
2978 auto get_children = [&](auto &p, auto &c) {
2980 CHKERR bit_mng->updateRangeByChildren(p, c);
2981 c = c.subset_by_dimension(SPACE_DIM);
2983 };
2984
2985 for (auto l = 1; l != nb_levels; ++l) {
2986 CHK_THROW_MESSAGE(get_children(ele_plus[l - 1], ele_plus[l]),
2987 "get children");
2988 CHK_THROW_MESSAGE(get_children(ele_minus[l - 1], ele_minus[l]),
2989 "get children");
2990 }
2991
2992 auto get_level = [&](auto &p, auto &m, auto z, auto bit, auto mask) {
2993 Range l;
2995 bit_mng->getEntitiesByDimAndRefLevel(bit, mask, SPACE_DIM, l),
2996 "can not get vertices on bit");
2997 l = subtract(l, p);
2998 l = subtract(l, m);
2999 for (auto f = 0; f != z; ++f) {
3000 Range conn;
3001 CHK_MOAB_THROW(moab.get_connectivity(l, conn, true), "");
3002 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(conn);
3003 l = get_elems(conn, bit, mask);
3004 }
3005 return l;
3006 };
3007
3008 std::vector<Range> vec_levels(nb_levels);
3009 for (auto l = nb_levels - 1; l >= 0; --l) {
3010 vec_levels[l] = get_level(ele_plus[l], ele_minus[l], 2 * overlap, bit(l),
3011 BitRefLevel().set());
3012 }
3013
3014 if constexpr (debug) {
3015 if (mField.get_comm_size() == 1) {
3016 for (auto l = 0; l != nb_levels; ++l) {
3017 std::string name = (boost::format("out_r%d.h5m") % l).str();
3018 CHK_THROW_MESSAGE(save_range(mField.get_moab(), name, vec_levels[l]),
3019 "save mesh");
3020 }
3021 }
3022 }
3023
3024 return vec_levels;
3025}
3026
3029
3030 auto bit_mng = mField.getInterface<BitRefManager>();
3031
3032 BitRefLevel start_mask;
3033 for (auto s = 0; s != get_start_bit(); ++s)
3034 start_mask[s] = true;
3035
3036 // store prev_level
3037 Range prev_level;
3038 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_current_bit()),
3039 BitRefLevel().set(), prev_level);
3040 Range prev_level_skin;
3041 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_parent_bit()),
3042 BitRefLevel().set(), prev_level_skin);
3043 // reset bit ref levels
3044 CHKERR bit_mng->lambdaBitRefLevel(
3045 [&](EntityHandle ent, BitRefLevel &bit) { bit &= start_mask; });
3046 CHKERR bit_mng->setNthBitRefLevel(prev_level, get_projection_bit(), true);
3047 CHKERR bit_mng->setNthBitRefLevel(prev_level_skin, get_skin_projection_bit(),
3048 true);
3049
3050 // set refinement levels
3051 auto set_levels = [&](auto &&
3052 vec_levels /*entities are refined on each level*/) {
3054
3055 // start with zero level, which is the coarsest mesh
3056 Range level0;
3057 CHKERR bit_mng->getEntitiesByRefLevel(bit(0), BitRefLevel().set(), level0);
3058 CHKERR bit_mng->setNthBitRefLevel(level0, get_start_bit(), true);
3059
3060 // get lower dimension entities
3061 auto get_adj = [&](auto ents) {
3062 Range conn;
3063 CHK_MOAB_THROW(mField.get_moab().get_connectivity(ents, conn, true),
3064 "get conn");
3065 for (auto d = 1; d != SPACE_DIM; ++d) {
3066 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
3067 ents.subset_by_dimension(SPACE_DIM), d, false, conn,
3068 moab::Interface::UNION),
3069 "get adj");
3070 }
3071 ents.merge(conn);
3072 return ents;
3073 };
3074
3075 // set bit levels
3076 for (auto l = 1; l != nb_levels; ++l) {
3077 Range level_prev;
3078 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(get_start_bit() + l - 1),
3079 BitRefLevel().set(),
3080 SPACE_DIM, level_prev);
3081 Range parents;
3082 CHKERR bit_mng->updateRangeByParent(vec_levels[l], parents);
3083 // subtract entities from previous level, which are refined, so should be
3084 // not there
3085 level_prev = subtract(level_prev, parents);
3086 // and instead add their children
3087 level_prev.merge(vec_levels[l]);
3088 // set bit to each level
3089 CHKERR bit_mng->setNthBitRefLevel(level_prev, get_start_bit() + l, true);
3090 }
3091
3092 // set bit levels to lower dimension entities
3093 for (auto l = 1; l != nb_levels; ++l) {
3094 Range level;
3095 CHKERR bit_mng->getEntitiesByDimAndRefLevel(
3096 bit(get_start_bit() + l), BitRefLevel().set(), SPACE_DIM, level);
3097 level = get_adj(level);
3098 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(level);
3099 CHKERR bit_mng->setNthBitRefLevel(level, get_start_bit() + l, true);
3100 }
3101
3103 };
3104
3105 // resolve skin between refined levels
3106 auto set_skins = [&]() {
3108
3109 moab::Skinner skinner(&mField.get_moab());
3110 ParallelComm *pcomm =
3111 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
3112
3113 // get skin of bit level
3114 auto get_bit_skin = [&](BitRefLevel bit, BitRefLevel mask) {
3115 Range bit_ents;
3118 bit, mask, SPACE_DIM, bit_ents),
3119 "can't get bit level");
3120 Range bit_skin;
3121 CHK_MOAB_THROW(skinner.find_skin(0, bit_ents, false, bit_skin),
3122 "can't get skin");
3123 return bit_skin;
3124 };
3125
3126 auto get_level_skin = [&]() {
3127 Range skin;
3128 BitRefLevel bit_prev;
3129 for (auto l = 1; l != nb_levels; ++l) {
3130 auto skin_level_mesh = get_bit_skin(bit(l), BitRefLevel().set());
3131 // filter (remove) all entities which are on partition borders
3132 CHKERR pcomm->filter_pstatus(skin_level_mesh,
3133 PSTATUS_SHARED | PSTATUS_MULTISHARED,
3134 PSTATUS_NOT, -1, nullptr);
3135 auto skin_level =
3136 get_bit_skin(bit(get_start_bit() + l), BitRefLevel().set());
3137 skin_level = subtract(skin_level,
3138 skin_level_mesh); // get only internal skins, not
3139 // on the body boundary
3140 // get lower dimension adjacencies (FIXME: add edges if 3D)
3141 Range skin_level_verts;
3142 CHKERR mField.get_moab().get_connectivity(skin_level, skin_level_verts,
3143 true);
3144 skin_level.merge(skin_level_verts);
3145
3146 // remove previous level
3147 bit_prev.set(l - 1);
3148 Range level_prev;
3149 CHKERR bit_mng->getEntitiesByRefLevel(bit_prev, BitRefLevel().set(),
3150 level_prev);
3151 skin.merge(subtract(skin_level, level_prev));
3152 }
3153
3154 return skin;
3155 };
3156
3157 auto resolve_shared = [&](auto &&skin) {
3158 Range tmp_skin = skin;
3159
3160 map<int, Range> map_procs;
3161 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
3162 tmp_skin, &map_procs);
3163
3164 Range from_other_procs; // entities which also exist on other processors
3165 for (auto &m : map_procs) {
3166 if (m.first != mField.get_comm_rank()) {
3167 from_other_procs.merge(m.second);
3168 }
3169 }
3170
3171 auto common = intersect(
3172 skin, from_other_procs); // entities which are on internal skin
3173 skin.merge(from_other_procs);
3174
3175 // entities which are on processor borders, and several processors are not
3176 // true skin.
3177 if (!common.empty()) {
3178 // skin is internal exist on other procs
3179 skin = subtract(skin, common);
3180 }
3181
3182 return skin;
3183 };
3184
3185 // get parents of entities
3186 auto get_parent_level_skin = [&](auto skin) {
3187 Range skin_parents;
3188 CHKERR bit_mng->updateRangeByParent(
3189 skin.subset_by_dimension(SPACE_DIM - 1), skin_parents);
3190 Range skin_parent_verts;
3191 CHKERR mField.get_moab().get_connectivity(skin_parents, skin_parent_verts,
3192 true);
3193 skin_parents.merge(skin_parent_verts);
3194 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
3195 skin_parents);
3196 return skin_parents;
3197 };
3198
3199 auto child_skin = resolve_shared(get_level_skin());
3200 auto parent_skin = get_parent_level_skin(child_skin);
3201
3202 child_skin = subtract(child_skin, parent_skin);
3203 CHKERR bit_mng->setNthBitRefLevel(child_skin, get_skin_child_bit(), true);
3204 CHKERR bit_mng->setNthBitRefLevel(parent_skin, get_skin_parent_bit(), true);
3205
3207 };
3208
3209 // take last level, remove childs on boarder, and set bit
3210 auto set_current = [&]() {
3212 Range last_level;
3213 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_start_bit() + nb_levels - 1),
3214 BitRefLevel().set(), last_level);
3215 Range skin_child;
3216 CHKERR bit_mng->getEntitiesByRefLevel(bit(get_skin_child_bit()),
3217 BitRefLevel().set(), skin_child);
3218
3219 last_level = subtract(last_level, skin_child);
3220 CHKERR bit_mng->setNthBitRefLevel(last_level, get_current_bit(), true);
3222 };
3223
3224 // set bits to levels
3225 CHKERR set_levels(findEntitiesCrossedByPhaseInterface(overlap));
3226 // set bits to skin
3227 CHKERR set_skins();
3228 // set current level bit
3229 CHKERR set_current();
3230
3231 if constexpr (debug) {
3232 if (mField.get_comm_size() == 1) {
3233 for (auto l = 0; l != nb_levels; ++l) {
3234 std::string name = (boost::format("out_level%d.h5m") % l).str();
3235 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_start_bit() + l),
3236 BitRefLevel().set(), name.c_str(), "MOAB",
3237 "PARALLEL=WRITE_PART");
3238 }
3239 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_current_bit()),
3240 BitRefLevel().set(), "current_bit.h5m",
3241 "MOAB", "PARALLEL=WRITE_PART");
3242 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_projection_bit()),
3243 BitRefLevel().set(), "projection_bit.h5m",
3244 "MOAB", "PARALLEL=WRITE_PART");
3245
3246 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_child_bit()),
3247 BitRefLevel().set(), "skin_child_bit.h5m",
3248 "MOAB", "PARALLEL=WRITE_PART");
3249 CHKERR bit_mng->writeBitLevel(BitRefLevel().set(get_skin_parent_bit()),
3250 BitRefLevel().set(), "skin_parent_bit.h5m",
3251 "MOAB", "PARALLEL=WRITE_PART");
3252 }
3253 }
3254
3256};
3257
3259 boost::shared_ptr<FEMethod> fe_top, std::string field_name,
3260 ForcesAndSourcesCore::UserDataOperator::OpType op,
3261 BitRefLevel child_ent_bit,
3262 boost::function<boost::shared_ptr<ForcesAndSourcesCore>()> get_elem,
3263 int verbosity, LogManager::SeverityLevel sev) {
3265
3266 /**
3267 * @brief Collect data from parent elements to child
3268 */
3269 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
3270 add_parent_level =
3271 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
3272 // Evaluate if not last parent element
3273 if (level > 0) {
3274
3275 // Create domain parent FE
3276 auto fe_ptr_current = get_elem();
3277
3278 // Call next level
3279 add_parent_level(
3280 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
3281 fe_ptr_current),
3282 level - 1);
3283
3284 // Add data to curent fe level
3285 if (op == ForcesAndSourcesCore::UserDataOperator::OPSPACE) {
3286
3287 // Only base
3288 parent_fe_pt->getOpPtrVector().push_back(
3289
3291
3292 H1, op, fe_ptr_current,
3293
3294 BitRefLevel().set(), BitRefLevel().set(),
3295
3296 child_ent_bit, BitRefLevel().set(),
3297
3298 verbosity, sev));
3299
3300 } else {
3301
3302 // Filed data
3303 parent_fe_pt->getOpPtrVector().push_back(
3304
3306
3307 field_name, op, fe_ptr_current,
3308
3309 BitRefLevel().set(), BitRefLevel().set(),
3310
3311 child_ent_bit, BitRefLevel().set(),
3312
3313 verbosity, sev));
3314 }
3315 }
3316 };
3317
3318 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
3319 nb_levels);
3320
3322}
3323
3326
3327 if (auto ptr = tsPrePostProc.lock()) {
3328
3329 /**
3330 * @brief cut-off values at nodes, i.e. abs("H") <= 1
3331 *
3332 */
3333 auto cut_off_dofs = [&]() {
3335
3336 auto &m_field = ptr->fsRawPtr->mField;
3337
3338 Range current_verts;
3339 auto bit_mng = m_field.getInterface<BitRefManager>();
3341 bit(get_current_bit()), BitRefLevel().set(), MBVERTEX, current_verts);
3342
3343 auto cut_off_verts = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
3345 for (auto &h : ent_ptr->getEntFieldData()) {
3346 h = cut_off(h);
3347 }
3349 };
3350
3351 auto field_blas = m_field.getInterface<FieldBlas>();
3352 CHKERR field_blas->fieldLambdaOnEntities(cut_off_verts, "H",
3353 &current_verts);
3355 };
3356
3357 CHKERR cut_off_dofs();
3358 }
3359
3360 if (auto ptr = tsPrePostProc.lock()) {
3361 MOFEM_LOG("FS", Sev::inform) << "Run step pre proc";
3362
3363 auto &m_field = ptr->fsRawPtr->mField;
3364 auto simple = m_field.getInterface<Simple>();
3365
3366 // get vector norm
3367 auto get_norm = [&](auto x) {
3368 double nrm;
3369 CHKERR VecNorm(x, NORM_2, &nrm);
3370 return nrm;
3371 };
3372
3373 // refine problem and project data, including theta data
3374 auto refine_problem = [&]() {
3376 MOFEM_LOG("FS", Sev::inform) << "Refine problem";
3377 CHKERR ptr->fsRawPtr->refineMesh(refine_overlap);
3378 CHKERR ptr->fsRawPtr->projectData();
3380 };
3381
3382 // set new jacobin operator, since problem and thus tangent matrix size has
3383 // changed
3384 auto set_jacobian_operators = [&]() {
3386 ptr->subB = createDMMatrix(ptr->solverSubDM);
3387 CHKERR KSPReset(ptr->subKSP);
3389 };
3390
3391 // set new solution
3392 auto set_solution = [&]() {
3394 MOFEM_LOG("FS", Sev::inform) << "Set solution";
3395
3396 PetscObjectState state;
3397
3398 // Record the state, and set it again. This is to fool PETSc that solution
3399 // vector is not updated. Otherwise PETSc will treat every step as a first
3400 // step.
3401
3402 // globSol is updated as result mesh refinement - this is not really set
3403 // a new solution.
3404
3405 CHKERR PetscObjectStateGet(getPetscObject(ptr->globSol.get()), &state);
3406 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), ptr->globSol,
3407 INSERT_VALUES, SCATTER_FORWARD);
3408 CHKERR PetscObjectStateSet(getPetscObject(ptr->globSol.get()), state);
3409 MOFEM_LOG("FS", Sev::verbose)
3410 << "Set solution, vector norm " << get_norm(ptr->globSol);
3412 };
3413
3414 PetscBool is_theta;
3415 PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &is_theta);
3416 if (is_theta) {
3417
3418 CHKERR refine_problem(); // refine problem
3419 CHKERR set_jacobian_operators(); // set new jacobian
3420 CHKERR set_solution(); // set solution
3421
3422 } else {
3423 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
3424 "Sorry, only TSTheta handling is implemented");
3425 }
3426
3427 // Need barriers, some functions in TS solver need are called collectively
3428 // and require the same state of variables
3429 PetscBarrier((PetscObject)ts);
3430
3431 MOFEM_LOG_CHANNEL("SYNC");
3432 MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PreProc done";
3433 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
3434 }
3435
3437}
3438
3440 if (auto ptr = tsPrePostProc.lock()) {
3441 auto &m_field = ptr->fsRawPtr->mField;
3442 MOFEM_LOG_CHANNEL("SYNC");
3443 MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "FS") << "PostProc done";
3444 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
3445 }
3446 return 0;
3447}
3448
3449MoFEMErrorCode TSPrePostProc::tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t,
3450 Vec f, void *ctx) {
3452 if (auto ptr = tsPrePostProc.lock()) {
3453 auto sub_u = ptr->getSubVector();
3454 auto sub_u_t = vectorDuplicate(sub_u);
3455 auto sub_f = vectorDuplicate(sub_u);
3456 auto scatter = ptr->getScatter(sub_u, u, R);
3457
3458 auto apply_scatter_and_update = [&](auto x, auto sub_x) {
3460 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3461 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3462 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3463 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3465 };
3466
3467 CHKERR apply_scatter_and_update(u, sub_u);
3468 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3469
3470 CHKERR TsSetIFunction(ts, t, sub_u, sub_u_t, sub_f, ptr->tsCtxPtr.get());
3471 CHKERR VecScatterBegin(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3472 CHKERR VecScatterEnd(scatter, sub_f, f, INSERT_VALUES, SCATTER_FORWARD);
3473 }
3475}
3476
3477MoFEMErrorCode TSPrePostProc::tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t,
3478 PetscReal a, Mat A, Mat B,
3479 void *ctx) {
3481 if (auto ptr = tsPrePostProc.lock()) {
3482 auto sub_u = ptr->getSubVector();
3483 auto sub_u_t = vectorDuplicate(sub_u);
3484 auto scatter = ptr->getScatter(sub_u, u, R);
3485
3486 auto apply_scatter_and_update = [&](auto x, auto sub_x) {
3488 CHKERR VecScatterBegin(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3489 CHKERR VecScatterEnd(scatter, x, sub_x, INSERT_VALUES, SCATTER_REVERSE);
3490 CHKERR VecGhostUpdateBegin(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3491 CHKERR VecGhostUpdateEnd(sub_x, INSERT_VALUES, SCATTER_FORWARD);
3493 };
3494
3495 CHKERR apply_scatter_and_update(u, sub_u);
3496 CHKERR apply_scatter_and_update(u_t, sub_u_t);
3497
3498 CHKERR TsSetIJacobian(ts, t, sub_u, sub_u_t, a, ptr->subB, ptr->subB,
3499 ptr->tsCtxPtr.get());
3500 }
3502}
3503
3504MoFEMErrorCode TSPrePostProc::tsMonitor(TS ts, PetscInt step, PetscReal t,
3505 Vec u, void *ctx) {
3507 if (auto ptr = tsPrePostProc.lock()) {
3508 auto get_norm = [&](auto x) {
3509 double nrm;
3510 CHKERR VecNorm(x, NORM_2, &nrm);
3511 return nrm;
3512 };
3513
3514 auto sub_u = ptr->getSubVector();
3515 auto scatter = ptr->getScatter(sub_u, u, R);
3516 CHKERR VecScatterBegin(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3517 CHKERR VecScatterEnd(scatter, u, sub_u, INSERT_VALUES, SCATTER_REVERSE);
3518 CHKERR VecGhostUpdateBegin(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3519 CHKERR VecGhostUpdateEnd(sub_u, INSERT_VALUES, SCATTER_FORWARD);
3520
3521 MOFEM_LOG("FS", Sev::verbose)
3522 << "u norm " << get_norm(u) << " u sub nom " << get_norm(sub_u);
3523
3524 CHKERR TsMonitorSet(ts, step, t, sub_u, ptr->tsCtxPtr.get());
3525 }
3527}
3528
3531 if (auto ptr = tsPrePostProc.lock()) {
3532 MOFEM_LOG("FS", Sev::verbose) << "SetUP sub PC";
3533 ptr->subKSP = createKSP(ptr->fsRawPtr->mField.get_comm());
3534 CHKERR KSPSetFromOptions(ptr->subKSP);
3535 }
3537};
3538
3539MoFEMErrorCode TSPrePostProc::pcApply(PC pc, Vec pc_f, Vec pc_x) {
3541 if (auto ptr = tsPrePostProc.lock()) {
3542 auto sub_x = ptr->getSubVector();
3543 auto sub_f = vectorDuplicate(sub_x);
3544 auto scatter = ptr->getScatter(sub_x, pc_x, R);
3545 CHKERR VecScatterBegin(scatter, pc_f, sub_f, INSERT_VALUES,
3546 SCATTER_REVERSE);
3547 CHKERR VecScatterEnd(scatter, pc_f, sub_f, INSERT_VALUES, SCATTER_REVERSE);
3548 CHKERR KSPSetOperators(ptr->subKSP, ptr->subB, ptr->subB);
3549 MOFEM_LOG("FS", Sev::verbose) << "PCShell solve";
3550 CHKERR KSPSolve(ptr->subKSP, sub_f, sub_x);
3551 MOFEM_LOG("FS", Sev::verbose) << "PCShell solve <- done";
3552 CHKERR VecScatterBegin(scatter, sub_x, pc_x, INSERT_VALUES,
3553 SCATTER_FORWARD);
3554 CHKERR VecScatterEnd(scatter, sub_x, pc_x, INSERT_VALUES, SCATTER_FORWARD);
3555 }
3557};
3558
3560 if (auto ptr = tsPrePostProc.lock()) {
3561 auto prb_ptr = ptr->fsRawPtr->mField.get_problem("SUB_SOLVER");
3562 if (auto sub_data = prb_ptr->getSubData()) {
3563 auto is = sub_data->getSmartColIs();
3564 VecScatter s;
3565 if (fr == R) {
3566 CHK_THROW_MESSAGE(VecScatterCreate(x, PETSC_NULLPTR, y, is, &s),
3567 "crate scatter");
3568 } else {
3569 CHK_THROW_MESSAGE(VecScatterCreate(x, is, y, PETSC_NULLPTR, &s),
3570 "crate scatter");
3571 }
3572 return SmartPetscObj<VecScatter>(s);
3573 }
3574 }
3577}
3578
3582
3584
3585 auto &m_field = fsRawPtr->mField;
3586 auto simple = m_field.getInterface<Simple>();
3587
3589
3590 auto dm = simple->getDM();
3591
3593 CHKERR TSSetIFunction(ts, globRes, tsSetIFunction, nullptr);
3594 CHKERR TSSetIJacobian(ts, PETSC_NULLPTR, PETSC_NULLPTR, tsSetIJacobian, nullptr);
3595 CHKERR TSMonitorSet(ts, tsMonitor, fsRawPtr, PETSC_NULLPTR);
3596
3597 SNES snes;
3598 CHKERR TSGetSNES(ts, &snes);
3599 auto snes_ctx_ptr = getDMSnesCtx(dm);
3600
3601 auto set_section_monitor = [&](auto snes) {
3603 CHKERR SNESMonitorSet(snes,
3604 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
3605 void *))MoFEMSNESMonitorFields,
3606 (void *)(snes_ctx_ptr.get()), nullptr);
3608 };
3609
3610 CHKERR set_section_monitor(snes);
3611
3612 auto ksp = createKSP(m_field.get_comm());
3613 CHKERR KSPSetType(ksp, KSPPREONLY); // Run KSP internally in ShellPC
3614 auto sub_pc = createPC(fsRawPtr->mField.get_comm());
3615 CHKERR PCSetType(sub_pc, PCSHELL);
3616 CHKERR PCShellSetSetUp(sub_pc, pcSetup);
3617 CHKERR PCShellSetApply(sub_pc, pcApply);
3618 CHKERR KSPSetPC(ksp, sub_pc);
3619 CHKERR SNESSetKSP(snes, ksp);
3620
3623
3624 CHKERR TSSetPreStep(ts, tsPreProc);
3625 CHKERR TSSetPostStep(ts, tsPostProc);
3626
3628}
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
constexpr int SPACE_DIM
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainParentEle DomainParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
BoundaryEle::UserDataOperator BoundaryEleOp
Kronecker Delta class symmetric.
@ QUIET
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define BITREFLEVEL_SIZE
max number of refinements
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_INVALID_DATA
Definition definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
@ LAST_COORDINATE_SYSTEM
@ CYLINDRICAL
@ CARTESIAN
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
auto cylindrical
Coordinate system scaling factor.
double kappa
auto marker
Create marker bit reference level
auto bubble_device
Bubble device initialization.
double mu_diff
ElementsAndOps< SPACE_DIM >::DomainParentEle DomainParentEle
auto save_range
Save range of entities to file.
auto get_M2
Degenerate mobility function M₂
constexpr int U_FIELD_DIM
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
FTensor::Index< 'j', SPACE_DIM > j
auto get_M2_dh
Derivative of degenerate mobility M₂
double tol
@ R
@ F
auto integration_rule
static char help[]
auto get_M_dh
auto kernel_eye
Eye-shaped interface initialization
double mu_m
double lambda
surface tension
FTensor::Index< 'k', SPACE_DIM > k
OpDomainMassH OpDomainMassG
auto get_M3_dh
Derivative of non-linear mobility M₃
FTensor::Index< 'i', SPACE_DIM > i
double rho_diff
auto init_h
Initialisation function.
constexpr auto t_kd
auto d_phase_function_h
Derivative of phase function with respect to h.
auto get_fe_bit
Get bit reference level from finite element.
constexpr int SPACE_DIM
double rho_ave
FTensor::Index< 'l', SPACE_DIM > l
int coord_type
double W
auto kernel_oscillation
Oscillating interface initialization.
double eta2
auto get_skin_parent_bit
auto get_M0
Constant mobility function M₀
auto get_f_dh
Derivative of double-well potential.
auto get_M3
Non-linear mobility function M₃
double rho_p
double mu_p
auto get_start_bit
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesVector< BASE_DIM, SPACE_DIM, 1 > OpDomainAssembleVector
OpDomainMassH OpDomainMassP
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpDomainMassH
constexpr int BASE_DIM
int nb_levels
auto wetting_angle
Wetting angle function (placeholder)
double h
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpDomainAssembleScalar
auto get_M
double mu_ave
double eps
auto get_dofs_ents_all
Get all entities with DOFs in the problem - used for debugging.
auto phase_function
Phase-dependent material property interpolation.
constexpr IntegrationType I
auto get_skin_projection_bit
auto get_global_size
Get global size across all processors.
auto get_dofs_ents_by_field_name
Get entities of DOFs by field name - used for debugging.
constexpr AssemblyType A
auto get_D
Create deviatoric stress tensor.
auto wetting_angle_sub_stepping
Wetting angle sub-stepping for gradual application.
int refine_overlap
auto get_M0_dh
Derivative of constant mobility.
auto my_max
Maximum function with smooth transition.
FormsIntegrators< BoundaryEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, 1 > OpBoundaryMassL
auto cut_off
Phase field cutoff function.
auto get_skin_child_bit
int order
approximation order
SideEle::UserDataOperator SideOp
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMixScalarTimesDiv< SPACE_DIM, COORD_TYPE > OpMixScalarTimesDiv
FormsIntegrators< BoundaryEleOp >::Assembly< A >::LinearForm< I >::OpBaseTimesScalar< BASE_DIM > OpBoundaryAssembleScalar
double a0
double rho_m
auto get_current_bit
dofs bit used to do calculations
BoundaryEle::UserDataOperator BoundaryEleOp
double eta
auto bit
Create bit reference level.
double md
auto capillary_tube
Capillary tube initialization.
auto get_f
Double-well potential function.
auto get_projection_bit
auto d_cut_off
Derivative of cutoff function.
auto my_min
Minimum function with smooth transition
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpMass< BASE_DIM, U_FIELD_DIM > OpDomainMassU
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition DMMoFEM.cpp:114
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
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:514
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
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:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1234
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition DMMoFEM.cpp:410
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1191
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:536
virtual const DofEntity_multiIndex * get_dofs() const =0
Get the dofs object.
MoFEMErrorCode getEntitiesByDimAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const int dim, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
MoFEMErrorCode getEntitiesByTypeAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const EntityType type, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
MoFEMErrorCode buildSubProblem(const std::string out_name, const std::vector< std::string > &fields_row, const std::vector< std::string > &fields_col, const std::string main_problem, const bool square_matrix=true, const map< std::string, boost::shared_ptr< Range > > *entityMapRow=nullptr, const map< std::string, boost::shared_ptr< Range > > *entityMapCol=nullptr, int verb=VERBOSE)
build sub problem
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
double D
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
auto createKSP(MPI_Comm comm)
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:1046
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition TsCtx.cpp:169
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition TsCtx.cpp:263
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1276
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition DMMoFEM.cpp:434
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
static const bool debug
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:592
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
Definition TsCtx.cpp:56
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVecScatter(Vec x, IS ix, Vec y, IS iy)
Create a Vec Scatter object.
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
Definition Schur.cpp:1296
auto createTS(MPI_Comm comm)
auto createPC(MPI_Comm comm)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscObject getPetscObject(T obj)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition DMMoFEM.hpp:1262
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition DMMoFEM.hpp:1179
int save_every_nth_step
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition plastic.cpp:61
constexpr double t
plate stiffness
Definition plate.cpp:58
constexpr auto field_name
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:56
constexpr double g
FTensor::Index< 'm', 3 > m
MoFEM::Interface & mField
Reference to MoFEM interface.
Definition plastic.cpp:226
Range findParentsToRefine(Range ents, BitRefLevel level, BitRefLevel mask)
Find parent entities that need refinement.
MoFEMErrorCode setupProblem()
Setup problem fields and parameters.
MoFEMErrorCode runProblem()
Main function to run the complete free surface simulation.
FreeSurface(MoFEM::Interface &m_field)
Constructor.
MoFEMErrorCode projectData()
Project solution data between mesh levels.
MoFEMErrorCode boundaryCondition()
Apply boundary conditions and initialize fields.
MoFEMErrorCode readMesh()
Read mesh from input file.
MoFEMErrorCode setParentDofs(boost::shared_ptr< FEMethod > fe_top, std::string field_name, ForcesAndSourcesCore::UserDataOperator::OpType op, BitRefLevel child_ent_bit, boost::function< boost::shared_ptr< ForcesAndSourcesCore >()> get_elem, int verbosity, LogManager::SeverityLevel sev)
Create hierarchy of elements run on parent levels.
std::vector< Range > findEntitiesCrossedByPhaseInterface(size_t overlap)
Find entities on refinement levels.
MoFEMErrorCode assembleSystem()
Assemble system operators and matrices.
MoFEM::Interface & mField
MoFEMErrorCode refineMesh(size_t overlap)
Perform adaptive mesh refinement.
MoFEMErrorCode makeRefProblem()
Create refined problem for mesh adaptation.
MoFEMErrorCode solveSystem()
Solve the time-dependent free surface problem.
Add operators pushing bases from local to physical configuration.
boost::weak_ptr< CacheTuple > getCacheWeakPtr() const
Get the cache weak pointer object.
Boundary condition manager for finite element problem setup.
Managing BitRefLevels.
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =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:118
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
static UId getLoFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
static UId getHiFieldEntityUId(const FieldBitNumber bit, const EntityHandle ent)
Data on single entity (This is passed as argument to DataOperator::doWork)
Structure for user loop methods on finite elements.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Shared pointer to finite element database structure.
Basic algebra on fields.
Definition FieldBlas.hpp:21
MoFEMErrorCode fieldLambdaOnEntities(OneFieldFunctionOnEntities lambda, const std::string field_name, Range *ents_ptr=nullptr)
field lambda
Definition FieldBlas.cpp:50
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
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.
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition Natural.hpp:65
Operator to project base functions from parent entity to child.
Calculate divergence of vector field at integration points.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision scalar field values calculation.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Specialization for double precision vector field values calculation.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Operator to execute finite element instance on parent element. This operator is typically used to pro...
Template struct for dimension-specific finite element types.
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800
bool & getParentAdjacencies()
Get the addParentAdjacencies flag.
Definition Simple.hpp:555
intrusive_ptr for managing petsc objects
PetscReal ts_t
Current time value.
PetscInt ts_step
Current time step number.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
[Push operators to pipeline]
boost::shared_ptr< moab::Core > postProcMesh
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
Post-processing function executed at loop completion.
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< moab::Core > post_proc_mesh, boost::shared_ptr< PostProcEleDomainCont > post_proc, boost::shared_ptr< PostProcEleBdyCont > post_proc_edge, std::pair< boost::shared_ptr< BoundaryEle >, boost::shared_ptr< VectorDouble > > p)
boost::shared_ptr< VectorDouble > liftVec
boost::shared_ptr< BoundaryEle > liftFE
boost::shared_ptr< PostProcEle > postProc
boost::shared_ptr< PostProcEleBdyCont > postProcEdge
boost::shared_ptr< PostProcEleDomainCont > postProc
Set of functions called by PETSc solver used to refine and update mesh.
boost::shared_ptr< TsCtx > tsCtxPtr
virtual ~TSPrePostProc()=default
SmartPetscObj< Vec > globRes
SmartPetscObj< Mat > subB
SmartPetscObj< Vec > globSol
SmartPetscObj< Vec > getSubVector()
Create sub-problem vector.
static MoFEMErrorCode tsMonitor(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Monitor solution during time stepping.
static MoFEMErrorCode pcApply(PC pc, Vec pc_f, Vec pc_x)
Apply preconditioner.
SmartPetscObj< VecScatter > getScatter(Vec x, Vec y, enum FR fr)
Get scatter context for vector operations.
SmartPetscObj< DM > solverSubDM
static MoFEMErrorCode tsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set implicit Jacobian for time stepping.
FreeSurface * fsRawPtr
boost::shared_ptr< SnesCtx > snesCtxPtr
static MoFEMErrorCode tsPreProc(TS ts)
Pre process time step.
TSPrePostProc()=default
static MoFEMErrorCode tsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec f, void *ctx)
Set implicit function for time stepping.
static MoFEMErrorCode pcSetup(PC pc)
Setup preconditioner.
static MoFEMErrorCode tsPostProc(TS ts)
Post process time step.
static MoFEMErrorCode tsPreStage(TS ts)
Pre-stage processing for time stepping.
MoFEMErrorCode tsSetUp(TS ts)
Used to setup TS solver.
SmartPetscObj< KSP > subKSP
constexpr CoordinateTypes COORD_TYPE