Compare commits

...
173 Commits
Author SHA1 Message Date
Sohail Reddy 9a12ddc9c1 Updated apps to use new architecture 2026-02-25 15:19:10 -08:00
Sohail Reddy 1d57520c3c Added field tracking by ID and generalized interface 2026-01-27 10:24:09 -07:00
Sohail Reddy b7753cb4db fixed cmake 2026-01-06 18:59:44 -08:00
Sohail Reddy 842fec483f Merge branch 'master' into multiapp 2026-01-06 18:53:05 -08:00
Sohail Reddy 4b0eafd047 Fixed const correctness 2026-01-06 18:49:38 -08:00
Tzanio Kolev 4a17243415 Merge pull request #5079 from mfem/ode-stage-solve
Update ODESolver to work with ImplicitSolve for stage-state, u, and stage-slope, du/dt
2026-01-02 11:41:38 -08:00
Tzanio Kolev 5731d2f22b Merge pull request #5123 from mfem/fix-nurbs-orientation
Fix nurbs orientations
2025-12-29 14:23:39 -08:00
Tzanio Kolev 150cbcf8c5 Merge pull request #5111 from mfem/sparse-smoothers-transpose
Add MultTranspose for DSmoother and GSSmoother
2025-12-29 13:36:17 -08:00
Will Pazner 1e35d9e49f Merge pull request #5092 from mfem/gpu-uuid
Added a way to get the GPU uuid
2025-12-29 12:04:11 -08:00
Tzanio Kolev ba2ebd318b Merge pull request #5061 from mfem/dg-lor-neumann
Add support for Neumann boundary conditions in the LOR DG preconditioner
2025-12-27 14:25:25 -08:00
Tzanio Kolev a9fd723a4e Merge pull request #5159 from mfem/table-readwrite-fix
Add missing HostReadWriteJ in Table
2025-12-27 14:24:21 -08:00
Sohail Reddy 0fa8a5c90d moved ComputeSlopeFromState function to public 2025-12-26 16:15:25 -08:00
Tzanio Kolev 8e3a7b90fc Merge pull request #5168 from mfem/remove-culib-from-magma
Fix HIP+MAGMA MFEM Link
2025-12-25 04:46:29 -08:00
Tzanio Kolev 2a6f07a2fb Merge pull request #5164 from mfem/particles-navier-dir-fix
Fix directory path in navier-bifurcation
2025-12-24 07:50:14 -08:00
Tzanio Kolev bab0f55127 Merge pull request #5166 from mfem/fix-1d-nurbs
Fix patch offset error in for 1D nurbs
2025-12-24 07:49:00 -08:00
Tzanio Kolev 821b7b31cf Merge pull request #5160 from mfem/task/chapman39/rename-serac-to-smith
Rename Serac to Smith
2025-12-24 07:48:25 -08:00
Tzanio Kolev 8bdb7f9a0d Merge pull request #5163 from mfem/fix-5151
mfem CMake C++ standard
2025-12-24 07:48:02 -08:00
Tzanio Kolev 7f314e25ce Merge pull request #5014 from mfem/nbeams/gko-updates-distrib
Add distributed Ginkgo support
2025-12-24 07:47:31 -08:00
Andrew Ho dee699eede switch to use the cuda runtime API for compatibility with MFEM_GPU_CHECK 2025-12-22 06:41:53 -08:00
Tom Stitt 6eb6a632c9 remove cublas and cusparse from MAGMA_LIB because it breaks hip 2025-12-18 12:58:49 -08:00
Justin Laughlin 5db0bffc2f Address feedback in review - merge my changes with Dylans 2025-12-17 14:55:01 -08:00
Dylan Copeland 3446841376 Name change. 2025-12-17 14:32:16 -08:00
Dylan Copeland c0d5cd54c5 Documentation. 2025-12-17 14:26:04 -08:00
Ido Akkerman eb606dd1b7 Fix patch offset error in for 1D nurbs 2025-12-17 13:27:21 +01:00
Mittal, Ketan 1282722ecf fix directory path due to change in navier location 2025-12-16 12:11:22 -08:00
Andrew Ho d1ca02310d Merge branch 'master' into nbeams/gko-updates-distrib 2025-12-16 12:05:14 -08:00
Andrew Ho b6285fb7eb ensure downstream libraries build with at least the same C++ standard as mfem 2025-12-16 09:31:07 -08:00
Ketan Mittal 81272cb468 Merge branch 'master' into gpu-uuid 2025-12-15 13:03:13 -08:00
Justin Laughlin 9ed202928e Remove repeated unit tests in CMakeLists.txt 2025-12-15 12:27:16 -08:00
Justin Laughlin 809eb996ca Merge branch 'master' into fix-nurbs-orientation 2025-12-15 12:21:40 -08:00
Justin Laughlin 62b27d3989 Add unit test for CorrectPatchTopoOrientations 2025-12-15 12:21:27 -08:00
Justin Laughlin b444c89ab0 Add test nurbs mesh - flipped edges 2025-12-15 11:39:39 -08:00
Alex Tyler Chapman eafb82943b Merge branch 'master' into task/chapman39/rename-serac-to-smith 2025-12-15 10:37:21 -08:00
chapman39 a2243413e6 rename serac to smith 2025-12-15 10:32:16 -08:00
Tzanio Kolev 8452c8cc7f Merge pull request #5153 from mfem/doc-fixes
ex41 and contact doc fixes
2025-12-15 08:39:46 -08:00
Will Pazner bd4df5fad0 Add missing HostReadWriteJ in Table 2025-12-14 11:31:57 -08:00
Socratis PetridesandJan Nikl 3d25fe2e3d Update doc/CodeDocumentation.dox
Co-authored-by: Jan Nikl <nikl1@llnl.gov>
2025-12-12 12:48:59 -08:00
Socratis Petrides 9093d35699 ex41 and contact doc fixes 2025-12-12 11:42:00 -08:00
Will Pazner b718da4997 Merge remote-tracking branch 'origin/master' into sparse-smoothers-transpose
# Conflicts:
#	tests/unit/CMakeLists.txt
2025-12-09 09:17:13 -08:00
Dylan Copeland 3c1f6d96f6 Fixed memory leaks. 2025-12-03 14:27:41 -08:00
Dylan Copeland 6e46adeae2 Unused variable. 2025-12-03 12:40:26 -08:00
Dylan Copeland a25c71d57e Merge branch 'master' of github.com:mfem/mfem into fix-nurbs-orientation 2025-12-03 11:34:59 -08:00
Dylan Copeland de99fd0caa Generalized algorithm for correcting knotvector orientations. Changed the reflector miniapp to output a NURBS mesh for NURBS input. 2025-12-03 11:34:19 -08:00
Andrew Ho 30e470ec12 Merge branch 'master' into gpu-uuid 2025-12-01 12:51:44 -08:00
Andrew Ho 15060b8da4 expanded documentation 2025-12-01 12:39:43 -08:00
Dylan Copeland 117fb91d34 Shadowed variable. 2025-11-28 11:23:44 -08:00
Dylan Copeland 6227904dbd Remove main.cpp. 2025-11-28 11:13:08 -08:00
Dylan Copeland 1b07ccd1b8 Remove testing drivers. 2025-11-28 11:12:04 -08:00
Dylan Copeland d8576eeb3b Set knotvector orientations by sweeping through patches. 2025-11-28 10:38:46 -08:00
Justin Laughlin 7579eb32c9 Cleanup 2025-11-23 12:05:42 -08:00
Justin Laughlin edc818f036 Cleanup 2025-11-23 11:58:11 -08:00
Justin Laughlin 3f4d577042 Style 2025-11-23 11:43:05 -08:00
Justin Laughlin df23247bc8 Generalize CheckPatchTopoOrientation to 3D 2025-11-23 11:42:23 -08:00
Justin Laughlin c9246e28f4 Add CheckPatchTopoOrientation 2025-11-23 11:31:07 -08:00
Justin Laughlin 4b2cc48687 Add 3patch mesh to miniapps/nurbs/meshes 2025-11-20 20:39:28 -08:00
Justin Laughlin 3c910d2651 Cleanup + WIP orientation flipper 2025-11-20 20:38:28 -08:00
Justin Laughlin 2f4f364a27 small fixes to sbox.cpp 2025-11-20 20:37:29 -08:00
Justin Laughlin 46714ce661 Make CheckPatches optional 2025-11-20 12:36:25 -08:00
Sohail Reddy f768ba1c75 Changed ImplicitVariable to ImplicitVariableType 2025-11-20 11:04:10 -08:00
Justin Laughlin d5969dfb94 Flip orientations - minimally working 2025-11-19 22:03:30 -08:00
Sohail Reddy e793c70058 Updated fsi app for constrained inner product on ess_dof 2025-11-19 14:28:38 -08:00
Sohail Reddy 1856912d63 Fixed style 2025-11-19 14:24:23 -08:00
Sohail Reddy 3b6817e36b Added custom inner product to fixed point relaxation methods 2025-11-19 14:24:04 -08:00
Sohail Reddy e4dcfb684f Merge branch 'ode-stage-solve' into multiapp 2025-11-18 23:26:55 -08:00
Sohail Reddy 7dcf77da54 Merge branch 'master' into multiapp 2025-11-18 23:25:23 -08:00
Sohail Reddy 3e823f1da7 updated morpher coupled bc 2025-11-18 23:16:06 -08:00
Will Pazner 5bf0ff54d0 Use relative tolerance in SparseSmoother::EnsureTranspose 2025-11-17 14:21:37 -08:00
Will Pazner 2def15c3f0 constexpr 2025-11-17 10:10:18 -08:00
Will Pazner ac257b0ee9 Add old constructors for backwards compatibility
There is no implicit conversion from int to enum types
2025-11-17 10:10:11 -08:00
Will Pazner 3bdcae8928 Improve Doxygen for DSmoother and GSSmoother 2025-11-17 09:36:52 -08:00
Will Pazner c95ab91de1 Special case for DSmoother::MultTranspose
Don't need to form the transpose when only the action of D (and not A) is
needed.
2025-11-17 09:36:34 -08:00
Andrew Ho 29cfcec850 Merge branch 'master' into gpu-uuid 2025-11-14 00:46:33 -08:00
Sohail Reddy 7c89e72c98 fixed style 2025-11-13 15:11:09 -08:00
Sohail Reddy 5c938034c8 Removed html tables, cleaned documentation 2025-11-13 14:13:30 -08:00
Will Pazner a5f68ad7a6 Add MultTranspose for DSmoother and GSSmoother
For non-symmetric matrices, this requires forming the transpose matrix.

For l1-Jacobi and lumped Jacobi, we don't implement the transpose action for
non-symmetric matrices.
2025-11-12 17:05:53 -08:00
Andrew Ho b993e42c24 Merge branch 'master' into gpu-uuid 2025-11-10 13:16:26 -08:00
Andrew HoandWill Pazner fc4f4fb5ed Update general/device.cpp
Co-authored-by: Will Pazner <11493037+pazner@users.noreply.github.com>
2025-10-27 14:39:39 -07:00
Andrew Ho b74f76bdc4 Added a way to get the GPU uuid 2025-10-27 13:32:49 -07:00
Sohail Reddy 0d7b7d0766 fixed typo in doc 2025-10-24 16:00:04 -07:00
Sohail Reddy a5835794e5 fixed examples to use new variable names 2025-10-24 14:31:21 -07:00
Sohail Reddy db8c62eaf5 Added documentation on switching between solving for stage state and slope, updated examples and variable names. 2025-10-24 14:22:12 -07:00
Sohail Reddy 83232d67d1 disabled Implicit stage-stage solve for SecondOrderODESolvers 2025-10-20 13:46:45 -07:00
Sohail Reddy 9a2fee8b6f moved FPI relaxation member functions to cpp 2025-10-19 21:44:25 -07:00
Sohail Reddy 5fd5b76e47 Updated examples ex9 and ex16 to demonstrate use of implicit stage-state vs stage-slope solve. 2025-10-19 16:51:07 -07:00
Sohail Reddy b57c18b664 Added capability to ImplicitSolve for either stage slope, k=du/dt, or the stage state. Defaults to stage-slope for compatibility with existing TimeDependentOperators 2025-10-19 16:42:33 -07:00
Sohail Reddy 4b38bd79e9 clean up fsi example 2025-10-16 18:05:41 -07:00
Sohail Reddy d1378b8a7b bug fix fpisolver; reset relaxation each call 2025-10-16 18:04:17 -07:00
Sohail Reddy 3972b24f3b fixed ale velocity calculation 2025-10-15 22:19:59 -07:00
Sohail Reddy 48bd66640d added command line option for relaxation factor and fixed mpi_c_bool 2025-10-15 22:19:35 -07:00
Sohail Reddy a230b87a20 Fixed pure virtual function overide if GSLib not built 2025-10-14 18:31:08 -07:00
Sohail Reddy 5f9f782ad2 fixed style 2025-10-14 17:01:51 -07:00
Sohail Reddy 159a259f30 fixed redeclared variable in FPISolver 2025-10-14 16:56:20 -07:00
Sohail Reddy 67da324c26 Merge branch 'master' into multiapp 2025-10-14 16:28:08 -07:00
Sohail Reddy b5cf20e76c added cmake for multiapp miniapps 2025-10-14 16:26:00 -07:00
Sohail Reddy 1510492dfd fixed style 2025-10-14 16:24:23 -07:00
Sohail Reddy 360e662cf8 FSI miniapp and updated elasticity and navier-stokes apps 2025-10-14 16:22:53 -07:00
Sohail Reddy c1ec6a8c72 bug fix in navier miniapp for Step() function 2025-10-14 13:27:04 -07:00
Sohail Reddy 805499d799 bugfix for updating grid functions for transfer 2025-10-14 13:26:15 -07:00
Sohail Reddy 7b0330ffca Removed repetitive code in cht miniapp 2025-10-14 13:23:55 -07:00
Sohail Reddy 68fd37de1d Added completed conjugate heat transfer test case 2025-10-11 17:44:41 -07:00
Sohail Reddy 9f1d2ac249 reverted changes to navier miniapp 2025-10-11 17:43:24 -07:00
Sohail Reddy a594ab5a3a fixed namespace std use 2025-10-11 17:42:57 -07:00
Sohail Reddy 343a8b473e moved fixed point solver to solvers.cpp 2025-10-11 17:13:31 -07:00
Sohail Reddy 911709c619 fixed use of MFEM_USE_MPI in multiapp 2025-10-11 16:10:17 -07:00
Sohail Reddy c303956b9f fixed source field ownership and cleaning 2025-10-11 01:29:50 -07:00
Sohail Reddy c031ec4466 Fixed const correctness and added operator ownership in multiapp 2025-10-10 17:37:12 -07:00
Sohail Reddy 58c56ec9dc added LinkedFieldsCollection to manage linked fields for each application 2025-10-10 00:01:12 -07:00
Will Pazner 9338a6ea61 Add support for Neumann boundary conditions in the LOR DG preconditioner
Don't add boundary penalty terms where markers are unset or if no boundary face
integrator has been added to the form.
2025-10-09 18:17:38 -07:00
nbeamsandAndrew Ho 65b2525ad7 Improve HYPRE/Ginkgo integer type interop
Co-authored-by: Andrew Ho <ho37@llnl.gov>
2025-10-09 02:26:00 +00:00
Andrew Ho 67311104a5 Merge branch 'master' into nbeams/gko-updates-distrib 2025-10-06 16:40:50 -07:00
Sohail Reddy 7db740ede3 Added documentation, updated class names for coupling schemes 2025-09-28 23:36:06 -07:00
Sohail Reddy 0436adfb1c updated fixed point solver 2025-09-28 23:31:33 -07:00
Andrew Ho addb6db73e Merge branch 'master' into nbeams/gko-updates-distrib 2025-09-25 15:24:25 -07:00
Tzanio Kolev d1333c2732 Merge branch 'master' into nbeams/gko-updates-distrib 2025-09-17 03:30:12 -07:00
nbeams 05b60c9c7b Remove sqrt from setting iterative solver tolerances in Ginkgo solvers 2025-09-10 22:21:32 +00:00
nbeams d630bf5250 Add distributed support for Ginkgo solvers and preconditioners 2025-09-10 22:21:27 +00:00
Sohail Reddy a6f2d8981a Merge branch 'master' into multiapp 2025-09-08 18:02:54 -07:00
Sohail Reddy 90edc198e3 update fsi app 2025-09-08 17:58:24 -07:00
Sohail Reddy c691c7c854 Formatted multiapp files and added documentation 2025-09-08 10:21:16 -07:00
Sohail Reddy 7c9105bd0e fixed style 2025-09-08 01:41:26 -07:00
Sohail Reddy 773bac0114 updated navier-stokes miniapp for ALE 2025-09-08 01:29:33 -07:00
Sohail Reddy 440db8189c Added option for ALE to VectorConvectionNLFIntegrator 2025-09-08 01:28:26 -07:00
Sohail Reddy 94efb4cd81 Updated multiapp interface for type-erased operator, field transfer and coupling scheme 2025-09-08 01:27:47 -07:00
Sohail Reddy ce3bdc0dfc Simplified mesh-morphing/diffusion app 2025-09-08 01:24:53 -07:00
Sohail Reddy aaf59d3a83 updated elasticity miniapp to use Schur complement 2025-09-08 01:22:39 -07:00
Sohail Reddy 3b3ffa7d7a updated native navier miniapp to use standard ODESolve::Step signature 2025-09-05 13:53:12 -07:00
Sohail Reddy 2898ec669c Updated native navier miniapp to own/not-own the boundary conditions 2025-09-04 17:19:45 -07:00
Sohail Reddy f34bfbbd8d Added field transfer based on FindPoints 2025-09-04 10:35:50 -07:00
Sohail Reddy 3cec2cba78 updated navier-stokes and elasticity examples 2025-09-03 16:34:45 -07:00
Sohail Reddy f769e2a7d7 updated morpher to diffuse displacement or velocity 2025-09-03 16:34:08 -07:00
Sohail Reddy 493674e650 updated navier-stokes miniapp to work with scaled pressure 2025-09-03 16:31:58 -07:00
Sohail Reddy 13632cde51 updated elasticity miniapp 2025-09-03 16:31:03 -07:00
Sohail Reddy ae63fda0fd Updated multiapp object for multi-stage coupling 2025-08-19 11:54:02 -07:00
Sohail Reddy 48b3e428fb clean up 2025-08-19 11:53:32 -07:00
Sohail Reddy 89dd191866 bug-fix: fixed point solver; updated relaxation schemes to clamp relaxation factors 2025-08-19 11:47:14 -07:00
Sohail Reddy 2a4f9184be Updated Navier-Stokes miniapp for multi-stage coupling and artificial compressibility 2025-08-19 11:40:47 -07:00
Sohail Reddy ff66dce92e updated elasticity and morpher miniapps for multi-stage coupling 2025-08-19 11:38:25 -07:00
Sohail Reddy 01ef886faf Added parallelized DistanceSquared for Vectors 2025-08-19 11:09:47 -07:00
Sohail Reddy 4be67907c4 Working FSI miniapp 2025-08-13 16:12:04 -07:00
Sohail Reddy 9443523dec Updated the multiapp (CoupledApplication) framework 2025-08-13 16:10:19 -07:00
Sohail Reddy fac6dfe283 Added mesh morphing app 2025-08-13 16:09:50 -07:00
Sohail Reddy 7776766ee7 updated FPISolver 2025-08-13 16:09:09 -07:00
Sohail Reddy 8e24224a18 updated cht example...unused atm 2025-08-13 16:08:33 -07:00
Sohail Reddy fa4c607622 updated Navier-Stokes miniapp for partial assembly 2025-08-13 16:07:09 -07:00
Sohail Reddy 29bd541a3e updated elasticity miniapp 2025-08-13 16:03:23 -07:00
Sohail Reddy 369f82adad Added function to fetch TimeDependantOperator from ODESolver 2025-08-13 16:01:39 -07:00
Sohail Reddy e2bc1aefd8 Working elasticity miniapp 2025-08-01 19:09:05 -07:00
Sohail Reddy 033814fe70 fixed code style 2025-07-31 10:14:48 -07:00
Sohail Reddy f4d4a5c4da Working fully implicit incompressible navier-stokes solver 2025-07-31 10:13:33 -07:00
Sohail Reddy 6b8d341478 Renamed some multiapp member functions 2025-07-31 09:56:27 -07:00
Sohail Reddy 7a02905d32 Updated constructors for some bilinear and nonlinearform integrator to now take integration rule as input 2025-07-31 09:53:16 -07:00
Sohail Reddy a4c2ad2f3c Updated type erasure to support non-mfem apps 2025-07-21 10:27:35 -07:00
Sohail Reddy 4a11684323 Working monolithic (linear and nonlinear) solves for multiapp 2025-07-17 18:22:41 -07:00
Sohail Reddy 7b4b2de092 Merge branch 'master' into multiapp 2025-07-15 16:37:35 -07:00
Sohail Reddy 285aa03049 conv-diff mulitapp example 2025-07-15 16:35:56 -07:00
Sohail Reddy 27af36adc5 Changed class naming and fixed point solver and relaxation schemes 2025-07-15 16:35:03 -07:00
Sohail Reddy de781c05b8 Fixed function pointer for relaxation method in FPISolver 2025-07-10 17:04:58 -07:00
Sohail Reddy 125f5c9e48 Added relaxation methods to fixed point solvers and updated multiapp for explicit coupling 2025-07-10 16:59:03 -07:00
Sohail Reddy 0d6c8fbbf1 Updated multiapp partition implicit solve to use fixed point Solver object 2025-07-09 19:48:04 -07:00
Sohail Reddy 34ebc0f158 Changed scheme naming convention and added/separated Addtive and Alternating Schwarz coupling schemes 2025-07-08 18:16:32 -07:00
Sohail Reddy 8db7813d86 Working paritioned implicit solves 2025-07-08 11:06:27 -07:00
Sohail Reddy d31e5c8af0 Merge branch 'master' into multiapp 2025-07-02 20:46:57 -07:00
Sohail Reddy d7e43cfc74 added test case for two-way partitioned coupled solve 2025-07-02 20:45:51 -07:00
Sohail Reddy 82193e0404 Merge branch 'master' into multiapp 2025-06-16 15:09:01 -07:00
Sohail Reddy 5796633307 updated example for implicit time-integration 2025-06-16 15:05:22 -07:00
Sohail Reddy 3d55c021e3 updated header copyright 2025-06-12 14:47:11 -07:00
Sohail Reddy a30303ecf9 fixed type checking with SFINAE 2025-06-12 14:31:05 -07:00
Sohail Reddy 9287488024 Merge branch 'master' into multiapp 2025-06-11 14:11:32 -07:00
Sohail Reddy 9aa2866b0d Merge branch 'master' into multiapp 2025-05-29 14:00:03 -07:00
Sohail Reddy ad16d2f5ae added field transfer, multiapp abstractions and working one-way coupling 2025-05-29 13:58:58 -07:00
Sohail Reddy c552925a54 Merge branch 'master' into multiapp 2025-03-13 15:41:44 -07:00
Sohail Reddy c3b9b62336 updated mfem_application to perform index controlled operators 2025-02-20 11:28:02 -08:00
Sohail Reddy c36e5facd1 initial multiapp skeleton with MFEMApplication class and type-erased multiapp 2024-10-02 17:11:14 -07:00
62 changed files with 9430 additions and 367 deletions
+1
View File
@@ -723,6 +723,7 @@ set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX})
# Declaring the library
mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
target_compile_features(mfem PUBLIC cxx_std_${CMAKE_CXX_STANDARD})
# message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES} ${TPL_TARGETS})
if (TPL_TARGETS)
+3 -1
View File
@@ -725,7 +725,9 @@ The specific libraries and their options are:
URL: https://ginkgo-project.github.io
Options: GINKGO_OPT, GINKGO_LIB, GINKGO_DIR, GINKGO_BUILD_TYPE (Release or
Debug).
Versions: Ginkgo >= 1.9.0.
Versions: Ginkgo >= 1.9.0. When building Ginkgo with distributed support, a
recent version of the "develop" branch is required (1.11 as defined
in include/ginkgo/config.hpp).
- AmgX (optional), used when MFEM_USE_AMGX = YES.
URL: https://github.com/NVIDIA/AMGX
+1 -1
View File
@@ -407,7 +407,7 @@ AMGX_LIB = -L$(AMGX_DIR)/lib -lamgx -lcusparse -lcusolver -lcublas -lnvToolsExt
# MAGMA library configuration
MAGMA_DIR = @MFEM_DIR@/../magma
MAGMA_OPT = -I$(MAGMA_DIR)/include
MAGMA_LIB = -L$(MAGMA_DIR)/lib -l:libmagma.a -lcublas -lcusparse $(LAPACK_LIB)
MAGMA_LIB = -L$(MAGMA_DIR)/lib -l:libmagma.a $(LAPACK_LIB)
# GnuTLS library configuration
GNUTLS_OPT =
+4 -1
View File
@@ -117,6 +117,8 @@ namespace mfem {
* - <a class="el" href="ex39p_8cpp_source.html">Example 39p</a>: parallel named mesh attributes
* - <a class="el" href="ex40_8cpp_source.html">Example 40</a>: eikonal equation
* - <a class="el" href="ex40p_8cpp_source.html">Example 40p</a>: parallel eikonal equation
* - <a class="el" href="ex41_8cpp_source.html">Example 41</a>: DG/CG IMEX time-dependent advection-diffusion
* - <a class="el" href="ex41p_8cpp_source.html">Example 41p</a>: parallel DG/CG IMEX time-dependent advection-diffusion
*
* <H4>AmgX Examples</H4>
* - Variants of Examples
@@ -234,7 +236,8 @@ namespace mfem {
* - <a class="el" href="miniapps_2performance_2ex1_8cpp_source.html">HPC Example 1</a>: high-performance nodal H1 FEM for the Poisson problem
* - <a class="el" href="miniapps_2performance_2ex1p_8cpp_source.html">HPC Example 1p</a>: high-performance parallel nodal H1 FEM for the Poisson problem
* - <a class="el" href="generate__random__field_8cpp_source.html">SPDE Solvers</a>: SPDE solver random field generation
* - <a class="el" href="contact-patch-test_8cpp_source.html">Contact</a>: mortar contact patch test for elasticity
* - <a class="el" href="contact-patch-test_8cpp_source.html">Tribol</a>: mortar contact patch test for elasticity
* - <a class="el" href="contact_8cpp_source.html">Contact</a>: Frictionless contact examples using <a class="el" href="classmfem_1_1IPSolver.html#details">IP optimization</a> and the <a class="el" href="classmfem_1_1AMGFSolver.html#details">AMGF solver</a>
* - <a class="el" href="multidomain_8cpp_source.html">Multidomain miniapp</a>: Multidomain and Submesh demonstration miniapp
* - <a class="el" href="pdiffusion_8cpp_source.html">DPG Diffusion example</a>: DPG formulation for the diffusion problem
* - <a class="el" href="pmaxwell_8cpp_source.html">DPG Maxwell example</a>: DPG formulation for the indefinite Maxwell problem
+29 -6
View File
@@ -105,6 +105,7 @@ int main(int argc, char *argv[])
bool visualization = true;
bool visit = false;
int vis_steps = 5;
bool solve_implicit_state = false;
int precision = 8;
cout.precision(precision);
@@ -126,6 +127,9 @@ int main(int argc, char *argv[])
"Alpha coefficient.");
args.AddOption(&kappa, "-k", "--kappa",
"Kappa coefficient offset.");
args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
"-imp-slope", "--implicit-slope",
"Implicitly solve for stage state or slope.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
@@ -179,6 +183,11 @@ int main(int argc, char *argv[])
// 7. Initialize the conduction operator and the visualization.
ConductionOperator oper(fespace, alpha, kappa, u);
using ImplicitVariableType = ConductionOperator::ImplicitVariableType;
ImplicitVariableType imp_var = solve_implicit_state ?
ImplicitVariableType::STATE
: ImplicitVariableType::SLOPE;
oper.SetImplicitVariableType(imp_var);
u_gf.SetFromTrueDofs(u);
{
@@ -316,11 +325,14 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
}
void ConductionOperator::ImplicitSolve(const real_t dt,
const Vector &u, Vector &du_dt)
const Vector &u, Vector &k)
{
// Solve the equation:
// du_dt = M^{-1}*[-K(u + dt*du_dt)]
// for du_dt, where K is linearized by using u from the previous timestep
// M*k = -K(u + dt*k) for k = du/dt, if solving for stage-slope
// or
// M*k = -dt*K(k) + M*u for k = u_s, if solving for stage-state
// where K is linearized by using u from the previous timestep, and
// the stage-state and slope relation: du/dt = (u_s - u)/dt.
if (!T)
{
T = Add(1.0, Mmat, dt, Kmat);
@@ -328,9 +340,20 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
T_solver.SetOperator(*T);
}
MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
Kmat.Mult(u, z);
z.Neg();
T_solver.Mult(z, du_dt);
// Construct current right-hand side for stage state vs. slope solve
if (ImplicitVarTypeIsState())
{
// k, on return, is the stage value u_s
Mmat.Mult(u, z);
}
else
{
// k, on return, is the stage slope du/dt
Kmat.Mult(u, z);
z.Neg();
}
T_solver.Mult(z, k);
}
void ConductionOperator::SetParameters(const Vector &u)
+29 -6
View File
@@ -115,6 +115,7 @@ int main(int argc, char *argv[])
bool visit = false;
int vis_steps = 5;
bool adios2 = false;
bool solve_implicit_state = false;
int precision = 8;
cout.precision(precision);
@@ -138,6 +139,9 @@ int main(int argc, char *argv[])
"Alpha coefficient.");
args.AddOption(&kappa, "-k", "--kappa",
"Kappa coefficient offset.");
args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
"-imp-slope", "--implicit-slope",
"Implicitly solve for stage state or slope.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
@@ -212,6 +216,11 @@ int main(int argc, char *argv[])
// 9. Initialize the conduction operator and the VisIt visualization.
ConductionOperator oper(fespace, alpha, kappa, u);
using ImplicitVariableType = ConductionOperator::ImplicitVariableType;
ImplicitVariableType imp_var = solve_implicit_state ?
ImplicitVariableType::STATE
: ImplicitVariableType::SLOPE;
oper.SetImplicitVariableType(imp_var);
u_gf.SetFromTrueDofs(u);
{
@@ -407,11 +416,14 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
}
void ConductionOperator::ImplicitSolve(const real_t dt,
const Vector &u, Vector &du_dt)
const Vector &u, Vector &k)
{
// Solve the equation:
// du_dt = M^{-1}*[-K(u + dt*du_dt)]
// for du_dt, where K is linearized by using u from the previous timestep
// M*k = -K(u + dt*k) for k = du/dt, if solving for stage-slope
// or
// M*k = -dt*K(k) + M*u for k = u_s, if solving for stage-state
// where K is linearized by using u from the previous timestep, and
// the stage-state and slope relation: du/dt = (u_s - u)/dt.
if (!T)
{
T = Add(1.0, Mmat, dt, Kmat);
@@ -419,9 +431,20 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
T_solver.SetOperator(*T);
}
MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
Kmat.Mult(u, z);
z.Neg();
T_solver.Mult(z, du_dt);
// Construct current right-hand side for stage state vs. slope solve
if (ImplicitVarTypeIsState())
{
// k, on return, is the stage value u
Mmat.Mult(u, z);
}
else
{
// k, on return, is the stage slope du/dt
Kmat.Mult(u, z);
z.Neg();
}
T_solver.Mult(z, k);
}
void ConductionOperator::SetParameters(const Vector &u)
+20 -1
View File
@@ -160,6 +160,7 @@ int main(int argc, char *argv[])
bool paraview = false;
bool binary = false;
int vis_steps = 5;
bool solve_implicit_state = false;
int precision = 8;
cout.precision(precision);
@@ -187,6 +188,9 @@ int main(int argc, char *argv[])
"Final time; start time is 0.");
args.AddOption(&dt, "-dt", "--time-step",
"Time step.");
args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
"-imp-slope", "--implicit-slope",
"Implicitly solve for stage state or slope.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
@@ -366,6 +370,11 @@ int main(int argc, char *argv[])
// right-hand side, and perform time-integration (looping over the time
// iterations, ti, with a time-step dt).
FE_Evolution adv(m, k, b);
using ImplicitVariableType = FE_Evolution::ImplicitVariableType;
ImplicitVariableType imp_var = solve_implicit_state ?
ImplicitVariableType::STATE
: ImplicitVariableType::SLOPE;
adv.SetImplicitVariableType(imp_var);
real_t t = 0.0;
adv.SetTime(t);
@@ -459,7 +468,17 @@ void FE_Evolution::ImplicitSolve(const real_t dt, const Vector &x, Vector &k)
{
MFEM_VERIFY(dg_solver != NULL,
"Implicit time integration is not supported with partial assembly");
K.Mult(x, z);
// Construct current right-hand side for stage state vs. slope solve
if (ImplicitVarTypeIsState())
{
// k, on return, is the stage value u
M.Mult(x, z);
}
else
{
// k, on return, is the stage slope du/dt
K.Mult(x, z);
}
z += b;
dg_solver->SetTimeStep(dt);
dg_solver->Mult(z, k);
+20 -1
View File
@@ -257,6 +257,7 @@ int main(int argc, char *argv[])
bool adios2 = false;
bool binary = false;
int vis_steps = 5;
bool solve_implicit_state = false;
#if MFEM_HYPRE_VERSION >= 21800
PrecType prec_type = PrecType::AIR;
#else
@@ -290,6 +291,9 @@ int main(int argc, char *argv[])
"Final time; start time is 0.");
args.AddOption(&dt, "-dt", "--time-step",
"Time step.");
args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
"-imp-slope", "--implicit-slope",
"Implicitly solve for stage state or slope.");
args.AddOption((int *)&prec_type, "-pt", "--prec-type", "Preconditioner for "
"implicit solves. 0 for ILU, 1 for pAIR-AMG.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -536,6 +540,11 @@ int main(int argc, char *argv[])
// right-hand side, and perform time-integration (looping over the time
// iterations, ti, with a time-step dt).
FE_Evolution adv(*m, *k, *B, prec_type);
using ImplicitVariableType = FE_Evolution::ImplicitVariableType;
ImplicitVariableType imp_var = solve_implicit_state ?
ImplicitVariableType::STATE
: ImplicitVariableType::SLOPE;
adv.SetImplicitVariableType(imp_var);
real_t t = 0.0;
adv.SetTime(t);
@@ -676,7 +685,17 @@ FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
// (M - dt*K) d = K*u + b
void FE_Evolution::ImplicitSolve(const real_t dt, const Vector &x, Vector &k)
{
K->Mult(x, z);
// Construct current right-hand side for stage state vs. slope solve
if (ImplicitVarTypeIsState())
{
// k, on return, is the stage value u
M->Mult(x, z);
}
else
{
// k, on return, is the stage slope du/dt
K->Mult(x, z);
}
z += b;
dg_solver->SetTimeStep(dt);
dg_solver->Mult(z, k);
+6
View File
@@ -14,6 +14,12 @@ list(APPEND GINKGO_EXAMPLES_SRCS
ex1.cpp
)
if (MFEM_USE_MPI AND GINKGO_BUILD_MPI)
list(APPEND GINKGO_EXAMPLES_SRCS
ex1p.cpp
)
endif()
# Include the source directory where mfem.hpp and mfem-performance.hpp are.
include_directories(BEFORE ${PROJECT_BINARY_DIR})
+3 -3
View File
@@ -207,7 +207,7 @@ int main(int argc, char *argv[])
Ginkgo::IcPreconditioner ginkgo_precond(exec, "paric", 30);
Ginkgo::CGSolver ginkgo_solver(exec, ginkgo_precond);
ginkgo_solver.SetPrintLevel(print_lvl);
ginkgo_solver.SetRelTol(1e-12);
ginkgo_solver.SetRelTol(sqrt(1e-12));
ginkgo_solver.SetAbsTol(0.0);
ginkgo_solver.SetMaxIter(400);
ginkgo_solver.SetOperator(*(A.Ptr()));
@@ -225,7 +225,7 @@ int main(int argc, char *argv[])
Ginkgo::MFEMPreconditioner gko_M(exec, M);
Ginkgo::CGSolver ginkgo_solver(exec, gko_M);
ginkgo_solver.SetPrintLevel(print_lvl);
ginkgo_solver.SetRelTol(1e-12);
ginkgo_solver.SetRelTol(sqrt(1e-12));
ginkgo_solver.SetAbsTol(0.0);
ginkgo_solver.SetMaxIter(400);
ginkgo_solver.SetOperator(*(A.Ptr()));
@@ -283,7 +283,7 @@ int main(int argc, char *argv[])
Ginkgo::MFEMPreconditioner gko_M(exec, M);
Ginkgo::CGSolver ginkgo_solver(exec, gko_M);
ginkgo_solver.SetPrintLevel(print_lvl);
ginkgo_solver.SetRelTol(1e-12);
ginkgo_solver.SetRelTol(sqrt(1e-12));
ginkgo_solver.SetAbsTol(0.0);
ginkgo_solver.SetMaxIter(400);
ginkgo_solver.SetOperator(*(A.Ptr()));
+436
View File
@@ -0,0 +1,436 @@
// MFEM Example 1 - Parallel Version
// GINKGO Modification
//
// Compile with: make ex1p
//
// Sample runs: mpirun -np 4 ex1p -m ../data/square-disc.mesh
// mpirun -np 4 ex1p -m ../data/star.mesh
// mpirun -np 4 ex1p -m ../data/star-mixed.mesh
// mpirun -np 4 ex1p -m ../data/escher.mesh
// mpirun -np 4 ex1p -m ../data/fichera.mesh
// mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
// mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
// mpirun -np 4 ex1p -m ../data/octahedron.mesh -o 1
// mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
// mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
// mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
// mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
// mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
// mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
// mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
// mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
// mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
// mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
// mpirun -np 4 ex1p -m ../data/star-surf.mesh
// mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
// mpirun -np 4 ex1p -m ../data/inline-segment.mesh
// mpirun -np 4 ex1p -m ../data/amr-quad.mesh
// mpirun -np 4 ex1p -m ../data/amr-hex.mesh
// mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
// mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
//
// Device sample runs:
// mpirun -np 4 ex1p -pa -d cuda
// mpirun -np 4 ex1p -fa -d cuda
// mpirun -np 4 ex1p -pa -d occa-cuda
// mpirun -np 4 ex1p -pa -d raja-omp
// mpirun -np 4 ex1p -pa -d ceed-cpu
// mpirun -np 4 ex1p -pa -d ceed-cpu -o 4 -a
// mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/square-mixed.mesh
// mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/fichera-mixed.mesh
// * mpirun -np 4 ex1p -pa -d ceed-cuda
// * mpirun -np 4 ex1p -pa -d ceed-hip
// mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
// mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/square-mixed.mesh
// mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/fichera-mixed.mesh
// mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
//
// Description: This example code demonstrates the use of MFEM to define a
// simple finite element discretization of the Laplace problem
// -Delta u = 1 with homogeneous Dirichlet boundary conditions.
// Specifically, we discretize using a FE space of the specified
// order, or if order < 1 using an isoparametric/isogeometric
// space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
// NURBS mesh, etc.)
//
// The example highlights the use of mesh refinement, finite
// element grid functions, as well as linear and bilinear forms
// corresponding to the left-hand side and right-hand side of the
// discrete linear system. We also cover the explicit elimination
// of essential boundary conditions, static condensation, and the
// optional connection to the GLVis tool for visualization.
#include "mfem.hpp"
#include <fstream>
#include <iostream>
#ifndef MFEM_USE_GINKGO
#error This example requires that MFEM is built with MFEM_USE_GINKGO=YES
#endif
using namespace std;
using namespace mfem;
int main(int argc, char *argv[])
{
// 1. Initialize MPI and HYPRE.
Mpi::Init();
int num_procs = Mpi::WorldSize();
int myid = Mpi::WorldRank();
Hypre::Init();
// 2. Parse command-line options.
const char *mesh_file = "../../data/star.mesh";
int order = 1;
bool static_cond = false;
bool pa = false;
bool fa = false;
const char *device_config = "cpu";
bool visualization = true;
int solver_config = 0;
int print_lvl = 1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree) or -1 for"
" isoparametric space.");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
"--no-partial-assembly", "Enable Partial Assembly.");
args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
"--no-full-assembly", "Enable Full Assembly.");
args.AddOption(&device_config, "-d", "--device",
"Device configuration string, see Device::Configure().");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.AddOption(&solver_config, "-s", "--solver-config",
"Solver and preconditioner combination: \n\t"
" 0 - Ginkgo solver and Ginkgo preconditioner, \n\t"
" 1 - Ginkgo solver and MFEM preconditioner, \n\t"
" 2 - MFEM solver and Ginkgo preconditioner, \n\t"
" 3 - MFEM solver and MFEM preconditioner.");
args.AddOption(&print_lvl, "-pl", "--print-level",
"Print level for iterative solver (1 prints every iteration).");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
// 3. Enable hardware devices such as GPUs, and programming models such as
// CUDA, OCCA, RAJA and OpenMP based on command line options.
Device device(device_config);
device.SetGPUAwareMPI(true);
if (myid == 0) { device.Print(); }
// 4. Read the (serial) mesh from the given mesh file on all processors. We
// can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
// and volume meshes with the same code.
Mesh mesh(mesh_file, 1, 1);
int dim = mesh.Dimension();
// 5. Refine the serial mesh on all processors to increase the resolution. In
// this example we do 'ref_levels' of uniform refinement. We choose
// 'ref_levels' to be the largest number that gives a final mesh with no
// more than 10,000 elements.
{
int ref_levels =
(int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
for (int l = 0; l < ref_levels; l++)
{
mesh.UniformRefinement();
}
}
// 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
// this mesh further in parallel to increase the resolution. Once the
// parallel mesh is defined, the serial mesh can be deleted.
ParMesh pmesh(MPI_COMM_WORLD, mesh);
mesh.Clear();
{
int par_ref_levels = 2;
for (int l = 0; l < par_ref_levels; l++)
{
pmesh.UniformRefinement();
}
}
// 7. Define a parallel finite element space on the parallel mesh. Here we
// use continuous Lagrange finite elements of the specified order. If
// order < 1, we instead use an isoparametric/isogeometric space.
FiniteElementCollection *fec;
bool delete_fec;
if (order > 0)
{
fec = new H1_FECollection(order, dim);
delete_fec = true;
}
else if (pmesh.GetNodes())
{
fec = pmesh.GetNodes()->OwnFEC();
delete_fec = false;
if (myid == 0)
{
cout << "Using isoparametric FEs: " << fec->Name() << endl;
}
}
else
{
fec = new H1_FECollection(order = 1, dim);
delete_fec = true;
}
ParFiniteElementSpace fespace(&pmesh, fec);
HYPRE_BigInt size = fespace.GlobalTrueVSize();
if (myid == 0)
{
cout << "Number of finite element unknowns: " << size << endl;
}
// 8. Determine the list of true (i.e. parallel conforming) essential
// boundary dofs. In this example, the boundary conditions are defined
// by marking all the boundary attributes from the mesh as essential
// (Dirichlet) and converting them to a list of true dofs.
Array<int> ess_tdof_list;
if (pmesh.bdr_attributes.Size())
{
Array<int> ess_bdr(pmesh.bdr_attributes.Max());
ess_bdr = 1;
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}
// 9. Set up the parallel linear form b(.) which corresponds to the
// right-hand side of the FEM linear system, which in this case is
// (1,phi_i) where phi_i are the basis functions in fespace.
ParLinearForm b(&fespace);
ConstantCoefficient one(1.0);
b.AddDomainIntegrator(new DomainLFIntegrator(one));
b.Assemble();
// 10. Define the solution vector x as a parallel finite element grid
// function corresponding to fespace. Initialize x with initial guess of
// zero, which satisfies the boundary conditions.
ParGridFunction x(&fespace);
x = 0.0;
// 11. Set up the parallel bilinear form a(.,.) on the finite element space
// corresponding to the Laplacian operator -Delta, by adding the
// Diffusion domain integrator.
ParBilinearForm a(&fespace);
if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
if (fa)
{
a.SetAssemblyLevel(AssemblyLevel::FULL);
// Sort the matrix column indices when running on GPU or with OpenMP (i.e.
// when Device::IsEnabled() returns true). This makes the results
// bit-for-bit deterministic at the cost of somewhat longer run time.
a.EnableSparseMatrixSorting(Device::IsEnabled());
}
a.AddDomainIntegrator(new DiffusionIntegrator(one));
// 12. Assemble the parallel bilinear form and the corresponding linear
// system, applying any necessary transformations such as: parallel
// assembly, eliminating boundary conditions, applying conforming
// constraints for non-conforming AMR, static condensation, etc.
if (static_cond) { a.EnableStaticCondensation(); }
a.Assemble();
OperatorPtr A;
Vector B, X;
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
// 13. Solve the linear system A X = B.
if (!pa)
{
switch (solver_config)
{
// Solve the linear system with CG + Schwarz (with IC) from Ginkgo
case 0:
{
if (myid == 0) { cout << "Using Ginkgo solver + preconditioner...\n"; }
Ginkgo::GinkgoExecutor exec(device);
Ginkgo::IcPreconditioner local_solver(exec, "exact");
Ginkgo::SchwarzPreconditioner gko_M(exec, MPI_COMM_WORLD, local_solver);
Ginkgo::CGSolver ginkgo_solver(exec, MPI_COMM_WORLD, gko_M);
ginkgo_solver.SetPrintLevel(print_lvl);
ginkgo_solver.SetRelTol(sqrt(1e-12));
ginkgo_solver.SetAbsTol(0.0);
ginkgo_solver.SetMaxIter(400);
ginkgo_solver.SetOperator(*(A.Ptr()));
ginkgo_solver.Mult(B, X);
break;
}
// Solve the linear system with CG from Ginkgo + MFEM preconditioner
case 1:
{
if (myid == 0) { cout << "Using Ginkgo solver + MFEM preconditioner...\n"; }
Ginkgo::GinkgoExecutor exec(device);
//Create MFEM preconditioner and wrap it for Ginkgo's use.
HypreBoomerAMG M((HypreParMatrix&)(*A));
Ginkgo::MFEMPreconditioner gko_M(exec, M, MPI_COMM_WORLD);
Ginkgo::CGSolver ginkgo_solver(exec, MPI_COMM_WORLD, gko_M);
ginkgo_solver.SetPrintLevel(print_lvl);
ginkgo_solver.SetRelTol(sqrt(1e-12));
ginkgo_solver.SetAbsTol(0.0);
ginkgo_solver.SetMaxIter(400);
ginkgo_solver.SetOperator(*(A.Ptr()));
ginkgo_solver.Mult(B, X);
break;
}
// Ginkgo Schwarz preconditioner (local ParIC) + MFEM CG solver
case 2:
{
if (myid == 0) { cout << "Using MFEM solver + Ginkgo preconditioner...\n"; }
Ginkgo::GinkgoExecutor exec(device);
Ginkgo::IcPreconditioner local_M(exec, "exact");
Ginkgo::SchwarzPreconditioner M(exec, MPI_COMM_WORLD, local_M);
M.SetOperator(*(A.Ptr())); // Generate the preconditioner for the matrix A.
CGSolver cg(MPI_COMM_WORLD);
cg.SetRelTol(sqrt(1e-12));
cg.SetMaxIter(400);
cg.SetPrintLevel(1);
cg.SetPreconditioner(M);
cg.SetOperator(*A);
cg.Mult(B, X);
break;
}
// MFEM solver + MFEM preconditioner
case 3:
{
if (myid == 0) { cout << "Using MFEM solver + MFEM preconditioner...\n"; }
HypreBoomerAMG M((HypreParMatrix&)(*A));
CGSolver cg(MPI_COMM_WORLD);
cg.SetRelTol(sqrt(1e-12));
cg.SetMaxIter(400);
cg.SetPrintLevel(1);
cg.SetPreconditioner(M);
cg.SetOperator(*A);
cg.Mult(B, X);
break;
}
} // End switch on solver_config
}
// Partial assembly mode. Cannot use Ginkgo preconditioners, but can use Ginkgo
// solvers.
else
{
if (UsesTensorBasis(fespace))
{
// Use Jacobi preconditioning in partial assembly mode.
OperatorJacobiSmoother M(a, ess_tdof_list);
switch (solver_config)
{
case 0:
{
if (myid == 0) { cout << "Using Ginkgo solver + preconditioner...\n"; }
MFEM_ABORT("Cannot use Ginkgo preconditioner in partial assembly mode.\n"
" Try -s 1 to test Ginkgo solver with an MFEM preconditioner.");
break;
}
// Use Ginkgo solver with MFEM preconditioner
case 1:
{
if (myid == 0) { cout << "Using Ginkgo solver + MFEM preconditioner...\n"; }
Ginkgo::GinkgoExecutor exec(device);
// Wrap MFEM preconditioner for Ginkgo's use.
Ginkgo::MFEMPreconditioner gko_M(exec, M, MPI_COMM_WORLD);
Ginkgo::CGSolver ginkgo_solver(exec, MPI_COMM_WORLD, gko_M);
ginkgo_solver.SetPrintLevel(print_lvl);
ginkgo_solver.SetRelTol(sqrt(1e-12));
ginkgo_solver.SetAbsTol(0.0);
ginkgo_solver.SetMaxIter(400);
ginkgo_solver.SetOperator(*(A.Ptr()));
ginkgo_solver.Mult(B, X);
break;
}
// No Ginkgo preconditioners work with matrix-free; error
case 2:
{
if (myid == 0) { cout << "Using Ginkgo solver + preconditioner...\n"; }
MFEM_ABORT("Cannot use Ginkgo preconditioner in partial assembly mode.\n"
" Try -s 1 to test Ginkgo solver with an MFEM preconditioner.");
break;
}
// Use MFEM solver and preconditioner
case 3:
{
if (myid == 0) { cout << "Using MFEM solver + MFEM preconditioner...\n"; }
CGSolver cg(MPI_COMM_WORLD);
cg.SetRelTol(sqrt(1e-12));
cg.SetMaxIter(400);
cg.SetPrintLevel(1);
cg.SetPreconditioner(M);
cg.SetOperator(*A);
cg.Mult(B, X);
break;
}
} // End switch on solver_config
}
else // CG with no preconditioning
{
if (myid == 0) { cout << "Using MFEM solver + no preconditioner...\n"; }
CGSolver cg(MPI_COMM_WORLD);
cg.SetRelTol(sqrt(1e-12));
cg.SetMaxIter(400);
cg.SetPrintLevel(1);
cg.SetOperator(*A);
cg.Mult(B, X);
}
}
// 14. Recover the parallel grid function corresponding to X. This is the
// local finite element solution on each processor.
a.RecoverFEMSolution(X, b, x);
// 15. Save the refined mesh and the solution in parallel. This output can
// be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
{
ostringstream mesh_name, sol_name;
mesh_name << "mesh." << setfill('0') << setw(6) << myid;
sol_name << "sol." << setfill('0') << setw(6) << myid;
ofstream mesh_ofs(mesh_name.str().c_str());
mesh_ofs.precision(8);
pmesh.Print(mesh_ofs);
ofstream sol_ofs(sol_name.str().c_str());
sol_ofs.precision(8);
x.Save(sol_ofs);
}
// 16. Send the solution by socket to a GLVis server.
if (visualization)
{
char vishost[] = "localhost";
int visport = 19916;
socketstream sol_sock(vishost, visport);
sol_sock << "parallel " << num_procs << " " << myid << "\n";
sol_sock.precision(8);
sol_sock << "solution\n" << pmesh << x << flush;
}
// 17. Free the used memory.
if (delete_fec)
{
delete fec;
}
return 0;
}
+1 -2
View File
@@ -20,9 +20,8 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
MFEM_LIB_FILE = mfem_is_not_built
-include $(CONFIG_MK)
# Currently there are only serial Ginkgo examples
SEQ_EXAMPLES = ex1
PAR_EXAMPLES =
PAR_EXAMPLES = ex1p
ifeq ($(MFEM_USE_MPI),NO)
EXAMPLES = $(SEQ_EXAMPLES)
else
+21 -14
View File
@@ -2137,14 +2137,14 @@ private:
int trial_dofs1D, test_dofs1D, quad1D;
public:
GradientIntegrator() :
Q{NULL}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
GradientIntegrator(const IntegrationRule *ir = NULL) :
BilinearFormIntegrator(ir), Q{NULL}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
{ }
GradientIntegrator(Coefficient *q_) :
Q{q_}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
GradientIntegrator(Coefficient *q_, const IntegrationRule *ir = NULL) :
BilinearFormIntegrator(ir), Q{q_}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
{ }
GradientIntegrator(Coefficient &q) :
Q{&q}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
GradientIntegrator(Coefficient &q, const IntegrationRule *ir = NULL) :
BilinearFormIntegrator(ir), Q{&q}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
{ }
void AssembleElementMatrix2(const FiniteElement &trial_fe,
@@ -2612,6 +2612,11 @@ protected:
Vector pa_data;
public:
VectorMassIntegrator(const IntegrationRule *ir)
: BilinearFormIntegrator(ir), vdim(-1), Q_order(0), Q(NULL), VQ(NULL),
MQ(NULL) { }
/// Construct an integrator with coefficient 1.0
VectorMassIntegrator() = default;
@@ -3038,14 +3043,15 @@ private:
int trial_dofs1D, test_dofs1D, quad1D;
public:
VectorDivergenceIntegrator() :
Q(NULL), trial_maps(NULL), test_maps(NULL), geom(NULL)
VectorDivergenceIntegrator(const IntegrationRule *ir = NULL) :
BilinearFormIntegrator(ir), Q(NULL), trial_maps(NULL), test_maps(NULL),
geom(NULL)
{ }
VectorDivergenceIntegrator(Coefficient *q_) :
Q(q_), trial_maps(NULL), test_maps(NULL), geom(NULL)
VectorDivergenceIntegrator(Coefficient *q_, const IntegrationRule *ir = NULL) :
BilinearFormIntegrator(ir), Q(q_), trial_maps(NULL), test_maps(NULL), geom(NULL)
{ }
VectorDivergenceIntegrator(Coefficient &q) :
Q(&q), trial_maps(NULL), test_maps(NULL), geom(NULL)
VectorDivergenceIntegrator(Coefficient &q, const IntegrationRule *ir = NULL) :
BilinearFormIntegrator(ir), Q(&q), trial_maps(NULL), test_maps(NULL), geom(NULL)
{ }
void AssembleElementMatrix2(const FiniteElement &trial_fe,
@@ -3263,8 +3269,9 @@ private:
void SetUpQuadratureSpaceAndCoefficients(const FiniteElementSpace &fes);
public:
ElasticityIntegrator(Coefficient &l, Coefficient &m)
{ lambda = &l; mu = &m; }
ElasticityIntegrator(Coefficient &l, Coefficient &m,
const IntegrationRule *ir = NULL) :
BilinearFormIntegrator(ir) { lambda = &l; mu = &m; }
/** With this constructor $\lambda = q_l m$ and $\mu = q_m m$
if $dim q_l + 2 q_m = 0$ then $tr(\sigma) = 0$. */
ElasticityIntegrator(Coefficient &m, real_t q_l, real_t q_m)
+1 -1
View File
@@ -10,7 +10,7 @@
// CONTRIBUTING.md for details.
#pragma once
// This is serac's tuple implementation
// This is smith's tuple implementation
#include <ostream>
#include "../../config/config.hpp"
+22
View File
@@ -652,6 +652,16 @@ void forall(func_t f,
class FDJacobian : public Operator
{
public:
FDJacobian(const Operator &op, real_t fixed_eps = 0.0) :
Operator(op.Height(), op.Width()),
op(op),
x(op.Width()),
fixed_eps(fixed_eps)
{
f.SetSize(Height());
xpev.SetSize(Width());
}
FDJacobian(const Operator &op, const Vector &x, real_t fixed_eps = 0.0) :
Operator(op.Height(), op.Width()),
op(op),
@@ -671,6 +681,18 @@ public:
MPI_COMM_WORLD);
}
void Update(const Vector &x_new)
{
x = x_new;
f.UseDevice(x.UseDevice());
xpev.UseDevice(x.UseDevice());
op.Mult(x, f);
const real_t xnorm_local = x.Norml2();
MPI_Allreduce(&xnorm_local, &xnorm, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
MPI_COMM_WORLD);
}
void Mult(const Vector &v, Vector &y) const override
{
// See [1] for choice of eps.
+16 -2
View File
@@ -23,6 +23,8 @@ class BatchedLOR_DG : BatchedLORKernel
{
IntegrationRule ir_face; ///< Collocated Gauss-Lobatto face quadrature rule.
real_t kappa; ///< DG penalty parameter.
bool has_bdr_integ; ///< Is there a boundary integrator?
const Array<int> *bdr_markers; ///< Boundary integrator markers.
public:
template <int ORDER, int SDIM> void Assemble2D();
template <int ORDER> void Assemble3D();
@@ -38,8 +40,7 @@ public:
ProjectLORCoefficient<MassIntegrator>(a, c1);
ProjectLORCoefficient<DiffusionIntegrator>(a, c2);
auto *integ = GetInteriorFaceIntegrator<DGDiffusionIntegrator>(a);
if (integ)
if (auto *integ = GetInteriorFaceIntegrator<DGDiffusionIntegrator>(a))
{
kappa = integ->GetPenaltyParameter();
}
@@ -47,6 +48,19 @@ public:
{
kappa = 0.0;
}
has_bdr_integ = false;
auto *bdr_face_integs = a.GetBFBFI();
for (int i = 0; i < bdr_face_integs->Size(); ++i)
{
if (auto *integ = dynamic_cast<DGDiffusionIntegrator*>((*bdr_face_integs)[i]))
{
kappa = integ->GetPenaltyParameter();
bdr_markers = (*a.GetBFBFI_Marker())[i];
has_bdr_integ = true;
break;
}
}
}
/// @brief Compute and return the face info array.
+18
View File
@@ -22,9 +22,13 @@ namespace mfem
Array<int> BatchedLOR_DG::GetFaceInfo() const
{
Mesh &mesh = *fes_ho.GetMesh();
const Array<int> &bdr_face_attrs = mesh.GetBdrFaceAttributes();
const int nf = mesh.GetNumFaces();
Array<int> face_info(nf * 6); // (e0, f0, o0, e1, f1, o1)
auto h_face_info = Reshape(face_info.HostWrite(), 6, nf);
int bdr_face_counter = 0;
for (int f = 0; f < nf; ++f)
{
auto finfo = mesh.GetFaceInformation(f);
@@ -43,6 +47,19 @@ Array<int> BatchedLOR_DG::GetFaceInfo() const
h_face_info(4, f) = -1;
h_face_info(5, f) = -1;
}
if (finfo.IsBoundary())
{
// Check if Neumann boundary; skip these when adding boundary penalties
const int bdr_attr = bdr_face_attrs[bdr_face_counter];
if (!has_bdr_integ || (bdr_markers && !(*bdr_markers)[bdr_attr - 1]))
{
h_face_info(0, f) = -1;
h_face_info(1, f) = -1;
h_face_info(2, f) = -1;
}
bdr_face_counter += 1;
}
}
return face_info;
}
@@ -144,6 +161,7 @@ void BatchedLOR_DG::AssembleFaceTerms()
{
const int f_0 = d_face_info(1, f);
const int f_1 = d_face_info(4, f);
if (f_0 < 0) { return; } // Skip Neumann boundary faces
const int nsides = (f_1 >= 0) ? 2 : 1;
for (int el_i = 0; el_i < nsides; ++el_i)
{
+13 -2
View File
@@ -758,7 +758,7 @@ void VectorConvectionNLFIntegrator::AssembleElementVector(
EF.UseExternalData(elfun.GetData(), nd, dim);
ELV.UseExternalData(elvect.GetData(), nd, dim);
Vector vec1(dim), vec2(dim);
Vector vec1(dim), vec2(dim), vecvq(dim);
const IntegrationRule *ir = GetIntegrationRule(el, T);
ELV = 0.0;
for (int i = 0; i < ir->GetNPoints(); i++)
@@ -772,6 +772,11 @@ void VectorConvectionNLFIntegrator::AssembleElementVector(
MultAtB(EF, dshape, gradEF);
EF.MultTranspose(shape, vec1);
if (VQ)
{
VQ->Eval(vecvq, T, ip);
vec1 += vecvq;
}
gradEF.Mult(vec1, vec2);
vec2 *= w;
AddMultVWt(shape, vec2, ELV);
@@ -797,7 +802,7 @@ void VectorConvectionNLFIntegrator::AssembleElementGrad(
EF.UseExternalData(elfun.GetData(), nd, dim);
real_t w;
Vector vec1(dim), vec2(dim), vec3(nd);
Vector vec1(dim), vec2(dim), vec3(nd), vecvq(dim);
const IntegrationRule *ir = GetIntegrationRule(el, trans);
@@ -822,6 +827,12 @@ void VectorConvectionNLFIntegrator::AssembleElementGrad(
MultAtB(EF, dshapex, gradEF);
EF.MultTranspose(shape, vec1);
if (VQ)
{
VQ->Eval(vecvq, trans, ip);
vec1 += vecvq;
}
trans.AdjugateJacobian().Mult(vec1, vec2);
vec2 *= w;
+8 -1
View File
@@ -381,6 +381,7 @@ class VectorConvectionNLFIntegrator : public NonlinearFormIntegrator
{
private:
Coefficient *Q{};
VectorCoefficient *VQ{};
DenseMatrix dshape, dshapex, EF, gradEF, ELV, elmat_comp;
Vector shape;
// PA extension
@@ -390,7 +391,13 @@ private:
int dim, ne, nq;
public:
VectorConvectionNLFIntegrator(Coefficient &q): Q(&q) { }
VectorConvectionNLFIntegrator(Coefficient &q,
const IntegrationRule *ir = NULL) :
NonlinearFormIntegrator(ir), Q(&q) { }
VectorConvectionNLFIntegrator(Coefficient &q, VectorCoefficient &vq,
const IntegrationRule *ir = NULL) :
NonlinearFormIntegrator(ir), Q(&q), VQ(&vq) { }
VectorConvectionNLFIntegrator() = default;
+25
View File
@@ -22,6 +22,8 @@
#include <unordered_map>
#include <map>
#include <sstream>
#include <iomanip>
namespace mfem
{
@@ -716,6 +718,29 @@ void Device::DeviceMem(size_t *free, size_t *total)
#endif
}
std::string Device::GetUUID(const int device_id)
{
std::stringstream res;
#if defined(MFEM_USE_CUDA)
cudaDeviceProp prop;
MFEM_GPU_CHECK(cudaGetDeviceProperties(&prop, device_id));
for (int i = 0; i < 16; ++i)
{
res << std::setfill('0') << std::setw(2) << std::hex
<< static_cast<unsigned>(prop.uuid.bytes[i]);
}
#elif defined(MFEM_USE_HIP)
hipUUID uuid;
MFEM_GPU_CHECK(hipDeviceGetUuid(&uuid, device_id));
for (int i = 0; i < 16; ++i)
{
res << std::setfill('0') << std::setw(2) << std::hex
<< static_cast<unsigned>(uuid.bytes[i]);
}
#endif
return res.str();
}
int Device::NumMultiprocessors(int dev)
{
#if defined(MFEM_USE_CUDA)
+4
View File
@@ -255,6 +255,10 @@ public:
/// Get the number of available devices (may be called before configuration).
static int GetDeviceCount();
/// Gets a string representation of the GPU UUID.
/// 0 <= @a device_id < GetDeviceCount()
static std::string GetUUID(const int device_id = 0);
/** @brief Return true if any of the backends in the backend mask, @a b_mask,
are allowed. */
/** This method can be used with any of the Backend::Id constants, the
+1 -1
View File
@@ -146,7 +146,7 @@ public:
int *ReadWriteJ(bool on_dev = true) { return J.ReadWrite(on_dev); }
const int *HostReadJ() const { return J.HostRead(); }
int *HostWriteJ() { return J.HostWrite(); }
int *ReadWriteJ() { return J.HostReadWrite(); }
int *HostReadWriteJ() { return J.HostReadWrite(); }
/// Sort the column (TYPE II) indices in each row.
void SortRows();
+823 -146
View File
File diff suppressed because it is too large Load Diff
+755 -80
View File
File diff suppressed because it is too large Load Diff
+7
View File
@@ -1681,6 +1681,13 @@ void HypreParMatrix::GetOffd(SparseMatrix &offd, HYPRE_BigInt* &cmap) const
cmap = A->col_map_offd;
}
void HypreParMatrix::GetOffdColMap(HYPRE_BigInt* &cmap,
HYPRE_Int &num_cols) const
{
cmap = A->col_map_offd;
num_cols = hypre_CSRMatrixNumCols(A->offd);
}
void HypreParMatrix::MergeDiagAndOffd(SparseMatrix &merged)
{
HostRead();
+10
View File
@@ -665,6 +665,8 @@ public:
void GetDiag(SparseMatrix &diag) const;
/// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
void GetOffd(SparseMatrix &offd, HYPRE_BigInt* &cmap) const;
/// Get the global column mapping for the local off-diagonal block.
void GetOffdColMap(HYPRE_BigInt* &cmap, HYPRE_Int &num_cols) const;
/** @brief Get a single SparseMatrix containing all rows from this processor,
merged from the diagonal and off-diagonal blocks stored by the
HypreParMatrix. */
@@ -959,6 +961,14 @@ public:
const Memory<HYPRE_Int> &GetDiagMemoryJ() const { return mem_diag.J; }
const Memory<real_t> &GetDiagMemoryData() const { return mem_diag.data; }
Memory<HYPRE_Int> &GetOffdMemoryI() { return mem_offd.I; }
Memory<HYPRE_Int> &GetOffdMemoryJ() { return mem_offd.J; }
Memory<real_t> &GetOffdMemoryData() { return mem_offd.data; }
const Memory<HYPRE_Int> &GetOffdMemoryI() const { return mem_offd.I; }
const Memory<HYPRE_Int> &GetOffdMemoryJ() const { return mem_offd.J; }
const Memory<real_t> &GetOffdMemoryData() const { return mem_offd.data; }
/// @brief Prints the locally owned rows in parallel. The resulting files can
/// be read with Read_IJMatrix().
void Print(const std::string &fname, HYPRE_Int offi = 0,
+101 -6
View File
@@ -10,6 +10,7 @@
// CONTRIBUTING.md for details.
#include "../general/communication.hpp"
#include "../general/forall.hpp"
#include "operator.hpp"
#include "ode.hpp"
@@ -184,6 +185,23 @@ void ODESolver::Init(TimeDependentOperator &f_)
mem_type = GetMemoryType(f_.GetMemoryClass());
}
void ODESolver::ComputeSlopeFromState(const real_t dt, const Vector &u,
Vector &k)
{
// k currently holds state u(t+dt),
// convert to slope k = du/dt ~= (u(t+dt)-u(t))/dt
const int usz = u.Size();
real_t fac = 1.0/dt;
auto d_u = u.Read();
auto d_k = k.ReadWrite();
mfem::forall(usz, [=] MFEM_HOST_DEVICE (int i)
{
d_k[i] -= d_u[i];
d_k[i] *= fac;
});
}
void ForwardEulerSolver::Init(TimeDependentOperator &f_)
{
ODESolver::Init(f_);
@@ -629,6 +647,10 @@ void AdamsMoultonSolver::Step(Vector &x, real_t &t, real_t &dt)
}
state.ShiftStages();
f->ImplicitSolve(a[0]*dt, x, state[0]);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a[0]*dt, x, state[0]);
}
x.Add(a[0]*dt, state[0]);
t += dt;
}
@@ -661,7 +683,15 @@ void BackwardEulerSolver::Step(Vector &x, real_t &t, real_t &dt)
{
f->SetTime(t + dt);
f->ImplicitSolve(dt, x, k); // solve for k: k = f(x + dt*k, t + dt)
x.Add(dt, k);
if (f->ImplicitVarTypeIsState())
{
x = k; // x = u_{i+1}
}
else
{
x.Add(dt, k);
}
t += dt;
}
@@ -676,7 +706,16 @@ void ImplicitMidpointSolver::Step(Vector &x, real_t &t, real_t &dt)
{
f->SetTime(t + dt/2);
f->ImplicitSolve(dt/2, x, k);
x.Add(dt, k);
if (f->ImplicitVarTypeIsState())
{
x.Neg();
x.Add(2.0, k);
}
else
{
x.Add(dt, k);
}
t += dt;
}
@@ -718,11 +757,19 @@ void SDIRK23Solver::Step(Vector &x, real_t &t, real_t &dt)
// note: with gamma_opt=3, both solve are outside [t,t+dt] since a>1
f->SetTime(t + gamma*dt);
f->ImplicitSolve(gamma*dt, x, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(gamma*dt, x, k);
}
add(x, (1.-2.*gamma)*dt, k, y); // y = x + (1-2*gamma)*dt*k
x.Add(dt/2, k);
f->SetTime(t + (1.-gamma)*dt);
f->ImplicitSolve(gamma*dt, y, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(gamma*dt, y, k);
}
x.Add(dt/2, k);
t += dt;
}
@@ -749,17 +796,29 @@ void SDIRK34Solver::Step(Vector &x, real_t &t, real_t &dt)
f->SetTime(t + a*dt);
f->ImplicitSolve(a*dt, x, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, x, k);
}
add(x, (0.5-a)*dt, k, y);
add(x, (2.*a)*dt, k, z);
x.Add(b*dt, k);
f->SetTime(t + dt/2);
f->ImplicitSolve(a*dt, y, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, y, k);
}
z.Add((1.-4.*a)*dt, k);
x.Add((1.-2.*b)*dt, k);
f->SetTime(t + (1.-a)*dt);
f->ImplicitSolve(a*dt, z, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, z, k);
}
x.Add(b*dt, k);
t += dt;
}
@@ -785,15 +844,27 @@ void SDIRK33Solver::Step(Vector &x, real_t &t, real_t &dt)
f->SetTime(t + a*dt);
f->ImplicitSolve(a*dt, x, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, x, k);
}
add(x, (c-a)*dt, k, y);
x.Add(b*dt, k);
f->SetTime(t + c*dt);
f->ImplicitSolve(a*dt, y, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, y, k);
}
x.Add((1.0-a-b)*dt, k);
f->SetTime(t + dt);
f->ImplicitSolve(a*dt, x, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, x, k);
}
x.Add(a*dt, k);
t += dt;
}
@@ -818,6 +889,10 @@ void TrapezoidalRuleSolver::Step(Vector &x, real_t &t, real_t &dt)
f->SetTime(t + dt);
f->ImplicitSolve(dt/2.0, y, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(0.5*dt, y, k);
}
x.Add(dt/2.0, k);
t += dt;
}
@@ -848,11 +923,19 @@ void ESDIRK32Solver::Step(Vector &x, real_t &t, real_t &dt)
f->SetTime(t + (2.0*a)*dt);
f->ImplicitSolve(a*dt, y, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, y, k);
}
z.Add(b*dt, k);
x.Add(b*dt, k);
f->SetTime(t + dt);
f->ImplicitSolve(a*dt, z, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, z, k);
}
x.Add(a*dt, k);
t += dt;
}
@@ -885,11 +968,19 @@ void ESDIRK33Solver::Step(Vector &x, real_t &t, real_t &dt)
f->SetTime(t + (2.0*a)*dt);
f->ImplicitSolve(a*dt, y, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, y, k);
}
z.Add(b*dt, k);
x.Add(b_2*dt, k);
f->SetTime(t + dt);
f->ImplicitSolve(a*dt, z, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(a*dt, z, k);
}
x.Add(b_3*dt, k);
t += dt;
}
@@ -955,6 +1046,10 @@ void GeneralizedAlphaSolver::Step(Vector &x, real_t &t, real_t &dt)
real_t dt_eff = (gamma*alpha_f/alpha_m)*dt;
f->SetTime(t + alpha_f*dt);
f->ImplicitSolve(dt_eff, y, k);
if (f->ImplicitVarTypeIsState())
{
ComputeSlopeFromState(dt_eff, y, k);
}
// Update x and xdot
x.Add((1.0 - (gamma/alpha_m))*dt, state[0]);
@@ -1116,8 +1211,8 @@ void SecondOrderODESolver::EulerStep(Vector &x, Vector &dxdt, real_t &t,
f->SetTime(t + dt);
f->ImplicitSolve(0.5*dt*dt, dt, x, dxdt, state[0]);
x .Add(0.5*dt*dt, state[0]);
dxdt.Add(dt, state[0]);
x.Add(0.5*dt*dt, state[0]);
dxdt.Add(dt, state[0]);
t += dt;
}
@@ -1203,8 +1298,8 @@ void NewmarkSolver::Step(Vector &x, Vector &dxdt, real_t &t, real_t &dt)
f->SetTime(t + dt);
f->ImplicitSolve(fac3*dt*dt, fac4*dt, x, dxdt, state[0]);
x .Add(fac3*dt*dt, state[0]);
dxdt.Add(fac4*dt, state[0]);
x.Add(fac3*dt*dt, state[0]);
dxdt.Add(fac4*dt, state[0]);
t += dt;
}
+79
View File
@@ -120,6 +120,7 @@ public:
class ODESolver
{
protected:
using ImplicitVariableType = TimeDependentOperator::ImplicitVariableType;
/// Pointer to the associated TimeDependentOperator.
TimeDependentOperator *f; // f(.,t) : R^n --> R^n
MemoryType mem_type;
@@ -192,6 +193,25 @@ public:
/// Returns how many State vectors the ODE requires
virtual int GetStateSize() { return 0; };
/// Returns the associated TimeDependentOperator
TimeDependentOperator* GetTimeDependentOperator() { return f; }
///@brief Returns @a true if the ODESolver supports the given
/// #ImplicitVariableType, @a var, and returns @a false otherwise.
///@note Should be overriden in ODESolver that calls TimeDependentOperator::ImplicitSolve().
virtual bool SupportsImplicitVariableType(ImplicitVariableType var) const
{ return false; };
/** @brief Compute the finite-difference slope, @a $\frac{du}{dt} \approx \frac{u(t+dt)-u(t)}{dt}$,
* and store it in @a k.
* @param [in] dt Finite difference step size.
* @param [in] u state vector, @a u(t).
* @param [in,out] k On input, @a k contains the state vector, @a u( @a t+ @a dt).
* On output, @a k contains the computed slope, @a du/dt.
* */
virtual void ComputeSlopeFromState(const real_t dt, const Vector &u,
Vector &k);
// Help info for ODESolver options
static MFEM_EXPORT std::string ExplicitTypes;
static MFEM_EXPORT std::string ImplicitTypes;
@@ -361,6 +381,12 @@ public:
void Init(TimeDependentOperator &f_) override;
void Step(Vector &x, real_t &t, real_t &dt) override;
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
@@ -374,6 +400,12 @@ public:
void Init(TimeDependentOperator &f_) override;
void Step(Vector &x, real_t &t, real_t &dt) override;
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
@@ -395,6 +427,12 @@ public:
void Init(TimeDependentOperator &f_) override;
void Step(Vector &x, real_t &t, real_t &dt) override;
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
@@ -409,6 +447,12 @@ public:
void Init(TimeDependentOperator &f_) override;
void Step(Vector &x, real_t &t, real_t &dt) override;
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
@@ -423,6 +467,12 @@ public:
void Init(TimeDependentOperator &f_) override;
void Step(Vector &x, real_t &t, real_t &dt) override;
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
@@ -437,6 +487,12 @@ public:
void Init(TimeDependentOperator &f_) override;
void Step(Vector &x, real_t &t, real_t &dt) override;
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
@@ -451,6 +507,12 @@ public:
void Init(TimeDependentOperator &f_) override;
void Step(Vector &x, real_t &t, real_t &dt) override;
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
@@ -465,6 +527,12 @@ public:
void Init(TimeDependentOperator &f_) override;
void Step(Vector &x, real_t &t, real_t &dt) override;
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
@@ -490,6 +558,12 @@ public:
ODEStateData& GetState() override { return state; }
const ODEStateData& GetState() const override { return state; }
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
@@ -606,6 +680,11 @@ public:
ODEStateData& GetState() override { return state; }
const ODEStateData& GetState() const override { return state; }
bool SupportsImplicitVariableType(ImplicitVariableType var) const override
{
return (var == ImplicitVariableType::STATE ||
var == ImplicitVariableType::SLOPE);
}
};
/** A 1-stage, 2nd order AM method. */
+37 -3
View File
@@ -381,11 +381,24 @@ public:
ADDITIVE_TERM_2
};
/** Used to specify the variable being returned by ImplicitSolve(). This can
* be queried by ODESolver to identify the variable being solved for.
* @warning Not all ODESolver may support all options. See ODESolver::SupportsImplicitVariableType() */
enum ImplicitVariableType
{
SLOPE, ///< stage slope, $k = \frac{du}{dt}$.
STATE ///< stage state, $k = u$.
};
protected:
real_t t; ///< Current time.
Type type; /**< @brief Describes the form of the TimeDependentOperator, see
the documentation of #Type. */
EvalMode eval_mode; ///< Current evaluation mode.
ImplicitVariableType implicit_variable_type =
ImplicitVariableType::SLOPE; /**< @brief
Return variable for
ImplicitSolve()*/
public:
/** @brief Construct a "square" TimeDependentOperator (u,t) -> k(u,t), where
@@ -429,6 +442,24 @@ public:
virtual void SetEvalMode(const EvalMode new_eval_mode)
{ eval_mode = new_eval_mode; }
/** @brief Sets the #ImplicitVariableType for ImplicitSolve()*/
virtual void SetImplicitVariableType(const ImplicitVariableType variable_type)
{ implicit_variable_type = variable_type; }
/** @brief Returns the #ImplicitVariableType for ImplicitSolve(). */
virtual ImplicitVariableType GetImplicitVariableType() const
{ return implicit_variable_type; }
/** @brief Returns @a true if implicit variable is #STATE and @a false otherwise.
* Used by ODESolver to identify the stage variable returned by ImplicitSolve() */
virtual bool ImplicitVarTypeIsState() const
{ return (implicit_variable_type == ImplicitVariableType::STATE); }
/** @brief Returns @a true if implicit variable is #SLOPE and @a false otherwise.
* Used by ODESolver to identify the stage variable returned by ImplicitSolve() */
virtual bool ImplicitVarTypeIsSlope() const
{ return (implicit_variable_type == ImplicitVariableType::SLOPE); }
/** @brief Perform the action of the explicit part of the operator, G:
@a v = G(@a u, t) where t is the current time.
@@ -462,7 +493,8 @@ public:
/** @brief Solve for the unknown @a k, at the current time t, the following
equation:
F(@a u + @a gamma @a k, @a k, t) = G(@a u + @a gamma @a k, t).
1. $F( u + \gamma k, k, t) = G( u + \gamma k, t)$, if solving for stage-slope (default)
2. $F( u , \frac{k-u}{\gamma}, t) = G(k, t)$, if solving for stage-state
For solving an ordinary differential equation of the form
$ M \frac{dy}{dt} = g(y,t) $, recall that F and G can be defined in
@@ -472,8 +504,9 @@ public:
2. F(u,k,t) = M k and G(u,t) = g(u,t)
3. F(u,k,t) = M k - g(u,t) and G(u,t) = 0
Regardless of the choice of F and G, this function should solve for @a k
in M @a k = g(@a u + @a gamma @a k, t).
Regardless of the choice of F and G, this function should solve for @a k:
- $~Mk = g( u + \gamma k, t)~$, if solving for stage-slope.
- $~Mk = \gamma g(k, t) + Mu~$, if solving for stage-state
To see how @a k can be useful, consider the backward Euler method defined
by $ y(t + \Delta t) = y(t) + \Delta t k_0 $ where
@@ -491,6 +524,7 @@ public:
$ y(t) + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j $ and @a gamma set to
$ a_{ii} \Delta t $, for $ k_i $. For example, see class SDIRK33Solver.
See SetImplicitVariableType() to switch between different variable modes.
If not re-implemented, this method simply generates an error. */
virtual void ImplicitSolve(const real_t gamma, const Vector &u, Vector &k);
+210
View File
@@ -851,6 +851,216 @@ void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
sli.Mult(b, x);
}
real_t FPIRelaxation::Dot(const Vector &x, const Vector &y) const
{
if (dot_oper) { return dot_oper->Eval(x,y); } // Use custom inner product (if provided)
#ifndef MFEM_USE_MPI
return (x * y);
#else
return InnerProduct(comm, x, y);
#endif
}
real_t FPIRelaxation::SquaredDistance(const Vector &x, const Vector &y) const
{
#ifndef MFEM_USE_MPI
return x.DistanceSquaredTo(y);
#else
return DistanceSquared(comm, x, y);
#endif
}
void AitkenRelaxation::Init()
{
FPIRelaxation::Init();
MemoryType mt = GetMemoryType(oper->GetMemoryClass());
rold.SetSize(oper->Width(), mt);
rold.UseDevice(true);
rold = 0.0;
rold_nsq = 0.0;
}
real_t AitkenRelaxation::Eval(const Vector &state, const Vector &residual,
real_t res_norm, real_t rfactor)
{
real_t rdot = Dot(rold, residual);
real_t rnsq = res_norm * res_norm;
real_t denom = rnsq - 2.0*rdot + rold_nsq; // ||rold - residual||^2
real_t ratio = (rdot - rold_nsq) / denom;
rold_nsq = rnsq;
rold = residual;
if (ratio == 0.0) { return rfactor; } // Avoid factor = 0.0 at first call
return Clamp(-rfactor * ratio);
}
void SteepestDescentRelaxation::Init()
{
FPIRelaxation::Init();
MemoryType mt = GetMemoryType(oper->GetMemoryClass());
z.SetSize(oper->Width(), mt);
z.UseDevice(true);
}
real_t SteepestDescentRelaxation::Eval(const Vector &state,
const Vector &residual,
real_t res_norm, real_t rfactor)
{
MFEM_VERIFY(oper,"Operator not set; set using SetOperator(Operator&)")
Operator *J = &oper->GetGradient(state);
real_t num = res_norm * res_norm;
J->Mult(residual, z); // rold = F'(x) * rnew;
real_t denom = Dot(z, residual);
return Clamp(num/denom);
}
/// Fixed point iteration solver: x <- f(x)
void FPISolver::UpdateVectors()
{
MemoryType mt = GetMemoryType(oper->GetMemoryClass());
r.SetSize(width, mt);
r.UseDevice(true);
z.SetSize(width, mt);
z.UseDevice(true);
}
void FPISolver::SetOperator(const Operator &op)
{
IterativeSolver::SetOperator(op);
UpdateVectors();
if (!relax_method)
{
relax_method = new FPIRelaxation(); // Default relaxation strategy
relax_owned = true;
}
relax_method->SetOperator(*oper);
}
void FPISolver::SetRelaxation(real_t rfactor, FPIRelaxation *relaxation,
bool own)
{
relax_factor = rfactor;
if (relaxation)
{
if (relax_method && relax_owned) { delete relax_method; }
relax_method = relaxation;
relax_owned = own;
}
#ifdef MFEM_USE_MPI
if (relax_method) { relax_method->SetComm(this->GetComm()); }
#endif
}
/// Iterative solution of the (non)linear system using Fixed Point Iteration
void FPISolver::Mult(const Vector &b, Vector &x) const
{
int i;
real_t factor = relax_factor;
real_t r0, nom, nom0, nomold = 1, cf, fac_old = factor;
if (iterative_mode)
{
oper->Mult(x, z); // z = F(x)
subtract(z, x, r); // r = z - x
}
else
{
x = 0.0;
oper->Mult(x, r); // r = F(x)
}
nom0 = nom = sqrt(Dot(r, r));
initial_norm = nom0;
if (print_options.iterations | print_options.first_and_last)
{
mfem::out << " Iteration : " << setw(3) << right << 0 << " ||Br|| = "
<< nom << (print_options.first_and_last ? " ..." : "") << '\n';
}
r0 = std::max(nom*rel_tol, abs_tol);
if (nom <= r0)
{
converged = true;
final_iter = 0;
final_norm = nom;
return;
}
// start iteration
converged = false;
final_iter = max_iter;
relax_method->Init();
for (i = 1; true; )
{
factor = relax_method->Eval(x,r,nom,fac_old);
x.Add(factor, r); // x = x + factor * r
oper->Mult(x, z); // z = F(x)
subtract(z, x, r); // r = z - x
nom = sqrt(Dot(r, r));
cf = nom/nomold;
nomold = nom;
fac_old = factor;
bool done = false;
if (nom < r0)
{
converged = true;
final_iter = i;
done = true;
}
if (++i > max_iter)
{
done = true;
}
if (print_options.iterations || (done && print_options.first_and_last))
{
mfem::out << " Iteration : " << setw(3) << right << (i-1)
<< " ||r|| = " << setw(11) << left << nom
<< "\trlx. fac.: " << fac_old << '\n';
}
if (done) { break; }
}
if (print_options.summary || (print_options.warnings && !converged))
{
const auto rf = pow (nom/nom0, 1.0/final_iter);
mfem::out << "FPI: Number of iterations: " << final_iter << '\n'
<< "Conv. rate: " << cf << '\n'
<< "Average reduction factor: "<< rf << '\n';
}
if (print_options.warnings && !converged)
{
mfem::out << "FPI: No convergence!" << '\n';
}
final_norm = nom;
}
void FPI(const Operator &A, const Vector &b, Vector &x,
int print_iter, int max_num_iter,
real_t RTOLERANCE, real_t ATOLERANCE,
real_t relax_factor, FPIRelaxation *relax_method)
{
MFEM_PERF_FUNCTION;
FPISolver fpi;
fpi.SetPrintLevel(print_iter);
fpi.SetMaxIter(max_num_iter);
fpi.SetRelTol(sqrt(RTOLERANCE));
fpi.SetAbsTol(sqrt(ATOLERANCE));
fpi.SetOperator(A);
fpi.SetRelaxation(relax_factor,relax_method,false);
fpi.Mult(b, x);
}
void CGSolver::UpdateVectors()
{
+178
View File
@@ -622,6 +622,184 @@ void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
real_t RTOLERANCE = 1e-12, real_t ATOLERANCE = 1e-24);
/**
@brief A class to handle fixed point iteration relaxation methods.
This class provides a base for implementing various relaxation strategies
for fixed point iteration solvers.
*/
class FPIRelaxation
{
#ifdef MFEM_USE_MPI
private:
MPI_Comm comm = MPI_COMM_NULL;
#endif
protected:
const Operator *oper = nullptr;
real_t lbnd = std::numeric_limits<real_t>::lowest();
real_t ubnd = std::numeric_limits<real_t>::max();
real_t abs_lbnd = 0.0;
InnerProductOperator *dot_oper = nullptr;
public:
FPIRelaxation() = default;
#ifdef MFEM_USE_MPI
FPIRelaxation(MPI_Comm comm_) {comm = comm_;}
void SetComm(MPI_Comm comm_) {comm = comm_;}
#endif
/// @brief Set the operator for the relaxation method.
virtual void SetOperator(const Operator &op) {oper = &op;}
/// @brief Initialize the relaxation method.
virtual void Init() {}
/// @brief Set the lower bound for the relaxation factor.
virtual void SetLowerBound(real_t lb) { lbnd = lb; }
/// @brief Set the upper bound for the relaxation factor.
virtual void SetUpperBound(real_t ub) { ubnd = ub; }
/// @brief Set both lower and upper bounds for the relaxation factor.
virtual void SetBounds(real_t lb, real_t ub)
{ lbnd = lb; ubnd = ub;}
/// @brief Set the absolute bound for the relaxation factor.
virtual void SetAbsoluteLowerBound(real_t abs_lb)
{ abs_lbnd = std::abs(abs_lb); }
/// @brief Set a user-defined inner product operator (not owned)
void SetInnerProduct(InnerProductOperator *ipo) { dot_oper = ipo; }
/// @brief Clamp the relaxation factor to the specified range.
virtual real_t Clamp(real_t factor) const
{
real_t afac = (abs_lbnd == 0.0) ? factor :
std::copysign(std::max(std::abs(factor), abs_lbnd), factor);
return std::max(std::min(afac, ubnd), lbnd);
}
/**
@brief Compute the relaxation factor for Fixed Point Iteration.
@param state Current state vector
@param residual Current residual vector
@param res_norm Norm of the current residual vector
@param rfactor Current relaxation factor
@return real_t The computed relaxation factor
*/
virtual real_t Eval(const Vector &state, const Vector &residual,
real_t res_norm, real_t rfactor)
{
return rfactor; // Default implementation returns the fixed factor
};
real_t Dot(const Vector &x, const Vector &y) const;
real_t SquaredDistance(const Vector &x, const Vector &y) const;
virtual ~FPIRelaxation() {}
};
/**
@brief A class to implement Aitken relaxation for Fixed Point Iteration
to accelerate the convergence of fixed point iterations.
*/
class AitkenRelaxation : public FPIRelaxation
{
protected:
Vector rold; // Old residual vector
real_t rold_nsq = 0.0;
public:
AitkenRelaxation() = default;
#ifdef MFEM_USE_MPI
AitkenRelaxation(MPI_Comm comm_) : FPIRelaxation(comm_) {}
#endif
void Init() override;
/**
@brief Compute the Aitken relaxation factor for Fixed Point Iteration.
@param state Current state vector
@param residual Current residual vector
@param res_norm Norm of the current residual vector
@param rfactor Current relaxation factor
@return The computed relaxation factor
*/
real_t Eval(const Vector &state, const Vector &residual, real_t res_norm,
real_t rfactor) override;
};
/**
@brief A class to implement steepest descent relaxation for Fixed Point Iteration
to accelerate the convergence of fixed point iterations.
*/
class SteepestDescentRelaxation : public FPIRelaxation
{
protected:
Vector z;
public:
SteepestDescentRelaxation() = default;
#ifdef MFEM_USE_MPI
SteepestDescentRelaxation(MPI_Comm comm_) : FPIRelaxation(comm_) {}
#endif
void Init() override;
/**
@brief Compute the steepest descent relaxation factor for Fixed Point Iteration.
@param state Current state vector
@param residual Current residual vector
@param res_norm Norm of the current residual vector
@param rfactor Current relaxation factor
@return The computed relaxation factor
*/
real_t Eval(const Vector &state, const Vector &residual, real_t res_norm,
real_t rfactor) override;
};
/// Fixed point iteration solver: x <- f(x)
class FPISolver : public IterativeSolver
{
protected:
FPIRelaxation *relax_method = nullptr; ///< Relaxation strategy for FPI
bool relax_owned = false;
mutable real_t relax_factor = 1.0;
mutable Vector r, z;
void UpdateVectors();
public:
FPISolver() { }
#ifdef MFEM_USE_MPI
FPISolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
#endif
virtual void SetOperator(const Operator &op) override;
void SetRelaxation(real_t rfactor, FPIRelaxation *relaxation = nullptr,
bool own = false);
/// Iterative solution of the (non)linear system using Fixed Point Iteration
/// b is not used
void Mult(const Vector &b, Vector &x) const override;
};
/// Fixed point iteration. (tolerances are squared)
void FPI(const Operator &A, const Vector &b, Vector &x,
int print_iter = 0, int max_num_iter = 1000,
real_t RTOLERANCE = 1e-12, real_t ATOLERANCE = 1e-24,
real_t relax_factor = 1.0, FPIRelaxation *relax_method = nullptr);
/// Conjugate gradient method
class CGSolver : public IterativeSolver
{
+65 -18
View File
@@ -23,15 +23,31 @@ namespace mfem
void SparseSmoother::SetOperator(const Operator &a)
{
oper = dynamic_cast<const SparseMatrix*>(&a);
if (oper == NULL)
{
mfem_error("SparseSmoother::SetOperator : not a SparseMatrix!");
}
MFEM_VERIFY(oper != nullptr, "Operator must be a SparseMatrix");
height = oper->Height();
width = oper->Width();
At.reset();
oper_T = nullptr;
}
void SparseSmoother::EnsureTranspose() const
{
if (oper_T) { return; }
const real_t tol = 1e-14;
if (oper->IsSymmetric() > tol * oper->MaxNorm())
{
At.reset(Transpose(*oper));
oper_T = At.get();
}
else
{
At.reset();
oper_T = oper;
}
}
/// Matrix vector multiplication with GS Smoother.
void GSSmoother::Mult(const Vector &x, Vector &y) const
{
if (!iterative_mode)
@@ -51,21 +67,33 @@ void GSSmoother::Mult(const Vector &x, Vector &y) const
}
}
/// Create the Jacobi smoother.
DSmoother::DSmoother(const SparseMatrix &a, int t, real_t s, int it)
: SparseSmoother(a)
void GSSmoother::MultTranspose(const Vector &x, Vector &y) const
{
type = t;
scale = s;
iterations = it;
EnsureTranspose();
if (!iterative_mode)
{
y = 0.0;
}
for (int i = 0; i < iterations; i++)
{
if (type != 1)
{
oper_T->Gauss_Seidel_forw(x, y);
}
if (type != 2)
{
oper_T->Gauss_Seidel_back(x, y);
}
}
}
/// Matrix vector multiplication with Jacobi smoother.
void DSmoother::Mult(const Vector &x, Vector &y) const
void DSmoother::Mult_(const SparseMatrix &A, const Vector &x, Vector &y) const
{
if (!iterative_mode && type == 0 && iterations == 1)
{
oper->DiagScale(x, y, scale, use_abs_diag);
A.DiagScale(x, y, scale, use_abs_diag);
return;
}
@@ -90,22 +118,41 @@ void DSmoother::Mult(const Vector &x, Vector &y) const
{
if (type == 0)
{
oper->Jacobi(x, *p, *r, scale, use_abs_diag);
A.Jacobi(x, *p, *r, scale, use_abs_diag);
}
else if (type == 1)
{
oper->Jacobi2(x, *p, *r, scale);
A.Jacobi2(x, *p, *r, scale);
}
else if (type == 2)
{
oper->Jacobi3(x, *p, *r, scale);
A.Jacobi3(x, *p, *r, scale);
}
else
{
mfem_error("DSmoother::Mult wrong type");
MFEM_ABORT("Invalid type.");
}
Swap<Vector*>(r, p);
}
}
void DSmoother::Mult(const Vector &x, Vector &y) const
{
Mult_(*oper, x, y);
}
void DSmoother::MultTranspose(const Vector &x, Vector &y) const
{
if (iterations == 1 && !iterative_mode)
{
Mult_(*oper, x, y);
return;
}
EnsureTranspose();
MFEM_VERIFY(type == 0 || !At, "l1 or lumped Jacobi transpose not implemented"
" for non-symmetric matrices");
Mult_(*oper_T, x, y);
}
}
+117 -25
View File
@@ -15,67 +15,159 @@
#include "../config/config.hpp"
#include "sparsemat.hpp"
#include <memory>
namespace mfem
{
/// Abstract base class for smoothers created from a SparseMatrix.
class SparseSmoother : public MatrixInverse
{
protected:
const SparseMatrix *oper;
const SparseMatrix *oper = nullptr; ///< The underlying matrix.
/// Pointer to the transpose of the underlying matrix. If the matrix is
/// symmetric, this will be the same as @a oper. If the matrix is not
/// symmetric, the transpose will be formed and stored in @a At. The
/// transpose will only be formed if MultTranspose() is called.
mutable const SparseMatrix *oper_T = nullptr;
mutable std::unique_ptr<SparseMatrix> At; ///< Transpose of A, if needed.
void EnsureTranspose() const; ///< Ensure that the transpose is set.
public:
SparseSmoother() { oper = NULL; }
SparseSmoother() = default;
SparseSmoother(const SparseMatrix &a)
: MatrixInverse(a) { oper = &a; }
SparseSmoother(const SparseMatrix &a) { SetOperator(a); }
/// Sets the underlying matrix. @a a must be a SparseMatrix.
void SetOperator(const Operator &a) override;
};
/// Data type for Gauss-Seidel smoother of sparse matrix
/// Gauss-Seidel smoother of a sparse matrix.
class GSSmoother : public SparseSmoother
{
public:
enum GSType
{
SYMMETRIC, ///< Forward Gauss-Seidel, then backward.
FORWARD, ///< Forward Gauss-Seidel ($L^{-1}$).
BACKWARD ///< Backward Gauss-Seidel ($U^{-1}$).
};
protected:
int type; // 0, 1, 2 - symmetric, forward, backward
int iterations;
GSType type; ///< Type of Gauss-Seidel, see GSSmoother::GSType.
int iterations; ///< Number of stationary iterations.
public:
/// Create GSSmoother.
GSSmoother(int t = 0, int it = 1) { type = t; iterations = it; }
/// @brief Create a Gauss-Seidel smoother. SetOperator() will need to be
/// called with a SparseMatrix before first use.
///
/// @param[in] t Type of GS smoother (see GSSmoother::GSType)
/// @param[in] it Number of stationary iterations to perform
GSSmoother(GSType t = SYMMETRIC, int it = 1) { type = t; iterations = it; }
/// Create GSSmoother.
GSSmoother(const SparseMatrix &a, int t = 0, int it = 1)
: SparseSmoother(a) { type = t; iterations = it; }
/// @brief Create a Jacobi smoother using the SparseMatrix @a a.
///
/// @param[in] a The underlying SparseMatrix
/// @param[in] t Type of GS smoother (see GSSmoother::GSType)
/// @param[in] it Number of stationary iterations to perform
GSSmoother(const SparseMatrix &a, GSType t = SYMMETRIC, int it = 1)
: GSSmoother(t, it) { SetOperator(a); }
/// Matrix vector multiplication with GS Smoother.
/// Same as GSSmoother(GSType,int), for backwards compatibility.
GSSmoother(int t, int it = 1) : GSSmoother(GSType(t), it) { }
/// @brief Same as GSSmoother(const SparseMatrix&,GSType,int), for
/// backwards compatibility.
GSSmoother(const SparseMatrix &a, int t, int it = 1)
: GSSmoother(a, GSType(t), it) { }
/// @brief Application of the Gauss-Seidel smoother.
///
/// Applies a stationary Gauss-Seidel iteration. If Solver::iterative_mode is
/// true, then @a y is used as the initial guess, and Gauss-Seidel is applied
/// to the residual $x - Ay$.
void Mult(const Vector &x, Vector &y) const override;
/// Application of the transpose of the Gauss-Seidel smoother.
void MultTranspose(const Vector &x, Vector &y) const override;
};
/// Data type for scaled Jacobi-type smoother of sparse matrix
/// Jacobi-type diagonal smoother of a sparse matrix.
class DSmoother : public SparseSmoother
{
public:
enum JacobiType
{
JACOBI, ///< Scale by the diagonal of the matrix.
L1_JACOBI, ///< Scale by the l1-norm of the rows.
LUMPED_JACOBI ///< Scale by the sum of the rows.
};
protected:
int type; // 0, 1, 2 - scaled Jacobi, scaled l1-Jacobi, scaled lumped-Jacobi
real_t scale;
int iterations;
/// Uses abs values of the diagonal entries. Relevant only when type = 0.
JacobiType type; ///< Type of diagonal scaling, see DSmoother::JacobiType.
real_t scale; ///< Scaling (damping) factor.
int iterations; ///< Number of stationary iterations to perform.
/// @brief Uses abs values of the diagonal entries. Relevant only with type
/// JacobiType::JACOBI.
bool use_abs_diag = false;
mutable Vector z;
mutable Vector z; ///< Temporary work vector.
/// Apply the Jacobi smoother (used internally by Mult() and MultTranspose())
void Mult_(const SparseMatrix &A, const Vector &x, Vector &y) const;
public:
/// Create Jacobi smoother.
DSmoother(int t = 0, real_t s = 1., int it = 1)
/// @brief Create a Jacobi smoother. SetOperator() will need to be called
/// with a SparseMatrix before first use.
///
/// @param[in] t Type of Jacobi smoother (see DSmoother::JacobiType)
/// @param[in] s Scaling factor
/// @param[in] it Number of stationary iterations to perform
DSmoother(JacobiType t = JACOBI, real_t s = 1., int it = 1)
{ type = t; scale = s; iterations = it; }
/// Create Jacobi smoother.
DSmoother(const SparseMatrix &a, int t = 0, real_t s = 1., int it = 1);
/// @brief Create a Jacobi smoother using the SparseMatrix @a a.
///
/// @param[in] a The underlying SparseMatrix
/// @param[in] t Type of Jacobi smoother (see DSmoother::JacobiType)
/// @param[in] s Scaling factor
/// @param[in] it Number of stationary iterations to perform
DSmoother(const SparseMatrix &a, JacobiType t = JACOBI, real_t s = 1.,
int it = 1) : DSmoother(t, s, it) { SetOperator(a); }
/// Replace diag entries with their abs values. Relevant only when type = 0.
/// @brief Same as DSmoother(JacobiType,real_t,int), for backwards compatbility.
DSmoother(int t, real_t s = 1., int it = 1)
: DSmoother(JacobiType(t), s, it) { }
/// @brief Same as DSmoother(const SparseMatrix&,JacobiType,real_t,int), for
/// backwards compatbility.
DSmoother(const SparseMatrix &a, int t, real_t s = 1., int it = 1)
: DSmoother(a, JacobiType(t), s, it) { }
/// @brief Replace diagonal entries with their absolute values. Relevant only
/// with JacobiType::JACOBI.
void SetPositiveDiagonal(bool pos_diag = true) { use_abs_diag = pos_diag; }
/// Matrix vector multiplication with Jacobi smoother.
/// @brief Apply the Jacobi smoother.
///
/// Applies a stationary iteration with diagonal scaling. If
/// Solver::iterative_mode is true, then @a y is used as the initial guess
/// (and the diagonal scaling is applied to the residual $x - Ay$, giving
/// $D^{-1}(x - Ay)$).
///
/// By default, Solver::iterative_mode is false and only one iteration is
/// performed, corresponding to $y = D^{-1}x$.
void Mult(const Vector &x, Vector &y) const override;
/// @brief Apply the transpose of the Jacobi smoother.
///
/// If the underlying matrix is symmetric, or if only one iteration is
/// performed with zero initial guess (Solver::iterative_mode is false), then
/// this is the same as Mult(). For non-symmetric matrices with iteration
/// count greater than one, only JacobiType::JACOBI is supported.
void MultTranspose(const Vector &x, Vector &y) const override;
};
}
+23
View File
@@ -777,8 +777,31 @@ inline real_t InnerProduct(MPI_Comm comm, const Vector &x, const Vector &y)
MPI_Allreduce(&loc_prod, &glb_prod, 1, MFEM_MPI_REAL_T, MPI_SUM, comm);
return glb_prod;
}
/// Returns the square of the Euclidean distance between two vectors in parallel
/** In parallel this computes the square of the Euclidean distance of the local
vectors, producing identical results on each MPI rank.*/
inline real_t DistanceSquared(MPI_Comm comm, const real_t *x, const real_t *y,
const int n)
{
real_t d = DistanceSquared(x, y, n);
real_t glb_d;
MPI_Allreduce(&d, &glb_d, 1, MFEM_MPI_REAL_T, MPI_SUM, comm);
return glb_d;
}
inline real_t DistanceSquared(MPI_Comm comm, const Vector &x, const Vector &y)
{
MFEM_ASSERT(x.Size() == y.Size(), "Incompatible vector sizes.");
real_t d = x.DistanceSquaredTo(y);
real_t glb_d;
MPI_Allreduce(&d, &glb_d, 1, MFEM_MPI_REAL_T, MPI_SUM, comm);
return glb_d;
}
#endif
} // namespace mfem
#endif
+259 -10
View File
@@ -36,6 +36,7 @@
#include <numeric>
#include <unordered_map>
#include <unordered_set>
#include <list>
// Include the METIS header, if using version 5. If using METIS 4, the needed
// declarations are inlined below, i.e. no header is needed.
@@ -4772,12 +4773,12 @@ Mesh::Mesh(real_t *vertices_, int num_vertices,
FinalizeTopology();
}
Mesh::Mesh( const NURBSExtension& ext )
Mesh::Mesh(const NURBSExtension& ext)
: attribute_sets(attributes), bdr_attribute_sets(bdr_attributes)
{
SetEmpty();
/// make an internal copy of the NURBSExtension
NURBSext = new NURBSExtension( ext );
NURBSext = new NURBSExtension(ext);
Dim = NURBSext->Dimension();
NumOfVertices = NURBSext->GetNV();
@@ -6537,6 +6538,8 @@ void Mesh::LoadPatchTopo(std::istream &input, Array<int> &edge_to_ukv)
Array<int> ukv_to_rpkv;
GetEdgeToUniqueKnotvector(edge_to_ukv, ukv_to_rpkv);
}
CorrectPatchTopoOrientations(edge_to_ukv);
}
void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
@@ -6547,9 +6550,9 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
const int NPKV = NP * dim; // number of patch knotvectors
constexpr int notset = -9999999;
// Sign convention
auto sign = [](int i) { return -1 - i; };
auto unsign = [](int i) { return (i < 0) ? -1 - i : i; };
// Edge index -> dimension convention
auto flipSign = [](int i) { return -1 - i; };
auto unSign = [](int i) { return (i < 0) ? -1 - i : i; };
// Local edge index -> dimension convention
auto edge_to_dim = [](int i) { return (i < 8) ? ((i & 1) ? 1 : 0) : 2; };
Array<int> v(2); // vertices of an edge
@@ -6564,7 +6567,7 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
{
GetElementVertices(i, v);
// Sign is based on the edge's vertex indices
edge_to_ukv[i] = (v[1] > v[0]) ? i : sign(i);
edge_to_ukv[i] = (v[1] > v[0]) ? i : flipSign(i);
ukv_to_rpkv[i] = i;
}
return;
@@ -6614,14 +6617,14 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
// We've set this edge already - link this index to it
if (edge_to_pkv[edge] != notset)
{
const int pkv_other = unsign(edge_to_pkv[edge]);
const int pkv_other = unSign(edge_to_pkv[edge]);
unite(pkv, pkv_other);
}
else
{
GetEdgeVertices(edge, v);
// Sign is based on the edge's vertex indices
edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : sign(pkv);
edge_to_pkv[edge] = (v[1] > v[0]) ? pkv : flipSign(pkv);
}
}
}
@@ -6648,11 +6651,255 @@ void Mesh::GetEdgeToUniqueKnotvector(Array<int> &edge_to_ukv,
edge_to_ukv.SetSize(NumOfEdges);
for (int i = 0; i < NumOfEdges; i++)
{
const int pkv = unsign(edge_to_pkv[i]);
const int pkv = unSign(edge_to_pkv[i]);
const int rpkv = pkv_to_rpkv[pkv];
const int ukv = rpkv_to_ukv[rpkv];
edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? sign(ukv) : ukv;
edge_to_ukv[i] = (edge_to_pkv[i] < 0) ? flipSign(ukv) : ukv;
}
CorrectPatchTopoOrientations(edge_to_ukv);
}
void Mesh::CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const
{
const int dim = Dimension(); // Topological (not physical) dimension
if (dim == 1) { return; }
// Sign convention
auto flipSign = [](int i) { return -1 - i; };
const Table *face2elem = GetFaceToElementTable();
Array<int> pfaces, orient;
Array<int> fe, feo;
// Finds elements sharing a face containing knotvector kv.
auto faceNeighbors = [&](int p, int kv, std::unordered_set<int> &nghb)
{
if (dim == 2) { GetElementEdges(p, pfaces, orient); }
else { GetElementFaces(p, pfaces, orient); }
for (auto face : pfaces)
{
// Check whether this face contains kv.
GetFaceEdges(face, fe, feo);
bool hasKV = false;
for (auto e : fe)
{
const int skv = edge_to_ukv[e];
if (skv == kv || flipSign(skv) == kv) { hasKV = true; }
}
if (hasKV)
{
Array<int> row;
face2elem->GetRow(face, row);
for (auto elem : row) { nghb.insert(elem); }
}
}
};
std::vector<std::vector<int>> dir_edges;
if (dim == 2)
{
dir_edges =
{
{0,2},
{1,3}
};
}
else
{
dir_edges =
{
{0,2,4,6},
{1,3,5,7},
{8,9,10,11}
};
}
Array<int> ukvs((dim==2) ? 4 : 12);
Array<int> pe, oe;
bool initKV = false;
auto setPatchDirections = [&](int p, int kv, Array<bool> &edgeSet,
std::unordered_set<int> &visited)
{
// Edges and orientations for this patch
GetElementEdges(p, pe, oe);
// Get the signed unique knot vector indices
for (int i = 0; i < pe.Size(); i++)
{
ukvs[i] = edge_to_ukv[pe[i]];
ukvs[i] = (oe[i] < 0) ? flipSign(ukvs[i]) : ukvs[i];
}
// Find the direction with this kv.
int thisDir = -1;
for (int d=0; d<dim; ++d) // Loop over directions.
{
const int skv = edge_to_ukv[pe[dir_edges[d][0]]];
if (skv == kv || flipSign(skv) == kv)
{
thisDir = d;
}
}
MFEM_VERIFY(thisDir >= 0, "");
// For this direction, find any edge already set. If no edge is set, we
// arbitrarily take the first.
int ref_edge0 = dir_edges[thisDir][0];
for (auto ref_edge : dir_edges[thisDir])
{
const int edge = pe[ref_edge];
if (edgeSet[edge])
{
ref_edge0 = ref_edge;
}
}
if (initKV && !edgeSet[pe[ref_edge0]])
{
visited.erase(p);
return false; // There is no set edge in this direction on this patch.
}
initKV = true;
// Use ref_edge0 to set other edges in this direction.
edgeSet[pe[ref_edge0]] = true;
for (auto i : dir_edges[thisDir])
{
if (i == ref_edge0)
{
continue;
}
const int edge = pe[i];
if ((dim == 2 && ukvs[i] != flipSign(ukvs[ref_edge0])) ||
(dim == 3 && ukvs[i] == flipSign(ukvs[ref_edge0])))
{
// Flip the sign of this edge
MFEM_VERIFY(!edgeSet[edge], "");
edge_to_ukv[edge] = flipSign(edge_to_ukv[edge]);
}
edgeSet[edge] = true;
}
return true;
};
Array<bool> edgeSet(NumOfEdges); // Whether edge has orientation set
edgeSet = false;
std::unordered_set<int> unset; // Patches with an unset edge
for (int i=0; i<NumOfElements; ++i) { unset.insert(i); }
const int max_iter = 3 * NumOfElements;
for (int iter=0; iter<max_iter; ++iter)
{
// Iteratively choose an unset patch (meaning not all edges have
// orientation set), choose a knotvector index for which the corresponding
// edges on this patch are not set, and sweep over all patches containing
// this knotvector. The patch sweep is ordered, by maintaining an ordered
// list `nextPatches` set by finding face-neighbor patches of visited
// patches, where the common face contains the knotvector. When each patch
// is visited, the edge orientations are set consistently. This iteration
// terminates when all edges have been set on all patches.
std::list<int> nextPatches; // Next patches to visit, ordered
std::unordered_set<int> nextSet; // nextPatches as a set
std::unordered_set<int> visited; // Visit each patch only once
if (unset.size() == 0)
{
break;
}
const int p0 = *unset.begin();
nextPatches.push_back(p0); // Start from arbitrary unset patch
nextSet.insert(p0);
// Choose an arbitrary unset direction for the first patch.
GetElementEdges(p0, pe, oe);
int unsetDim = -1;
for (int d=0; d<dim; ++d) // Loop over dimensions.
{
if (!edgeSet[pe[dir_edges[d][0]]])
{
unsetDim = d;
}
}
if (unsetDim == -1)
{
unset.erase(p0);
continue;
}
const int kv_signed = edge_to_ukv[pe[dir_edges[unsetDim][0]]];
const int kv = kv_signed < 0 ? flipSign(kv_signed) : kv_signed;
MFEM_VERIFY(!edgeSet[pe[dir_edges[unsetDim][0]]], "");
initKV = false;
while (nextPatches.size() > 0)
{
const int p = nextPatches.front();
nextPatches.pop_front();
nextSet.erase(p);
visited.insert(p);
const bool somethingSet = setPatchDirections(p, kv, edgeSet, visited);
if (!somethingSet)
{
continue;
}
// Find neighbors of patch p sharing a conforming face, via face2elem.
std::unordered_set<int> neighbors;
faceNeighbors(p, kv, neighbors);
bool allSet = true;
GetElementEdges(p, pe, oe);
for (auto edge : pe)
{
if (!edgeSet[edge])
{
allSet = false;
}
}
if (allSet)
{
unset.erase(p);
}
// Add neighbors not done to nextPatches.
for (auto n : neighbors)
{
if (n != p && visited.count(n) == 0 && unset.count(n) > 0)
{
if (nextSet.count(n) == 0)
{
nextPatches.push_back(n);
nextSet.insert(n);
}
}
}
}
}
bool allSet = true;
for (auto eset : edgeSet)
{
if (!eset)
{
allSet = false;
}
}
MFEM_VERIFY(allSet && unset.size() == 0, "Some edge is not set");
delete face2elem;
}
void Mesh::LoadNonconformingPatchTopo(std::istream &input,
@@ -9563,6 +9810,8 @@ void Mesh::GetVertices(Vector &vert_coord) const
void Mesh::SetVertices(const Vector &vert_coord)
{
MFEM_VERIFY(vert_coord.Size() == spaceDim * NumOfVertices, "");
vertices.SetSize(NumOfVertices);
for (int i = 0, nv = vertices.Size(); i < nv; i++)
for (int j = 0; j < spaceDim; j++)
{
+11 -5
View File
@@ -527,6 +527,9 @@ protected:
void PrintTopoEdges(std::ostream &out, const Array<int> &e_to_k,
bool vmap = false) const;
/// Set signs to ensure knotvectors are pointed in the same direction.
void CorrectPatchTopoOrientations(Array<int> &edge_to_ukv) const;
/// Used in GetFaceElementTransformations (...)
void GetLocalPtToSegTransformation(IsoparametricTransformation &,
int i) const;
@@ -984,8 +987,8 @@ public:
///@}
/// Construct a Mesh from a NURBSExtension
explicit Mesh( const NURBSExtension& ext );
/// Construct a Mesh from a NURBSExtension, which is deep-copied.
explicit Mesh(const NURBSExtension& ext);
/** @anchor mfem_Mesh_construction
@name Methods for piecewise Mesh construction.
@@ -2538,13 +2541,16 @@ public:
changing the mesh file itself. Examples in miniapps/nurbs/meshes. */
void RefineNURBSFromFile(std::string ref_file);
/// For NURBS meshes, insert the new knots in @a kv, for each direction.
/// For NURBS meshes, insert the new knots in @a kv, for each KnotVector.
/// The size of @a kv should be the number of KnotVectors in NURBSExtension.
void KnotInsert(Array<KnotVector*> &kv);
/// For NURBS meshes, insert the knots in @a kv, for each direction.
/// For NURBS meshes, insert the knots in @a kv, for each KnotVector.
/// The size of @a kv should be the number of KnotVectors in NURBSExtension.
void KnotInsert(Array<Vector*> &kv);
/// For NURBS meshes, remove the knots in @a kv, for each direction.
/// For NURBS meshes, remove the knots in @a kv, for each KnotVector.
/// The size of @a kv should be the number of KnotVectors in NURBSExtension.
void KnotRemove(Array<Vector*> &kv);
/* For each knot vector:
+2
View File
@@ -2827,6 +2827,8 @@ void NCNURBSExtension::PropagateFactorsForKV(int rf_default)
}
}
delete face2elem;
// For any unset entries of kvf, set to default refinement factor rf_default.
for (size_t i=0; i<kvf.size(); ++i)
{
+23 -10
View File
@@ -2278,8 +2278,9 @@ void NURBSExtension::Load(std::istream &input, bool spacing)
{
own_topo = true;
CheckPatches();
// CheckBdrPatches();
MFEM_VERIFY(CheckPatches(),
"NURBSExtension::CheckPatch"
"\n Inconsistent edge-to-knotvector mapping!");
skip_comment_lines(input, '#');
@@ -2676,7 +2677,9 @@ NURBSExtension::NURBSExtension(const Mesh *patch_topology,
patchTopo->GetEdgeToUniqueKnotvector(edge_to_ukv, ukv_to_rpkv);
own_topo = true;
CheckPatches(); // This is checking the edge_to_ukv mapping
MFEM_VERIFY(CheckPatches(),
"NURBSExtension::CheckPatch"
"\n Inconsistent edge-to-knotvector mapping!");
// Set number of unique (not comprehensive) knot vectors
NumOfKnotVectors = ukv_to_rpkv.Size();
@@ -3246,9 +3249,9 @@ void NURBSExtension::MergeGridFunctions(
}
}
void NURBSExtension::CheckPatches()
bool NURBSExtension::CheckPatches()
{
if (Dimension() == 1 ) { return; }
if (Dimension() == 1 ) { return true; }
Array<int> edges, oedge;
@@ -3256,6 +3259,7 @@ void NURBSExtension::CheckPatches()
{
patchTopo->GetElementEdges(p, edges, oedge);
// Convert to ukv and apply sign-flip
for (int i = 0; i < edges.Size(); i++)
{
edges[i] = edge_to_ukv[edges[i]];
@@ -3265,6 +3269,16 @@ void NURBSExtension::CheckPatches()
}
}
// In 2d - opposite edges must be same knotvector with opposite sign.
// In 3d - opposite edges must be same knotvector with same sign.
// This logic is the result of Mesh::GetElementEdges setting orientation
// for edges based on ascending vertex indices, using reference vertex
// ordering
// {0, 1}, {1, 2}, {2, 3}, {3, 0} for Geometry::SQUARE in 2D
// and
// {0, 1}, {1, 2}, {3, 2}, {0, 3}, {4, 5}, {5, 6},
// {7, 6}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7} for Geometry::CUBE in 3D
// See fem/geom.cpp for these definitions.
if ((Dimension() == 2 &&
(edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
@@ -3275,11 +3289,10 @@ void NURBSExtension::CheckPatches()
edges[8] != edges[9] || edges[8] != edges[10] ||
edges[8] != edges[11])))
{
mfem::err << "NURBSExtension::CheckPatch (patch = " << p
<< ")\n Inconsistent edge-to-knotvector mapping!";
mfem_error();
return false;
}
}
return true;
}
void NURBSExtension::CheckBdrPatches()
@@ -3757,8 +3770,8 @@ void NURBSExtension::GetPatchOffsets(int &meshCounter, int &spaceCounter)
if (dim == 1)
{
meshCounter += KnotVec(0)->GetNE() - 1;
spaceCounter += KnotVec(0)->GetNCP() - 2;
meshCounter += KnotVec(p)->GetNE() - 1;
spaceCounter += KnotVec(p)->GetNCP() - 2;
}
else if (dim == 2)
{
+4 -4
View File
@@ -596,11 +596,8 @@ protected:
if the KnotVector index associated with edge @a edge is negative. */
inline const KnotVector *KnotVec(int edge, int oedge, int *okv) const;
/// Throw an error if any patch has an inconsistent edge_to_ukv mapping.
void CheckPatches();
/// Throw an error if any boundary patch has invalid KnotVector orientation.
void CheckBdrPatches();
MFEM_DEPRECATED void CheckBdrPatches();
/** @brief Return the directions in @a kvdir of the KnotVectors in patch @a p
based on the patch edge orientations. Each entry of @a kvdir is -1 if the
@@ -794,6 +791,9 @@ public:
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces,
GridFunction &merged);
/// Returns false if any patch has an inconsistent edge_to_ukv mapping.
bool CheckPatches();
/// Destroy a NURBSExtension.
virtual ~NURBSExtension();
+1
View File
@@ -22,6 +22,7 @@ add_subdirectory(common)
add_subdirectory(contact)
add_subdirectory(dfem)
add_subdirectory(diag-smoothers)
add_subdirectory(multiapp)
add_subdirectory(dpg)
add_subdirectory(electromagnetics)
add_subdirectory(fluids/navier)
@@ -127,7 +127,7 @@ int main(int argc, char *argv[])
}
// Load mesh + complete any serial refinements
Mesh mesh("../../data/channel-bifurcation-2d.mesh");
Mesh mesh("../../../data/channel-bifurcation-2d.mesh");
for (int lev = 0; lev < ctx.rs_levels; lev++)
{
mesh.UniformRefinement();
+16
View File
@@ -362,6 +362,22 @@ void NavierSolver::UpdateTimestepHistory(real_t dt)
un_gf.SetFromTrueDofs(un);
}
void NavierSolver::Step(Vector &up, real_t &t, real_t &dt)
{
Array<int> offsets({0, vfes->GetTrueVSize(), pfes->GetTrueVSize()});
offsets.PartialSum();
BlockVector upb(up.GetData(), offsets);
un_gf.SetFromTrueDofs(upb.GetBlock(0));
un_next_gf.SetFromTrueDofs(upb.GetBlock(0));
pn_gf.SetFromTrueDofs(upb.GetBlock(1));
Step(t, dt, (int) t/dt,false);
un_gf.GetTrueDofs(upb.GetBlock(0));
pn_gf.GetTrueDofs(upb.GetBlock(1));
}
void NavierSolver::Step(real_t &time, real_t dt, int current_step,
bool provisional)
{
+2
View File
@@ -182,6 +182,8 @@ public:
*/
void Step(real_t &time, real_t dt, int cur_step, bool provisional = false);
void Step(Vector &up, real_t &t, real_t &dt);
/// Return a pointer to the provisional velocity ParGridFunction.
ParGridFunction *GetProvisionalVelocity() { return &un_next_gf; }
+234 -12
View File
@@ -33,8 +33,7 @@
using namespace std;
using namespace mfem;
void ReflectPoint(Vector & p, Vector const& origin, Vector const& normal)
void ReflectPoint(Vector &p, const Vector &origin, const Vector &normal)
{
Vector diff(3);
Vector proj(3);
@@ -57,12 +56,12 @@ private:
// Map from reflected to original mesh elements
std::vector<int> *r2o;
std::vector<std::vector<int>> *perm;
std::vector<std::array<int, 8>> *perm;
public:
ReflectedCoefficient(VectorCoefficient &A, Vector const& origin_,
Vector const& normal_, std::vector<int> *r2o_,
Mesh *mesh, std::vector<std::vector<int>> *refPerm) :
Mesh *mesh, std::vector<std::array<int, 8>> *refPerm) :
VectorCoefficient(3), a(&A), origin(origin_), normal(normal_),
meshOrig(mesh), r2o(r2o_), perm(refPerm)
{ }
@@ -109,7 +108,7 @@ void ReflectedCoefficient::Eval(Vector &V, ElementTransformation &T,
// give the columns of A.
// Permutation p is such that hex_reflected[i] = hex_init[p[i]]
const std::vector<int>& p = (*perm)[elem];
const std::array<int, 8>& p = (*perm)[elem];
// ip is on reflected hex. We map from the reflected hex to the initial
// hex, in reference space. Thus we use y = Ax + b, where x is in the
@@ -164,7 +163,7 @@ void ReflectedCoefficient::Eval(Vector &V, ElementTransformation &T,
// Find perm such that h1[i] = h2[perm[i]]
void GetHexPermutation(Array<int> const& h1, Array<int> const& h2,
std::vector<int> & perm)
std::array<int, 8> &perm)
{
std::map<int, int> h2inv;
const int n = perm.size();
@@ -236,7 +235,7 @@ public:
int AddElement(Array<int> const& vertices, const bool reorder);
Mesh *mesh;
std::vector<std::vector<int>> refPerm;
std::vector<std::array<int, 8>> refPerm;
private:
std::vector<std::vector<int>> faces;
@@ -279,13 +278,13 @@ int HexMeshBuilder::AddElement(Array<int> const& vertices, const bool reorder)
}
while (reordered);
std::vector<int> perm_e(8);
std::array<int, 8> perm_e;
GetHexPermutation(rvert, vertices, perm_e);
refPerm.push_back(perm_e);
}
else
{
refPerm.push_back(std::vector<int> {0, 1, 2, 3, 4, 5, 6, 7});
refPerm.push_back(std::array<int, 8> {0, 1, 2, 3, 4, 5, 6, 7});
}
SaveHexFaces(mesh->GetNE(), rvert);
@@ -762,7 +761,10 @@ bool GetMeshElementOrder(Mesh const& mesh, Vector const& origin,
return true;
}
Mesh* ReflectHighOrderMesh(Mesh & mesh, Vector origin, Vector normal)
Mesh* ReflectHighOrderMesh(Mesh &mesh,
const Vector &origin, const Vector &normal,
std::vector<std::array<int, 8>> &hexPerm,
std::vector<int> &elOrder)
{
MFEM_VERIFY(mesh.Dimension() == 3, "Only 3D meshes can be reflected");
@@ -837,7 +839,6 @@ Mesh* ReflectHighOrderMesh(Mesh & mesh, Vector origin, Vector normal)
}
}
std::vector<int> elOrder;
const bool onPlane = GetMeshElementOrder(mesh, origin, normal, elOrder);
for (int eidx=0; eidx<mesh.GetNE(); eidx++)
@@ -1006,6 +1007,133 @@ Mesh* ReflectHighOrderMesh(Mesh & mesh, Vector origin, Vector normal)
*reflected_nodes = newReflectedNodes;
}
hexPerm = builder.refPerm;
return reflected;
}
void ReorderHexArray(const std::array<int, 3> &dim,
const array<int, 8> &hexperm,
std::array<int, 3> &dir, std::array<int, 3> &dims,
Array3D<int> &permArray);
NURBSPatch* ReflectPatch(NURBSPatch *patch, int nx, int ny, int nz,
const Vector &origin, const Vector &normal,
const std::array<int, 8> &hexPerm)
{
// The hexahedral element for this patch in the reflected patch topology mesh
// is the reflection of an original patch topology mesh element, with
// reference vertices permuted according to hexPerm. The original grid of
// (nx + 1) x (ny + 1) x (nz + 1)
// control points has a new size and ordering, depending on hexPerm. Now,
// ReorderHexArray finds the new dimensions of this grid in `dims`, maps the
// directions in `dir`, and sets the permutation of grid indices as triples
// in `permArray`.
std::array<int, 3> dims, dir;
Array3D<int> permArray;
ReorderHexArray({nx+1, ny+1, nz+1}, hexPerm, dir, dims, permArray);
const KnotVector *kv0 = patch->GetKV(dir[0]);
const KnotVector *kv1 = patch->GetKV(dir[1]);
const KnotVector *kv2 = patch->GetKV(dir[2]);
NURBSPatch *rpatch = new NURBSPatch(kv0, kv1, kv2, 4);
// Reflect the control points in this reflected patch `rpatch`.
Vector vr(3);
for (int i=0; i<dims[0]; ++i)
{
for (int j=0; j<dims[1]; ++j)
{
for (int k=0; k<dims[2]; ++k)
{
const int old = permArray(i,j,k);
const int i0 = old / ((ny + 1) * (nz + 1));
const int j0 = (old - (i0 * (ny + 1) * (nz + 1))) / (nz + 1);
const int k0 = old - (i0 * (ny + 1) * (nz + 1)) - (j0 * (nz+1));
const real_t w = (*patch)(i0,j0,k0,3); // Weight
for (int l=0; l<3; ++l) { vr[l] = (*patch)(i0,j0,k0,l) / w; }
ReflectPoint(vr, origin, normal);
for (int l=0; l<3; ++l) { (*rpatch)(i,j,k,l) = vr[l] * w; }
(*rpatch)(i,j,k,3) = w;
}
}
}
return rpatch;
}
Mesh* ReflectNURBSMesh(Mesh &mesh, const Vector &origin, const Vector &normal)
{
MFEM_VERIFY(mesh.NURBSext && mesh.Dimension() == 3,
"Only 3D NURBS meshes can be reflected");
Mesh patchTopo = mesh.NURBSext->GetPatchTopology(); // Deep copy
Array<NURBSPatch*> patchesOriginal, patches;
mesh.GetNURBSPatches(patchesOriginal); // Deep copy
NURBSPatchMap p2g(mesh.NURBSext);
const KnotVector *kv[3];
const int pnv = patchTopo.GetNV();
Vector vert_coord(3 * patchTopo.GetNV());
for (int p=0; p<patchesOriginal.Size(); ++p)
{
p2g.SetPatchDofMap(p, kv);
const int nx = p2g.nx();
const int ny = p2g.ny();
const int nz = p2g.nz();
Array<int> vert;
patchTopo.GetElementVertices(p, vert);
for (int l=0; l<3; ++l)
{
const int os = l * pnv;
vert_coord[vert[0] + os] = (*patchesOriginal[p])(0,0,0,l);
vert_coord[vert[1] + os] = (*patchesOriginal[p])(nx,0,0,l);
vert_coord[vert[2] + os] = (*patchesOriginal[p])(nx,ny,0,l);
vert_coord[vert[3] + os] = (*patchesOriginal[p])(0,ny,0,l);
vert_coord[vert[4] + os] = (*patchesOriginal[p])(0,0,nz,l);
vert_coord[vert[5] + os] = (*patchesOriginal[p])(nx,0,nz,l);
vert_coord[vert[6] + os] = (*patchesOriginal[p])(nx,ny,nz,l);
vert_coord[vert[7] + os] = (*patchesOriginal[p])(0,ny,nz,l);
}
}
patchTopo.SetVertices(vert_coord);
std::vector<std::array<int, 8>> hexPerm;
std::vector<int> elOrder;
Mesh *reflectedPatchTopo = ReflectHighOrderMesh(patchTopo, origin, normal,
hexPerm, elOrder);
// Construct reflected patches. Note that reflectedPatchTopo has patch
// ordering depending on patchTopo.
for (int p=0; p<patchesOriginal.Size(); ++p)
{
const int p_orig = elOrder[p]; // TODO: use r2o instead?
p2g.SetPatchDofMap(p_orig, kv);
const int nx = p2g.nx();
const int ny = p2g.ny();
const int nz = p2g.nz();
patches.Append(patchesOriginal[p_orig]);
patches.Append(ReflectPatch(patchesOriginal[p_orig], nx, ny, nz,
origin, normal, hexPerm[(2 * p) + 1]));
}
NURBSExtension *ne = new NURBSExtension(reflectedPatchTopo, patches);
delete reflectedPatchTopo;
for (auto patch : patches) { delete patch; }
Mesh *reflected = new Mesh(*ne);
delete ne;
return reflected;
}
@@ -1045,7 +1173,19 @@ int main(int argc, char *argv[])
Mesh mesh(mesh_file, 0, 0);
Mesh *reflected = ReflectHighOrderMesh(mesh, origin, normal);
Mesh *reflected{nullptr};
//if (mesh.IsNURBS()) // TODO: available in PR 4936
if (mesh.NURBSext)
{
reflected = ReflectNURBSMesh(mesh, origin, normal);
}
else
{
std::vector<std::array<int, 8>> hexPerm;
std::vector<int> elOrder;
reflected = ReflectHighOrderMesh(mesh, origin, normal, hexPerm, elOrder);
}
// Save the final mesh
ofstream mesh_ofs("reflected.mesh");
@@ -1065,3 +1205,85 @@ int main(int argc, char *argv[])
return 0;
}
void HexVertexIJK(const int idx, std::array<int, 3>& ijk)
{
ijk[2] = idx / 4;
const int id2d = idx - (4 * ijk[2]);
ijk[1] = id2d / 2;
ijk[0] = (ijk[1] == 0) ? id2d : 3 - id2d;
}
void ReorderHexArray(const std::array<int, 3> &dim,
const array<int, 8> &hexperm,
std::array<int, 3> &dir, std::array<int, 3> &dims,
Array3D<int> &permArray)
{
int prinV[4] = {0, 1, 3, 4}; // Vertices in principal directions (after 0)
int newPrinV[4];
// newVertices[i] = oldVertices[hexperm[i]]
// Hence newPrinV[0] = hexperm[0] is the index
// of new vertex 0 in the old hex.
std::array<int, 3> newIJK[4];
for (int i = 0; i < 4; ++i)
{
newPrinV[i] = hexperm[prinV[i]];
HexVertexIJK(newPrinV[i], newIJK[i]);
}
// For direction i in the new hex, dir[i] is the direction in the old hex.
Array<bool> rev(3);
for (int i = 0; i < 3; ++i)
{
bool iset = false;
for (int j = 0; j < 3; ++j)
{
const int d = newIJK[i + 1][j] - newIJK[0][j];
if (d != 0)
{
MFEM_VERIFY(!iset, "");
MFEM_VERIFY(d == 1 || d == -1, "");
dir[i] = j;
rev[i] = (d == -1);
iset = true;
}
}
MFEM_VERIFY(iset, "");
dims[i] = dim[dir[i]];
}
MFEM_VERIFY(dir[0] + dir[1] + dir[2] == 3, "");
permArray.SetSize(dims[0], dims[1], dims[2]);
Array<int> old_ijk(3);
Array<int> new_ijk(3);
for (int i = 0; i < dims[0]; ++i)
for (int j = 0; j < dims[1]; ++j)
for (int k = 0; k < dims[2]; ++k)
{
new_ijk[0] = i;
new_ijk[1] = j;
new_ijk[2] = k;
for (int m = 0; m < 3; ++m)
{
const int d = dir[m]; // Old hex direction
if (rev[m])
{
old_ijk[d] = dim[d] - 1 - new_ijk[m];
}
else
{
old_ijk[d] = new_ijk[m];
}
}
permArray(i, j, k) =
old_ijk[2] + (old_ijk[1] * dim[2]) + (old_ijk[0] * dim[1] * dim[2]);
}
}
+58
View File
@@ -0,0 +1,58 @@
# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
# LICENSE and NOTICE for details. LLNL-CODE-806117.
#
# This file is part of the MFEM library. For more information and source code
# availability visit https://mfem.org.
#
# MFEM is free software; you can redistribute it and/or modify it under the
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
set(MESH_FILES
backward-facing-step.msh
channel-cylinder.msh
)
# Add a target to copy the mesh files from the source directory; used by sample
# runs.
set(SRC_MESH_FILES)
foreach(MESH_FILE ${MESH_FILES})
list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/${MESH_FILE})
endforeach()
add_custom_command(OUTPUT data_is_copied
COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} .
COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
COMMENT "Copying multiapp miniapps data files ...")
add_custom_target(copy_miniapps_multiapp_data DEPENDS data_is_copied)
list(APPEND MULTIAPP_COMMON_SOURCES
multiapp.cpp)
list(APPEND MULTIAPP_COMMON_HEADERS
multiapp.hpp)
set(MULTIAPP_COMMON_FILES
EXTRA_SOURCES ${MULTIAPP_COMMON_SOURCES}
EXTRA_HEADERS ${MULTIAPP_COMMON_HEADERS})
# Parallel apps.
if (MFEM_USE_MPI)
add_mfem_miniapp(cht-BFS
MAIN cht-BFS.cpp
${MFEM_MINIAPPS_COMMON_HEADERS} ${MULTIAPP_COMMON_FILES} ../fluids/navier/navier_solver.hpp ../fluids/navier/navier_solver.cpp
LIBRARIES mfem-common)
add_dependencies(cht-BFS copy_miniapps_multiapp_data)
# add_mfem_miniapp(fsi
# MAIN fsi.cpp
# ${MFEM_MINIAPPS_COMMON_HEADERS} ${MULTIAPP_COMMON_FILES}
# LIBRARIES mfem-common)
# add_dependencies(fsi copy_miniapps_multiapp_data)
# Add parallel tests.
# if (MFEM_ENABLE_TESTING)
# endif()
endif()
+323
View File
@@ -0,0 +1,323 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#ifndef MFEM_ELASTICITY_HPP
#define MFEM_ELASTICITY_HPP
#include "mfem.hpp"
#include "../multiapp.hpp"
using namespace mfem;
/**
* @brief Elasticity time dependent operator
*
* ρdv/dt = (-σ + f)
* du/dt = v
*
* where σ is the stress tensor, f is the body force,
* and ρ is the density.
*/
class Elasticity : public Application
{
public:
// Mesh and finite element space
ParMesh &mesh;
ParFiniteElementSpace &fes;
Array<int> offsets;
/// Essential and natural dof array.
Array<int> ess_attr, nat_attr;
Array<int> ess_tdofs, nat_tdofs;
/// Material properties
ConstantCoefficient density, mu, lambda;
/// Grid functions for the displacement, velocity, and traction
mutable ParGridFunction x_gf; ///< displacement
mutable ParGridFunction u_gf; ///< velocity (dx/dt)
mutable ParGridFunction stress_gf; ///< traction
mutable ParGridFunction bc_send_gf; ///< Grid functions for transfering BCs
/// Mass and Stiffness forms
mutable ParBilinearForm Mform, Kform, Kform_e;
/// RHS form
mutable ParLinearForm bform;
mutable Vector b;
/// Mass and Stiffness operators
mutable HypreParMatrix Mmat, Kmat, Kmat_e, Mrhomat, Mmat_e;
mutable HypreParMatrix *T = nullptr;
/// Force
VectorGridFunctionCoefficient fcoeff;
ScalarVectorProductCoefficient scaled_fcoeff;
/// Mass matrix and implicit solver
GMRESSolver implicit_solver;
CGSolver M_solver;
/// Preconditioner
HypreSmoother M_prec;
mutable HypreBoomerAMG *amg = nullptr;
mutable Solver *pc = nullptr;
/// Auxiliary variables
mutable Vector z;
real_t current_dt = -1.0;
bool updated = false;
public:
Elasticity(ParFiniteElementSpace &fes_,
Array<int> ess_attr_,
Array<int> nat_attr_,
real_t density_ = 1.0,
real_t mu_ = 1.0,
real_t lambda_ = 1.0) :
Application(2*fes_.GetTrueVSize()),
mesh(*fes_.GetParMesh()),
fes(fes_),
ess_attr(ess_attr_),
nat_attr(nat_attr_),
density(density_),
mu(mu_),
lambda(lambda_),
x_gf(&fes), u_gf(&fes), stress_gf(&fes),
bc_send_gf(&fes),
Mform(&fes), Kform(&fes),
Kform_e(&fes), bform(&fes),
fcoeff(&stress_gf), scaled_fcoeff(-1.0, fcoeff),
implicit_solver(mesh.GetComm()),
M_solver(mesh.GetComm())
{
fes.GetEssentialTrueDofs(ess_attr, ess_tdofs);
fes.GetEssentialTrueDofs(nat_attr, nat_tdofs);
x_gf = 0.0;
u_gf = 0.0;
stress_gf = 0.0;
bc_send_gf = 0.0;
// Setup field collection for output and transfer
field_collection.SetName("Elasticity");
field_collection.AddField("Displacement", &x_gf);
field_collection.AddField("Velocity", &u_gf);
field_collection.AddField("Traction", &stress_gf);
field_collection.AddSourceField("Displacement_BC", &bc_send_gf);
field_collection.AddSourceField("Velocity_BC", &bc_send_gf);
offsets = Array<int>({0, fes.GetTrueVSize(), fes.GetTrueVSize()});
offsets.PartialSum();
Mform.AddDomainIntegrator(new VectorMassIntegrator(density));
Kform.AddDomainIntegrator(new ElasticityIntegrator(lambda, mu));
Kform_e.AddDomainIntegrator(new ElasticityIntegrator(lambda, mu));
if(nat_attr.Size() > 0)
{
bform.AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(scaled_fcoeff), nat_attr);
}
b.SetSize(fes.GetTrueVSize());
z.SetSize(fes.GetTrueVSize());
Assemble();
BuildSolvers();
}
/// Assemble linear and bilinear forms; called if mesh is updated
void Assemble() override
{
AssembleLinearForms();
AssembleBilinearForms();
}
/// Assemble linear forms for traction
void AssembleLinearForms()
{
bform.Assemble();
b.SetSize(fes.GetTrueVSize());
bform.ParallelAssemble(b);
}
/// Assemble bilinear forms for mass and stiffness
void AssembleBilinearForms()
{
Mform.Assemble();
Kform.Assemble();
Kform_e.Assemble();
Array<int> empty;
Mform.FormSystemMatrix(ess_tdofs, Mrhomat);
Kform_e.FormSystemMatrix(ess_tdofs, Kmat_e);
Kform.FormSystemMatrix(empty, Kmat);
}
/// Update finite element space and re-assemble forms
/// if the mesh has changed
void Update() override
{
fes.Update();
u_gf.Update();
x_gf.Update();
stress_gf.Update();
bc_send_gf.Update();
Mform.Update();
Kform.Update();
Kform_e.Update();
bform.Update();
Assemble();
updated = true;
}
/// Build implicit solver and AMG preconditioner
void BuildSolvers()
{
M_solver.iterative_mode = false;
M_solver.SetRelTol(1e-8);
M_solver.SetAbsTol(1e-8);
M_solver.SetMaxIter(1000);
M_solver.SetPrintLevel(0);
M_solver.SetOperator(Mrhomat);
M_prec.SetType(HypreSmoother::Jacobi);
M_prec.SetOperator(Mrhomat);
M_solver.SetPreconditioner(M_prec);
if(amg) delete amg;
amg = new HypreBoomerAMG;
HYPRE_BoomerAMGSetSmoothType(*amg, 5);
amg->SetOperator(Kmat_e);
amg->SetSystemsOptions(2, true);
amg->SetElasticityOptions(&fes);
amg->SetPrintLevel(0);
pc = amg;
implicit_solver.iterative_mode = false;
implicit_solver.SetRelTol(1e-4);
implicit_solver.SetAbsTol(0.0);
implicit_solver.SetMaxIter(500);
implicit_solver.SetKDim(200);
implicit_solver.SetPrintLevel(0);
implicit_solver.SetPreconditioner(*pc);
}
/// Apply operator
void Mult(const Vector &u, Vector &k) const override
{
BlockVector ub(u.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
Vector &vel = ub.GetBlock(0);
Vector &pos = ub.GetBlock(1);
Vector &ku = kb.GetBlock(0);
Vector &kx = kb.GetBlock(1);
Kmat.Mult(pos,z);
z.Neg();
z.Add(1.0, b);
kx = vel; // dx/dt = u
M_solver.Mult(z, ku);
ku.SetSubVector(ess_tdofs, 0.0);
kx.SetSubVector(ess_tdofs, 0.0);
}
/// Solve implicit system in Schur complement form
void ImplicitSolve(const real_t dt, const Vector &u, Vector &k)
{
BlockVector ub(u.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
Vector &vel = ub.GetBlock(0);
Vector &pos = ub.GetBlock(1);
Vector &ku = kb.GetBlock(0);
Vector &kx = kb.GetBlock(1);
AssembleLinearForms();
Kmat.Mult(pos,z);
z.Neg();
z.Add(1.0, b);
Kmat.AddMult(vel,z,-dt); // z = z - dt*K*vel
if((current_dt != dt) || updated)
{
if (T) { delete T; }
current_dt = dt;
T = Add(1.0, Mrhomat, dt*dt, Kmat_e);
implicit_solver.SetOperator(*T);
amg->SetOperator(*T);
updated = false;
}
implicit_solver.Mult(z, ku);
add(1.0, vel, dt, ku, kx);
ku.SetSubVector(ess_tdofs, 0.0);
kx.SetSubVector(ess_tdofs, 0.0);
}
/// Computes residual of the elasticity equations
void ImplicitMult(const Vector &x, const Vector &k, Vector &v ) const override
{
BlockVector xb(x.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
BlockVector vb(v.GetData(), offsets);
// Apply the implicit operator
Mrhomat.Mult(kb.GetBlock(0), vb.GetBlock(0)); // v = A(x) + M*k
vb.GetBlock(0).Add(-1.0, b);
Kmat.AddMult(xb.GetBlock(1), vb.GetBlock(0), 1.0); // v = K*x
vb.GetBlock(1) = kb.GetBlock(1); // v = u
vb.GetBlock(1) -= xb.GetBlock(0); // v = k
vb.GetBlock(0).SetSubVector(ess_tdofs, 0.0);
vb.GetBlock(1).SetSubVector(ess_tdofs, 0.0);
}
void PreProcess(Vector &x) override {}
void PostProcess(Vector &x) override {}
void Transfer(const Vector &x) override
{
BlockVector xb(x.GetData(), offsets);
field_collection.Transfer("Velocity_BC", xb.GetBlock(0));
field_collection.Transfer("Displacement_BC", xb.GetBlock(1));
}
void Transfer(const Vector &u, const Vector &k, real_t dt = 0.0) override
{
BlockVector kb(k.GetData(), offsets);
field_collection.Transfer("Velocity_BC", kb.GetBlock(0));
field_collection.Transfer("Displacement_BC", kb.GetBlock(1));
}
~Elasticity() override
{
if (amg) { delete amg; }
if (T) { delete T; }
}
};
#endif
+235
View File
@@ -0,0 +1,235 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#ifndef MFEM_MESH_MORPHER_HPP
#define MFEM_MESH_MORPHER_HPP
#include "mfem.hpp"
#include "../multiapp.hpp"
using namespace mfem;
/**
* @brief Mesh morphing is modeled as a time dependent vector diffusion equation
*
* dx/dt = κΔx
*
* with vector diffusion coefficient, κ.
*/
class MeshDiffusion : public Application
{
public:
// Mesh and finite element space
ParMesh &mesh;
ParFiniteElementSpace &fes;
/// Essential true dof array.
Array<int> ess_attr, ess_tdofs;
/// Diffusion coefficient
ConstantCoefficient kappa;
/// Grid function for the mesh displacement variable
mutable ParGridFunction x_gf;
mutable ParGridFunction x_gf_bc;
mutable ParGridFunction u_gf_bc;
mutable ParGridFunction bc_send_gf; ///< Grid functions for transfering BCs
/// Mass and Stiffness forms
mutable ParBilinearForm Mform, Kform, Kform_e;
/// Mass and Stiffness operators
OperatorHandle M, K;
mutable HypreParMatrix Mmat, Kmat, Mmat_e, Kmat_e;
/// Mass matrix and implicit solver
mutable CGSolver M_solver;
mutable GMRESSolver implicit_solver;
HypreParMatrix *T = nullptr; // T = M + dt K
/// Mass matrix preconditioner
HypreSmoother M_prec;
HypreSmoother T_prec;
real_t current_dt = -1.0;
/// Auxiliary vectors
mutable Vector z, zv;
bool updated = false;
public:
MeshDiffusion(ParFiniteElementSpace &fes_,
Array<int> ess_attr_,
real_t kappa_ = 1.0e0) :
Application(fes_.GetTrueVSize()),
mesh(*fes_.GetParMesh()),
fes(fes_),
ess_attr(ess_attr_),
kappa(kappa_),
x_gf(&fes),
x_gf_bc(&fes), u_gf_bc(&fes),
bc_send_gf(&fes),
Mform(&fes), Kform(&fes),Kform_e(&fes),
M_solver(mesh.GetComm()),
implicit_solver(mesh.GetComm()),
z(fes_.GetTrueVSize())
{
fes.GetEssentialTrueDofs(ess_attr, ess_tdofs);
x_gf = 0.0;
x_gf_bc = 0.0;
u_gf_bc = 0.0;
bc_send_gf = 0.0;
// Setup field collection for output and transfer
field_collection.SetName("Mesh-Diffusion");
field_collection.AddSourceField("Displacement", &x_gf);
field_collection.AddSourceField("dxdt", &bc_send_gf);
field_collection.AddField("Displacement_BC", &x_gf_bc);
field_collection.AddField("Velocity_BC", &u_gf_bc);
Mform.AddDomainIntegrator(new VectorMassIntegrator);
Kform.AddDomainIntegrator(new VectorDiffusionIntegrator(kappa));
Kform_e.AddDomainIntegrator(new VectorDiffusionIntegrator(kappa));
Assemble();
BuildSolvers();
}
void Assemble() override
{
Mform.Assemble();
Kform.Assemble();
Kform_e.Assemble();
Array<int> empty;
Mform.FormSystemMatrix(ess_tdofs, Mmat_e);
Kform_e.FormSystemMatrix(ess_tdofs, Kmat_e);
Kform.FormSystemMatrix(empty, Kmat);
}
void Update() override
{
fes.Update();
x_gf.Update();
x_gf_bc.Update();
u_gf_bc.Update();
bc_send_gf.Update();
Mform.Update();
Kform.Update();
Kform_e.Update();
Assemble();
updated = true;
}
void BuildSolvers()
{
M_solver.iterative_mode = false;
M_solver.SetRelTol(1e-8);
M_solver.SetAbsTol(0.0);
M_solver.SetMaxIter(100);
M_solver.SetPrintLevel(0);
M_solver.SetOperator(Mmat_e);
M_prec.SetType(HypreSmoother::Jacobi);
M_solver.SetPreconditioner(M_prec);
implicit_solver.iterative_mode = false;
implicit_solver.SetRelTol(1e-8);
implicit_solver.SetAbsTol(1e-5);
implicit_solver.SetMaxIter(500);
implicit_solver.SetPrintLevel(0);
T_prec.SetType(HypreSmoother::Jacobi);
implicit_solver.SetPreconditioner(T_prec);
}
void Mult(const Vector &x, Vector &k) const override
{
Kmat.Mult(x, z);
z.Neg();
M_solver.Mult(z, k);
k.SetSubVector(ess_tdofs, 0.0);
}
void ImplicitSolve(const real_t dt, const Vector &x, Vector &k)
{
if((current_dt != dt) || updated)
{
if (T) delete T;
T = Add(1.0, Mmat_e, dt, Kmat_e);
implicit_solver.SetOperator(*T);
current_dt = dt;
updated = false;
}
Kmat.Mult(x, z);
z.Neg();
implicit_solver.Mult(z, k);
if(IsCoupled())
{
x_gf_bc.GetTrueDofs(z); // contains the BC in terms of kx (dx/dt)
for (int i = 0; i < ess_tdofs.Size(); i++)
{
int idx = ess_tdofs[i];
k(idx) = z(idx);
}
}
else
{
k.SetSubVector(ess_tdofs, 0.0);
}
}
void ImplicitMult(const Vector &x, const Vector &k, Vector &v ) const override
{}
void PreProcess(Vector &x) override
{}
void PostProcess(Vector &x) override {}
void Transfer(const Vector &x) override
{
field_collection.Transfer("Displacement", x);
zv.SetSize(fes.GetTrueVSize());
ParGridFunction *dxdt = field_collection.GetField("Velocity_BC");
dxdt->GetTrueDofs(zv);
bc_send_gf.GetTrueDofs(z);
for (int i = 0; i < ess_tdofs.Size(); i++)
{
int idx = ess_tdofs[i];
z(idx) = zv(idx);
}
field_collection.Transfer("dxdt", z);
}
void Transfer(const Vector &u, const Vector &k, real_t dt = 0.0) override
{
field_collection.Transfer("dxdt", k);
}
~MeshDiffusion() override
{
if (T) delete T;
}
};
#endif
+503
View File
@@ -0,0 +1,503 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#ifndef MFEM_NAVIER_STOKES_HPP
#define MFEM_NAVIER_STOKES_HPP
#include "mfem.hpp"
#include "../multiapp.hpp"
using namespace mfem;
/// @brief Compute deviatoric stress on a boundary
class DeviatoricStressCoefficient : public VectorCoefficient
{
protected:
ParGridFunction *p_gf = nullptr; ///< Grid function for pressure
ParGridFunction *u_gf = nullptr; ///< Grid function for velocity
Coefficient *viscosity = nullptr; ///< Coefficient for the kinematic viscosity
Coefficient *density = nullptr; ///< Coefficient for the density
bool scaled_pressure = false; ///< True if pressure is density-scaled
DenseMatrix dudx, tau; ///< Velocity gradient tensor
Vector normal;
public:
// For velocity and pressure case
DeviatoricStressCoefficient(ParGridFunction *p_gf_,
ParGridFunction *u_gf_,
Coefficient *viscosity_,
Coefficient *density_,
bool scaled_pressure_ = false) :
VectorCoefficient(p_gf_->FESpace()->GetMesh()->Dimension()),
p_gf(p_gf_), u_gf(u_gf_),
viscosity(viscosity_), density(density_), scaled_pressure(scaled_pressure_)
{
int ndim = p_gf->FESpace()->GetMesh()->Dimension();
normal.SetSize(ndim);
if(u_gf)
{
dudx.SetSize(ndim, vdim);
tau.SetSize(ndim, vdim);
}
}
// For pressure-only case
DeviatoricStressCoefficient(ParGridFunction *p_gf_) :
DeviatoricStressCoefficient(p_gf_,NULL,NULL,NULL) {}
/// Evaluate the deviatoric stress tensor at the given point.
void Eval(Vector &v, ElementTransformation &T,
const IntegrationPoint &ip) override
{
v.SetSize(vdim);
v=0.0;
const DenseMatrix &jacobian = T.Jacobian();
MFEM_ASSERT( (jacobian.Height() - 1 == jacobian.Width()),
"Incorrect Jacobian dimension. Coefficient only supported for boundary elements.");
CalcOrtho(jacobian, normal);
const double scale = normal.Norml2();
normal /= scale;
real_t rho = density ? density->Eval(T, ip) : 1.0;
if(p_gf)
{ // Add pressure term to the deviatoric stress tensor
real_t pscale = scaled_pressure ? rho : 1.0;
for (int i = 0; i < vdim; ++i)
{
v(i) = -pscale*p_gf->GetValue(T, ip) * normal(i);
}
}
if(u_gf)
{
real_t nu = viscosity ? viscosity->Eval(T, ip) : 1.0;
u_gf->GetVectorGradient(T, dudx);
for (int i = 0; i < vdim; i++)
{
for (int j = 0; j < vdim; j++)
{
tau(i,j) = rho*nu*(dudx(i,j) + dudx(j,i));
}
}
tau.AddMult(normal, v);
}
}
};
/**
* @brief Navier-stokes time dependent operator
* du/dt = - div(u u) + nu * laplacian(u) - grad(p) + f
* div(u) = 0
*
* M du/dt = [ -N(u) + L, -G][u]
* 0 = [ D , 0][p]
*/
class NavierStokes : public Application
{
protected:
class NewtonResidual : public TimeDependentOperator
{
protected:
NavierStokes *app;
const Vector *u = nullptr; ///< Pointer to the input vector, used in Mult() methods
mutable future::FDJacobian grad;
real_t dt= 0.0; ///< Time step size, used for time-dependent applications
mutable Vector upk; ///< Temporary vector for u + dt*k
public:
NewtonResidual(NavierStokes *app_) : TimeDependentOperator(app_->Width()),
app(app_), grad(*this,1e-6), upk(Width()) {}
virtual void SetTimeStep(real_t dt_){ dt = dt_;}
virtual void SetState(const Vector *u_){u = u_;}
virtual void Mult(const Vector &k, Vector &y) const override
{
add(1.0,*u, dt, k, upk); // upk = u + dt*k
app->ImplicitMult(upk, k, y); // y = f(upk,k,t)
}
Operator& GetGradient(const Vector &k) const override {
// grad.Update(k);
// return const_cast<future::FDJacobian&>(grad);
add(1.0,*u, dt, k, upk); // upk = u + dt*k
return app->GetGradient(upk);
}
};
public:
// Mesh and finite element spaces
ParMesh &mesh;
ParFiniteElementSpace &u_fes, &p_fes;
Array<int> offsets; ///< Offsets for the velocity and pressure spaces
/// Essential true dof array. Relevant for eliminating boundary conditions
Array<int> u_ess_attr, p_ess_attr;
Array<int> u_ess_tdofs, p_ess_tdofs;
ConstantCoefficient viscosity, density, inv_density, compressibility; // Fluid properties
mutable ParGridFunction u_gf, p_gf; ///< Grid function for the velocity and pressure
mutable ParGridFunction u_gf_bc, p_gf_bc; ///< Grid functions for enforcing/transfering BCs
mutable ParBilinearForm Mpform, Muform, Kuform; ///< Mass and Stiffness forms
mutable ParMixedBilinearForm Gform, Dform; ///< Divergence and Gradient forms
mutable ParNonlinearForm Nform; ///< Nonlinear form for the Navier-Stokes equations
mutable ParNonlinearForm Nform_e; ///< Gradient for nonlinear form for the Navier-Stokes equations
/// Mass matrix and implicit solver
CGSolver M_solver;
FGMRESSolver linear_solver;
NewtonSolver newton_solver;
mutable OperatorHandle Mumat, Mpmat;
mutable OperatorHandle Kumat, Dmat, Gmat;
mutable OperatorHandle Kumat_e, Dmat_e, Gmat_e, Mpmat_e;
BlockOperator *implicit_op = nullptr; ///< Block operator for the Navier-Stokes system
BlockOperator *implicit_grad = nullptr; ///< Block linear operator for the Navier-Stokes gradient
NewtonResidual *newton_residual = nullptr; ///< Newton residual operator
/// Mass matrix preconditioner
HypreSmoother M_prec;
mutable BlockLowerTriangularPreconditioner *pc = nullptr;
mutable Solver *Nu_pc = nullptr;
mutable Vector z, pz, uz;
ConstantCoefficient zero, one; ///< Zero coefficient for boundary conditions
ConstantCoefficient inv_dtc;
IntegrationRules intrules;
IntegrationRule ir, ir_nl;
DeviatoricStressCoefficient stress_coeff;
ParGridFunction stress_gf;
VectorCoefficient *ale_velocity; ///< ALE velocity for moving mesh problems
real_t current_dt = -1.0; ///< Current time in the simulation
bool updated = false;
public:
NavierStokes(ParFiniteElementSpace &u_fes_,
ParFiniteElementSpace &p_fes_,
Array<int> u_ess_attr,
Array<int> p_ess_attr,
real_t density_ = 1.0,
real_t viscosity_ = 1.0,
real_t compressibility_ = 1.0,
bool scaled_pressure = true,
VectorCoefficient *ale_velocity_ = nullptr) :
Application(u_fes_.GetTrueVSize()+p_fes_.GetTrueVSize()),
mesh(*p_fes_.GetParMesh()),
u_fes(u_fes_), p_fes(p_fes_),
u_ess_attr(u_ess_attr), p_ess_attr(p_ess_attr),
viscosity(viscosity_), density(density_),
inv_density(scaled_pressure ? 1.0 : 1.0/density_),
compressibility(compressibility_),
u_gf(&u_fes), p_gf(&p_fes),
u_gf_bc(&u_fes), p_gf_bc(&p_fes),
Mpform(&p_fes), Muform(&u_fes),
Kuform(&u_fes), Gform(&p_fes, &u_fes),
Dform(&u_fes, &p_fes),
Nform(&u_fes), Nform_e(&u_fes),
linear_solver(mesh.GetComm()),
newton_solver(mesh.GetComm()),
zero(0.0), one(1.0), inv_dtc(1.0),
stress_coeff(&p_gf, &u_gf, &viscosity, &density, scaled_pressure),
stress_gf(&u_fes), ale_velocity(ale_velocity_)
{
u_fes.GetEssentialTrueDofs(u_ess_attr, u_ess_tdofs);
p_fes.GetEssentialTrueDofs(p_ess_attr, p_ess_tdofs);
u_gf = 0.0;
p_gf = 0.0;
stress_gf = 0.0;
u_gf_bc = 0.0;
p_gf_bc = 0.0;
// Setup field collection for output and transfer
field_collection.SetName("Navier-Stokes");
field_collection.AddField("Velocity", &u_gf);
field_collection.AddField("Pressure", &p_gf);
field_collection.AddSourceField("Stress", &stress_gf);
field_collection.AddField("Velocity_BC", &u_gf_bc);
field_collection.AddField("Pressure_BC", &p_gf_bc);
offsets = Array<int>({0, u_fes.GetTrueVSize(), p_fes.GetTrueVSize()});
offsets.PartialSum();
auto geom_type = u_fes.GetFE(0)->GetGeomType();
ir = intrules.Get(geom_type, (int)(2*(u_fes.GetOrder(0)+1) - 3));
ir_nl = intrules.Get(geom_type, (int)(ceil(1.5 * 2*(u_fes.GetOrder(0)+1) - 3)));
Mpform.AddDomainIntegrator(new MassIntegrator(compressibility,&ir));
Muform.AddDomainIntegrator(new VectorMassIntegrator(&ir));
// Kuform.AddDomainIntegrator(new VectorDiffusionIntegrator(viscosity, &ir));
Gform.AddDomainIntegrator(new GradientIntegrator(inv_density, &ir));
Dform.AddDomainIntegrator(new VectorDivergenceIntegrator(&ir));
Nform.AddDomainIntegrator(new VectorConvectionNLFIntegrator(one, *ale_velocity, &ir_nl));
Nform.AddDomainIntegrator(new VectorDiffusionIntegrator(viscosity, &ir));
Nform_e.AddDomainIntegrator(new VectorConvectionNLFIntegrator(one, *ale_velocity, &ir_nl));
Nform_e.AddDomainIntegrator(new VectorDiffusionIntegrator(viscosity, &ir));
Nform_e.AddDomainIntegrator(new VectorMassIntegrator(inv_dtc,&ir));
// Muform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
Mpform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
Gform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
Dform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
Nform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
Assemble();
BuildSolvers();
}
void Assemble() override
{
AssembleLinearForms();
AssembleNonlinearForm();
AssembleBilinearForms();
}
void AssembleLinearForms() {}
void AssembleBilinearForms()
{
Mpform.Assemble();
Muform.Assemble();
// Kuform.Assemble();
Dform.Assemble();
Gform.Assemble();
// Muform.FormSystemMatrix(u_ess_tdofs, Mumat);
Mpform.FormSystemMatrix(p_ess_tdofs, Mpmat_e);
Gform.FormRectangularSystemMatrix(p_ess_tdofs, u_ess_tdofs, Gmat_e);
Dform.FormRectangularSystemMatrix(u_ess_tdofs, p_ess_tdofs, Dmat_e);
Array<int> empty;
Mpform.FormSystemMatrix(empty, Mpmat);
Muform.FormSystemMatrix(empty, Mumat);
Gform.FormRectangularSystemMatrix(empty, empty, Gmat);
Dform.FormRectangularSystemMatrix(empty, empty, Dmat);
if(!implicit_op) implicit_op = new BlockOperator(offsets);
implicit_op->SetBlock(0, 0, &Nform);
implicit_op->SetBlock(0, 1, Gmat.Ptr());
implicit_op->SetBlock(1, 0, Dmat.Ptr());
if(!implicit_grad) implicit_grad = new BlockOperator(offsets);
implicit_grad->SetBlock(0, 1, Gmat_e.Ptr());
implicit_grad->SetBlock(1, 0, Dmat_e.Ptr());
implicit_grad->SetBlock(1, 1, Mpmat_e.Ptr());
}
void AssembleNonlinearForm()
{
Nform.Setup();
Nform_e.SetEssentialTrueDofs(u_ess_tdofs);
Nform_e.Setup();
}
void Update() override
{
u_fes.Update();
p_fes.Update();
u_gf.Update();
p_gf.Update();
stress_gf.Update();
p_gf_bc.Update();
u_gf_bc.Update();
Mpform.Update();
Muform.Update();
// Kuform.Update();
Dform.Update();
Gform.Update();
Nform.Update();
Nform_e.Update();
Assemble();
updated = true;
}
void BuildSolvers()
{
M_solver.iterative_mode = false;
M_solver.SetRelTol(1e-8);
M_solver.SetMaxIter(100);
M_prec.SetType(HypreSmoother::Jacobi);
M_solver.SetPreconditioner(M_prec);
M_solver.SetOperator(*Mumat.Ptr());
linear_solver.iterative_mode = false;
linear_solver.SetRelTol(1e-4);
linear_solver.SetAbsTol(1e-4);
linear_solver.SetMaxIter(500);
linear_solver.SetKDim(300);
linear_solver.SetPrintLevel(0);
// pc = new BlockLowerTriangularPreconditioner(offsets);
// pc->SetBlock(1, 0, Dmat_e.Ptr());
// linear_solver.SetPreconditioner(*pc);
if(newton_residual) delete newton_residual;
newton_residual = new NewtonResidual(this);
newton_solver.iterative_mode = true;
newton_solver.SetRelTol(0.0);
newton_solver.SetAbsTol(1e-4);
newton_solver.SetMaxIter(30);
newton_solver.SetPrintLevel(0);
newton_solver.SetOperator(*newton_residual);
newton_solver.SetPreconditioner(linear_solver);
}
void UpdatePreconditioner() const
{
if(Nu_pc) delete Nu_pc;
auto amg = new HypreBoomerAMG;
HYPRE_BoomerAMGSetSmoothType(*amg, 5);
amg->SetOperator(*static_cast<HypreParMatrix*>(&implicit_grad->GetBlock(0, 0)));
amg->SetSystemsOptions(2, true);
amg->SetPrintLevel(0);
Nu_pc = amg;
pc->SetBlock(0, 0, Nu_pc);
}
Operator& GetGradient(const Vector &x) const override
{
BlockVector xb(x.GetData(), offsets);
implicit_grad->SetBlock(0, 0, &Nform_e.GetGradient(xb.GetBlock(0)));
implicit_grad->SetBlockCoef(0,0,current_dt);
implicit_grad->SetBlockCoef(0,1,current_dt);
implicit_grad->SetBlockCoef(1,0,current_dt);
if(pc) UpdatePreconditioner();
return *implicit_grad;
}
void Mult(const Vector &x, Vector &k) const override
{
BlockVector xb(x.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
implicit_op->Mult(xb, kb);
M_solver.Mult(kb.GetBlock(0), z);
kb.GetBlock(0) = z;
kb.GetBlock(0).SetSubVector(u_ess_tdofs, 0.0);
kb.GetBlock(1).SetSubVector(p_ess_tdofs, 0.0);
}
void ImplicitSolve(const real_t dt, const Vector &u, Vector &k) override
{
BlockVector ub(u.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
if((current_dt != dt) || updated)
{
inv_dtc.constant = 1.0/dt;
current_dt = dt;
AssembleNonlinearForm();
updated = false;
}
if(IsCoupled())
{ // Enforce velocity BCs through initial guess
Vector &ku = kb.GetBlock(0);
u_gf_bc.GetTrueDofs(z);
for (int i = 0; i < u_ess_tdofs.Size(); i++)
{
int idx = u_ess_tdofs[i];
ku(idx) = z(idx);
}
}
else
{
kb.GetBlock(0) = 0.0;
kb.GetBlock(1) = 0.0;
}
Vector zero_vec;
newton_residual->SetTimeStep(dt);
newton_residual->SetState(&u);
newton_solver.Mult(zero_vec, k); // Solve the nonlinear system
}
// Computes residual of the Navier-Stokes equations
void ImplicitMult(const Vector &x, const Vector &k, Vector &v ) const override
{
BlockVector xb(x.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
BlockVector vb(v.GetData(), offsets);
// Apply the implicit operator
implicit_op->Mult(xb, vb); // v = A(x)
Mumat->AddMult(kb.GetBlock(0), vb.GetBlock(0), 1.0); // v = A(x) + M*k
Mpmat->AddMult(kb.GetBlock(1), vb.GetBlock(1), 1.0); // v = A(x) + M*k
vb.GetBlock(0).SetSubVector(u_ess_tdofs, 0.0);
vb.GetBlock(1).SetSubVector(p_ess_tdofs, 0.0);
}
void PreProcess(Vector &x) override {}
void PostProcess(Vector &x) override {}
void Transfer(const Vector &x) override
{
BlockVector xb(x.GetData(), offsets);
u_gf.SetFromTrueDofs(xb.GetBlock(0));
p_gf.SetFromTrueDofs(xb.GetBlock(1));
stress_gf.ProjectBdrCoefficient(stress_coeff,u_ess_attr);
field_collection.Transfer("Stress");
}
void Transfer(const Vector &u, const Vector &k, real_t dt = 0.0) override
{
BlockVector ub(u.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
uz.SetSize(u_fes.GetTrueVSize());
pz.SetSize(p_fes.GetTrueVSize());
// compute stage-updated velocity and pressure
add(1.0,ub.GetBlock(0), dt, kb.GetBlock(0), uz);
add(1.0,ub.GetBlock(1), dt, kb.GetBlock(1), pz);
u_gf.SetFromTrueDofs(uz);
p_gf.SetFromTrueDofs(pz);
stress_gf.ProjectBdrCoefficient(stress_coeff,u_ess_attr);
field_collection.Transfer("Stress");
}
~NavierStokes()
{
if(implicit_op) delete implicit_op;
if(implicit_grad) delete implicit_grad;
if(newton_residual) delete newton_residual;
if(Nu_pc) delete Nu_pc;
if(pc) delete pc;
}
};
#endif
@@ -0,0 +1,57 @@
SetFactory("OpenCASCADE");
Mesh.ElementOrder = 2;
Mesh.SecondOrderLinear = 0;
lc = 1.0;
// channel height
H = 1.0;
// step height
h = H / 2.0;
// slab thickness
b = 4.0 * h;
// channel length
L = 30.0 * H;
Point(1) = {0, 0, 0, lc};
Point(2) = {L, 0, 0, lc};
Point(3) = {L, b, 0, lc};
Point(4) = {L, b+h, 0, lc};
Point(5) = {L, b+H, 0, lc};
Point(6) = {0, b+H, 0, lc};
Point(7) = {0, b+h, 0, lc};
Point(8) = {0, b, 0, lc};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 1};
Line(9) = {8, 3};
Line(10) = {7, 4};
Curve Loop(1) = {10, 4, 5, 6};
Plane Surface(1) = {1};
Curve Loop(2) = {9, 3, -10, 7};
Plane Surface(2) = {2};
Curve Loop(3) = {9, -2, -1, -8};
Plane Surface(3) = {3};
Physical Surface("fluid", 1) = {1, 2};
Physical Surface("solid", 2) = {3};
Physical Curve("inlet", 1) = {6};
Physical Curve("outlet", 2) = {3, 4};
Physical Curve("interface", 3) = {9};
Physical Curve("wall fluid", 4) = {7, 5};
Physical Curve("wall side", 5) = {8, 2};
Physical Curve("wall bottom", 6) = {1};
Transfinite Surface {1:3};
Recombine Surface {1:3};
Transfinite Curve {1, -5, 10, 9} = 8 Using Progression 1;
Transfinite Curve {6, 4, 7, 3} = 1 Using Progression 1;
Transfinite Curve {8, 2} = 2 Using Progression 1;
+173
View File
@@ -0,0 +1,173 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
8
1 1 "inlet"
1 2 "outlet"
1 3 "interface"
1 4 "wall fluid"
1 5 "wall side"
1 6 "wall bottom"
2 1 "fluid"
2 2 "solid"
$EndPhysicalNames
$Nodes
105
1 0 0 0
2 30 0 0
3 30 2 0
4 30 2.5 0
5 30 3 0
6 0 3 0
7 0 2.5 0
8 0 2 0
9 4.285714285714286 0 0
10 8.571428571428571 0 0
11 12.85714285714286 0 0
12 17.14285714285714 0 0
13 21.42857142857143 0 0
14 25.71428571428571 0 0
15 2.142857142857143 0 0
16 6.428571428571429 0 0
17 10.71428571428571 0 0
18 15 0 0
19 19.28571428571428 0 0
20 23.57142857142857 0 0
21 27.85714285714285 0 0
22 30 0.9999999999999999 0
23 30 2.25 0
24 30 2.75 0
25 25.71428571428572 3 0
26 21.42857142857143 3 0
27 17.14285714285715 3 0
28 12.85714285714286 3 0
29 8.571428571428573 3 0
30 4.285714285714288 3 0
31 27.85714285714286 3 0
32 23.57142857142857 3 0
33 19.28571428571428 3 0
34 15 3 0
35 10.71428571428572 3 0
36 6.428571428571431 3 0
37 2.142857142857146 3 0
38 0 2.75 0
39 0 2.25 0
40 0 1 0
41 4.285714285714286 2 0
42 8.571428571428571 2 0
43 12.85714285714286 2 0
44 17.14285714285714 2 0
45 21.42857142857143 2 0
46 25.71428571428571 2 0
47 2.142857142857143 2 0
48 6.428571428571429 2 0
49 10.71428571428571 2 0
50 15 2 0
51 19.28571428571428 2 0
52 23.57142857142857 2 0
53 27.85714285714285 2 0
54 4.285714285714286 2.5 0
55 8.571428571428571 2.5 0
56 12.85714285714286 2.5 0
57 17.14285714285714 2.5 0
58 21.42857142857143 2.5 0
59 25.71428571428571 2.5 0
60 2.142857142857143 2.5 0
61 6.428571428571429 2.5 0
62 10.71428571428571 2.5 0
63 15 2.5 0
64 19.28571428571428 2.5 0
65 23.57142857142857 2.5 0
66 27.85714285714285 2.5 0
67 4.285714285714286 2.75 0
68 2.142857142857146 2.75 0
69 8.571428571428573 2.75 0
70 6.428571428571431 2.75 0
71 12.85714285714286 2.75 0
72 10.71428571428572 2.75 0
73 17.14285714285714 2.75 0
74 15 2.75 0
75 21.42857142857143 2.75 0
76 19.28571428571428 2.75 0
77 25.71428571428572 2.75 0
78 23.57142857142857 2.75 0
79 27.85714285714286 2.75 0
80 4.285714285714285 2.25 0
81 2.142857142857142 2.25 0
82 8.571428571428571 2.25 0
83 6.428571428571431 2.25 0
84 12.85714285714286 2.25 0
85 10.71428571428571 2.25 0
86 17.14285714285714 2.25 0
87 15 2.25 0
88 21.42857142857143 2.25 0
89 19.28571428571428 2.25 0
90 25.71428571428571 2.25 0
91 23.57142857142857 2.25 0
92 27.85714285714285 2.25 0
93 4.285714285714285 1 0
94 2.142857142857142 1 0
95 8.571428571428571 1 0
96 6.428571428571431 1 0
97 12.85714285714286 1 0
98 10.71428571428571 1 0
99 17.14285714285714 1 0
100 15 1 0
101 21.42857142857143 1 0
102 19.28571428571428 1 0
103 25.71428571428571 1 0
104 23.57142857142857 1 0
105 27.85714285714285 1 0
$EndNodes
$Elements
48
1 8 2 6 1 1 9 15
2 8 2 6 1 9 10 16
3 8 2 6 1 10 11 17
4 8 2 6 1 11 12 18
5 8 2 6 1 12 13 19
6 8 2 6 1 13 14 20
7 8 2 6 1 14 2 21
8 8 2 5 2 2 3 22
9 8 2 2 3 3 4 23
10 8 2 2 4 4 5 24
11 8 2 4 5 5 25 31
12 8 2 4 5 25 26 32
13 8 2 4 5 26 27 33
14 8 2 4 5 27 28 34
15 8 2 4 5 28 29 35
16 8 2 4 5 29 30 36
17 8 2 4 5 30 6 37
18 8 2 1 6 6 7 38
19 8 2 4 7 7 8 39
20 8 2 5 8 8 1 40
21 8 2 3 9 8 41 47
22 8 2 3 9 41 42 48
23 8 2 3 9 42 43 49
24 8 2 3 9 43 44 50
25 8 2 3 9 44 45 51
26 8 2 3 9 45 46 52
27 8 2 3 9 46 3 53
28 10 2 1 1 7 54 30 6 60 67 37 38 68
29 10 2 1 1 54 55 29 30 61 69 36 67 70
30 10 2 1 1 55 56 28 29 62 71 35 69 72
31 10 2 1 1 56 57 27 28 63 73 34 71 74
32 10 2 1 1 57 58 26 27 64 75 33 73 76
33 10 2 1 1 58 59 25 26 65 77 32 75 78
34 10 2 1 1 59 4 5 25 66 24 31 77 79
35 10 2 1 2 8 41 54 7 47 80 60 39 81
36 10 2 1 2 41 42 55 54 48 82 61 80 83
37 10 2 1 2 42 43 56 55 49 84 62 82 85
38 10 2 1 2 43 44 57 56 50 86 63 84 87
39 10 2 1 2 44 45 58 57 51 88 64 86 89
40 10 2 1 2 45 46 59 58 52 90 65 88 91
41 10 2 1 2 46 3 4 59 53 23 66 90 92
42 10 2 2 3 8 41 9 1 47 93 15 40 94
43 10 2 2 3 41 42 10 9 48 95 16 93 96
44 10 2 2 3 42 43 11 10 49 97 17 95 98
45 10 2 2 3 43 44 12 11 50 99 18 97 100
46 10 2 2 3 44 45 13 12 51 101 19 99 102
47 10 2 2 3 45 46 14 13 52 103 20 101 104
48 10 2 2 3 46 3 2 14 53 22 21 103 105
$EndElements
+395
View File
@@ -0,0 +1,395 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
9
1 1 "inlet"
1 2 "outlet"
1 3 "top wall"
1 4 "bottom wall"
1 5 "beam wall"
1 6 "cylinder"
1 7 "beam cylinder curve"
2 1 "fluid"
2 2 "beam"
$EndPhysicalNames
$Nodes
268
1 0 0 0
2 2.2 0 0
3 2.2 0.41 0
4 0 0.41 0
5 0.2353553390593274 0.2353553390593274 0
6 0.2353553390593274 0.1646446609406726 0
7 0.1646446609406726 0.2353553390593274 0
8 0.1646446609406726 0.1646446609406726 0
9 0.2 0.25 0
10 0.15 0.2 0
11 0.2 0.15 0
12 0.2489897948556636 0.19 0
13 0.2489897948556636 0.21 0
14 0.3310660171779822 0.3060660171779822 0
15 0.3310660171779822 0.09393398282201787 0
16 0.06893398282201785 0.3060660171779822 0
17 0.06893398282201785 0.09393398282201787 0
18 0.2 0.325 0
19 0.04749999999999999 0.2 0
20 0.2 0.07499999999999998 0
21 0.6 0.21 0
22 0.6 0.19 0
23 0.3525 0.21 0
24 0.3525 0.19 0
25 0.2 0.41 0
26 0.2 0 0
27 0.06893398282201785 0.41 0
28 0.06893398282201785 0 0
29 0.3310660171779822 0 0
30 0.3310660171779822 0.41 0
31 0.6 0.41 0
32 0.6 0 0
33 0.6 0.3075 0
34 0.6 0.1025 0
35 2.2 0.3075 0
36 2.2 0.1025 0
37 2.2 0.21 0
38 2.2 0.19 0
39 0 0.2 0
40 0 0.3060660171779822 0
41 0 0.09393398282201787 0
42 0.1808658283817455 0.2461939766255643 0
43 0.1538060233744357 0.2191341716182545 0
44 0.1538060233744357 0.1808658283817456 0
45 0.1808658283817455 0.1538060233744357 0
46 0.2191341716182545 0.1538060233744357 0
47 0.2440369016676673 0.1763198122576657 0
48 0.2440369016676674 0.2236801877423343 0
49 0.2191341716182545 0.2461939766255643 0
50 0.2 0.2875 0
51 0.2 0.3675 0
52 0.1344669914110089 0.41 0
53 0.06893398282201785 0.3580330085889911 0
54 0.1167893218813452 0.2707106781186548 0
55 0.1344669914110089 0.3155330085889911 0
56 0.2655330085889911 0.3155330085889911 0
57 0.3310660171779822 0.3580330085889911 0
58 0.2655330085889911 0.41 0
59 0.2832106781186548 0.2707106781186548 0
60 0.5103553390593271 0.41 0
61 0.4207106781186545 0.41 0
62 0.5551776695296635 0.41 0
63 0.4655330085889908 0.41 0
64 0.3758883476483184 0.41 0
65 0.03446699141100892 0.41 0
66 0 0.3580330085889911 0
67 0.03446699141100892 0.3060660171779822 0
68 0.05821699141100892 0.2530330085889911 0
69 0.09875 0.2 0
70 0.02374999999999999 0.2 0
71 0 0.2530330085889911 0
72 0 0.1469669914110089 0
73 0.05821699141100892 0.1469669914110089 0
74 0.03446699141100892 0.09393398282201787 0
75 0 0.04696699141100893 0
76 0.06893398282201785 0.04696699141100893 0
77 0.03446699141100892 0 0
78 0.1344669914110089 0 0
79 0.2 0.03749999999999999 0
80 0.1344669914110089 0.08446699141100893 0
81 0.1167893218813452 0.1292893218813453 0
82 0.2 0.1125 0
83 0.2832106781186548 0.1292893218813453 0
84 0.2655330085889911 0.08446699141100893 0
85 0.2655330085889911 0 0
86 0.3310660171779822 0.04696699141100893 0
87 0.3417830085889911 0.1419669914110089 0
88 0.3007448974278318 0.19 0
89 0.3007448974278318 0.21 0
90 0.3417830085889911 0.2580330085889911 0
91 0.4207106781186546 0.3065440114519881 0
92 0.5103553390593272 0.307022005725994 0
93 0.3758883476483184 0.3063050143149851 0
94 0.4655330085889909 0.306783008588991 0
95 0.5551776695296636 0.307261002862997 0
96 0.6 0.35875 0
97 0.6 0.25875 0
98 0.5174999999999998 0.21 0
99 0.4349999999999998 0.21 0
100 0.5587499999999999 0.21 0
101 0.4762499999999998 0.21 0
102 0.3937499999999999 0.21 0
103 0.4350000000000002 0.19 0
104 0.5175000000000003 0.19 0
105 0.3937500000000001 0.19 0
106 0.4762500000000002 0.19 0
107 0.5587500000000001 0.19 0
108 0.6 0.2 0
109 0.6 0.14625 0
110 0.5103553390593273 0.09964466094067262 0
111 0.4207106781186548 0.09678932188134524 0
112 0.5551776695296636 0.1010723304703363 0
113 0.465533008588991 0.09821699141100892 0
114 0.3758883476483185 0.09536165235168155 0
115 0.420710678118655 0 0
116 0.5103553390593276 0 0
117 0.3758883476483186 0 0
118 0.4655330085889914 0 0
119 0.5551776695296639 0 0
120 0.6 0.05125 0
121 0.7769304775334696 0.41 0
122 1.00694009683966 0.41 0
123 1.305952603331686 0.41 0
124 1.69466886165541 0.41 0
125 0.6884652387667348 0.41 0
126 0.8919352871865647 0.41 0
127 1.156446350085673 0.41 0
128 1.500310732493548 0.41 0
129 1.947334430827705 0.41 0
130 2.2 0.35875 0
131 2.2 0.05125 0
132 0.7769304775334696 0 0
133 1.00694009683966 0 0
134 1.305952603331686 0 0
135 1.69466886165541 0 0
136 0.6884652387667348 0 0
137 0.8919352871865647 0 0
138 1.156446350085673 0 0
139 1.500310732493548 0 0
140 1.947334430827705 0 0
141 0.7769304775334696 0.1025 0
142 1.00694009683966 0.1025 0
143 1.305952603331686 0.1025 0
144 1.69466886165541 0.1025 0
145 0.6884652387667348 0.1025 0
146 0.8919352871865647 0.1025 0
147 1.156446350085673 0.1025 0
148 1.500310732493548 0.1025 0
149 1.947334430827705 0.1025 0
150 0.7769304775334696 0.3075 0
151 1.00694009683966 0.3075 0
152 1.305952603331686 0.3075 0
153 1.69466886165541 0.3075 0
154 0.6884652387667348 0.3075 0
155 0.8919352871865647 0.3075 0
156 1.156446350085673 0.3075 0
157 1.500310732493548 0.3075 0
158 1.947334430827705 0.3075 0
159 0.7769304775334696 0.21 0
160 1.00694009683966 0.21 0
161 1.305952603331686 0.21 0
162 1.69466886165541 0.21 0
163 0.6884652387667348 0.21 0
164 0.8919352871865647 0.21 0
165 1.156446350085673 0.21 0
166 1.500310732493548 0.21 0
167 1.947334430827705 0.21 0
168 0.7769304775334696 0.19 0
169 1.00694009683966 0.19 0
170 1.305952603331686 0.19 0
171 1.69466886165541 0.19 0
172 0.6884652387667348 0.19 0
173 0.8919352871865647 0.19 0
174 1.156446350085673 0.19 0
175 1.500310732493548 0.19 0
176 1.947334430827705 0.19 0
177 0.3525 0.2 0
178 2.2 0.2 0
179 2.2 0.25875 0
180 2.2 0.14625 0
181 0.25 0.2 0
182 0.30125 0.2 0
183 0.2929099551283292 0.2408565981656627 0
184 0.4278553390593272 0.2582720057259941 0
185 0.3848191738241591 0.2581525071574926 0
186 0.5139276695296635 0.2585110028629971 0
187 0.4708915042944953 0.2583915042944955 0
188 0.5569638347648318 0.2586305014314985 0
189 0.5175000000000001 0.2 0
190 0.55875 0.2 0
191 0.435 0.2 0
192 0.47625 0.2 0
193 0.39375 0.2 0
194 0.2929099551283292 0.1591434018343373 0
195 0.2423335901036228 0.2808634926072777 0
196 0.2655330085889911 0.3627665042944955 0
197 0.5103553390593272 0.358511002862997 0
198 0.5551776695296635 0.3586305014314985 0
199 0.4207106781186545 0.358272005725994 0
200 0.4655330085889908 0.3583915042944955 0
201 0.3758883476483184 0.3581525071574925 0
202 0.1344669914110089 0.3627665042944955 0
203 0.1576664098963772 0.2808634926072777 0
204 0.03446699141100892 0.3580330085889911 0
205 0.02910849570550446 0.2530330085889911 0
206 0.1060115073927223 0.2360835901036228 0
207 0.1060115073927223 0.1639164098963772 0
208 0.03446699141100892 0.04696699141100893 0
209 0.1576664098963772 0.1191365073927223 0
210 0.1344669914110089 0.04223349570550446 0
211 0.2423335901036228 0.1191365073927223 0
212 0.2655330085889911 0.04223349570550446 0
213 0.4278553390593275 0.1433946609406726 0
214 0.3848191738241593 0.1426808261758408 0
215 0.5139276695296637 0.1448223304703363 0
216 0.4708915042944956 0.1441084957055045 0
217 0.5569638347648319 0.1455361652351682 0
218 0.5103553390593275 0.04982233047033631 0
219 0.5551776695296637 0.05053616523516816 0
220 0.4207106781186549 0.04839466094067262 0
221 0.4655330085889912 0.04910849570550446 0
222 0.3758883476483185 0.04768082617584078 0
223 0.7769304775334696 0.2 0
224 0.6884652387667349 0.2 0
225 1.00694009683966 0.2 0
226 0.8919352871865647 0.2 0
227 1.305952603331686 0.2 0
228 1.156446350085673 0.2 0
229 1.69466886165541 0.2 0
230 1.500310732493548 0.2 0
231 1.947334430827705 0.2 0
232 1.69466886165541 0.25875 0
233 1.947334430827705 0.25875 0
234 1.305952603331686 0.25875 0
235 1.500310732493548 0.25875 0
236 1.00694009683966 0.25875 0
237 1.156446350085673 0.25875 0
238 0.7769304775334696 0.25875 0
239 0.8919352871865647 0.25875 0
240 0.6884652387667347 0.25875 0
241 0.7769304775334696 0.35875 0
242 0.6884652387667349 0.35875 0
243 1.00694009683966 0.35875 0
244 0.8919352871865647 0.35875 0
245 1.305952603331686 0.35875 0
246 1.156446350085673 0.35875 0
247 1.69466886165541 0.35875 0
248 1.500310732493548 0.35875 0
249 1.947334430827705 0.35875 0
250 1.69466886165541 0.14625 0
251 1.947334430827705 0.14625 0
252 1.305952603331686 0.14625 0
253 1.500310732493548 0.14625 0
254 1.00694009683966 0.14625 0
255 1.156446350085673 0.14625 0
256 0.7769304775334696 0.14625 0
257 0.8919352871865647 0.14625 0
258 0.6884652387667349 0.14625 0
259 0.7769304775334696 0.05125 0
260 0.6884652387667349 0.05125 0
261 1.00694009683966 0.05125 0
262 0.8919352871865647 0.05125 0
263 1.305952603331686 0.05125 0
264 1.156446350085673 0.05125 0
265 1.69466886165541 0.05125 0
266 1.500310732493548 0.05125 0
267 1.947334430827705 0.05125 0
268 0.02910849570550446 0.1469669914110089 0
$EndNodes
$Elements
106
1 8 2 6 1 9 7 42
2 8 2 6 2 7 10 43
3 8 2 6 3 10 8 44
4 8 2 6 4 8 11 45
5 8 2 6 5 11 6 46
6 8 2 6 6 6 12 47
7 8 2 6 9 13 5 48
8 8 2 6 10 5 9 49
9 8 2 3 13 25 27 52
10 8 2 3 19 30 25 58
11 8 2 3 21 31 60 62
12 8 2 3 21 60 61 63
13 8 2 3 21 61 30 64
14 8 2 3 22 27 4 65
15 8 2 1 23 4 40 66
16 8 2 1 28 39 40 71
17 8 2 1 29 41 39 72
18 8 2 1 32 41 1 75
19 8 2 4 34 28 1 77
20 8 2 4 35 28 26 78
21 8 2 4 42 26 29 85
22 8 2 5 45 24 12 88
23 8 2 5 46 13 23 89
24 8 2 5 51 21 98 100
25 8 2 5 51 98 99 101
26 8 2 5 51 99 23 102
27 8 2 5 52 24 103 105
28 8 2 5 52 103 104 106
29 8 2 5 52 104 22 107
30 8 2 5 53 22 21 108
31 8 2 4 56 29 115 117
32 8 2 4 56 115 116 118
33 8 2 4 56 116 32 119
34 8 2 3 58 31 121 125
35 8 2 3 58 121 122 126
36 8 2 3 58 122 123 127
37 8 2 3 58 123 124 128
38 8 2 3 58 124 3 129
39 8 2 2 59 3 35 130
40 8 2 2 61 36 2 131
41 8 2 4 62 32 132 136
42 8 2 4 62 132 133 137
43 8 2 4 62 133 134 138
44 8 2 4 62 134 135 139
45 8 2 4 62 135 2 140
46 8 2 2 68 37 38 178
47 8 2 2 69 35 37 179
48 8 2 2 70 38 36 180
49 8 2 7 71 13 12 181
50 10 2 2 1 13 23 24 12 89 177 88 181 182
51 10 2 1 2 14 5 13 23 59 48 89 90 183
52 10 2 1 3 23 14 91 99 90 93 184 102 185
53 10 2 1 3 99 91 92 98 184 94 186 101 187
54 10 2 1 3 98 92 33 21 186 95 97 100 188
55 10 2 2 4 21 98 104 22 100 189 107 108 190
56 10 2 2 4 98 99 103 104 101 191 106 189 192
57 10 2 2 4 99 23 24 103 102 177 105 191 193
58 10 2 1 5 15 24 12 6 87 88 47 83 194
59 10 2 1 6 9 18 14 5 50 56 59 49 195
60 10 2 1 7 18 25 30 14 51 58 57 56 196
61 10 2 1 8 31 60 92 33 62 197 95 96 198
62 10 2 1 8 60 61 91 92 63 199 94 197 200
63 10 2 1 8 61 30 14 91 64 57 93 199 201
64 10 2 1 9 18 25 27 16 51 52 53 55 202
65 10 2 1 10 16 7 9 18 54 42 50 55 203
66 10 2 1 11 27 16 40 4 53 67 66 65 204
67 10 2 1 12 16 19 39 40 68 70 71 67 205
68 10 2 1 13 19 10 7 16 69 43 54 68 206
69 10 2 1 14 10 8 17 19 44 81 73 69 207
70 10 2 1 16 41 1 28 17 75 77 76 74 208
71 10 2 1 17 17 8 11 20 81 45 82 80 209
72 10 2 1 18 20 17 28 26 80 76 78 79 210
73 10 2 1 19 15 20 11 6 84 82 46 83 211
74 10 2 1 20 15 20 26 29 84 79 85 86 212
75 10 2 1 21 15 24 103 111 87 105 213 114 214
76 10 2 1 21 111 103 104 110 213 106 215 113 216
77 10 2 1 21 110 104 22 34 215 107 109 112 217
78 10 2 1 22 34 32 116 110 120 119 218 112 219
79 10 2 1 22 110 116 115 111 218 118 220 113 221
80 10 2 1 22 111 115 29 15 220 117 86 114 222
81 10 2 1 23 21 159 168 22 163 223 172 108 224
82 10 2 1 23 159 160 169 168 164 225 173 223 226
83 10 2 1 23 160 161 170 169 165 227 174 225 228
84 10 2 1 23 161 162 171 170 166 229 175 227 230
85 10 2 1 23 162 37 38 171 167 178 176 229 231
86 10 2 1 24 35 37 162 153 179 167 232 158 233
87 10 2 1 24 153 162 161 152 232 166 234 157 235
88 10 2 1 24 152 161 160 151 234 165 236 156 237
89 10 2 1 24 151 160 159 150 236 164 238 155 239
90 10 2 1 24 150 159 21 33 238 163 97 154 240
91 10 2 1 25 33 31 121 150 96 125 241 154 242
92 10 2 1 25 150 121 122 151 241 126 243 155 244
93 10 2 1 25 151 122 123 152 243 127 245 156 246
94 10 2 1 25 152 123 124 153 245 128 247 157 248
95 10 2 1 25 153 124 3 35 247 129 130 158 249
96 10 2 1 26 38 36 144 171 180 149 250 176 251
97 10 2 1 26 171 144 143 170 250 148 252 175 253
98 10 2 1 26 170 143 142 169 252 147 254 174 255
99 10 2 1 26 169 142 141 168 254 146 256 173 257
100 10 2 1 26 168 141 34 22 256 145 109 172 258
101 10 2 1 27 32 132 141 34 136 259 145 120 260
102 10 2 1 27 132 133 142 141 137 261 146 259 262
103 10 2 1 27 133 134 143 142 138 263 147 261 264
104 10 2 1 27 134 135 144 143 139 265 148 263 266
105 10 2 1 27 135 2 36 144 140 131 149 265 267
106 10 2 1 28 17 19 39 41 73 70 72 74 268
$EndElements
+960
View File
@@ -0,0 +1,960 @@
/**
* Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
* at the Lawrence Livermore National Laboratory. All Rights reserved. See files
* LICENSE and NOTICE for details. LLNL-CODE-806117.
*
* This file is part of the MFEM library. For more information and source code
* availability visit https://mfem.org.
*
* MFEM is free software; you can redistribute it and/or modify it under the
* terms of the BSD-3 license. We welcome feedback and contributions, see file
* CONTRIBUTING.md for details.
*
* --------------------------------------------------
* Conjugate heat transfer miniapp
* --------------------------------------------------
*
* This is a miniapp demonstrates conjugate heat transfer by coupling different
* physics in different domains:
* 1) Incompressible Navier-Stokes equations in a fluid domain
* 2) Heat equation in fluid
* 3) Heat equation solid domains
*
* The test case is a benchmark Backward Facing Step (BFS) with a heated base.
* The following boundary conditions are applied:
* 1) Fluid inlet (attribute 1): parabolic velocity profile, T = 0.0
* 2) Fluid outlet (attribute 2): zero-pressure, -kappa grad(T)n = 0.0
* 3) Fluid walls (attribute 3 & 4): no-slip, -kappa grad(T)n = 0.0
* 4) Solid base (attribute 6): T = 1.0
* 5) Solid walls (attribute 5): -kappa grad(T)n = 0.0
*
* 4
* -----------------------------------------------------------------
* 1 | fluid | 2
* 4 | | 2
* ---------------------------- 3 ----------------------------------
* 5 | | 5
* 5 | solid | 5
* -----------------------------------------------------------------
* 6
*
* This example demonstrates nested coupling with the fluid flow and heat
* transfer solvers coupling where the heat transfer solver is itself a
* coupled solver with the fluid and solid heat transfer solvers.
*
* The fluid flow and heat transfer solvers are flow maps (i.e ODE Solvers) and
* coupled with one-way, serial or parallel, coupling. The fluid and solid heat
* transfer solvers can be partitioned-coupled or monolithically-coupled.
*
* For partitioned-coupling, the following boundary conditions are applied
* at the fluid-solid interface:
* 1) Dirichlet: T = T_f on attr 3 in solid domain
* 2) Neumann: Q = -kappa grad(T_s)n on attr 3 in fluid domain
*
* For monolithic-coupling, the following conditions are imposed at the
* fluid-solid interface:
* 1) Continuity of temperature: T_f = T_s
* 2) Continuity of heat flux: -k_f grad(T_f)n = k_s grad(T_s)n
*
* Sample run:
* mpirun -np 6 cht-BFS -vs 500 -dt 1e-3 -tf 1000 -rs 2 -o 3 -ode 21 -scheme -1 -cht
*/
#include "mfem.hpp"
#include "multiapp.hpp"
#include "../fluids/navier/navier_solver.hpp"
using namespace mfem;
using namespace navier;
using namespace std;
struct BFSContext
{
int ser_ref = 1; // Serial mesh refinement
int order = 3; // Finite element order
int ode_solver = 21; // ODE solver
real_t dt = 1e-2; // Time step size
real_t t_final = 3.0; // Final time
int vis_steps = 100; // Visualization steps
int couple_scheme = -1; // Coupling scheme
// -1: Monolithic, 0: Alternating Schwarz,
// >0: Additive Schwarz with number of iterations
bool ht_only = false; // Conjugate heat transfer on/off
bool visualization = true;// Visualization on/off
real_t Re = 800.0; // Reynolds number
real_t Pr = 0.71; // Prandtl number
real_t density = 1.0; // Density
real_t kappa_ratio = 1e2; // Conductivity ratio solid/fluid
bool checkres = false; // Check results
#if defined(MFEM_USE_DOUBLE)
real_t tol_T = 1e-4;
real_t tol_Q = 1e-4;
#elif defined(MFEM_USE_SINGLE)
real_t tol_T = 1e-3;
real_t tol_Q = 1e-3;
#else
#error "Only single and double precision are supported!"
real_t tol_T = 0;
real_t tol_Q = 0;
#endif
} ctx;
void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol , int max_it,
int print_level, bool iterative_mode);
/// Temperature profile
double temp_profile(const Vector& x)
{
return x(1) == 0.0 ? 1.0 : 0.0;
}
/// Parabolic velocity profile for channel inlet
void velocity_profile(const Vector &x, Vector &u)
{
double xi = x(0), yi = x(1) - 2.5;
u = 0.0;
u(0) = 24.0*yi*(0.5-yi);
}
/**
* @brief Coefficient to compute normal heat flux Q = -kappa grad(T) n
*/
class HeatFluxCoefficient : public Coefficient
{
protected:
ParGridFunction *T_gf = nullptr;
Coefficient *conductivity = nullptr;
Vector grad_T, normal;
public:
HeatFluxCoefficient(ParGridFunction *T_gf_,
Coefficient *conductivity_) :
Coefficient(), T_gf(T_gf_),
conductivity(conductivity_)
{
int dim = T_gf->FESpace()->GetMesh()->Dimension();
grad_T.SetSize(dim);
normal.SetSize(dim);
}
real_t Eval(ElementTransformation &T,
const IntegrationPoint &ip) override
{
const DenseMatrix &jacobian = T.Jacobian();
MFEM_ASSERT( (jacobian.Height() - 1 == jacobian.Width()),
"Incorrect Jacobian dimension. Coefficient only "
"supported for boundary elements.");
CalcOrtho(jacobian, normal);
const double scale = normal.Norml2();
normal /= scale;
real_t kappa = conductivity ? conductivity->Eval(T, ip) : 1.0;
T_gf->GetGradient(T,grad_T);
real_t flux = -kappa*(grad_T * normal);
return flux;
}
};
/**
* @brief Convection-diffusion time dependent operator
*
* dT/dt = κΔT - αTu
*
* Can also be used to create a diffusion or convection
* only operator by setting α or κ to zero.
*/
class ConvectionDiffusion : public Application
{
public:
// Mesh and finite element space
ParMesh &mesh;
ParFiniteElementSpace &fes;
/// Essential and natural dof array.
Array<int> ess_attr, nat_attr;
Array<int> ess_tdofs, nat_tdofs;
/// Material properties
ConstantCoefficient diffusivity, kappa, alpha;
/// Grid functions for the temperature and heat flux
mutable ParGridFunction T_gf, Q_gf;
/// Used to store boundary condition data
mutable ParGridFunction T_gf_bc, Q_gf_bc;
mutable GridFunctionCoefficient Q_fgc;
mutable HeatFluxCoefficient Q_coeff;
/// Mass form and Stiffness form. Might include
/// diffusion, convection or both.
mutable ParBilinearForm Mform, Kform, Mform_e, Kform_e;
/// RHS form
mutable ParLinearForm bform;
mutable Vector b;
/// Mass and Stiffness operators
mutable HypreParMatrix Mmat, Kmat, Kmat_e, Mmat_e;
/// Mass matrix solver
CGSolver M_solver;
GMRESSolver implicit_solver;
HypreParMatrix *T = nullptr; // T = M + dt K
/// Preconditioners
HypreSmoother M_prec, T_prec;
/// Velocity coefficient
VectorCoefficient *velocity_coeff = nullptr;
/// Auxiliary variables
real_t current_dt = -1.0;
bool updated = false;
mutable Vector z, q;
public:
ConvectionDiffusion(ParFiniteElementSpace &fes_,
Array<int> ess_attr_,
Array<int> nat_attr_,
real_t diffusivity_ = 1.0,
real_t kappa_ = 1.0,
VectorCoefficient *velocity_coeff_ = nullptr,
real_t alpha_ = 1.0)
: Application(fes_.GetTrueVSize()),
mesh(*fes_.GetParMesh()),
fes(fes_),
ess_attr(ess_attr_),
nat_attr(nat_attr_),
diffusivity(diffusivity_),
kappa(kappa_), alpha(-alpha_),
T_gf(&fes), Q_gf(&fes),
T_gf_bc(&fes), Q_gf_bc(&fes),
Q_fgc(&Q_gf_bc), Q_coeff(&T_gf,&kappa),
Mform(&fes), Kform(&fes),
Mform_e(&fes), Kform_e(&fes),
bform(&fes),
M_solver(mesh.GetComm()),
implicit_solver(mesh.GetComm()),
velocity_coeff(velocity_coeff_)
{
fes.GetEssentialTrueDofs(ess_attr, ess_tdofs);
fes.GetEssentialTrueDofs(nat_attr, nat_tdofs);
T_gf = 0.0;
Q_gf = 0.0;
T_gf_bc = 0.0;
Q_gf_bc = 0.0;
field_collection.AddSourceField("Temperature",&T_gf);
field_collection.AddSourceField("Flux",&Q_gf);
field_collection.AddField("Temperature_BC",&T_gf_bc);
field_collection.AddField("Flux_BC",&Q_gf_bc);
Mform.AddDomainIntegrator(new MassIntegrator);
Mform_e.AddDomainIntegrator(new MassIntegrator);
Kform.AddDomainIntegrator(new DiffusionIntegrator(diffusivity));
Kform_e.AddDomainIntegrator(new DiffusionIntegrator(diffusivity));
if(velocity_coeff)
{
Kform.AddDomainIntegrator(new ConvectionIntegrator(*velocity_coeff, alpha.constant));
Kform_e.AddDomainIntegrator(new ConvectionIntegrator(*velocity_coeff, alpha.constant));
}
if(nat_attr.Max() > 0)
{
bform.AddBoundaryIntegrator(new BoundaryLFIntegrator(Q_fgc),nat_attr);
}
z.SetSize(fes.GetTrueVSize());
q.SetSize(fes.GetTrueVSize());
Assemble();
BuildSolvers();
}
/// Assemble linear and bilinear forms; called if mesh is updated
void Assemble()
{
AssembleLinearForms();
AssembleBilinearForms();
}
void AssembleBilinearForms()
{
Mform.Assemble();
Kform.Assemble();
Mform_e.Assemble();
Kform_e.Assemble();
Array<int> empty;
Mform.FormSystemMatrix(ess_tdofs, Mmat);
Kform.FormSystemMatrix(ess_tdofs, Kmat);
Mform_e.FormSystemMatrix(empty, Mmat_e);
Kform_e.FormSystemMatrix(empty, Kmat_e);
}
void AssembleLinearForms()
{
b.SetSize(fes.GetTrueVSize());
b = 0.0;
bform.Assemble();
bform.ParallelAssemble(b);
}
/// Update finite element space and re-assemble forms
/// if the mesh has changed
void Update() override
{
fes.Update();
T_gf.Update();
T_gf_bc.Update();
Q_gf.Update();
Q_gf_bc.Update();
Mform.Update();
Kform.Update();
Mform_e.Update();
Kform_e.Update();
bform.Update();
Assemble();
updated = true;
}
void BuildSolvers()
{
M_solver.iterative_mode = false;
M_solver.SetRelTol(1e-8);
M_solver.SetAbsTol(0.0);
M_solver.SetMaxIter(1000);
M_solver.SetPrintLevel(0);
M_prec.SetType(HypreSmoother::Jacobi);
M_solver.SetPreconditioner(M_prec);
M_solver.SetOperator(Mmat);
implicit_solver.iterative_mode = false;
implicit_solver.SetRelTol(1e-8);
implicit_solver.SetAbsTol(0.0);
implicit_solver.SetMaxIter(500);
implicit_solver.SetPrintLevel(0);
T_prec.SetType(HypreSmoother::Jacobi);
implicit_solver.SetPreconditioner(T_prec);
}
/// For explict time integration: k = M^{-1} ( -K u + b )
/// Used for partitioned coupling of explicit methods
void Mult(const Vector &u, Vector &k) const override
{
Kmat.Mult(u, z);
z.Neg();
z.Add(1.0, b);
M_solver.Mult(z, k);
k.SetSubVector(ess_tdofs, 0.0);
}
/// For implicit time integration: k solves (M + dt K) k = -K u + b
/// Used for partitioned coupling of implicit methods
void ImplicitSolve(const real_t dt, const Vector &u, Vector &k)
{
AssembleLinearForms();
if((current_dt != dt) || updated)
{
if (T) delete T;
T = Add(1.0, Mmat, dt, Kmat);
implicit_solver.SetOperator(*T);
updated = false;
current_dt = dt;
}
Kmat_e.Mult(u, z);
z.Neg();
z.Add(1.0, b);
implicit_solver.Mult(z, k);
if(IsCoupled() && nat_attr.Max() == 0)
{ // Apply interface conditions on temperature
// Condition imposed by prescribing k = dT/dt
// T_gf_bc contains dT/dt from other domain
T_gf_bc.GetTrueDofs(z);
for (int i = 0; i < ess_tdofs.Size(); i++)
{
int idx = ess_tdofs[i];
k(idx) = z(idx);
}
}
else
{ // if uncoupled, apply standard essential BCs
k.SetSubVector(ess_tdofs, 0.0);
}
}
/// Computes the residual v = M k + K u - b + (Tf - Ts)_int + (Qf + Qs)_int
/// where ()_int represents the temperature and flux interface conditions
/// Used for monolithic coupling
void ImplicitMult(const Vector &u, const Vector &k, Vector &v ) const override
{
// v = M*k + K*u
Mmat_e.Mult(k, v);
Kmat_e.AddMult(u, v);
bform.Assemble();
bform.ParallelAssemble(b);
v.Add(-1.0,b); // v -= b
v.SetSubVector(ess_tdofs, 0.0); // Residual at uncoupled dofs is zero
if(IsCoupled()) // Residual from interface conditions
{ // nat_dofs represent the coupled dofs
T_gf_bc.GetTrueDofs(z); // T from other domain
for (int i = 0; i < nat_tdofs.Size(); i++)
{
int idx = nat_tdofs[i];
v(idx) += (u(idx) - z(idx)); // T_f = T_s
}
Q_gf.GetTrueDofs(q); // Flux from T in this domain
Q_gf_bc.GetTrueDofs(z); // Flux from the other domain
for (int i = 0; i < nat_tdofs.Size(); i++)
{
int idx = nat_tdofs[i];
v(idx) += (q(idx) + z(idx)); // Q_f = -Q_s
}
}
}
/// Transfer temperature and flux boundary data
/// Called if application is monolithically coupled
void Transfer(const Vector &x) override
{
T_gf.SetFromTrueDofs(x);
Q_gf.ProjectBdrCoefficient(Q_coeff,ess_attr);
Q_gf.ProjectBdrCoefficient(Q_coeff,nat_attr);
Application::Transfer(); // Transfer all source fields
}
/// Transfer temperature and flux boundary data
/// Called if application is partitioned coupled
void Transfer(const Vector &x, const Vector &k, real_t dt = 0.0) override
{
if(nat_attr.Max() == 0)
{ // This domain sends flux
add(1.0,x,dt,k,z);
T_gf.SetFromTrueDofs(z);
Q_gf.ProjectBdrCoefficient(Q_coeff,ess_attr);
field_collection.Transfer("Flux"); // Only transfer flux
}
else
{ // This domain is sending temperature
T_gf.SetFromTrueDofs(k);
field_collection.Transfer("Temperature"); // Only transfer temperature
}
}
~ConvectionDiffusion() override
{
if(T) delete T;
}
};
int main(int argc, char *argv[])
{
Mpi::Init();
Hypre::Init();
OptionsParser args(argc, argv);
args.AddOption(&ctx.order, "-o", "--order",
"Finite element order (polynomial degree).");
args.AddOption(&ctx.t_final, "-tf", "--t-final",
"Final time; start time is 0.");
args.AddOption(&ctx.dt, "-dt", "--time-step",
"Time step.");
args.AddOption(&ctx.visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.AddOption(&ctx.vis_steps, "-vs", "--visualization-steps",
"Visualize every n-th timestep.");
args.AddOption(&ctx.ode_solver, "-ode", "--ode-solver-type",
"ODESolver id.");
args.AddOption(&ctx.ser_ref, "-rs", "--serial-refine",
"Number of times to refine the mesh in serial.");
args.AddOption(&ctx.ht_only, "-ht", "--heat-transfer-only",
"-cht", "--conjugate-heat-transfer",
"Conjugate heat transfer or heat transfer only.");
args.AddOption(&ctx.couple_scheme, "-scheme", "--coupling-scheme",
"Coupling scheme: -1 = Monolithic; 0 = Alt. Schw.; >0 = Add. Schw.");
args.AddOption(&ctx.checkres, "-cr", "--checkresult", "-no-cr", "--no-checkresult",
"Enable or disable checking of the result. Returns -1 on failure.");
args.ParseCheck();
int order = ctx.order;
int ode_solver = ctx.ode_solver;
Mesh *serial_mesh = new Mesh("backward-facing-step.msh");
int dim = serial_mesh->Dimension();
for (int i = 0; i < ctx.ser_ref; ++i) { serial_mesh->UniformRefinement(); }
serial_mesh->SetCurvature(order, false, dim, Ordering::byNODES);
ParMesh parent_mesh = ParMesh(MPI_COMM_WORLD, *serial_mesh);
delete serial_mesh;
parent_mesh.UniformRefinement();
// Create the sub-domains and accompanying Finite Element spaces
Array<int> domain_attributes(1);
domain_attributes[0] = 1;
auto fluid_mesh = ParSubMesh::CreateFromDomain(parent_mesh, domain_attributes);
fluid_mesh.SetAttributes();
fluid_mesh.EnsureNodes();
fluid_mesh.Finalize();
domain_attributes[0] = 2;
auto solid_mesh = ParSubMesh::CreateFromDomain(parent_mesh,domain_attributes);
solid_mesh.SetAttributes();
solid_mesh.EnsureNodes();
solid_mesh.Finalize();
// Set essential and natural boundary conditions attributes
Array<int> u_ess_attr, p_ess_attr, noslip_ess_attr;
Array<int> Ts_ess_attr, Ts_nat_attr;
Array<int> Tf_ess_attr, Tf_nat_attr;
if (solid_mesh.bdr_attributes.Size() > 0)
{
Ts_ess_attr.SetSize(solid_mesh.bdr_attributes.Max());
Ts_ess_attr = 0;
Ts_ess_attr[2] = 1; // fluid-solid interface (T_fluid -> T_solid)
Ts_ess_attr[5] = 1; // bottom wall
Ts_nat_attr.SetSize(solid_mesh.bdr_attributes.Max());
Ts_nat_attr = 0;
if(ctx.couple_scheme < 0)
{ // If fully coupled, set nat. bc on interface for equality condition
Ts_nat_attr[2] = 1;
}
}
if (fluid_mesh.bdr_attributes.Size() > 0)
{
u_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
u_ess_attr = 0;
u_ess_attr[0] = 1; // inlet
p_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
p_ess_attr = 0;
p_ess_attr[1] = 1; // outlet
noslip_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
noslip_ess_attr = 1;
noslip_ess_attr[0] = 0; // inlet
noslip_ess_attr[1] = 0; // outlet
Tf_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
Tf_ess_attr = 0;
Tf_ess_attr[0] = 1; // inlet
Tf_nat_attr.SetSize(fluid_mesh.bdr_attributes.Max());
Tf_nat_attr = 0;
Tf_nat_attr[2] = 1; // fluid-solid interface (Qs -> Qf)
}
// Finite element spaces for solid and fluid domains
H1_FECollection ufec(order, dim); // Velocity field (fluid domain)
H1_FECollection pfec(order-1, dim); // Pressure field (fluid domain)
H1_FECollection Tfec(order, dim);
ParFiniteElementSpace u_fes(&fluid_mesh, &ufec, dim, Ordering::byNODES);
ParFiniteElementSpace p_fes(&fluid_mesh, &pfec);
ParFiniteElementSpace Tf_fes(&fluid_mesh, &Tfec);
ParFiniteElementSpace Tuf_fes(&fluid_mesh, &Tfec, dim, Ordering::byNODES);
ParFiniteElementSpace Ts_fes(&solid_mesh, &Tfec);
// Set material properties
// From Backward Facing Step (BFS) Benchmark
real_t Re = ctx.Re;
real_t fluid_density = ctx.density;
real_t viscosity = fluid_density/Re;
real_t Pr = ctx.Pr ;
real_t fluid_alpha = 1.0;
real_t fluid_diffusivity = 1/(Re*Pr);
real_t fluid_kappa = fluid_diffusivity;
real_t solid_alpha = 0.0e0;
real_t kappa_ratio = ctx.kappa_ratio;
real_t solid_diffusivity = kappa_ratio*fluid_diffusivity;
real_t solid_kappa = fluid_kappa*kappa_ratio;
if(Mpi::Root())
{
std::cout << "Fluid Density: " << fluid_density << std::endl;
std::cout << "Viscosity: " << viscosity << std::endl;
std::cout << "Prandtl Number: " << Pr << std::endl;
std::cout << "Reynolds Number: " << Re << std::endl;
std::cout << "Solid Conductivity: " << solid_kappa << std::endl;
std::cout << "Fluid Diffusivity: " << fluid_diffusivity << std::endl;
}
// Navier miniapp
NavierSolver nse_miniapp(&fluid_mesh, order, viscosity);
int max_bdf_order = 3;
nse_miniapp.EnablePA(true);
nse_miniapp.SetMaxBDFOrder(max_bdf_order);
ParGridFunction &uf_gf = *nse_miniapp.GetCurrentVelocity();
ParGridFunction &p_gf = *nse_miniapp.GetCurrentPressure();
Vector vzero(dim); vzero = 0.0;
ConstantCoefficient one_coeff(1.0);
Coefficient *zero_coeff = new ConstantCoefficient(0.0);
VectorCoefficient *zerovec = new VectorConstantCoefficient(vzero);
VectorCoefficient *u_coeff = new VectorFunctionCoefficient(dim, velocity_profile);
// Set initial conditions in fluid
p_gf.ProjectCoefficient(*zero_coeff);
uf_gf.ProjectCoefficient(*zerovec);
uf_gf.ProjectBdrCoefficient(*u_coeff,u_ess_attr);
nse_miniapp.AddVelDirichletBC(u_coeff, u_ess_attr);
nse_miniapp.AddVelDirichletBC(zerovec, noslip_ess_attr);
nse_miniapp.AddPresDirichletBC(zero_coeff, p_ess_attr);
nse_miniapp.Setup(ctx.dt);
/// Create navier block vector
Array<int> nse_offsets({0,uf_gf.ParFESpace()->GetTrueVSize(), p_gf.ParFESpace()->GetTrueVSize()});
nse_offsets.PartialSum();
BlockVector up(nse_offsets);
uf_gf.GetTrueDofs(up.GetBlock(0));
p_gf.GetTrueDofs(up.GetBlock(1));
// Fluid Heat Transfer
ParGridFunction Tuf_gf(&Tuf_fes); Tuf_gf = 0.0;
VectorGridFunctionCoefficient fluid_velocity(&Tuf_gf);
ConvectionDiffusion fluid_ht(Tf_fes, Tf_ess_attr, Tf_nat_attr, fluid_diffusivity,
fluid_kappa, &fluid_velocity, fluid_alpha);
std::unique_ptr<ODESolver> Tf_odesolver = ODESolver::Select(ode_solver);
Tf_odesolver->Init(fluid_ht);
ParGridFunction &Tf_gf = (ParGridFunction&)(*fluid_ht.Fields()["Temperature"]->GetField());
ParGridFunction &Qf_gf = (ParGridFunction&)(*fluid_ht.Fields()["Flux"]->GetField());
ParGridFunction &Tf_gf_bc = (ParGridFunction&)(*fluid_ht.Fields()["Temperature_BC"]->GetField());
ParGridFunction &Qf_gf_bc = (ParGridFunction&)(*fluid_ht.Fields()["Flux_BC"]->GetField());
// Set initial conditions in fluid
Tf_gf.ProjectCoefficient(*zero_coeff);
Tf_gf.ProjectBdrCoefficient(*zero_coeff,Tf_ess_attr);
// Solid Heat Transfer
ConvectionDiffusion solid_ht(Ts_fes, Ts_ess_attr, Ts_nat_attr, solid_diffusivity,
solid_kappa, nullptr, solid_alpha);
std::unique_ptr<ODESolver> Ts_odesolver = ODESolver::Select(ode_solver);
Ts_odesolver->Init(solid_ht);
ParGridFunction &Ts_gf = (ParGridFunction&)(*solid_ht.Fields()["Temperature"]->GetField());
ParGridFunction &Qs_gf = (ParGridFunction&)(*solid_ht.Fields()["Flux"]->GetField());
ParGridFunction &Ts_gf_bc = (ParGridFunction&)(*solid_ht.Fields()["Temperature_BC"]->GetField());
ParGridFunction &Qs_gf_bc = (ParGridFunction&)(*solid_ht.Fields()["Flux_BC"]->GetField());
// Set initial conditions in solid
FunctionCoefficient temp_coeff(temp_profile);
Ts_gf.ProjectCoefficient(temp_coeff);
// Set up the coupled heat transfer multiapp
CoupledOperator ht_operator(2); // two coupled applications: fluid and solid heat transfer
Application* fl_ht_app = ht_operator.AddOperator(&fluid_ht);
Application* sl_ht_app = ht_operator.AddOperator(&solid_ht);
fl_ht_app->SetCoupled(true);
sl_ht_app->SetCoupled(true);
// ODE solver for the coupled operator
std::unique_ptr<ODESolver> ht_odesolver = ODESolver::Select(ode_solver);
ht_odesolver->Init(ht_operator);
Array<int> ht_offsets({0,Tf_fes.GetTrueVSize(), Ts_fes.GetTrueVSize()});
ht_offsets.PartialSum();
BlockVector Tfsv(ht_offsets);
Tf_gf.GetTrueDofs(Tfsv.GetBlock( fl_ht_app->GetOperatorIndex() ));
Ts_gf.GetTrueDofs(Tfsv.GetBlock( sl_ht_app->GetOperatorIndex() ));
// Set up conjugate heat transfer app
CoupledOperator cht_app(2); // two coupled applications: navier, fluid-solid heat transfer
Application* nse_app = cht_app.AddOperator(&nse_miniapp,up.Size()); // (type-erased) Navier miniapp
Application* ht_app = cht_app.AddOperator(ht_odesolver.get()); // Coupled fluid-solid heat transfer ODESolver;
// Set up field transfers
SubMeshTransfer Tf_Ts_map(&Tf_fes, &Ts_fes); // default map if none provided
SubMeshTransfer Qs_Qf_map(&Ts_fes, &Tf_fes); // default map if none provided
SubMeshTransfer uf_Tuf_map(uf_gf.ParFESpace(), Tuf_gf.ParFESpace()); // default map if none provided
// GSLibTransfer Tf_Ts_map(Tf_fes, Ts_fes);
/// Create link between fields in different apps
LinkedFields uf_to_Tuf_lf(&uf_gf, &Tuf_gf, &uf_Tuf_map); // Navier velocity to fluid-heat convection velocity
LinkedFields Tf_to_Ts_lf(&Tf_gf, &Ts_gf_bc, &Tf_Ts_map); // Fluid temperature to solid temperature
LinkedFields Qs_to_Qf_lf(&Qs_gf, &Qf_gf_bc, &Qs_Qf_map); // Solid heat flux to fluid heat flux
/// Different methods for adding the linked field to their source apps
nse_app->AddLinkedFields("Velocity",&uf_to_Tuf_lf);
fl_ht_app->AddLinkedFields("Temperature", &Tf_to_Ts_lf);
fl_ht_app->Fields().AddTargetField("Flux", &Qs_gf_bc, &Tf_Ts_map);
sl_ht_app->Fields().AddLinkedFields("Flux", &Qs_to_Qf_lf);
sl_ht_app->Fields().AddTargetField("Temperature", &Tf_gf_bc,&Qs_Qf_map);
// Solvers
// Select solver for partitioned coupling
FPISolver fp_solver(MPI_COMM_WORLD); // For partitioned solves
AitkenRelaxation fp_relax;
SetSolverParameters(&fp_solver, 0.0, 5e-4, 500, 1, false);
fp_relax.SetBounds(0.0,1.0e-1);
fp_solver.SetRelaxation(5e-1, nullptr); // Use default relaxation method
// fp_solver.SetRelaxation(1e-1, &fp_relax);
// Select solver for monolithic/full coupling
GMRESSolver gmres_solver(MPI_COMM_WORLD);
SetSolverParameters(&gmres_solver, 1e-7, 1e-7, 500, 0, false);
gmres_solver.SetKDim(300);
NewtonSolver newton_solver(MPI_COMM_WORLD);
SetSolverParameters(&newton_solver, 0.0, 1e-7, 100, 0, false);
newton_solver.SetSolver(gmres_solver);
/// Set coupling scheme and corresponding solvers
/// The Navier miniapp and coupled heat transfer ODESolver (both flow maps) are coupled
/// can be coupled in parallel (Additive Schwarz) or serial (Alternating Schwarz) but not
/// monolithically. The fluid and solid heat transfer Applications in the coupled heat
/// transfer ODESolver can also be coupled monolithically.
if(ctx.couple_scheme == -1)
{
cht_app.SetCouplingScheme(CoupledOperator::Scheme::ALTERNATING_SCHWARZ);
ht_operator.SetCouplingScheme(CoupledOperator::Scheme::MONOLITHIC);
ht_operator.SetSolver(&newton_solver);
}
else if(ctx.couple_scheme == 0)
{
cht_app.SetCouplingScheme(CoupledOperator::Scheme::ALTERNATING_SCHWARZ);
ht_operator.SetCouplingScheme(CoupledOperator::Scheme::ALTERNATING_SCHWARZ);
ht_operator.SetSolver(&fp_solver);
}
else
{
cht_app.SetCouplingScheme(CoupledOperator::Scheme::ADDITIVE_SCHWARZ);
ht_operator.SetCouplingScheme(CoupledOperator::Scheme::ADDITIVE_SCHWARZ);
ht_operator.SetSolver(&fp_solver);
}
ht_operator.Assemble(false);
cht_app.Assemble(false);
auto nse_preprocess = [&nse_offsets, &uf_gf, &p_gf](Vector &x) mutable {
BlockVector up(x.GetData(), nse_offsets);
uf_gf.GetTrueDofs(up.GetBlock(0));
p_gf.GetTrueDofs(up.GetBlock(1));
};
auto nse_postprocess = [&nse_offsets, &uf_gf, &p_gf](Vector &x) mutable {
BlockVector up(x.GetData(), nse_offsets);
uf_gf.SetFromTrueDofs(up.GetBlock(0));
p_gf.SetFromTrueDofs(up.GetBlock(1));
};
auto ht_preprocess = [&ht_offsets, &Tf_gf, &Ts_gf](Vector &x) mutable {
BlockVector Tb(x.GetData(), ht_offsets);
Tf_gf.GetTrueDofs(Tb.GetBlock(0));
Ts_gf.GetTrueDofs(Tb.GetBlock(1));
};
auto ht_postprocess = [&ht_offsets, &Tf_gf, &Ts_gf](Vector &x) mutable {
BlockVector Tb(x.GetData(), ht_offsets);
Tf_gf.SetFromTrueDofs(Tb.GetBlock(0));
Ts_gf.SetFromTrueDofs(Tb.GetBlock(1));
};
/// Set pre/post processing lambdas to corresponding apps
ht_app->SetPreProcessFunction(ht_preprocess);
ht_app->SetPostProcessFunction(ht_postprocess);
/// Not strictly necessary since Navier owns and updates GridFunctions
/// internally but included here for completeness
nse_app->SetPreProcessFunction(nse_preprocess);
nse_app->SetPostProcessFunction(nse_postprocess);
// Set up the initial conditions in block vector for the
// coupled application in the correct order
int fl_id = nse_app->GetOperatorIndex();
int ht_id = ht_app->GetOperatorIndex();
BlockVector xb(cht_app.GetBlockOffsets());
xb.GetBlock(fl_id) = up; // Fluid velocity and pressure
xb.GetBlock(ht_id) = Tfsv;
// Set up visualization
ParaViewDataCollection *fluid_pv = nullptr;
ParaViewDataCollection *solid_pv = nullptr;
if(ctx.visualization)
{
fluid_pv = new ParaViewDataCollection("cht-BFS-fluid", &fluid_mesh);
solid_pv = new ParaViewDataCollection("cht-BFS-solid", &solid_mesh);
fluid_pv->SetLevelsOfDetail(order);
fluid_pv->SetDataFormat(VTKFormat::BINARY);
fluid_pv->SetHighOrderOutput(true);
fluid_pv->RegisterField("pressure",&p_gf);
fluid_pv->RegisterField("velocity",&uf_gf);
fluid_pv->RegisterField("convection",&Tuf_gf);
fluid_pv->RegisterField("Temperature",&Tf_gf);
fluid_pv->RegisterField("Flux",&Qf_gf);
solid_pv->SetLevelsOfDetail(order);
solid_pv->SetDataFormat(VTKFormat::BINARY);
solid_pv->SetHighOrderOutput(true);
solid_pv->RegisterField("Temperature",&Ts_gf);
solid_pv->RegisterField("Flux",&Qs_gf);
}
auto save_callback = [&](int cycle, double t)
{
if(fluid_pv)
{
fluid_pv->SetCycle(cycle);
fluid_pv->SetTime(t);
fluid_pv->Save();
}
if(solid_pv)
{
solid_pv->SetCycle(cycle);
solid_pv->SetTime(t);
solid_pv->Save();
}
};
if (Mpi::Root()) {
out << "Starting time integration..." << std::endl;
}
StopWatch timer;
timer.Start();
real_t t = 0.0;
bool last_step = false;
int tindex = 1;
save_callback(0, t);
last_step = false;
for (; !last_step; tindex++)
{
if (t + ctx.dt >= ctx.t_final - ctx.dt/2){ last_step = true; }
if(ctx.ht_only)
{
ht_odesolver->Step(Tfsv,t,ctx.dt);
Tf_gf.SetFromTrueDofs(Tfsv.GetBlock(0));
Ts_gf.SetFromTrueDofs(Tfsv.GetBlock(1));
}
else
{
cht_app.Step(xb, t, ctx.dt);
}
if (last_step || (tindex % ctx.vis_steps) == 0){
if (Mpi::Root()) { out << "step " << tindex << ", t = " << t << std::endl;}
save_callback(tindex, t);
}
}
timer.Stop();
if (Mpi::Root()){
out << "Total time: " << timer.RealTime() << " seconds." << std::endl;
}
/// Compute interface error
if(ctx.checkres)
{
Array<int> fl_int_attr(fluid_mesh.bdr_attributes.Max());
fl_int_attr[2] = 1; // fluid-solid interface
/// Create submesh and FE space for the interface
ParSubMesh int_mesh = ParSubMesh::CreateFromBoundary(fluid_mesh, fl_int_attr);
ParFiniteElementSpace int_fes(&int_mesh, &Tfec);
/// Receiving grid functions on the interface
ParGridFunction fl_int_gf(&int_fes);
ParGridFunction sl_int_gf(&int_fes);
/// Maps and linked fields to transfer domain grid functions to the interface
SubMeshTransfer fl_to_int_submesh(&Tf_fes, &int_fes);
LinkedFields TQf_to_Tint_lf(&Tf_gf, &fl_int_gf, &fl_to_int_submesh);
LinkedFields TQs_to_Tint_lf(&Tf_gf_bc, &sl_int_gf, &fl_to_int_submesh);
cht_app.Transfer(xb); // Transfer all fields *_gf to their target fields *_gf_bc
fl_int_gf = 0.0; sl_int_gf = 0.0;
TQf_to_Tint_lf.Transfer(); // Transfer Tf to interface
TQs_to_Tint_lf.Transfer(); // Transfer Ts (in Tf_gf_bc) to interface
real_t err_T = sqrt(DistanceSquared(int_mesh.GetComm(), fl_int_gf, sl_int_gf));
// Update sources in existing linked fields (can also create new ones)
TQf_to_Tint_lf.SetSource(&Qf_gf);
TQs_to_Tint_lf.SetSource(&Qf_gf_bc);
fl_int_gf = 0.0; sl_int_gf = 0.0;
TQf_to_Tint_lf.Transfer(); // Transfer Qf to interface
TQs_to_Tint_lf.Transfer(); // Transfer Qs (in Qf_gf_bc) to interface
real_t err_Q = sqrt(DistanceSquared(int_mesh.GetComm(), fl_int_gf, sl_int_gf));
if (sqrt(err_T) > ctx.tol_T || sqrt(err_Q) > ctx.tol_Q)
{
if (Mpi::Root())
{
mfem::out << "Result has a larger error than expected."
<< "T Error = " << sqrt(err_T)
<< ", Q Error = " << sqrt(err_Q)
<< std::endl;
}
return -1;
}
}
if(fluid_pv) delete fluid_pv;
if(solid_pv) delete solid_pv;
delete zero_coeff;
delete zerovec;
delete u_coeff;
return 0;
}
void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol , int max_it,
int print_level, bool iterative_mode)
{
solver->SetRelTol(rtol);
solver->SetAbsTol(atol);
solver->SetMaxIter(max_it);
solver->SetPrintLevel(print_level);
solver->iterative_mode = iterative_mode;
}
+616
View File
@@ -0,0 +1,616 @@
/**
* Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
* at the Lawrence Livermore National Laboratory. All Rights reserved. See files
* LICENSE and NOTICE for details. LLNL-CODE-806117.
*
* This file is part of the MFEM library. For more information and source code
* availability visit https://mfem.org.
*
* MFEM is free software; you can redistribute it and/or modify it under the
* terms of the BSD-3 license. We welcome feedback and contributions, see file
* CONTRIBUTING.md for details.
*
* --------------------------------------------------
* Fluid-Structure Interaction miniapp
* --------------------------------------------------
*
* This miniapp simulates fluid-structure interaction (FSI) problems described in
* the paper:
*
* with the incompressible Navier-Stokes equations, in Arbitrary Lagrangian-Eulerian (ALE)
* formulation, coupled with linear elasticity equations in the solid domain.
* The coupling is done with a partitioned approach using the alternating or additive
* Schwarz method. The fluid mesh motion is handled with a mesh displacement diffusion
* approach. The geometry is a channel with a cylinder and a flexible beam attached to
* the cylinder downstream.
*
* The following boundary conditions are applied:
* 1) Channel inlet (attribute 1): parabolic velocity profile
* 2) Channel outlet (attribute 2): zero-pressure
* 3) Channel walls: no-slip
*
* with the following interface conditions at the fluid-structure interface:
* 1) Continuity of velocity: u_f = u_s
* 2) Continuity of traction: -*pI + mu(grad(u_f)+grad(u_f)^T))n = sigma_sn
*
* The velocity continuity is imposed using the stage-slope in the implicit multistage method
* ku_f = du_f/dt = du_s/dt = ku_s, where ku_f and ku_s are the fluid and solid stage-slopes
*
* The mesh morphing is modeled as a displacement diffusion equation,
* dx/dt = κΔx, with the ALE mesh velocity, w = dx/dt.
*
* Sample run:
* mpirun -np 6 ./fsi -vs 5 -dt 5e-3 -tf 10 -o 3 -rs 2 -ode 21 -U 1.0 -cs 1 -idir fsi-turek
*/
#include "mfem.hpp"
#include "multiapp.hpp"
#include "apps/navier_stokes.hpp"
#include "apps/elasticity.hpp"
#include "apps/mesh_morpher.hpp"
#include <filesystem>
using namespace mfem;
using namespace std;
// mpirun -np 6 ./fsi -vs 5 -dt 1e-2 -tf 5 -o 2 -rs 2 -ode 21
// mpirun -np 6 ./fsi -vs 5 -dt 5e-3 -tf 10 -o 3 -rs 2 -ode 21 -U 1.0 -cs 0 -init
void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol , int max_it,
int print_level, bool iterative_mode);
bool ReadGridFunctionFromFile(const string &dirname, const string &gf_name,
ParMesh &mesh, ParGridFunction &gf);
void CollectArrays(std::vector<Array<int>*> &dof_arrays, Array<int> &tdof_array);
int main(int argc, char *argv[])
{
Mpi::Init();
Hypre::Init();
int num_procs = Mpi::WorldSize();
int myid = Mpi::WorldRank();
int order = 2;
int ser_ref = 0;
int ode_solver_type = 21; // (21) BackwardEulerSolver
// (22) SDIRK23Solver
// (23) SDIRK33Solver
// (34) SDIRK34Solver
real_t Uavg = 2.0;
real_t t_dev = 2.0;
real_t t_final = 1.0;
real_t dt = 1.0e-3;
real_t relax_factor = 1.0;
int couple_scheme = 1;
bool init = false;
bool lsave = false;
bool visualization = true;
int vis_steps = 10;
std::string init_dir= "";
OptionsParser args(argc, argv);
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree).");
args.AddOption(&t_final, "-tf", "--t-final",
"Final time; start time is 0.");
args.AddOption(&dt, "-dt", "--time-step",
"Time step.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.AddOption(&vis_steps, "-vs", "--visualization-steps",
"Visualize every n-th timestep.");
args.AddOption(&ode_solver_type, "-ode", "--ode-solver-type",
"ODESolver id.");
args.AddOption(&ser_ref, "-rs", "--serial-refine",
"Number of times to refine the mesh in serial.");
args.AddOption(&Uavg, "-U", "--velocity", "Mean velocity.");
args.AddOption(&t_dev, "-t_dev", "--transistion", "Developed flow transition time.");
args.AddOption(&couple_scheme, "-cs", "--coupling-scheme",
"Coupling scheme: -1 = Monolithic; 0 = Add. Schw.; >0 = Alt. Schw.");
args.AddOption(&relax_factor, "-rf", "--relaxation-factor",
"Initial relaxation factor for the FPI solver.");
args.AddOption(&init_dir, "-idir", "--init-directory",
"Directory containing intialization files. If doesn't exist, used to write init files.");
args.AddOption(&lsave, "-save", "--save-init", "-no-save",
"--no-save", "Enable or disable saving initialization files.");
args.ParseCheck();
init = !init_dir.empty(); // if init_dir is provided, then init = true
std::string mesh_file = "channel-cylinder.msh";
Mesh *serial_mesh = new Mesh(mesh_file);
int dim = serial_mesh->Dimension();
for (int i = 0; i < ser_ref; ++i) { serial_mesh->UniformRefinement(); }
serial_mesh->SetCurvature(order, false, dim, Ordering::byNODES);
serial_mesh->EnsureNCMesh();
ParMesh parent_mesh = ParMesh(MPI_COMM_WORLD, *serial_mesh);
delete serial_mesh;
// Create mesh (sub)domains
Array<int> domain_attributes(1);
// Create submesh for solid
domain_attributes[0] = 2;
auto solid_mesh = ParSubMesh::CreateFromDomain(parent_mesh, domain_attributes);
solid_mesh.SetAttributes();
solid_mesh.EnsureNodes();
// Create submesh for fluid
domain_attributes[0] = 1;
auto fluid_mesh = ParSubMesh::CreateFromDomain(parent_mesh, domain_attributes);
fluid_mesh.SetAttributes();
fluid_mesh.EnsureNodes();
// Set essential and natural boundary conditions attributes
Array<int> u_ess_attr, p_ess_attr, noslip_attr;
Array<int> xs_ess_attr, xs_nat_attr;
Array<int> xf_ess_attr;
Array<int> empty;
if (solid_mesh.bdr_attributes.Size() > 0)
{
xs_ess_attr.SetSize(solid_mesh.bdr_attributes.Max());
xs_nat_attr.SetSize(solid_mesh.bdr_attributes.Max());
xs_ess_attr = 0; xs_nat_attr = 0;
xs_nat_attr[4] = 1; // beam wall
xs_ess_attr[6] = 1; // beam cylinder curve
}
if (fluid_mesh.bdr_attributes.Size() > 0)
{
u_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
p_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
noslip_attr.SetSize(fluid_mesh.bdr_attributes.Max());
xf_ess_attr.SetSize(fluid_mesh.bdr_attributes.Max());
xf_ess_attr = 1;
u_ess_attr = 1;
u_ess_attr[1] = 0; // outlet
p_ess_attr = 0;
p_ess_attr[1] = 1; // outlet
noslip_attr = 1;
noslip_attr[0] = 0; // inlet
noslip_attr[1] = 0; // outlet
}
// Finite element spaces for solid and fluid domains
H1_FECollection xfec(order, dim); // Displacement field (solid and fluid domains)
H1_FECollection ufec(order, dim); // Velocity field (fluid domain)
H1_FECollection pfec(order-1, dim); // Pressure field (fluid domain)
ParFiniteElementSpace xs_fes(&solid_mesh, &xfec, dim, Ordering::byNODES);
ParFiniteElementSpace xf_fes(&fluid_mesh, &xfec, dim, Ordering::byNODES);
ParFiniteElementSpace u_fes(&fluid_mesh, &ufec, dim, Ordering::byNODES);
ParFiniteElementSpace p_fes(&fluid_mesh, &pfec);
Vector vzero(dim); vzero = 0.0;
VectorConstantCoefficient zerovec(vzero);
// Inlet velocity boundary condition
auto velocity_profile = [&Uavg, &t_dev](const Vector &x, double t, Vector &u) mutable
{
double xi = x(0), yi = x(1);
double U = Uavg;
double ramp_time = 2.0;
u = 0.0;
if (xi == 0.0)
{
u(0) = 1.5 * 4.0 * U * yi * (0.41 - yi) / (pow(0.41, 2.0));
u(1) = 0.0;
}
if(t < ramp_time) u(0) *= 0.5*(1.0 - cos(0.5*M_PI*t));
};
// Set material properties
// real_t E = 5.6e6; // Young's modulus
real_t E = 1.4e6; // Young's modulus
real_t nu = 4.0e-1; // Poisson's ratio
real_t solid_density = 1.0e4;
real_t viscosity = 1.0e-3;
real_t fluid_density = 1.0e3;
real_t mesh_diffusion = 1.0e0;
real_t lame_mu = E / (2.0 * (1.0 + nu));
real_t lame_lambda = 2.0*lame_mu*nu/(1.0-2.0*nu);
real_t compressibility = 0e-4; // Artificial Compressibility
bool scaled_pressure = true;
if(Mpi::Root())
{
std::cout << "Young's Modulus: " << E << std::endl;
std::cout << "Poisson's Ratio: " << nu << std::endl;
std::cout << "Lame's Mu: " << lame_mu << std::endl;
std::cout << "Lame's Lambda: " << lame_lambda << std::endl;
std::cout << "Solid Density: " << solid_density << std::endl;
std::cout << "Fluid Density: " << fluid_density << std::endl;
std::cout << "Viscosity: " << viscosity << std::endl;
std::cout << "Mesh Diffusion: " << mesh_diffusion << std::endl;
std::cout << "Artificial Compressibility: " << compressibility << std::endl;
std::cout << "Mean Velocity: " << Uavg << std::endl;
std::cout << "Reynolds Number: " << (Uavg*0.1)/viscosity << std::endl;
std::cout << "AE Number: " << (E)/(Uavg*Uavg*fluid_density) << std::endl;
std::cout << "Beta Number: " << solid_density/fluid_density << std::endl;
}
// Build individual applications (morpher, elasticity, navier-stokes)
// Mesh Morpher
MeshDiffusion morpher(xf_fes, xf_ess_attr, mesh_diffusion);
// Morphing grid functions
ParGridFunction &xf_gf = *morpher.Fields().GetField("Displacement");
ParGridFunction &dxf_gf = *morpher.Fields().GetField("dxdt");
ParGridFunction &xf_gf_bc = *morpher.Fields().GetField("Displacement_BC");
ParGridFunction &dxf_gf_bc = *morpher.Fields().GetField("Velocity_BC");
/// Morphing solution
Vector xf(xf_fes.GetTrueVSize());
xf_gf.GetTrueDofs(xf);
ParGridFunction mesh_disp(&xf_fes), mesh_vel(&xf_fes);
mesh_disp = 0.0; mesh_vel = 0.0;
// Elasticity
Elasticity elasticity(xs_fes, xs_ess_attr, xs_nat_attr, solid_density, lame_mu, lame_lambda);
// Elasticity grid functions
ParGridFunction &xs_gf = *elasticity.Fields().GetField("Displacement");
ParGridFunction &us_gf = *elasticity.Fields().GetField("Velocity");
ParGridFunction &stress_gf = *elasticity.Fields().GetField("Traction");
/// Solid solution
Array<int> solid_offsets({0,xs_fes.GetTrueVSize(), xs_fes.GetTrueVSize()});
solid_offsets.PartialSum();
BlockVector xu(solid_offsets); // Solid displacement and velocity
us_gf.GetTrueDofs(xu.GetBlock(0));
xs_gf.GetTrueDofs(xu.GetBlock(1));
// Navier-Stokes
VectorGridFunctionCoefficient ale_uf(&dxf_gf);
ScalarVectorProductCoefficient neg_ale(-1.0, ale_uf);
VectorCoefficient *ale_velocity = &neg_ale;
NavierStokes nse(u_fes, p_fes, u_ess_attr, p_ess_attr,
fluid_density, viscosity, compressibility,
scaled_pressure, ale_velocity);
std::unique_ptr<ODESolver> nse_solver = ODESolver::Select(ode_solver_type);
nse_solver->Init(nse);
// Navier-Stokes grid functions
ParGridFunction &p_gf = *nse.Fields().GetField("Pressure");
ParGridFunction &uf_gf = *nse.Fields().GetField("Velocity");
ParGridFunction &tau_gf = *nse.Fields().GetField("Stress");
ParGridFunction &uf_gf_bc = *nse.Fields().GetField("Velocity_BC");
ConstantCoefficient p_coeff(0.0);
VectorFunctionCoefficient u_coeff(dim, velocity_profile);
// Set initial conditions in fluid
p_gf.ProjectCoefficient(p_coeff);
uf_gf.ProjectCoefficient(u_coeff);
uf_gf.ProjectBdrCoefficient(zerovec,noslip_attr);
if(init) // Initialize from file
{
ReadGridFunctionFromFile(init_dir, "p-init.gf", fluid_mesh, p_gf);
ReadGridFunctionFromFile(init_dir, "u-init.gf", fluid_mesh, uf_gf);
u_coeff.SetTime(t_dev+2.0); // time beyond ramp-up
uf_gf.ProjectBdrCoefficient(u_coeff,u_ess_attr);
}
tau_gf.ProjectBdrCoefficient(nse.stress_coeff,u_ess_attr);
/// Fluid solution
Array<int> nse_offsets({0,u_fes.GetTrueVSize(), p_fes.GetTrueVSize()});
nse_offsets.PartialSum();
BlockVector up(nse_offsets);
uf_gf.GetTrueDofs(up.GetBlock(0));
p_gf.GetTrueDofs(up.GetBlock(1));
// Three coupled applications: navier, elasticity, and morpher
CoupledOperator multiapp(3);
std::unique_ptr<ODESolver> coupled_solver = ODESolver::Select(ode_solver_type);
Application* nse_app = multiapp.AddOperator(&nse);
Application* elasticity_app = multiapp.AddOperator(&elasticity);
Application* morpher_app = multiapp.AddOperator(&morpher);
// Set up field transfer
// NativeTransfer strsf_to_strss_map(u_fes, xs_fes); // default map if none provided
// GSLibTransfer strsf_to_strss_map(u_fes, xs_fes, xs_nat_attr);
nse_app->Fields().AddTargetField("Stress", &stress_gf);
elasticity_app->Fields().AddTargetField("Velocity_BC", &uf_gf_bc);
elasticity_app->Fields().AddTargetField("Displacement_BC", &xf_gf_bc);
elasticity_app->Fields().AddTargetField("Velocity_BC", &dxf_gf_bc);
// Assemble the true (offseted) dofs array for the coupled FSI system
Array<int> tdofs;
std::vector<Array<int>*> state_dofs = {&nse.u_ess_tdofs, &nse.p_ess_tdofs,
&elasticity.ess_tdofs, &elasticity.ess_tdofs,
&morpher.ess_tdofs};
CollectArrays(state_dofs, tdofs);
// Set up coupling parameters (schemes and solvers)
FPISolver fp_solver(MPI_COMM_WORLD); // For partitioned solves
AitkenRelaxation fp_relax;
ConstrainedInnerProduct constr_ipo(MPI_COMM_WORLD , tdofs);
SetSolverParameters(&fp_solver, 0.0, 5e-4, 100, 1, false);
fp_solver.SetRelaxation(relax_factor, &fp_relax);
// fp_solver.SetRelaxation(relax_factor, nullptr); // Use default relaxation method
// fp_relax.SetBounds(-1.0e0,1.0e0);
fp_relax.SetAbsoluteLowerBound(1e-1);
fp_solver.SetInnerProduct(&constr_ipo);
fp_relax.SetInnerProduct(&constr_ipo);
NewtonSolver newton_solver(MPI_COMM_WORLD); // For fully coupled
GMRESSolver gmres_solver(MPI_COMM_WORLD); // For fully coupled
if(couple_scheme == -1)
{
MFEM_ABORT("Monolithic coupling not supported for FSI.")
multiapp.SetCouplingScheme(CoupledOperator::Scheme::MONOLITHIC);
SetSolverParameters(&gmres_solver, 1e-3, 1e-3, 500, 1, false);
gmres_solver.SetKDim(300);
SetSolverParameters(&newton_solver, 0.0, 1e-4, 30, 1, false);
newton_solver.SetSolver(gmres_solver);
multiapp.SetSolver(&newton_solver); // Set the solver for the multiapp
}
else if(couple_scheme == 0)
{
multiapp.SetCouplingScheme(CoupledOperator::Scheme::ADDITIVE_SCHWARZ);
multiapp.SetSolver(&fp_solver); // Set the solver for the multiapp
}
else
{
multiapp.SetCouplingScheme(CoupledOperator::Scheme::ALTERNATING_SCHWARZ);
multiapp.SetSolver(&fp_solver); // Set the solver for the multiapp
}
multiapp.Assemble(false); // Assemble the multiapp (build OperatorCoupler)
multiapp.Finalize(false); // Finalize the multiapp (perform checks)
coupled_solver->Init(multiapp);
// Set up the initial conditions in block vector for the
// coupled application in the correct order
int fl_id = nse_app->GetOperatorIndex();
int el_id = elasticity_app->GetOperatorIndex();
int morph_id = morpher_app->GetOperatorIndex();
BlockVector xb(multiapp.GetBlockOffsets());
xb.GetBlock(el_id) = xu; // Solid displacement and velocity
xb.GetBlock(fl_id) = up; // Fluid velocity and pressure
xb.GetBlock(morph_id) = xf; // Fluid displacement
auto update_grid_functions = [&](Vector &x) mutable
{
BlockVector xb(x.GetData(), multiapp.GetBlockOffsets());
BlockVector elas_x(xb.GetBlock(el_id).GetData(), solid_offsets);
BlockVector nse_x(xb.GetBlock(fl_id).GetData(), nse_offsets);
Vector morph_x(xb.GetBlock(morph_id).GetData(), xb.BlockSize(morph_id));
xf_gf.SetFromTrueDofs(morph_x);
us_gf.SetFromTrueDofs(elas_x.GetBlock(0));
xs_gf.SetFromTrueDofs(elas_x.GetBlock(1));
uf_gf.SetFromTrueDofs(nse_x.GetBlock(0));
p_gf.SetFromTrueDofs(nse_x.GetBlock(1));
tau_gf.ProjectBdrCoefficient(nse.stress_coeff,u_ess_attr);
};
auto update_nse_grid_functions = [&](Vector &x) mutable
{
BlockVector nse_x(x.GetData(), nse_offsets);
uf_gf.SetFromTrueDofs(nse_x.GetBlock(0));
p_gf.SetFromTrueDofs(nse_x.GetBlock(1));
tau_gf.ProjectBdrCoefficient(nse.stress_coeff,u_ess_attr);
};
// Set up visualization
ParaViewDataCollection *fluid_pv = nullptr;
ParaViewDataCollection *solid_pv = nullptr;
if(visualization)
{
fluid_pv = new ParaViewDataCollection("fsi-fluid", &fluid_mesh);
solid_pv = new ParaViewDataCollection("fsi-solid", &solid_mesh);
fluid_pv->SetLevelsOfDetail(order);
fluid_pv->SetDataFormat(VTKFormat::BINARY);
fluid_pv->SetHighOrderOutput(true);
fluid_pv->RegisterField("displacement",&mesh_disp);
fluid_pv->RegisterField("dxdt",&mesh_vel);
fluid_pv->RegisterField("pressure",&p_gf);
fluid_pv->RegisterField("velocity",&uf_gf);
fluid_pv->RegisterField("stress",&tau_gf);
fluid_pv->RegisterField("ale_velocity",&dxf_gf);
solid_pv->SetLevelsOfDetail(order);
solid_pv->SetDataFormat(VTKFormat::BINARY);
solid_pv->SetHighOrderOutput(true);
solid_pv->RegisterField("displacement",&xs_gf);
solid_pv->RegisterField("velocity",&us_gf);
solid_pv->RegisterField("stress",&stress_gf);
}
auto save_callback = [&](int cycle, double t)
{
if(fluid_pv)
{
fluid_pv->SetCycle(cycle);
fluid_pv->SetTime(t);
fluid_pv->Save();
}
if(solid_pv)
{
solid_pv->SetCycle(cycle);
solid_pv->SetTime(t);
solid_pv->Save();
}
};
if (Mpi::Root()) out << "Starting time integration..." << std::endl;
StopWatch timer;
timer.Start();
real_t t = 0.0;
bool last_step = false;
int tindex = 1;
save_callback(0, t);
// Solve the Navier-Stokes equations to fully developed flow time, t_dev
if(t_dev > 0.0 && !init)
{
nse_app->SetOperationID(Application::OperationID::STEP);
for (; !last_step; tindex++)
{
if (t + dt >= t_dev - dt/2){ last_step = true; }
u_coeff.SetTime(t); // Slowly ramp-up inlet velocity
uf_gf.ProjectBdrCoefficient(u_coeff,u_ess_attr);
nse_solver->Step(up,t,dt);
if (last_step || (tindex % vis_steps) == 0){
if (Mpi::Root()) { out << "step " << tindex << ", t = " << t << std::endl;}
update_nse_grid_functions(up);
save_callback(tindex, t);
}
}
if ((myid==0) && lsave)
{
std::filesystem::path dir_path = init_dir;
if(!std::filesystem::is_directory(dir_path)) std::filesystem::create_directory(dir_path);
}
if(lsave){
p_gf.Save((init_dir+"/p-init.gf").c_str());
uf_gf.Save((init_dir+"/u-init.gf").c_str());
}
xb.GetBlock(fl_id) = up; // Update nse block
}
vis_steps = 1;
last_step = false;
nse_app->SetCoupled(true);
morpher_app->SetCoupled(true);
elasticity_app->SetCoupled(true);
multiapp.Transfer(xb);
// Store original fluid nodes; morphing is done w.r.t. original configuration
GridFunction fluid_nodes_orig = *(fluid_mesh.GetNodes());
for (; !last_step; tindex++)
{
if (t + dt >= t_final - dt/2){ last_step = true; }
coupled_solver->Step(xb, t, dt);
// multiapp.Transfer(xb);
update_grid_functions(xb);
// Morph fluid mesh
mesh_disp.SetFromTrueDofs(xb.GetBlock(morph_id));
mesh_vel = dxf_gf;
GridFunction *fluid_nodes = fluid_mesh.GetNodes();
*fluid_nodes = fluid_nodes_orig;
*fluid_nodes += mesh_disp;
fluid_mesh.DeleteGeometricFactors();
// Update FE spaces, grid functions and forms
// after fluid mesh update
nse_app->Update();
morpher_app->Update();
if (last_step || (tindex % vis_steps) == 0){
if (Mpi::Root()){ out << "step " << tindex << ", t = " << t << std::endl;}
save_callback(tindex, t);
}
}
timer.Stop();
if (Mpi::Root()){
out << "Total time: " << timer.RealTime() << " seconds." << std::endl;
}
return 0;
}
void SetSolverParameters(IterativeSolver *solver, real_t rtol, real_t atol , int max_it,
int print_level, bool iterative_mode)
{
solver->SetRelTol(rtol);
solver->SetAbsTol(atol);
solver->SetMaxIter(max_it);
solver->SetPrintLevel(print_level);
solver->iterative_mode = iterative_mode;
}
bool ReadGridFunctionFromFile(const string &dirname, const string &gf_name,
ParMesh &mesh, ParGridFunction &gf)
{
int myid = Mpi::WorldRank();
std::string mpirank = std::to_string(myid);
std::string filename = dirname+"/"+gf_name+"."+mpirank.insert(0,6-mpirank.length(),'0');
bool sucess = false;
if (std::filesystem::exists(filename))
{
istream *ifile;
ifile = new ifstream(filename);
gf = ParGridFunction(&mesh,*ifile);
delete ifile;
sucess = true;
}
return sucess;
}
void CollectArrays(std::vector<Array<int>*> &dof_arrays, Array<int> &tdof_array)
{
int total_dofs = 0;
for (const auto& arr : dof_arrays)
{
total_dofs += arr->Size();
}
int n = 0;
tdof_array.SetSize(total_dofs);
for (auto arr : dof_arrays)
{
std::transform(arr->begin(), arr->end(), tdof_array.begin() + n,
[&](int i) { return i+n; });
n += arr->Size();
}
}
+548
View File
@@ -0,0 +1,548 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "multiapp.hpp"
#ifdef MFEM_USE_MPI
namespace mfem
{
GridFunctionTransfer* GridFunctionTransfer::Select(Type type,
ParFiniteElementSpace *src,
ParFiniteElementSpace *tar)
{
switch (type)
{
case Type::SUBMESH:
return new SubMeshTransfer(src, tar);
// case Type::GSLIB:
// return new GSLibTransfer(src, tar);
default:
MFEM_ABORT("Unknown GridFunctionTransfer scheme: " << static_cast<int>(type));
}
}
void SubMeshTransfer::Transfer(const Field &src, Field &tar)
{
// Use GetField() to get the underlying Vector and ParGridFunction objects
const ParGridFunction &src_gf = dynamic_cast<const ParGridFunction&>(*src.GetField());
ParGridFunction &tar_gf = dynamic_cast<ParGridFunction&>(*tar.GetField());
transfer_map->Transfer(src_gf, tar_gf);
}
OperatorCoupler* OperatorCoupler::Select(CoupledOperator *op,
Scheme scheme)
{
switch (scheme)
{
case Scheme::MONOLITHIC:
{
real_t fd_eps = 1e-6;
return new JacobianFreeFullCoupler(op, fd_eps);
}
case Scheme::ADDITIVE_SCHWARZ:
return new AdditiveSchwarzCoupler(op);
case Scheme::ALTERNATING_SCHWARZ:
return new AlternatingSchwarzCoupler(op);
case Scheme::NONE:
return nullptr;
default:
MFEM_ABORT("Unknown coupling scheme: " << static_cast<int>(scheme));
}
}
CoupledOperator::~CoupledOperator()
{
if(solver && own_solver) delete solver;
if(op_coupler && own_op_coupler) delete op_coupler;
for(int i=0; i < nops; i++)
{
if(operators_owned[i] && operators[i]) delete operators[i];
}
}
void CoupledOperator::SetOperatorCoupler(OperatorCoupler* op, bool own)
{
if(op_coupler && own_op_coupler) delete op_coupler;
op_coupler = op;
own_op_coupler = own;
coupler_type = op_coupler->GetType();
}
void CoupledOperator::Initialize(bool do_ops)
{
if (do_ops)
{
for (auto &op : operators)
{
op->Initialize();
}
}
}
void CoupledOperator::Assemble(bool do_ops)
{
if (do_ops)
{
for (auto &op : operators)
{
op->Assemble();
}
}
// Check block offsets against operator size
Array<int> true_offsets(Size()+1);
bool offset_consistent = true;
true_offsets = 0;
int max_size = 0;
for (int i=0; i < nops; i++)
{
auto op = GetOperator(i);
int block_size = offsets[i+1]-offsets[i];
true_offsets[i+1] = true_offsets[i] + op->Width();
if (block_size != op->Width())
{
offset_consistent = false;
}
}
if (!offset_consistent)
{
MFEM_WARNING("Block offsets inconsistent with operator sizes."
"Using default offsets.");
offsets = true_offsets;
max_op_size = max_size;
}
if(op_coupler && own_op_coupler) delete op_coupler;
op_coupler = OperatorCoupler::Select(this, coupler_type);
if(solver) solver->SetOperator(*op_coupler);
}
void CoupledOperator::Finalize(bool do_ops)
{
if (do_ops)
{
for (auto &op : operators)
{
op->Finalize();
}
}
}
void CoupledOperator::PreProcess(Vector &x, bool do_ops)
{
if (do_ops)
{
BlockVector xb(x.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
operators[i]->PreProcess(xi);
}
}
}
void CoupledOperator::PostProcess(Vector &x, bool do_ops)
{
if (do_ops)
{
BlockVector xb(x.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
operators[i]->PostProcess(xi);
}
}
}
void CoupledOperator::SetOperationID(OperationID id, bool do_ops)
{
Application::SetOperationID(id);
if (do_ops)
{
for (auto &op : operators)
{
op->SetOperationID(id);
}
}
}
void CoupledOperator::SetTime(const real_t t_)
{
TimeDependentOperator::SetTime(t_);
if(op_coupler) op_coupler->SetTime(t_);
for (auto &op : operators)
{
op->SetTime(t_);
}
}
void CoupledOperator::Transfer(const Vector &x)
{
BlockVector xb(x.GetData(), offsets);
for (int i=0; i < nops; i++)
{
operators[i]->Transfer(xb.GetBlock(i));
}
}
void CoupledOperator::Transfer(const Vector &u, const Vector &k, real_t dt)
{
BlockVector ub(u.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
for (int i=0; i < nops; i++)
{
operators[i]->Transfer(ub.GetBlock(i), kb.GetBlock(i), dt);
}
}
void CoupledOperator::Mult(const Vector &x, Vector &y) const
{
if(op_coupler && coupler_type != Scheme::NONE)
{
if(solver) {
op_coupler->SetOperationID(OperationID::MULT);
op_coupler->SetInput(&x);
solver->Mult(b,y);
}
else {
op_coupler->Mult(x,y);
}
}
else
{
BlockVector xb(x.GetData(), offsets);
BlockVector yb(y.GetData(), offsets);
for (int i=0; i < nops; i++)
{
operators[i]->SetOperationID(OperationID::MULT);
operators[i]->Mult(xb.GetBlock(i), yb.GetBlock(i));
}
}
}
void CoupledOperator::ImplicitSolve(const real_t dt, const Vector &x, Vector &k ){
if(op_coupler && coupler_type != Scheme::NONE)
{
if(solver) {
op_coupler->SetOperationID(OperationID::IMPLICIT_SOLVE); ///< OperatorCoupler::Mult() -> OperatorCoupler::ImplicitSolve()
op_coupler->SetTimeStep(dt);
op_coupler->SetInput(&x);
solver->Mult(b,k);
}
else {
op_coupler->ImplicitSolve(dt,x,k);
}
}
else
{
BlockVector xb(x.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
Vector &ki = kb.GetBlock(i);
operators[i]->SetOperationID(OperationID::IMPLICIT_SOLVE);
operators[i]->ImplicitSolve(dt,xi,ki); ///< Solve the implicit system for the application
}
}
}
void CoupledOperator::Step(Vector &x, real_t &t, real_t &dt)
{
if(op_coupler && coupler_type != Scheme::NONE)
{
if(solver) {
op_coupler->SetOperationID(OperationID::STEP); ///< OperatorCoupler::Mult() -> OperatorCoupler::Mult()
op_coupler->SetTimeStep(dt); ///< Set the time step for the ODE Solver
op_coupler->SetTime(t);
op_coupler->SetInput(&x);
solver->Mult(b,x);
}
else {
op_coupler->Step(x,t,dt);
}
}
else
{
BlockVector xb(x.GetData(), offsets);
for (int i=0; i < nops; i++)
{
real_t t0 = t; ///< Store the current time
real_t dt0 = dt; ///< Store the current time step
Vector &xi = xb.GetBlock(i);
operators[i]->SetOperationID(OperationID::STEP);
operators[i]->Step(xi,t0,dt0); ///< Advance the time step for application
}
t += dt; ///< Update the time after all applications have been stepped forward
///< NOTE: does not work for adaptive time-stepping
}
}
void CoupledOperator::ImplicitMult(const Vector &u, const Vector &k, Vector &v) const
{
BlockVector ub(u.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
BlockVector vb(v.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &ui = ub.GetBlock(i);
Vector &ki = kb.GetBlock(i);
Vector &vi = vb.GetBlock(i);
operators[i]->SetOperationID(OperationID::IMPLICIT_MULT);
operators[i]->ImplicitMult(ui,ki,vi); ///< Solve the implicit system for the application
}
}
void CoupledOperator::ExplicitMult(const Vector &u, Vector &v) const
{
BlockVector ub(u.GetData(), offsets);
BlockVector vb(v.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &ui = ub.GetBlock(i);
Vector &vi = vb.GetBlock(i);
operators[i]->SetOperationID(OperationID::EXPLICIT_MULT);
operators[i]->ExplicitMult(ui,vi); ///< Solve the implicit system for the application
}
}
// AdditiveSchwarzCoupler methods
void AdditiveSchwarzCoupler::Mult(const Vector &x, Vector &y) const
{
/// This is use to call either ImplicitSolve or Step when Solver::Mult()
/// calls Solver.Operator::Mult()
if(GetOperationID() == OperationID::IMPLICIT_SOLVE)
{
y=x; // input vector passed as initial guess for k in ImpliicitSolve
ImplicitSolve(timestep,*input,y);
return;
}
else if(GetOperationID() == OperationID::STEP)
{
y=x; // input vector passed as initial condition in Step
real_t t_ = t, dt = timestep;
Step(y,t_,dt);
return;
}
int nops = coupled_op->Size();
const Array<int> offsets = coupled_op->GetBlockOffsets();
BlockVector xb(x.GetData(), offsets);
BlockVector yb(y.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
Vector &yi = yb.GetBlock(i);
auto op = coupled_op->GetOperator(i);
op->Transfer(xi,yi,0.0);
}
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
Vector &yi = yb.GetBlock(i);
auto op = coupled_op->GetOperator(i);
op->SetOperationID(OperationID::MULT);
op->PreProcess(xi); ///< Postprocess the data for the application
op->Mult(xi,yi);
op->PostProcess(yi); ///< Postprocess the data for the application
}
}
void AdditiveSchwarzCoupler::ImplicitSolve(const real_t dt, const Vector &x, Vector &k ) const
{
int nops = coupled_op->Size();
const Array<int> offsets = coupled_op->GetBlockOffsets();
BlockVector xb(x.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
Vector &ki = kb.GetBlock(i);
auto op = coupled_op->GetOperator(i);
op->Transfer(xi,ki,dt);
}
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
Vector &ki = kb.GetBlock(i);
auto op = coupled_op->GetOperator(i);
op->SetOperationID(OperationID::IMPLICIT_SOLVE);
op->PreProcess(xi);
op->ImplicitSolve(dt,xi,ki);
op->PostProcess(ki);
}
}
void AdditiveSchwarzCoupler::Step(Vector &x, real_t &t_, real_t &dt) const
{
int nops = coupled_op->Size();
const Array<int> offsets = coupled_op->GetBlockOffsets();
BlockVector xb(x.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
auto op = coupled_op->GetOperator(i);
op->Transfer(xi);
}
// TODO: Add time-interpolation to enable different time step for each operator;
// currently, all operators are stepped forward with the same time step
for (int i=0; i < nops; i++)
{
real_t ti = t_; ///< Store the current time
real_t dti = dt; ///< Store the current time step
Vector &xi = xb.GetBlock(i);
auto op = coupled_op->GetOperator(i);
op->SetOperationID(OperationID::STEP);
op->PreProcess(xi);
op->Step(xi,ti,dti);
op->PostProcess(xi);
}
t_ += dt; ///< Update the time after all applications have been stepped forward
///< NOTE: does not work for adaptive time-stepping
}
// AlternatingSchwarzCoupler methods
void AlternatingSchwarzCoupler::Mult(const Vector &x, Vector &y) const
{
/// This is use to call either ImplicitSolve or Step when Solver::Mult()
/// calls Solver.Operator::Mult()
if(GetOperationID() == OperationID::IMPLICIT_SOLVE)
{
y=x; // input vector passed as initial guess for k in ImpliicitSolve
ImplicitSolve(timestep,*input,y);
return;
}
else if(GetOperationID() == OperationID::STEP)
{
y=x; // input vector passed as initial condition in Step
real_t t_ = t, dt = timestep;
Step(y,t_,dt);
return;
}
int nops = coupled_op->Size();
const Array<int> offsets = coupled_op->GetBlockOffsets();
BlockVector xb(x.GetData(), offsets);
BlockVector yb(y.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
Vector &yi = yb.GetBlock(i);
auto op = coupled_op->GetOperator(i);
op->SetOperationID(OperationID::MULT);
op->PreProcess(xi);
op->Mult(xi,yi);
op->PostProcess(yi);
op->Transfer(xi,yi,0.0);
}
}
void AlternatingSchwarzCoupler::ImplicitSolve(const real_t dt, const Vector &x, Vector &k ) const
{
int nops = coupled_op->Size();
const Array<int> offsets = coupled_op->GetBlockOffsets();
BlockVector xb(x.GetData(), offsets);
BlockVector kb(k.GetData(), offsets);
for (int i=0; i < nops; i++)
{
Vector &xi = xb.GetBlock(i);
Vector &ki = kb.GetBlock(i);
auto op = coupled_op->GetOperator(i);
op->SetOperationID(OperationID::IMPLICIT_SOLVE);
op->PreProcess(xi);
op->ImplicitSolve(dt,xi,ki);
op->PostProcess(ki);
op->Transfer(xi,ki,dt);
}
}
void AlternatingSchwarzCoupler::Step(Vector &x, real_t &t_, real_t &dt) const
{
int nops = coupled_op->Size();
const Array<int> offsets = coupled_op->GetBlockOffsets();
BlockVector xb(x.GetData(), offsets);
for (int i=0; i < nops; i++)
{
real_t ti = t_; ///< Store the current time
real_t dti = dt; ///< Store the current time step
Vector &xi = xb.GetBlock(i);
auto op = coupled_op->GetOperator(i);
op->SetOperationID(OperationID::STEP);
op->PreProcess(xi);
op->Step(xi,ti,dti);
op->PostProcess(xi);
op->Transfer(xi);
}
t_ += dt; ///< Update the time after all applications have been stepped forward
///< NOTE: does not work for adaptive time-stepping
}
// JacobianFreeFullCoupler methods
void JacobianFreeFullCoupler::Mult(const Vector &k, Vector &y) const
{
add(1.0,*input,timestep,k,u); // u = u + dt*k
coupled_op->Transfer(u);
coupled_op->ImplicitMult(u,k,y); //compute residual y = f(u,k,t)
}
Operator& JacobianFreeFullCoupler::GetGradient(const Vector &k) const
{
grad.Update(k);
return const_cast<future::FDJacobian&>(grad);
}
}
#endif // MFEM_USE_MPI
File diff suppressed because it is too large Load Diff
+2
View File
@@ -21,6 +21,8 @@ set(MESH_FILES
two-squares-nurbs-autoedge.mesh
two-squares-nurbs-rot.mesh
two-squares-nurbs.mesh
3patch-nurbs.mesh
3patch-nurbs-flipedge.mesh
)
# Add a target to copy the mesh files from the source directory; used by sample
# runs.
@@ -0,0 +1,76 @@
MFEM NURBS mesh v1.0
# Same as 3patch-nurbs.mesh but with some flipped edges
# This will fail to load without CorrectPatchTopoOrientations
dimension
2
elements
3
1 3 0 1 4 3
1 3 1 2 5 4
1 3 5 6 3 4
boundary
0
edges
9
2 1 0
1 1 4
2 4 3
1 0 3
0 1 2
1 5 2
0 5 4
2 5 6
0 3 6
vertices
7
patches
knotvectors
2
1 2 0 0 1 1
1 2 0 0 1 1
dimension
2
controlpoints
0 0 1
1 0 1
0 1 1
1 1 1
knotvectors
2
1 2 0 0 1 1
1 2 0 0 1 1
dimension
2
controlpoints
1 0 1
2 0 1
1 1 1
2 2 1
knotvectors
2
1 2 0 0 1 1
1 2 0 0 1 1
dimension
2
controlpoints
2 2 1
1 2 1
1 1 1
0 1 1
+72
View File
@@ -0,0 +1,72 @@
MFEM NURBS mesh v1.0
dimension
2
elements
3
1 3 0 1 4 3
1 3 1 2 5 4
1 3 5 6 3 4
boundary
0
edges
9
2 0 1
1 1 4
2 3 4
1 0 3
0 1 2
1 2 5
0 4 5
2 6 5
0 3 6
vertices
7
patches
knotvectors
2
1 2 0 0 1 1
1 2 0 0 1 1
dimension
2
controlpoints
0 0 1
1 0 1
0 1 1
1 1 1
knotvectors
2
1 2 0 0 1 1
1 2 0 0 1 1
dimension
2
controlpoints
1 0 1
2 0 1
1 1 1
2 2 1
knotvectors
2
1 2 0 0 1 1
1 2 0 0 1 1
dimension
2
controlpoints
2 2 1
1 2 1
1 1 1
0 1 1
+2 -2
View File
@@ -27,8 +27,8 @@ accelerations, so the relationship between forces/contact pressures and
deformations/contact gaps is linear and, therefore, the problem can be solved
exactly with a single linear solve. The mortar implementation is based on [Puso
and Laursen (2004)](https://doi.org/10.1016/j.cma.2003.10.010). A description of
the Tribol implementation is available in [Serac
documentation](https://serac.readthedocs.io/en/latest/sphinx/theory_reference/solid.html#contact-mechanics).
the Tribol implementation is available in [smith
documentation](https://llnlsmith.readthedocs.io/en/latest/sphinx/theory_reference/solid.html#contact-mechanics).
Lagrange multipliers are used to solve for the pressure required to prevent
violation of the contact constraints.
+1 -2
View File
@@ -53,6 +53,7 @@ set(UNIT_TESTS_SRCS
linalg/test_ode2.cpp
linalg/test_operator.cpp
linalg/test_particlevector.cpp
linalg/test_sparsesmoothers.cpp
linalg/test_vector.cpp
mesh/mesh_test_utils.cpp
mesh/test_exodus_reader.cpp
@@ -68,8 +69,6 @@ set(UNIT_TESTS_SRCS
mesh/test_psubmesh.cpp
mesh/test_submesh.cpp
mesh/test_vtu.cpp
mesh/test_nurbs.cpp
mesh/test_exodus_writer.cpp
fem/make_permuted_mesh.cpp
fem/test_1d_bilininteg.cpp
fem/test_2d_bilininteg.cpp
@@ -0,0 +1,76 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "mfem.hpp"
#include "unit_tests.hpp"
using namespace mfem;
static void TestTranspose(const Operator &A)
{
DenseMatrix A_dense(A.Height(), A.Width());
Vector e(A.Width());
e = 0.0;
for (int i = 0; i < A.Width(); ++i)
{
e[i] = 1.0;
Vector Ae(A.Height());
A.Mult(e, Ae);
A_dense.SetCol(i, Ae);
e[i] = 0.0;
}
Vector v(A.Height());
v.Randomize();
Vector w1(A.Width()), w2(A.Width());
A.MultTranspose(v, w1);
A_dense.MultTranspose(v, w2);
w1 -= w2;
REQUIRE(w1.Normlinf() == MFEM_Approx(0.0));
}
TEST_CASE("Sparse Smoothers Transposed", "[DSmoother][GSSmoother]")
{
const bool sym = GENERATE(true, false);
constexpr int n = 10;
SparseMatrix A(n, n);
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < n; ++j)
{
const real_t val = rand_real();
A.Set(i, j, val);
if (sym) { A.Set(j, i, val); }
}
A.Add(i, i, 10.0);
}
A.Finalize();
constexpr int nit = 2; // Number of smoother iterations
TestTranspose(DSmoother(A, 0, 1.0, nit)); // scaled
if (sym)
{
TestTranspose(DSmoother(A, 1, 1.0, nit)); // l1-Jacobi
TestTranspose(DSmoother(A, 2, 1.0, nit)); // lumped Jacobi
}
TestTranspose(GSSmoother(A, 0, nit)); // symmetric
TestTranspose(GSSmoother(A, 1, nit)); // forward
TestTranspose(GSSmoother(A, 2, nit)); // backward
}
+8
View File
@@ -132,6 +132,14 @@ TEST_CASE("NURBS mesh reconstruction", "[NURBS]")
for (auto *p : patches) { delete p; }
}
TEST_CASE("NURBS knotvector orientation", "[NURBS]")
{
// This will fail to load without CorrectPatchTopoOrientations
auto mesh_fname = "../../miniapps/nurbs/meshes/3patch-nurbs-flipedge.mesh";
Mesh mesh(mesh_fname, 1, 1);
REQUIRE(mesh.NURBSext->CheckPatches());
}
TEST_CASE("NURBS NC-patch mesh loading", "[NURBS]")
{
auto mesh_fname = GENERATE("../../data/nc3-nurbs.mesh",