updates workflow script (for macOS), Makefiles and SIEM mesh routines; adds explicit type conversions

.github/workflows/CI.yml

@@ -104,16 +104,27 @@ jobs:
         echo ""
         ## avoids sed -i '' issue on MacOS, using gnu sed to have the same sed command lines as in linux
         brew install gnu-sed
-        echo "PATH=/usr/local/opt/gnu-sed/libexec/gnubin:$PATH" >> $GITHUB_ENV
+        echo "PATH=/opt/homebrew/opt/gnu-sed/libexec/gnubin:$PATH" >> $GITHUB_ENV
         ## OpenMP
         echo "OMP_NUM_THREADS=2" >> $GITHUB_ENV
+        ## MPI
+        echo
+        echo "MPI environment:"
+        echo "mpif90 on path: $(which mpif90)"
+        echo
+        # add OpenMPI path to have ./configure detect mpi.h
+        echo "MPI_INC=/opt/homebrew/include" >> $GITHUB_ENV
         ## avoids MPI issue with number of slots
         echo "OMPI_MCA_rmaps_base_oversubscribe=1" >> $GITHUB_ENV
         echo "OMPI_MCA_rmaps_base_inherit=1" >> $GITHUB_ENV
         ## avoids MPI issue when running in MacOS
         echo "OMPI_MCA_btl=self,tcp" >> $GITHUB_ENV
+        # newer OpenMPI version (5.x)
+        echo "PRTE_MCA_rmaps_default_mapping_policy=:oversubscribe" >> $GITHUB_ENV
         # exports for xterm output (for make tests)
         echo "TERM=xterm" >> $GITHUB_ENV
+        # warning: (arm64)  could not find object file symbol for symbol ___mulsc3
+        echo "LDFLAGS=-lgcc" >> $GITHUB_ENV
         echo ""
         echo "exports:"
         export

Makefile.in

@@ -54,6 +54,10 @@ CXXFLAGS = -I${SETUP} @CXXFLAGS@
 
 FCLINK = $(MPIFCCOMPILE_CHECK)
 
+# all linker flags
+LDFLAGS = @LDFLAGS@
+MPILIBS += $(LDFLAGS) @LIBS@
+
 #######################################
 ####
 #### MPI

src/auxiliaries/create_movie_AVS_DX.f90

@@ -894,11 +894,11 @@
             ipoin = ipoin + 1
             ireorder(ibool_number) = ipoin
             ! point value
-            total_dat(ipoin) = field_display(ilocnum+ieoff)
+            total_dat(ipoin) = real(field_display(ilocnum+ieoff),kind=CUSTOM_REAL)
             ! point location
-            total_dat_xyz(1,ipoin) = xp_save(ilocnum+ieoff)
-            total_dat_xyz(2,ipoin) = yp_save(ilocnum+ieoff)
-            total_dat_xyz(3,ipoin) = zp_save(ilocnum+ieoff)
+            total_dat_xyz(1,ipoin) = real(xp_save(ilocnum+ieoff),kind=CUSTOM_REAL)
+            total_dat_xyz(2,ipoin) = real(yp_save(ilocnum+ieoff),kind=CUSTOM_REAL)
+            total_dat_xyz(3,ipoin) = real(zp_save(ilocnum+ieoff),kind=CUSTOM_REAL)
           endif
           mask_point(ibool_number) = .true.
         enddo

src/auxiliaries/create_movie_GMT_global.f90

@@ -358,9 +358,9 @@
       close(IIN)
     endif
     ! effective half duration in movie runs
-    hdur = sqrt( cmt_hdur**2 + HDUR_MOVIE**2)
+    hdur = real(sqrt( cmt_hdur**2 + HDUR_MOVIE**2),kind=CUSTOM_REAL)
     ! start time of simulation
-    t0 = - 1.5d0*( cmt_t_shift - hdur )
+    t0 = real(-1.5d0*( cmt_t_shift - hdur ),kind=CUSTOM_REAL)
 
     ! becomes time (s) from hypocenter to reach surface (using average 8 km/s s-wave speed)
     ! note: especially for deep sources, this helps determine a better starttime to mute
@@ -383,15 +383,15 @@
 
     ! converts values into radians
     ! colatitude [0, PI]
-    LAT_SOURCE = (90.0 - LAT_SOURCE)*DEGREES_TO_RADIANS
+    LAT_SOURCE = real((90.d0 - LAT_SOURCE)*DEGREES_TO_RADIANS,kind=CUSTOM_REAL)
 
     ! longitude [-PI, PI]
     if (LON_SOURCE < -180.0 ) LON_SOURCE = LON_SOURCE + 360.0
     if (LON_SOURCE > 180.0 ) LON_SOURCE = LON_SOURCE - 360.0
-    LON_SOURCE = LON_SOURCE * DEGREES_TO_RADIANS
+    LON_SOURCE = real(LON_SOURCE * DEGREES_TO_RADIANS,kind=CUSTOM_REAL)
 
     ! mute radius in rad
-    RADIUS_TO_MUTE = RADIUS_TO_MUTE * DEGREES_TO_RADIANS
+    RADIUS_TO_MUTE = real(RADIUS_TO_MUTE * DEGREES_TO_RADIANS,kind=CUSTOM_REAL)
   endif
 
   print *,'--------'
@@ -468,7 +468,7 @@
 
         ! approximate wavefront travel distance in degrees
         ! (~3.5 km/s wave speed for surface waves)
-        distance = SURFACE_WAVE_VELOCITY * ((it-1)*DT-t0) / (R_PLANET/1000.d0) * RADIANS_TO_DEGREES
+        distance = real(SURFACE_WAVE_VELOCITY * ((it-1)*DT-t0) / (R_PLANET/1000.d0) * RADIANS_TO_DEGREES,kind=CUSTOM_REAL)
 
         print *,'distance approximate: ',distance,'(degrees)'
 
@@ -498,7 +498,7 @@
         print *,'muting radius: ',0.7 * distance,'(degrees)'
 
         ! new radius of mute area (in rad)
-        RADIUS_TO_MUTE = 0.7 * distance * DEGREES_TO_RADIANS
+        RADIUS_TO_MUTE = real(0.7d0 * distance * DEGREES_TO_RADIANS,kind=CUSTOM_REAL)
       else
         ! mute_factor used at the beginning for scaling displacement values
         if (STARTTIME_TO_MUTE > TINYVAL) then
@@ -507,9 +507,9 @@
           ! linear scaling between [0,1]:
           ! from 0 (simulation time equal to zero )
           ! to 1 (simulation time equals starttime_to_mute)
-          mute_factor = 1.0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0)
+          mute_factor = real(1.d0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0),kind=CUSTOM_REAL)
           ! threshold value for mute_factor
-          if (mute_factor < TINYVAL ) mute_factor = TINYVAL
+          if (mute_factor < TINYVAL ) mute_factor = real(TINYVAL,kind=CUSTOM_REAL)
           if (mute_factor > 1.0 ) mute_factor = 1.0
         endif
       endif
@@ -611,17 +611,17 @@
                   ! checks source longitude range
                   if (LON_SOURCE - RADIUS_TO_MUTE < -PI .or. LON_SOURCE + RADIUS_TO_MUTE > PI) then
                     ! source close to 180. longitudes, shifts range to [0, 2PI]
-                    if (phival < 0.0 ) phival = phival + TWO_PI
+                    if (phival < 0.0 ) phival = phival + real(TWO_PI,kind=CUSTOM_REAL)
                     if (LON_SOURCE < 0.0) then
-                      dist_lon = phival - (LON_SOURCE + TWO_PI)
+                      dist_lon = phival - real((LON_SOURCE + TWO_PI),kind=CUSTOM_REAL)
                     else
                       dist_lon = phival - LON_SOURCE
                     endif
                   else
                     ! source well between range to [-PI, PI]
                     ! shifts phival to be like LON_SOURCE between [-PI,PI]
-                    if (phival > PI ) phival = phival - TWO_PI
-                    if (phival < -PI ) phival = phival + TWO_PI
+                    if (phival > PI ) phival = phival - real(TWO_PI,kind=CUSTOM_REAL)
+                    if (phival < -PI ) phival = phival + real(TWO_PI,kind=CUSTOM_REAL)
 
                     dist_lon = phival - LON_SOURCE
                   endif
@@ -634,7 +634,7 @@
                     if ((it-1)*DT-t0 > STARTTIME_TO_MUTE) then
                       ! wavefield will be tapered to mask out noise in source area
                       ! factor from 0 to 1
-                      mute_factor = ( 0.5*(1.0 - cos(distance/RADIUS_TO_MUTE*PI)) )**6
+                      mute_factor = real( ( 0.d5*(1.d0 - cos(distance/RADIUS_TO_MUTE*PI)) )**6,kind=CUSTOM_REAL)
                       ! factor from 0.01 to 1
                       mute_factor = mute_factor * 0.99 + 0.01
                       displn(i,j) = displn(i,j) * mute_factor
@@ -721,11 +721,11 @@
 
                 ! determines North / South pole index for stamping maximum values
                 if (USE_AVERAGED_MAXIMUM .and. AVERAGE_NORMALIZE_VALUES) then
-                  xmesh = xp(ieoff)
-                  ymesh = yp(ieoff)
-                  zmesh = zp(ieoff)
-                  if (zmesh > -SMALL_VAL_ANGLE .and. zmesh <= ZERO) zmesh = -SMALL_VAL_ANGLE
-                  if (zmesh < SMALL_VAL_ANGLE .and. zmesh >= ZERO) zmesh = SMALL_VAL_ANGLE
+                  xmesh = real(xp(ieoff),kind=CUSTOM_REAL)
+                  ymesh = real(yp(ieoff),kind=CUSTOM_REAL)
+                  zmesh = real(zp(ieoff),kind=CUSTOM_REAL)
+                  if (zmesh > -SMALL_VAL_ANGLE .and. zmesh <= ZERO) zmesh = - real(SMALL_VAL_ANGLE,kind=CUSTOM_REAL)
+                  if (zmesh < SMALL_VAL_ANGLE .and. zmesh >= ZERO) zmesh = real(SMALL_VAL_ANGLE,kind=CUSTOM_REAL)
                   thetaval = atan2(sqrt(xmesh*xmesh+ymesh*ymesh),zmesh)
                   ! thetaval between 0 and PI / 2
                   !print *,'thetaval:',thetaval * 180. / PI
@@ -812,7 +812,7 @@
 
         if (max_absol < max_average) then
           ! distance (in degree) of surface waves travelled
-          distance = SURFACE_WAVE_VELOCITY * ((it-1)*DT-t0) / (R_PLANET/1000.d0) * RADIANS_TO_DEGREES
+          distance = real(SURFACE_WAVE_VELOCITY * ((it-1)*DT-t0) / (R_PLANET/1000.d0) * RADIANS_TO_DEGREES,kind=CUSTOM_REAL)
           if (distance > 10.0 .and. distance <= 20.0) then
             ! smooth transition between 10 and 20 degrees
             ! sets positive and negative maximum
@@ -866,15 +866,15 @@
               ! linear scaling between [0,1]:
               ! from 0 (simulation time equal to -t0 )
               ! to 1 (simulation time equals starttime_to_mute)
-              mute_factor = 1.0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0)
+              mute_factor = real(1.d0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0),kind=CUSTOM_REAL)
               ! takes complement and shifts scale to (1,100)
               ! thus, mute factor is 100 at simulation start and 1.0 at starttime_to_mute
               mute_factor = abs(1.0 - mute_factor) * 99.0 + 1.0
               ! positive and negative maximum reach average when wavefield appears
-              val = mute_factor * max_average
+              val = real(mute_factor * max_average,kind=CUSTOM_REAL)
             else
               ! uses a constant factor
-              val = 100.0 * max_average
+              val = real(100.d0 * max_average,kind=CUSTOM_REAL)
             endif
             ! positive and negative maximum
             field_display(istamp1) = + val

src/auxiliaries/rules.mk

@@ -216,7 +216,7 @@ xcombine_paraview_strain_data_SHARED_OBJECTS = \
 	$(EMPTY_MACRO)
 
 ${E}/xcombine_paraview_strain_data: $(xcombine_paraview_strain_data_OBJECTS) $(xcombine_paraview_strain_data_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_paraview_strain_data.auxsolver.o: $O/specfem3D_par.solverstatic_module.o
@@ -240,7 +240,7 @@ xcombine_surf_data_SHARED_OBJECTS = \
 	$(EMPTY_MACRO)
 
 ${E}/xcombine_surf_data: $(xcombine_surf_data_OBJECTS) $(xcombine_surf_data_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_surf_data.auxsolver.o: $O/specfem3D_par.solverstatic_module.o
@@ -267,7 +267,7 @@ xcombine_surf_data_vtk_SHARED_OBJECTS = \
 	$(EMPTY_MACRO)
 
 ${E}/xcombine_surf_data_vtk: $(xcombine_surf_data_vtk_OBJECTS) $(xcombine_surf_data_vtk_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_surf_data.auxsolver_vtk.o: $O/specfem3D_par.solverstatic_module.o
@@ -282,7 +282,7 @@ xcombine_surf_data_vtu_OBJECTS = \
 xcombine_surf_data_vtu_SHARED_OBJECTS = $(xcombine_surf_data_vtk_SHARED_OBJECTS)
 
 ${E}/xcombine_surf_data_vtu: $(xcombine_surf_data_vtu_OBJECTS) $(xcombine_surf_data_vtu_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_surf_data.auxsolver_vtu.o: $O/specfem3D_par.solverstatic_module.o
@@ -317,7 +317,7 @@ xcombine_vol_data_SHARED_OBJECTS = \
 	$(EMPTY_MACRO)
 
 ${E}/xcombine_vol_data: $(xcombine_vol_data_OBJECTS) $(xcombine_vol_data_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_vol_data.auxsolver.o: $O/specfem3D_par.solverstatic_module.o
@@ -360,7 +360,7 @@ xcombine_vol_data_adios_SHARED_OBJECTS = \
 $O/combine_vol_data.auxadios.o: $O/combine_vol_data_adios_impl.auxadios.o
 
 ${E}/xcombine_vol_data_adios: $(xcombine_vol_data_adios_OBJECTS) $(xcombine_vol_data_adios_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_vol_data.auxadios.o: $O/specfem3D_par.solverstatic_module.o
@@ -396,7 +396,7 @@ xcombine_vol_data_vtk_SHARED_OBJECTS = \
 	$(EMPTY_MACRO)
 
 ${E}/xcombine_vol_data_vtk: $(xcombine_vol_data_vtk_OBJECTS) $(xcombine_vol_data_vtk_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_vol_data.auxsolver_vtk.o: $O/specfem3D_par.solverstatic_module.o
@@ -438,7 +438,7 @@ xcombine_vol_data_vtk_adios_SHARED_OBJECTS = \
 	$(EMPTY_MACRO)
 
 ${E}/xcombine_vol_data_vtk_adios: $(xcombine_vol_data_vtk_adios_OBJECTS) $(xcombine_vol_data_vtk_adios_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_vol_data.auxadios_vtk.o: $O/combine_vol_data_adios_impl.auxadios.o
@@ -454,7 +454,7 @@ xcombine_vol_data_vtu_OBJECTS = \
 xcombine_vol_data_vtu_SHARED_OBJECTS = $(xcombine_vol_data_vtk_SHARED_OBJECTS)
 
 ${E}/xcombine_vol_data_vtu: $(xcombine_vol_data_vtu_OBJECTS) $(xcombine_vol_data_vtu_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_vol_data.auxsolver_vtu.o: $O/specfem3D_par.solverstatic_module.o
@@ -469,7 +469,7 @@ xcombine_vol_data_vtu_adios_OBJECTS = \
 xcombine_vol_data_vtu_adios_SHARED_OBJECTS = $(xcombine_vol_data_vtk_adios_SHARED_OBJECTS)
 
 ${E}/xcombine_vol_data_vtu_adios: $(xcombine_vol_data_vtu_adios_OBJECTS) $(xcombine_vol_data_vtu_adios_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/combine_vol_data.auxadios_vtu.o: $O/combine_vol_data_adios_impl.auxadios.o
@@ -561,7 +561,7 @@ xextract_database_SHARED_OBJECTS = \
 	$(EMPTY_MACRO)
 
 ${E}/xextract_database: $(xextract_database_OBJECTS) $(xextract_database_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ### additional dependencies
 $O/extract_database.aux.o: $O/specfem3D_par.solverstatic_module.o
@@ -704,7 +704,7 @@ endif
 
 
 ${E}/xwrite_profile: $(xwrite_profile_OBJECTS) $(xwrite_profile_SHARED_OBJECTS)
-	${MPIFCCOMPILE_CHECK} -o $@ $+ $(LDFLAGS) $(MPILIBS) $(LIBS)
+	${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
 
 ## additional module dependencies
 $O/write_profile.aux.o: $O/meshfem3D_models.check.o

src/auxiliaries/write_profile.f90

@@ -1211,8 +1211,8 @@
   integer,save :: slice_proc_old = 0
 
   ! converts to lat/lon
-  lat = 90.0 - theta_degrees
-  lon = phi_degrees
+  lat = real(90.d0 - theta_degrees,kind=CUSTOM_REAL)
+  lon = real(phi_degrees,kind=CUSTOM_REAL)
   ! limits value range
   if (lat < -90.0) lat = -90.0
   if (lat > 90.0) lat = 90.0
@@ -1226,8 +1226,8 @@
     xi_width = 90.0
     eta_width = 90.0
   else
-    xi_width = ANGULAR_WIDTH_XI_IN_DEGREES
-    eta_width = ANGULAR_WIDTH_ETA_IN_DEGREES
+    xi_width = real(ANGULAR_WIDTH_XI_IN_DEGREES,kind=CUSTOM_REAL)
+    eta_width = real(ANGULAR_WIDTH_ETA_IN_DEGREES,kind=CUSTOM_REAL)
   endif
 
   ! converts chunk width to radians
@@ -1306,8 +1306,8 @@
   integer :: k
 
   ! converts lat/lon from degrees to radians (colatitute/longitude)
-  theta = (90.0 - lat)/180.0 * PI
-  phi = lon/180.0 * PI
+  theta = real((90.d0 - lat)/180.d0 * PI,kind=CUSTOM_REAL)
+  phi = real(lon/180.d0 * PI,kind=CUSTOM_REAL)
 
   ! converts (r,theta,phi) to (x,y,z) on unit sphere
   r = 1.0_CUSTOM_REAL
@@ -1317,9 +1317,9 @@
     ! compute rotation matrix from Euler angles
     call euler_angles(rotation_matrix,CENTER_LONGITUDE_IN_DEGREES,CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH)
 
-    xn = x * rotation_matrix(1,1) + y * rotation_matrix(2,1) + z * rotation_matrix(3,1)
-    yn = x * rotation_matrix(1,2) + y * rotation_matrix(2,2) + z * rotation_matrix(3,2)
-    zn = x * rotation_matrix(1,3) + y * rotation_matrix(2,3) + z * rotation_matrix(3,3)
+    xn = real(x * rotation_matrix(1,1) + y * rotation_matrix(2,1) + z * rotation_matrix(3,1),kind=CUSTOM_REAL)
+    yn = real(x * rotation_matrix(1,2) + y * rotation_matrix(2,2) + z * rotation_matrix(3,2),kind=CUSTOM_REAL)
+    zn = real(x * rotation_matrix(1,3) + y * rotation_matrix(2,3) + z * rotation_matrix(3,3),kind=CUSTOM_REAL)
 
     x = xn; y = yn; z = zn
   endif
@@ -1513,8 +1513,8 @@
   integer :: proc_xi,proc_eta
 
   ! converts chunk width to radians
-  xi_width_rad = xi_width/180.0 * PI
-  eta_width_rad = eta_width/180.0 * PI
+  xi_width_rad = real(xi_width/180.d0 * PI,kind=CUSTOM_REAL)
+  eta_width_rad = real(eta_width/180.d0 * PI,kind=CUSTOM_REAL)
 
   ! gets process/slice number
   xi1 = xi / xi_width_rad * 2.0