Merge branch 'master' into libxsmm-jit

CMakeLists.txt

@@ -11,6 +11,13 @@ endif()
 
 project(SeisSol LANGUAGES C CXX Fortran)
 
+if (CMAKE_CXX_COMPILER_ID MATCHES "NVHPC|PGI")
+  set(NVHPC_REQUIRED_VERSION "22.1.0")
+  if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS NVHPC_REQUIRED_VERSION)
+    message(FATAL_ERROR "NVHPC version ${NVHPC_REQUIRED_VERSION} or higher is required")
+  endif()
+endif()
+
 # set hardware specific definition needed for seissol compilation
 # 'process_users_input' returns the following:
 #
@@ -307,13 +314,21 @@ target_compile_definitions(SeisSol-lib PUBLIC
 )
 
 # Fortran compliler settings
-set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -cpp")
 if ("${CMAKE_Fortran_COMPILER_ID}" STREQUAL "GNU")
-  target_compile_options(SeisSol-lib PUBLIC $<$<COMPILE_LANGUAGE:Fortran>:-ffree-line-length-none -fdefault-real-8 -Wno-unused-parameter>)
+  target_compile_options(SeisSol-lib PUBLIC $<$<COMPILE_LANGUAGE:Fortran>:-cpp -ffree-line-length-none -fdefault-real-8 -Wno-unused-parameter>)
 elseif ("${CMAKE_Fortran_COMPILER_ID}" STREQUAL "Intel")
   # todo intel, is needed:  -align -align array64byte
   # todo -r8 -WB is needed for intel (8 byte precision for reals)
   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -cpp -r8 -WB")
+elseif(CMAKE_Fortran_COMPILER_ID MATCHES "NVHPC|PGI")
+    # NOTE:
+    # 1. -Mpreprocess enables a C preprocessor macro
+    # 2. -Mfree specifies fortran free format
+    # 3. -r8 Interpret REAL variables as DOUBLE PRECISION
+    # 4. BG (see, Initializer/preProcessorMacros.fpp) around logDebug, logInfo, etc.
+  target_compile_options(SeisSol-lib PUBLIC $<$<COMPILE_LANGUAGE:Fortran>:-Mpreprocess -Mfree -r8 -DBG>)
+elseif ("${CMAKE_Fortran_COMPILER_ID}" STREQUAL "Flang")
+  target_compile_options(SeisSol-lib PUBLIC $<$<COMPILE_LANGUAGE:Fortran>:-cpp -ffree-form -fdefault-real-8 -DBG>)
 endif()
 
 # C++ compiler settings
@@ -328,8 +343,19 @@ elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
 
   # Activate interprocedual optimization.
   #set_property(TARGET SeisSol-lib PROPERTY INTERPROCEDURAL_OPTIMIZATION TRUE) 
-elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+elseif (CMAKE_CXX_COMPILER_ID MATCHES "Clang|IntelLLVM")
   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++17 -Wall -Wextra -pedantic -Wno-unused-parameter")
+
+elseif(CMAKE_CXX_COMPILER_ID MATCHES "NVHPC|PGI")
+    set(WARNINGS "--display_error_number --diag_suppress186")
+    if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS "22.1.1")
+      set(WARNINGS "${WARNINGS} --diag_suppress1")
+    endif()
+
+    # NOTE:
+    # 1. --pending_instantiations=0 allows an infinite recursive template instantiation
+    # 2. EIGEN_DONT_VECTORIZE=1 waiting for eigen3 support for nvhpc compiler collection w.r.t. vectorization
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++17 -Wc,--pending_instantiations=0 ${WARNINGS} -DEIGEN_DONT_VECTORIZE=0")
 endif()
 
 find_package(FILESYSTEM REQUIRED)

Documentation/generating-a-megathrust-geometry.rst

@@ -0,0 +1,112 @@
+Proposed workflow for generating a CAD model of a megathrust earthquake
+========================================================================
+
+Here is present a workflow for generating a CAD model of the Japan subduction, including subduction interface, to be used in dynamic rupture models or for imposing a kinematic model on the subduction interface.
+We use scripts from https://github.com/SeisSol/Meshing. To best follow this tutorial, we suggest adding the geometry script folder to the path:
+
+.. code-block:: bash
+
+    export PATH=$PATH:~/SeisSol/Meshing/creating_geometric_models
+
+
+topography and slab interface
+-----------------------------
+
+First, we download topography and bathymetry data from GEBCO
+(`http://www.gebco.net/ <http://www.gebco.net/>`__), and we triangulate it into a GoCad ts file.
+Note that we downsample the topography data by a factor 5 for dealing with a reasonable size dataset in this tutorial.
+Note also that we use a custom transverse Mercator roughly centered at the domain center. 
+
+.. code-block:: bash
+
+    #!/bin/sh
+    myproj='+proj=tmerc +datum=WGS84 +k=0.9996 +lon_0=143 +lat_0=39'
+    topofile='data/gebco_2021_n42.61596679687501_s34.16381791234017_w137.10571333765984_e146.759033203125.nc'
+    create_surface_from_rectilinear_grid.py $topofile tmp/topo.ts --proj "$myproj" --sub 5
+
+Next, we generate a mesh of the subduction interface.
+We use the Kamchatka-Kuril Islands-Japan Region model of Slab2.0 available here `<https://www.sciencebase.gov/catalog/item/5aa4060de4b0b1c392eaaee2>`__.
+We download kur_slab2_dep_02.24.18.grd and rename it as kur_slab2_dep_02.24.18.nc.
+Finally, we crop and triangulate it using the following command:
+
+.. code-block:: bash
+
+    create_surface_from_rectilinear_grid.py data/kur_slab2_dep_02.24.18.nc tmp/ur_slab2_dep_02.24.18.ts --crop 140 145 35.5 41 --proj "$myproj"
+
+The Slab2.0 data are sampled every 5 km, but a dynamic rupture simulation may require one order of magnitude smaller resolution, or even more.
+If we mesh at 500m a geometry sampled at 5 km, the obtained mesh won't be smooth but will have edges every 5 km.
+Therefore, we smooth and refine the mesh with:
+
+.. code-block:: bash
+
+    refine_and_smooth_mesh.py --N 0 --P 1 tmp/ur_slab2_dep_02.24.18.ts
+
+The help (-h option) offers further details about the script options.
+Here we use P=1 to avoid dealing with a too large dataset in the tutorial.
+
+Terminal splay
+--------------
+
+The shallowest part of the subduction interface is 5-10km below the sea surface. We might want to acknowledge the possibility of surface rupture by extendiing the interface to the surface at steeper angle.
+We therefore build a terminal splay extending from the trace of the know trench location.
+For that, we first preprocess a USGS trench data file into a fault trace. Then we extend the trace along a constant dip angle.
+
+.. code-block:: bash
+
+    generate_trace_terminal_splay.py --shift -0.05 0 --filter_trench 137.2 146.7 35.2 41.2 jap --usgs_trench_file data/trenches_usgs_2017_depths.csv --bathymetry $topofile --plot tmp/extractedtrenchTohoku.dat
+    create_fault_from_trace.py tmp/extractedtrenchTohoku.dat 0 30 --proj "$myproj" --maxdepth  20e3 --extend 8000 --extrudeDir data/strike90.dat --dd 5000 --smoothingParameter 1e7
+    refine_and_smooth_mesh.py --N 0 --P 1 extractedtrenchTohoku0.ts
+
+data/strike90.dat contains 2 lines, indicating that the extrude direction is towards West.
+
+.. code-block:: 
+
+    0 -90
+    1 -90
+
+
+Merging terminal splay and subduction interface
+-----------------------------------------------
+
+First, we generate 2 planes that we will use to trim slab interface and extension to the same dimension.
+
+.. code-block:: bash
+
+    python generate_cutting_plane.py  --center -65000 -367000 -11000  --normal 0. 1. 0.  --dims 1000e3 300e3 --mesh 10e3 tmp/cut1.stl
+    python generate_cutting_plane.py  --center 125000 210000 -11000  --normal 0. 1. 0.  --dims 1000e3 300e3 --mesh 10e3 tmp/cut2.stl
+
+Then we load both processed subduction interface and terminal splay into SimModeler.
+We intersect them, and trim what can be trimmed.
+Then we import the cutting plane, and intersect all parts.
+We finally remove all parts that are not faults.
+At this point we have 2 connected surfaces, subduction interface and terminal splay, and a sharp angle between them.
+This sharp angle can be smoothed by extracting the discrete mesh as inp file (see :ref:`Remeshing the topography`), converting to ts and applying the refine_and_smooth_mesh.py script.
+
+
+.. figure:: LatexFigures/Tohoku_building_fault_model.png
+   :alt: Discrete surface for generating Tohoku's fault model 
+   :width: 11.00000cm
+   :align: center
+
+   Discrete geometry model of the slab interface (red), a potential terminal splay fault (yellow), and the cutting planes (blue and green) used to trim laterally these surfaces, before interesection.
+
+
+Final steps
+-----------
+
+Finally, we create a box mesh box domain with pygmsh as follow:
+
+.. code-block:: 
+
+    generate_box.py --proj "$myproj" --rangeFromTopo $topofile tmp/box.stl --zdim " -500e3" 5e3 --shrink 0.9
+
+The final step consists in intersecting all objects (topography, faults and domain box) in the GUI of SimModeler, as presented in :doc:`simmodelerCAD-workflow`.
+
+.. figure:: LatexFigures/Tohoku_final_cut_view.png
+   :alt: Cut view of the final Tohoku's model
+   :width: 11.00000cm
+   :align: center
+
+   Cut view of the final Tohoku's model
+
+

Documentation/index.rst

@@ -98,6 +98,7 @@ Characteristics of the SeisSol simulation software are:
 
   simmodelerCAD-workflow
   generating-a-cad-model-using-gocad-basic-tutorial
+  generating-a-megathrust-geometry
   remeshing-the-topography
   adapting-the-cad-model-resolution-using-gocad
   manually-fixing-an-intersection-in-gocad

Documentation/remeshing-the-topography.rst

@@ -1,3 +1,5 @@
+.. _Remeshing the topography:
+
 Remeshing the topography
 ========================
 

Documentation/simmodelerCAD-workflow.rst

@@ -15,7 +15,7 @@ The dataset for generating the Palu structural model is available `here <https:/
 Creating the topographic layer
 ------------------------------
 
-We create the topography from a netcdf file downloaded from https://www.gebco.net/
+We create the topography from a netcdf file downloaded from https://www.gebco.net/.
 The domain range from longitude 118.9 to 121.7 and from latitude -2.4 to 1.0.
 We then project the data (see :ref:`On the use of projections` for the choice of a projection), triangulate it, and export it as stl (list of triangles) using:
 

cmake/cpu_arch_flags.cmake

@@ -13,7 +13,7 @@ function(get_arch_flags architecture compiler)
     elseif ("${HOST_ARCH}" STREQUAL "hsw")
         if (compiler STREQUAL "Intel")
             set(CPU_ARCH_FLAGS "-xCORE-AVX2" "-fma" PARENT_SCOPE)
-        elseif(compiler MATCHES "GNU|Clang")
+        elseif(compiler MATCHES "GNU|Clang|IntelLLVM")
             set(CPU_ARCH_FLAGS "-mavx2" "-mfma" PARENT_SCOPE)
         endif()
 
@@ -25,15 +25,15 @@ function(get_arch_flags architecture compiler)
     elseif ("${HOST_ARCH}" STREQUAL "knl")
         if (compiler STREQUAL "Intel")
             set(CPU_ARCH_FLAGS "-xMIC-AVX512" "-fma" PARENT_SCOPE)
-        elseif(compiler MATCHES "GNU|Clang")
+        elseif(compiler MATCHES "GNU|Clang|IntelLLVM")
             set(CPU_ARCH_FLAGS "-mavx512f" "-mavx512cd" "-mavx512pf" "-mavx512er" "-mfma" PARENT_SCOPE)
         endif()
 
     # Skylake cpu architecture
     elseif ("${HOST_ARCH}" STREQUAL "skx")
         if (compiler STREQUAL "Intel")
             set(CPU_ARCH_FLAGS "-xCORE-AVX512" "-fma" PARENT_SCOPE)
-        elseif(compiler MATCHES "GNU|Clang")
+        elseif(compiler MATCHES "GNU|Clang|IntelLLVM")
             set(CPU_ARCH_FLAGS "-march=skylake-avx512" PARENT_SCOPE)
         endif()
 
@@ -50,7 +50,7 @@ function(get_arch_flags architecture compiler)
     elseif ("${HOST_ARCH}" STREQUAL "rome")
         if (compiler STREQUAL "Intel")
             set(CPU_ARCH_FLAGS "-march=core-avx2" "-fma" PARENT_SCOPE)
-        elseif(compiler MATCHES "GNU|Clang")
+        elseif(compiler MATCHES "GNU|Clang|IntelLLVM")
             set(CPU_ARCH_FLAGS "-march=znver2" "-mtune=znver2" PARENT_SCOPE)
     endif()
 
@@ -61,4 +61,9 @@ function(get_arch_flags architecture compiler)
 
     endif()
 
+    if (compiler MATCHES "NVHPC|PGI")
+        #NOTE: PGI-based compiler does not have `-mno-red-zone` flag
+        set(HAS_REDZONE OFF PARENT_SCOPE)
+    endif()
+
 endfunction()

src/Geometry/MeshReaderCBinding.f90

@@ -138,7 +138,7 @@ subroutine read_mesh_fast(io, eqn, disc, mesh, bnd, mpi)
 
         write(str, *) mpi%nCPU
         if (io%meshgenerator .eq. 'Gambit3D-fast') then
-            call read_mesh_gambitfast_c(mpi%myRank, trim(io%MeshFile) // c_null_char, \
+            call read_mesh_gambitfast_c(mpi%myRank, trim(io%MeshFile) // c_null_char, &
                 trim(io%MetisFile) // '.epart.' // trim(adjustl(str)) // c_null_char, hasFault, MESH%Displacement(:), m_mesh%ScalingMatrix(:,:))
         elseif (io%meshgenerator .eq. 'Netcdf') then
 #ifdef PARALLEL

src/Initializer/dg_setup.f90

@@ -425,9 +425,9 @@ SUBROUTINE iniGalerkin3D_us_level2_new(OptionalFields,EQN,DISC,MESH,BND,IC,SOURC
     iSide = 0
 
     materialVal = OptionalFields%BackgroundValue(iElem,:)
-    call c_interoperability_setMaterial( i_elem = iElem,                                         \
-                                         i_side = iSide,                                         \
-                                         i_materialVal = materialVal,                            \
+    call c_interoperability_setMaterial( i_elem = iElem,                                         &
+                                         i_side = iSide,                                         &
+                                         i_materialVal = materialVal,                            &
                                          i_numMaterialVals = EQN%nBackgroundVar                  )
 
     do iSide = 1,4
@@ -444,9 +444,9 @@ SUBROUTINE iniGalerkin3D_us_level2_new(OptionalFields,EQN,DISC,MESH,BND,IC,SOURC
               materialVal = OptionalFields%BackgroundValue(iElem,:)
           END SELECT
       ENDIF
-      call c_interoperability_setMaterial( i_elem = iElem,                        \
-                                           i_side = iSide,                        \
-                                           i_materialVal = materialVal,       \
+      call c_interoperability_setMaterial( i_elem = iElem,                        &
+                                           i_side = iSide,                        &
+                                           i_materialVal = materialVal,           &
                                            i_numMaterialVals = EQN%nBackgroundVar )
 
     enddo
@@ -905,11 +905,11 @@ SUBROUTINE icGalerkin3D_us_new(EQN, DISC, MESH, IC, SOURCE, IO, MPI)
 
            ! initialize loading in C
            oneRankedShaped_iniloading = pack( l_initialLoading, .true. )
-           call c_interoperability_setInitialLoading( i_meshId = iElem, \
+           call c_interoperability_setInitialLoading( i_meshId = iElem, &
                                                       i_initialLoading = oneRankedShaped_iniloading)
 
            !initialize parameters in C
-           call c_interoperability_setPlasticParameters( i_meshId            = iElem, \
+           call c_interoperability_setPlasticParameters( i_meshId            = iElem, &
                                                          i_plasticParameters = l_plasticParameters )
 
        END IF