magnetotellurics/ModEM: 2023-11-19-rev755 - New default settings for SP/SP2 Solvers (BICG)

Note: This is a release on the classic branch (SVN stable branch). Although this revision is after release 2022-07-07-rev630 (that was released on the main branch), it does not include the GPU solvers.

This version is based on the latest stable code but includes IMPORTANT technical changes relative to stable -rev 710, namely:

FWD,FWD_SP,FWD_SP2/EMsolve3D.f90 - NEW DEFAULTS; unless the solver is specified in the forward solver configuration file, BICG is now the default for sparse-matrix formulations SP & SP2. QMR is still the default for the traditional matrix-free formulation.

MPI/Main_MPI.f90 - incorporated all the deallocate statements from Naser's version of the code to, hopefully, eliminate the memory leaks; also cleaned up the formatting in some obvious places. Note that this version does not include the initialization of grid elements l_E and S_F. This is by design - this is now correctly done in modelOperator3D.f90 for all flavors of the code (FWD, FWD_SP, FWD_SP2, FWD_SPETSc2).

MPI/Declaration_MPI.f90 - added subroutine gather_runtime and Naser's solver_name, period and solver_residual_vec variables.

3D_MT/Sub_MPI.f90 - added rFile_Model1D variable to allow for secondary field formulation in the near future (this version is not fully equipped with it).

FWD,FWD_SP,FWD_SP2,FWD_SPETSc2/modelOperator3D.f90 - added initialization and deallocation of metric elements (rvectors & cvectors in addition to arrays used previously by SP*). These are used repeatedly in EMfieldInterp.f90 and ModelMap.f90 for all, and for other things in FWD as well. Note that this duplicates the storage of 5 metric element arrays. Currently unavoidable (to support spherical coords) unless the code is rewritten. Previous versions used these arrays as well but they were allocated in the wrong places.

SP_Topology/MetricElements.f90 - added Vcell variable to support modelOperator3D.f90 metric element initialization (needed for ModelMap.f90 in SP variants of the code).

FWD/modelParam/modelMap.inc - removed metric elements (V_C and V_E) initialization, now correctly done in modelOperator3D.f90; incorporated Gary's changes that remove the ghost nodes which I no longer remember anything about and which messed up Jmult for Gary.

UTILS/elements.f90 - added the logic for boundary elements and, importantly, replaced integral approximations with finite differences for the area and volume elements. This has NOT been carefully verified to be better. Needs to be carefully tested. This module is not used by the traditional 3D MT but it will be needed for the spherical coords modeling, which is not yet part of this distribution since it requires additional development before it can be released to the public.

--

Technical notes:

This version at present does NOT pass the symmetry tests, specifically Jmult. Matrix-free is broken relative to how well it passed the test on 2022-11-17 and SP2 is substantially broken just like it used to be. All versions compile and run, tested both in serial and in parallel. Additional length, area and volume elements storage adds just below 7% memory usage to matrix-free and 10% to sparse matrix formulations, based on a typical realistic parallel forward modeling example, because they are currently stored both in matrix form and as vectors, for different purposes - can be improved upon with code rewrite. A much more substantial memory reduction has been achieved by removing some of the memory leaks from the MPI code, based on Naser Meqbel's and Paulo Werdt's analysis.

SP2 BICG formulation is be far the most efficient at present. For a realistic MPI example with an optimized Intel-compiled code, forward modeling completes in 14.46min with SP2 BiCG compared to 1h 49min with matrix-free QMR. However, matrix-free code requires approx 2.5 times less operational memory than SP2. SP2 QMR runs 3.5 - 10 times SLOWER than SP2 BICG. The matrix-free QMR works while the matrix-free BICG never converges and is truly broken.

All these observations equally apply to stable -rev 710 that this is based upon, so I'm committing to the svn with the intention of figuring out what has broken the symmetry in a separate release sometime later - not the first priority right now.

Both the current stable version and this version pass the symmetry test perfectly using the old (broken!) model and air layers, unit_testing. It will take some time to understand why the new, corrected, model and air layers - deeper lower boundary and better resolution in the air - does not pass the Jmult symmetry test as perfectly. After wasting a couple of days on this again, will postpone for now.

And of course, the latest fix introduced a stupid MPI bug, which is no longer as of this commit. The forward modeling is tested to run efficiently in parallel (FWD and SP2).