Add Modular fp

.github/workflows/ci.yml

@@ -86,17 +86,29 @@ jobs:
           ${{ matrix.intel-command }}
           /bin/bash mfc.sh build -j $(nproc) ${{ matrix.debug }}
 
-      - name: Test Suite (Debug)
+      - name: Test Suite (Debug, Double Precision)
         if: matrix.debug == '--debug'
         run: |
           ${{ matrix.intel-command }}
           /bin/bash mfc.sh test  -j $(nproc) --debug
 
-      - name: Test Suite (No Debug)
+      - name: Test Suite (Debug, Single Precision)
+        if: matrix.debug == '--debug'
+        run: |
+          ${{ matrix.intel-command }}
+          /bin/bash mfc.sh test  -j $(nproc) --debug --single
+
+      - name: Test Suite (No Debug, Double Precision)
         if: matrix.debug == '--no-debug'
         run: |
           ${{ matrix.intel-command }}
           /bin/bash mfc.sh test -j $(nproc) -a
+
+      - name: Test Suite (No Debug, Single Precision)
+        if: matrix.debug == '--no-debug'
+        run: |
+          ${{ matrix.intel-command }}
+          /bin/bash mfc.sh test -j $(nproc) -a --single  
 
 
   self-cpu-release:

CMakeLists.txt

@@ -13,6 +13,7 @@ option(MFC_SIMULATION    "Build simulation"    OFF)
 option(MFC_POST_PROCESS  "Build post_process"  OFF)
 option(MFC_DOCUMENTATION "Build documentation" OFF)
 option(MFC_ALL           "Build everything"    OFF)
+option(MFC_SINGLE_PRECISION "Build double precision" OFF)
 
 if (MFC_BUILD_ALL)
     set(MFC_PRE_PROCESS   ON FORCE)
@@ -100,7 +101,6 @@ elseif (CMAKE_Fortran_COMPILER_ID STREQUAL "Flang")
     add_compile_options(
         $<$<COMPILE_LANGUAGE:Fortran>:-Mfreeform>
         $<$<COMPILE_LANGUAGE:Fortran>:-Mpreprocess>
-        $<$<COMPILE_LANGUAGE:Fortran>:-fdefault-real-8>
     )
 elseif (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
     add_compile_options($<$<COMPILE_LANGUAGE:Fortran>:-free>)
@@ -157,6 +157,12 @@ if (CMAKE_BUILD_TYPE STREQUAL "Debug")
     add_compile_definitions(MFC_DEBUG)
 endif()
 
+if (MFC_SINGLE_PRECISION)
+    add_compile_definitions(MFC_SINGLE_PRECISION)
+else()
+    add_compile_definitions(MFC_DOUBLE_PRECISION)
+endif()
+
 ## === HANDLE_SOURCES
 # Gather F90 source files for a given target, including common/, and preprocessing .fpp -> .f90.
 # Outputs:

docs/documentation/case.md

@@ -76,8 +76,11 @@ Definition of the parameters is described in the following subsections.
 | Parameter        | Type           | Description                      |
 | ---:             |    :----:      |          :---                    |
 | `run_time_info`  | Logical        | Output run-time information      |
+| `working_precision` | Integer     | Simulation working precision     |
 
 - `run_time_info` generates a text file that includes run-time information including the CFL number(s) at each time-step.
+- `working_precision` sets the working precision to either (1) - single precision, or
+(2) - double precision.
 
 ### 2. Computational Domain
 

docs/documentation/running.md

@@ -164,7 +164,13 @@ $ ./mfc.sh run examples/1D_vacuum_restart/case.py -t post_process
 $ ./mfc.sh run examples/1D_vacuum_restart/restart_case.py -t post_process
 ```
 
-### Example Runs
+### Changing Working Precision
+
+MFC has the capability to run in both single and double precision. To enable single
+precision, run `./mfc.sh run ... --single`. To disable single precision run
+`./mfc.sh run ... --no-single`.
+
+## Example Runs
 
 - Oak Ridge National Laboratory's [Summit](https://www.olcf.ornl.gov/summit/):
 

misc/run-phoenix-release-cpu.sh

@@ -11,4 +11,5 @@
 cd $SLURM_SUBMIT_DIR                            # Change to working directory
 echo $(pwd)
 . ./mfc.sh load -c p -m g
-./mfc.sh test -j 12 -b mpirun -a
+./mfc.sh test -j 12 -b mpirun -a --no-single
+./mfc.sh test -j 12 -b mpirun -a --single

misc/run-phoenix-release-gpu.sh

@@ -11,4 +11,5 @@
 cd $SLURM_SUBMIT_DIR                            # Change to working directory
 echo $(pwd)
 . ./mfc.sh load -c p -m g
-./mfc.sh test -j 1 -b mpirun -a --gpu
+./mfc.sh test -j 1 -b mpirun -a --gpu --no-single
+./mfc.sh test -j 1 -b mpirun -a --gpu --single

src/common/include/inline_conversions.fpp

@@ -1,50 +1,50 @@
 #:def s_compute_speed_of_sound()
     subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c)
 
-        real(kind(0d0)), intent(IN) :: pres
-        real(kind(0d0)), intent(IN) :: rho, gamma, pi_inf
-        real(kind(0d0)), intent(IN) :: H
-        real(kind(0d0)), dimension(num_fluids), intent(IN) :: adv
-        real(kind(0d0)), intent(IN) :: vel_sum
-        real(kind(0d0)), intent(OUT) :: c
+        real(wp), intent(IN) :: pres
+        real(wp), intent(IN) :: rho, gamma, pi_inf
+        real(wp), intent(IN) :: H
+        real(wp), dimension(num_fluids), intent(IN) :: adv
+        real(wp), intent(IN) :: vel_sum
+        real(wp), intent(OUT) :: c
 
-        real(kind(0d0)) :: blkmod1, blkmod2
+        real(wp) :: blkmod1, blkmod2
 
         integer :: q
 
         if (alt_soundspeed) then 
-            blkmod1 = ((gammas(1) + 1d0)*pres + & 
+            blkmod1 = ((gammas(1) + 1._wp)*pres + & 
                         pi_infs(1))/gammas(1)
-            blkmod2 = ((gammas(2) + 1d0)*pres + & 
+            blkmod2 = ((gammas(2) + 1._wp)*pres + & 
                         pi_infs(2))/gammas(2) 
-            c = (1d0/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2))) 
+            c = (1._wp/(rho*(adv(1)/blkmod1 + adv(2)/blkmod2))) 
         elseif (model_eqns == 3) then 
-            c = 0d0 
+            c = 0._wp 
 !$acc loop seq 
             do q = 1, num_fluids 
-                c = c + adv(q)*(1d0/gammas(q) + 1d0)* & 
-                    (pres + pi_infs(q)/(gammas(q) + 1d0)) 
+                c = c + adv(q)*(1._wp/gammas(q) + 1._wp)* & 
+                    (pres + pi_infs(q)/(gammas(q) + 1._wp)) 
             end do 
             c = c/rho
 
         elseif (((model_eqns == 4) .or. (model_eqns == 2 .and. bubbles))) then
             ! Sound speed for bubble mmixture to order O(\alpha)
 
             if (mpp_lim .and. (num_fluids > 1)) then
-                c = (1d0/gamma + 1d0)* &
-                      (pres + pi_inf/(gamma+1d0))/rho
+                c = (1._wp/gamma + 1._wp)* &
+                      (pres + pi_inf/(gamma+1._wp))/rho
             else
                 c = &
-                    (1d0/gamma + 1d0)* &    
-                    (pres + pi_inf/(gamma+1d0))/ &
-                    (rho*(1d0 - adv(num_fluids)))
+                    (1._wp/gamma + 1._wp)* &    
+                    (pres + pi_inf/(gamma + 1._wp))/ &
+                    (rho*(1._wp - adv(num_fluids)))
             end if
         else 
-            c = ((H - 5d-1*vel_sum)/gamma) 
+            c = ((H - (5._wp * (10._wp ** -(1)))*vel_sum)/gamma) 
         end if 
 
-        if (mixture_err .and. c < 0d0) then
-            c = 100.d0*sgm_eps
+        if (mixture_err .and. c < 0._wp) then
+            c = 100._wp*sgm_eps
         else
             c = sqrt(c)
         end if

src/common/m_constants.fpp

@@ -4,13 +4,28 @@
 
 module m_constants
 
+    use mpi
+
+    integer, parameter :: sp = selected_real_kind(6, 37)
+    integer, parameter :: dp = selected_real_kind(15, 307)
+
+#ifdef MFC_DOUBLE_PRECISION
+    integer, parameter :: wp = dp
+    integer, parameter :: mpi_p = MPI_DOUBLE_PRECISION
+#endif
+
+#ifdef MFC_SINGLE_PRECISION
+    integer, parameter :: wp = sp
+    integer, parameter :: mpi_p = MPI_REAL
+#endif
+
     character, parameter :: dflt_char = ' ' !< Default string value
 
-    real(kind(0d0)), parameter :: dflt_real = -1d6                !< Default real value
-    real(kind(0d0)), parameter :: sgm_eps   = 1d-16               !< Segmentation tolerance
-    real(kind(0d0)), parameter :: small_alf = 1d-7                !< Small alf tolerance
-    real(kind(0d0)), parameter :: pi        = 3.141592653589793d0 !< Pi
-    real(kind(0d0)), parameter :: verysmall = 1.d-12              !< Very small number
+    real(wp), parameter :: dflt_real = (-1._wp * (10._wp ** 6))                !< Default real value
+    real(wp), parameter :: sgm_eps   = (1._wp * (10._wp ** -(16)))               !< Segmentation tolerance
+    real(wp), parameter :: small_alf = (1._wp * (10._wp ** -(7)))                !< Small alf tolerance
+    real(wp), parameter :: pi        = 3.141592653589793_wp !< Pi
+    real(wp), parameter :: verysmall = (1._wp * (10._wp ** -(12)))              !< Very small number
 
     integer, parameter :: num_stcls_min   = 5    !< Mininum # of stencils
     integer, parameter :: path_len        = 400  !< Maximum path length

src/common/m_derived_types.f90

@@ -8,16 +8,17 @@ module m_derived_types
 
     use m_constants !< Constants
 
+
     implicit none
 
     !> Derived type adding the field position (fp) as an attribute
     type field_position
-        real(kind(0d0)), allocatable, dimension(:, :, :) :: fp !< Field position
+        real(wp), allocatable, dimension(:, :, :) :: fp !< Field position
     end type field_position
 
     !> Derived type annexing a scalar field (SF)
     type scalar_field
-        real(kind(0d0)), pointer, dimension(:, :, :) :: sf => null()
+        real(wp), pointer, dimension(:, :, :) :: sf => null()
     end type scalar_field
 
     type mpi_io_var
@@ -38,8 +39,8 @@ module m_derived_types
 
     !> Derived type adding beginning (beg) and end bounds info as attributes
     type bounds_info
-        real(kind(0d0)) :: beg
-        real(kind(0d0)) :: end
+        real(wp) :: beg
+        real(wp) :: end
     end type bounds_info
 
     !> bounds for the bubble dynamic variables
@@ -62,23 +63,23 @@ module m_derived_types
 
         integer :: geometry !< Type of geometry for the patch
 
-        real(kind(0d0)) :: x_centroid, y_centroid, z_centroid !<
+        real(wp) :: x_centroid, y_centroid, z_centroid !<
         !! Location of the geometric center, i.e. the centroid, of the patch. It
         !! is specified through its x-, y- and z-coordinates, respectively.
 
-        real(kind(0d0)) :: length_x, length_y, length_z !< Dimensions of the patch. x,y,z Lengths.
-        real(kind(0d0)) :: radius !< Dimensions of the patch. radius.
+        real(wp) :: length_x, length_y, length_z !< Dimensions of the patch. x,y,z Lengths.
+        real(wp) :: radius !< Dimensions of the patch. radius.
 
-        real(kind(0d0)), dimension(3) :: radii !<
+        real(wp), dimension(3) :: radii !<
         !! Vector indicating the various radii for the elliptical and ellipsoidal
         !! patch geometries. It is specified through its x-, y-, and z-components
         !! respectively.
 
-        real(kind(0d0)) :: epsilon, beta !<
+        real(wp) :: epsilon, beta !<
         !! The isentropic vortex parameters administrating, respectively, both
         !! the amplitude of the disturbance as well as its domain of influence.
 
-        real(kind(0d0)), dimension(3) :: normal !<
+        real(wp), dimension(3) :: normal !<
         !! Normal vector indicating the orientation of the patch. It is specified
         !! through its x-, y- and z-components, respectively.
         logical, dimension(0:num_patches_max - 1) :: alter_patch !<
@@ -94,79 +95,79 @@ module m_derived_types
         integer :: smooth_patch_id !<
         !! Identity (id) of the patch with which current patch is to get smoothed
 
-        real(kind(0d0)) :: smooth_coeff !<
+        real(wp) :: smooth_coeff !<
         !! Smoothing coefficient (coeff) adminstrating the size of the stencil of
         !! cells across which boundaries of the current patch will be smeared out
 
-        real(kind(0d0)), dimension(num_fluids_max) :: alpha_rho
-        real(kind(0d0)) :: rho
-        real(kind(0d0)), dimension(3) :: vel
-        real(kind(0d0)) :: pres
-        real(kind(0d0)), dimension(num_fluids_max) :: alpha
-        real(kind(0d0)) :: gamma
-        real(kind(0d0)) :: pi_inf !<
+        real(wp), dimension(num_fluids_max) :: alpha_rho
+        real(wp) :: rho
+        real(wp), dimension(3) :: vel
+        real(wp) :: pres
+        real(wp), dimension(num_fluids_max) :: alpha
+        real(wp) :: gamma
+        real(wp) :: pi_inf !<
 
 
         !! Primitive variables associated with the patch. In order, these include
         !! the partial densities, density, velocity, pressure, volume fractions,
         !! specific heat ratio function and the liquid stiffness function.
 
-        real(kind(0d0)), dimension(6) :: tau_e
+        real(wp), dimension(6) :: tau_e
         !! Elastic stresses added to primitive variables if hypoelasticity = True
 
-        real(kind(0d0)) :: R0 !< Bubble size
-        real(kind(0d0)) :: V0 !< Bubble velocity
+        real(wp) :: R0 !< Bubble size
+        real(wp) :: V0 !< Bubble velocity
 
-        real(kind(0d0)) :: p0 !< Bubble size
-        real(kind(0d0)) :: m0 !< Bubble velocity
+        real(wp) :: p0 !< Bubble size
+        real(wp) :: m0 !< Bubble velocity
 
     end type ic_patch_parameters
 
     !> Derived type annexing the physical parameters (PP) of the fluids. These
     !! include the specific heat ratio function and liquid stiffness function.
     type physical_parameters
-        real(kind(0d0)) :: gamma   !< Sp. heat ratio
-        real(kind(0d0)) :: pi_inf  !< Liquid stiffness
-        real(kind(0d0)), dimension(2) :: Re      !< Reynolds number
-        real(kind(0d0)) :: mul0    !< Bubble viscosity
-        real(kind(0d0)) :: ss      !< Bubble surface tension
-        real(kind(0d0)) :: pv      !< Bubble vapour pressure
-        real(kind(0d0)) :: gamma_v !< Bubble constants (see Preston (2007), Ando (2010))
-        real(kind(0d0)) :: M_v     !< Bubble constants (see Preston (2007), Ando (2010))
-        real(kind(0d0)) :: mu_v    !< Bubble constants (see Preston (2007), Ando (2010))
-        real(kind(0d0)) :: k_v     !< Bubble constants (see Preston (2007), Ando (2010))
-        real(kind(0d0)) :: G
+        real(wp) :: gamma   !< Sp. heat ratio
+        real(wp) :: pi_inf  !< Liquid stiffness
+        real(wp), dimension(2) :: Re      !< Reynolds number
+        real(wp) :: mul0    !< Bubble viscosity
+        real(wp) :: ss      !< Bubble surface tension
+        real(wp) :: pv      !< Bubble vapour pressure
+        real(wp) :: gamma_v !< Bubble constants (see Preston (2007), Ando (2010))
+        real(wp) :: M_v     !< Bubble constants (see Preston (2007), Ando (2010))
+        real(wp) :: mu_v    !< Bubble constants (see Preston (2007), Ando (2010))
+        real(wp) :: k_v     !< Bubble constants (see Preston (2007), Ando (2010))
+        real(wp) :: G
     end type physical_parameters
 
     !> Derived type annexing the flow probe location
     type probe_parameters
-        real(kind(0d0)) :: x !< First coordinate location
-        real(kind(0d0)) :: y !< Second coordinate location
-        real(kind(0d0)) :: z !< Third coordinate location
+        real(wp) :: x !< First coordinate location
+        real(wp) :: y !< Second coordinate location
+        real(wp) :: z !< Third coordinate location
     end type probe_parameters
 
     !> Derived type annexing integral regions
     type integral_parameters
-        real(kind(0d0)) :: xmin !< Min. boundary first coordinate direction
-        real(kind(0d0)) :: xmax !< Max. boundary first coordinate direction
-        real(kind(0d0)) :: ymin !< Min. boundary second coordinate direction
-        real(kind(0d0)) :: ymax !< Max. boundary second coordinate direction
-        real(kind(0d0)) :: zmin !< Min. boundary third coordinate direction
-        real(kind(0d0)) :: zmax !< Max. boundary third coordinate direction
+        real(wp) :: xmin !< Min. boundary first coordinate direction
+        real(wp) :: xmax !< Max. boundary first coordinate direction
+        real(wp) :: ymin !< Min. boundary second coordinate direction
+        real(wp) :: ymax !< Max. boundary second coordinate direction
+        real(wp) :: zmin !< Min. boundary third coordinate direction
+        real(wp) :: zmax !< Max. boundary third coordinate direction
     end type integral_parameters
 
     !> Monopole acoustic source parameters
     type mono_parameters
-        real(kind(0d0)), dimension(3) :: loc !< Physical location of acoustic source
-        real(kind(0d0)) :: mag !< Magnitude
-        real(kind(0d0)) :: length !< Length of line source
-        real(kind(0d0)) :: npulse !< Number of cycles of pulse
-        real(kind(0d0)) :: dir !< Direction of pulse
-        real(kind(0d0)) :: delay !< Time-delay of pulse start
+        real(wp), dimension(3) :: loc !< Physical location of acoustic source
+        real(wp) :: mag !< Magnitude
+        real(wp) :: length !< Length of line source
+        real(wp) :: npulse !< Number of cycles of pulse
+        real(wp) :: dir !< Direction of pulse
+        real(wp) :: delay !< Time-delay of pulse start
         integer :: pulse
         integer :: support
-        real(kind(0d0)) :: aperture
-        real(kind(0d0)) :: foc_length
+        real(wp) :: aperture
+        real(wp) :: foc_length
     end type mono_parameters
 
 end module m_derived_types