Fisrt implementation for the 1D axisymmetric spherical coordinates

JRChreim · JRChreim · commit 143f6d30ecb7 · 2025-01-12T01:20:26.000-08:00
diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp
@@ -187,6 +187,7 @@ contains
         @:PROHIBIT(n < 0)
         @:PROHIBIT(p < 0)
         @:PROHIBIT(cyl_coord .and. p > 0 .and. mod(p, 2) /= 1, "p must be odd for cylindrical coordinates")
+        @:PROHIBIT(sph_coord .and. m /= 1 .and. p > 0, "m must not be different than 1 and p not greater than 0 for axysimmetric spherical coordinates")
         @:PROHIBIT(n == 0 .and. p > 0, "p must be 0 if n = 0")
     end subroutine s_check_inputs_simulation_domain
 
@@ -250,10 +251,14 @@ contains
 
         ! Check for y and z directions for cylindrical coordinates
         @: PROHIBIT(cyl_coord .and. n == 0, "n must be positive (2D or 3D) for cylindrical coordinates")
+        @: PROHIBIT(sph_coord .and. m /= 1, "m must be one (1D) for axysimmetric spherical coordinates")
         @: PROHIBIT(cyl_coord .and. p == 0 .and. bc_y%beg /= -2, "bc_y%beg must be -2 for 2D cylindrical coordinates (p = 0)")
+        @: PROHIBIT(sph_coord .and. p == 0 .and. bc_y%beg /= -2, "bc_y%beg must be -2 axysimmetric spherical coordinates (p = 0)")
         @: PROHIBIT(cyl_coord .and. p > 0 .and. bc_y%beg /= -14, "bc_y%beg must be -14 for 3D cylindrical coordinates (p > 0)")
         @: PROHIBIT(cyl_coord .and. (bc_y%end > -1 .or. bc_y%end < -16), "bc_y%end must be between -1 and -16")
+        @: PROHIBIT(sph_coord .and. (bc_y%end > -1 .or. bc_y%end < -16), "bc_y%end must be between -1 and -16")
         @: PROHIBIT(cyl_coord .and. bc_y%end == -14, "bc_y%end must not be -14")
+        @: PROHIBIT(sph_coord .and. (bc_y%end == -2 .or. bc_y%end == -14), "bc_y%end must not be -2 or -14")        
 
         ! Check for y and z directions for 3D cylindrical coordinates
         @: PROHIBIT(cyl_coord .and. p > 0 .and. (bc_z%beg /= -1 .and. bc_z%beg /= -2), &
diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp
@@ -43,6 +43,7 @@ module m_global_parameters
     !> @name Cylindrical coordinates (either axisymmetric or full 3D)
     !> @{
     logical :: cyl_coord
+    logical :: sph_coord
     integer :: grid_geometry
     !> @}
 
@@ -318,6 +319,7 @@ contains
         m = dflt_int; n = 0; p = 0
         m_root = dflt_int
         cyl_coord = .false.
+        sph_coord = .false.
 
         t_step_start = dflt_int
         t_step_stop = dflt_int
@@ -834,7 +836,7 @@ contains
 
         if (cyl_coord .neqv. .true.) then ! Cartesian grid
             grid_geometry = 1
-        elseif (cyl_coord .and. p == 0) then ! Axisymmetric cylindrical grid
+        elseif ((cyl_coord .or. sph_coord) .and. p == 0) then ! Axisymmetric cylindrical grid
             grid_geometry = 2
         else ! Fully 3D cylindrical grid
             grid_geometry = 3
diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp
@@ -168,7 +168,8 @@ contains
             & 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt', 'bubbles_euler', 'qbmm',   &
             & 'polytropic', 'polydisperse', 'file_per_process', 'relax', 'cf_wrt',     &
             & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt',          &
-            & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange', 'rkck_adap_dt', 'output_partial_domain']
+            & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange', 'rkck_adap_dt', 'output_partial_domain', &
+            & 'sph_coord']
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90
@@ -84,7 +84,7 @@ subroutine s_read_input_file
             relax_model, cf_wrt, sigma, adv_n, ib, num_ibs, &
             cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, &
             cfl_target, surface_tension, bubbles_lagrange, rkck_adap_dt, &
-            sim_data, hyperelasticity
+            sim_data, hyperelasticity, sph_coord
 
         ! Inquiring the status of the post_process.inp file
         file_loc = 'post_process.inp'
diff --git a/src/pre_process/m_check_patches.fpp b/src/pre_process/m_check_patches.fpp
@@ -158,6 +158,7 @@ contains
         @:PROHIBIT(patch_icpp(patch_id)%length_x <= 0._wp, "Line segment patch "//trim(iStr)//": length_x must be greater than zero")
         @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), "Line segment patch "//trim(iStr)//": x_centroid must be set")
         @:PROHIBIT(cyl_coord, "Line segment patch "//trim(iStr)//": cyl_coord is not supported")
+        @:PROHIBIT(sph_coord, "Line segment patch "//trim(iStr)//": sph_coord is not supported")
 
     end subroutine s_check_line_segment_patch_geometry
 
diff --git a/src/pre_process/m_checker.fpp b/src/pre_process/m_checker.fpp
@@ -94,6 +94,17 @@ contains
         @:PROHIBIT(cyl_coord .and. p > 0 .and. (z_domain%beg /= 0._wp .or. z_domain%end /= 2._wp*pi), &
             "z_domain%beg must be 0 and z_domain%end must be 2*pi for 3D cylindrical coordinates")
 
+        @:PROHIBIT(sph_coord .and. n == 0, &
+            "n must be positive for 1D axisymmetric spherical coordinates")
+        @:PROHIBIT(sph_coord .and. (f_is_default(y_domain%beg) .or. f_is_default(y_domain%end)), &
+            "y_domain%beg and y_domain%end must be set for n > 0 (1D axisymmetric spherical coordinates)")
+        @:PROHIBIT(sph_coord .and. (y_domain%beg /= 0._wp .or. y_domain%end <= 0._wp), &
+            "y_domain%beg must be 0 and y_domain%end must be positive for 1D axisymmetric spherical coordinates")
+        @:PROHIBIT(sph_coord .and. p == 0 .and. ((.not. f_is_default(z_domain%beg)) .or. (.not. f_is_default(z_domain%end))), &
+            "z_domain%beg and z_domain%end are not supported for p = 0 (1D spherical coordinates)")
+        @:PROHIBIT(sph_coord .and. p > 0 &
+            "only 1D axisymmetric (m = 1, n > 0, p = 0) is supported with spherical coordinates")
+                    
         @:PROHIBIT(num_patches < 0)
         @:PROHIBIT(num_patches == 0 .and. t_step_old == dflt_int, &
             "num_patches must be positive for the non-restart case")
diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp
@@ -47,6 +47,7 @@ module m_global_parameters
     integer :: num_dims !< Number of spatial dimensions
 
     logical :: cyl_coord
+    logical :: sph_coord    
     integer :: grid_geometry !< Cylindrical coordinates (either axisymmetric or full 3D)
 
     real(wp), allocatable, dimension(:) :: x_cc, y_cc, z_cc !<
@@ -275,6 +276,7 @@ contains
         m = dflt_int; n = 0; p = 0
 
         cyl_coord = .false.
+        sph_coord = .false.
 
         x_domain%beg = dflt_real
         x_domain%end = dflt_real
@@ -828,7 +830,7 @@ contains
 
         if (cyl_coord .neqv. .true.) then ! Cartesian grid
             grid_geometry = 1
-        elseif (cyl_coord .and. p == 0) then ! Axisymmetric cylindrical grid
+        elseif ((cyl_coord .or. sph_coord) .and. p == 0) then ! Axisymmetric cylindrical grid
             grid_geometry = 2
         else ! Fully 3D cylindrical grid
             grid_geometry = 3
diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp
@@ -55,7 +55,7 @@ contains
             & 'mixlayer_perturb', 'bubbles_euler', 'polytropic', 'polydisperse',&
             & 'qbmm', 'file_per_process', 'adv_n', 'ib' , 'cfl_adap_dt',       &
             & 'cfl_const_dt', 'cfl_dt', 'surface_tension',                     &
-            & 'hyperelasticity', 'pre_stress' ]
+            & 'hyperelasticity', 'pre_stress', 'sph_coord']
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
         call MPI_BCAST(fluid_rho(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp
@@ -140,7 +140,8 @@ contains
             file_per_process, relax, relax_model, &
             palpha_eps, ptgalpha_eps, ib, num_ibs, patch_ib, &
             sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, &
-            n_start_old, surface_tension, hyperelasticity, pre_stress, rkck_adap_dt
+            n_start_old, surface_tension, hyperelasticity, pre_stress, rkck_adap_dt, &
+            sph_coord
 
         ! Inquiring the status of the pre_process.inp file
         file_loc = 'pre_process.inp'
diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp
@@ -901,7 +901,7 @@ contains
                         end if
 
                         ! Be careful about the cylindrical coordinate!
-                        if (cyl_coord .and. cbc_dir == 2 .and. cbc_loc == 1) then
+                        if ((cyl_coord .or. sph_coord) .and. cbc_dir == 2 .and. cbc_loc == 1) then
                             dpres_dt = -5e-1_wp*(L(advxe) + L(1)) + rho*c*c*vel(dir_idx(1)) &
                                        /y_cc(n)
                         else
@@ -929,7 +929,7 @@ contains
                         end do
 
                         ! The treatment of void fraction source is unclear
-                        if (cyl_coord .and. cbc_dir == 2 .and. cbc_loc == 1) then
+                        if ((cyl_coord .or. sph_coord) .and. cbc_dir == 2 .and. cbc_loc == 1) then
                             !$acc loop seq
                             do i = 1, advxe - E_idx
                                 dadv_dt(i) = -L(momxe + i) !+ adv(i) * vel(dir_idx(1))/y_cc(n)
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
@@ -51,9 +51,10 @@ module m_global_parameters
     !> @name Cylindrical coordinates (either axisymmetric or full 3D)
     !> @{
     logical :: cyl_coord
+    logical :: sph_coord
     integer :: grid_geometry
     !> @}
-    !$acc declare create(cyl_coord, grid_geometry)
+    !$acc declare create(cyl_coord, grid_geometry, sph_coord)
 
     !> @name Cell-boundary (CB) locations in the x-, y- and z-directions, respectively
     !> @{
@@ -488,6 +489,7 @@ contains
         m = dflt_int; n = 0; p = 0
 
         cyl_coord = .false.
+        sph_coord = .false.
 
         dt = dflt_real
 
@@ -1158,7 +1160,7 @@ contains
 
         if (cyl_coord .neqv. .true.) then ! Cartesian grid
             grid_geometry = 1
-        elseif (cyl_coord .and. p == 0) then ! Axisymmetric cylindrical grid
+        elseif ((cyl_coord .or. sph_coord) .and. p == 0) then ! Axisymmetric cylindrical grid
             grid_geometry = 2
         else ! Fully 3D cylindrical grid
             grid_geometry = 3
@@ -1188,7 +1190,7 @@ contains
         !$acc update device(cfl_target, m, n, p)
 
         !$acc update device(alt_soundspeed, acoustic_source, num_source)
-        !$acc update device(dt, sys_size, buff_size, pref, rhoref, gamma_idx, pi_inf_idx, E_idx, alf_idx, stress_idx, mpp_lim, bubbles_euler, hypoelasticity, alt_soundspeed, avg_state, num_fluids, model_eqns, num_dims, mixture_err, grid_geometry, cyl_coord, mp_weno, weno_eps, teno_CT, hyperelasticity, hyper_model, elasticity, xi_idx, low_Mach)
+        !$acc update device(dt, sys_size, buff_size, pref, rhoref, gamma_idx, pi_inf_idx, E_idx, alf_idx, stress_idx, mpp_lim, bubbles_euler, hypoelasticity, alt_soundspeed, avg_state, num_fluids, model_eqns, num_dims, mixture_err, grid_geometry, cyl_coord, mp_weno, weno_eps, teno_CT, hyperelasticity, hyper_model, elasticity, xi_idx, low_Mach, sph_coord)
 
         #:if not MFC_CASE_OPTIMIZATION
             !$acc update device(wenojs, mapped_weno, wenoz, teno)
diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
@@ -195,7 +195,7 @@ contains
             & 'bc_z%grcbc_in', 'bc_z%grcbc_out', 'bc_z%grcbc_vel_out',          &
             & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', 'surface_tension',        &
             & 'viscous', 'shear_stress', 'bulk_stress', 'bubbles_lagrange',     &
-            & 'hyperelasticity', 'rkck_adap_dt' ]
+            & 'hyperelasticity', 'rkck_adap_dt', 'sph_coord' ]
             call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
@@ -1122,6 +1122,22 @@ contains
                             end do
                         end do
                     end do
+                elseif (sph_coord) then
+                    !$acc parallel loop collapse(4) gang vector default(present)
+                    do l = 0, p
+                        do k = 0, n
+                            do j = 0, m
+                                do i = 1, num_fluids
+                                    rhs_vf(i + intxb - 1)%sf(j, k, l) = &
+                                        rhs_vf(i + intxb - 1)%sf(j, k, l) - 1._wp/y_cc(k)* &
+                                        q_cons_vf%vf(i + advxb - 1)%sf(j, k, l)* &
+                                        q_prim_vf%vf(E_idx)%sf(j, k, l)* &
+                                        (flux_src_n(2)%vf(advxb)%sf(j, k, l) + &
+                                         flux_src_n(2)%vf(advxb)%sf(j, k - 1, l))
+                                end do
+                            end do
+                        end do
+                    end do
                 end if
             end if
 
@@ -1139,6 +1155,20 @@ contains
                         end do
                     end do
                 end do
+            elseif (sph_coord) then
+                !$acc parallel loop collapse(4) gang vector default(present)
+                do j = 1, sys_size
+                    do l = 0, p
+                        do k = 0, n
+                            do q = 0, m
+                                rhs_vf(j)%sf(q, k, l) = &
+                                    rhs_vf(j)%sf(q, k, l) - 1._wp/y_cc(k)* &
+                                    (flux_gsrc_n(2)%vf(j)%sf(q, k - 1, l) &
+                                     + flux_gsrc_n(2)%vf(j)%sf(q, k, l))
+                            end do
+                        end do
+                    end do
+                end do
             end if
 
             if (riemann_solver == 1) then
@@ -1186,6 +1216,19 @@ contains
                                         end do
                                     end do
                                 end do
+                            elseif (sph_coord) then
+                                !$acc parallel loop collapse(3) gang vector default(present)
+                                do l = 0, p
+                                    do k = 0, n
+                                        do q = 0, m
+                                            rhs_vf(j)%sf(q, k, l) = &
+                                                rhs_vf(j)%sf(q, k, l) - &
+                                                (Kterm(q, k, l)/1._wp/y_cc(k))* &
+                                                (flux_src_n(2)%vf(j)%sf(q, k, l) &
+                                                 + flux_src_n(2)%vf(j)%sf(q, k - 1, l))
+                                        end do
+                                    end do
+                                end do
                             end if
                         else if ((j == advxb) .and. (bubbles_euler .neqv. .true.)) then
                             !$acc parallel loop collapse(3) gang vector default(present)
@@ -1213,6 +1256,19 @@ contains
                                         end do
                                     end do
                                 end do
+                            elseif (sph_coord) then
+                                !$acc parallel loop collapse(3) gang vector default(present)
+                                do l = 0, p
+                                    do k = 0, n
+                                        do q = 0, m
+                                            rhs_vf(j)%sf(q, k, l) = &
+                                                rhs_vf(j)%sf(q, k, l) + &
+                                                (Kterm(q, k, l)/1._wp/y_cc(k))* &
+                                                (flux_src_n(2)%vf(j)%sf(q, k, l) &
+                                                 + flux_src_n(2)%vf(j)%sf(q, k - 1, l))
+                                        end do
+                                    end do
+                                end do
                             end if
                         end if
                     end do
@@ -1359,6 +1415,19 @@ contains
                                             end do
                                         end do
                                     end do
+                                elseif (sph_coord) then
+                                    !$acc parallel loop collapse(3) gang vector default(present)
+                                    do l = 0, p
+                                        do k = 0, n
+                                            do q = 0, m
+                                                rhs_vf(j)%sf(q, k, l) = &
+                                                    rhs_vf(j)%sf(q, k, l) - &
+                                                    (Kterm(q, k, l)/1._wp/y_cc(k))* &
+                                                    (flux_src_n(2)%vf(j)%sf(q, k, l) &
+                                                     + flux_src_n(2)%vf(j)%sf(q, k - 1, l))
+                                            end do
+                                        end do
+                                    end do
                                 end if
                             else if ((j == advxb) .and. (bubbles_euler .neqv. .true.)) then
                                 !$acc parallel loop collapse(3) gang vector default(present)
@@ -1386,6 +1455,19 @@ contains
                                             end do
                                         end do
                                     end do
+                                elseif (sph_coord) then
+                                    !$acc parallel loop collapse(3) gang vector default(present)
+                                    do l = 0, p
+                                        do k = 0, n
+                                            do q = 0, m
+                                                rhs_vf(j)%sf(q, k, l) = &
+                                                    rhs_vf(j)%sf(q, k, l) + &
+                                                    (Kterm(q, k, l)/1._wp/y_cc(k))* &
+                                                    (flux_src_n(2)%vf(j)%sf(q, k, l) &
+                                                     + flux_src_n(2)%vf(j)%sf(q, k - 1, l))
+                                            end do
+                                        end do
+                                    end do
                                 end if
                             end if
                         end do
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
@@ -1430,7 +1430,7 @@ contains
 
                                 ! Geometrical source flux for cylindrical coordinates
                                 #:if (NORM_DIR == 2)
-                                    if (cyl_coord) then
+                                    if (cyl_coord .or. sph_coord) then
                                         !Substituting the advective flux into the inviscid geometrical source flux
                                         !$acc loop seq
                                         do i = 1, E_idx
@@ -1670,7 +1670,7 @@ contains
                                 ! Geometrical source flux for cylindrical coordinates
 
                                 #:if (NORM_DIR == 2)
-                                    if (cyl_coord) then
+                                    if (cyl_coord .or. sph_coord) then
                                         ! Substituting the advective flux into the inviscid geometrical source flux
                                         !$acc loop seq
                                         do i = 1, E_idx
@@ -2134,7 +2134,7 @@ contains
 
                                 ! Geometrical source flux for cylindrical coordinates
                                 #:if (NORM_DIR == 2)
-                                    if (cyl_coord) then
+                                    if (cyl_coord .or. sph_coord) then
                                         ! Substituting the advective flux into the inviscid geometrical source flux
                                         !$acc loop seq
                                         do i = 1, E_idx
@@ -2633,7 +2633,7 @@ contains
 
                                 ! Geometrical source flux for cylindrical coordinates
                                 #:if (NORM_DIR == 2)
-                                    if (cyl_coord) then
+                                    if (cyl_coord .or. sph_coord) then
                                         !Substituting the advective flux into the inviscid geometrical source flux
                                         !$acc loop seq
                                         do i = 1, E_idx
@@ -4419,7 +4419,7 @@ contains
                 end do
             end do
 
-            if (cyl_coord) then
+            if (cyl_coord .or. sph_coord) then
                 !$acc parallel loop collapse(4) gang vector default(present)
                 do i = 1, sys_size
                     do l = is3%beg, is3%end
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
