October Changes from Julius, mainly time statistics

LIB/EQUATION/ACMnew/module_ACM.f90

@@ -83,7 +83,7 @@ module module_acm
     real(kind=rk) :: mask_volume=0.0_rk, meanflow_channel(1:3) = 0.0_rk
 
     logical :: use_passive_scalar = .false.
-    integer(kind=ik) :: N_scalars = 0, nsave_stats = 999999
+    integer(kind=ik) :: N_scalars = 0
     real(kind=rk), allocatable :: schmidt_numbers(:), x0source(:), y0source(:), &
     z0source(:), scalar_Ceta(:), widthsource(:)
     character(len=cshort), allocatable :: scalar_inicond(:), scalar_source_type(:)
@@ -192,9 +192,9 @@ subroutine READ_PARAMETERS_ACM( filename, N_mask_components, g )
 
     character(len=*), intent(in) :: filename
     integer(kind=ik) :: mpicode, nx_max, n_entries
-    real(kind=rk) :: dx_min, dt_min
+    real(kind=rk) :: dx_min, dt_min_c0, dt_min_vpm, dt_min_nu
     character(len=cshort) :: Bs_str, Bs_conc
-    character(len=16834) :: input_files
+    character(len=400) :: input_files
     character(len=12) :: timestamp
     character(:), allocatable :: Bs_short
     real(kind=rk), dimension(3) :: ddx
@@ -348,7 +348,6 @@ subroutine READ_PARAMETERS_ACM( filename, N_mask_components, g )
 
     call read_param_mpi(FILE, 'Time', 'time_max', params_acm%T_end, 1.0_rk   )
 
-    call read_param_mpi(FILE, 'Statistics', 'nsave_stats', params_acm%nsave_stats, 999999   )
     call read_param_mpi(FILE, 'FreeFlightSolver', 'use_free_flight_solver', params_acm%use_free_flight_solver, .false.   )
 
     call read_param_mpi(FILE, 'Blocks', 'max_treelevel', params_acm%Jmax, 1   )
@@ -428,11 +427,23 @@ subroutine READ_PARAMETERS_ACM( filename, N_mask_components, g )
     ! uniqueGrid modification
     nx_max = maxval( (params_acm%Bs) * 2**(params_acm%Jmax) )
 
+    ! print some time-step related information
     if (params_acm%c_0 > 0.0_rk) then
-        dt_min = params_acm%CFL*dx_min/params_acm%c_0
+        dt_min_c0 = params_acm%CFL*dx_min/params_acm%c_0
     else
-        dt_min = 0.0_rk
+        dt_min_c0 = 0.0_rk
     endif
+    if (params_acm%nu > 0.0_rk) then
+        dt_min_nu = params_acm%CFL_nu*dx_min**2/params_acm%nu
+    else
+        dt_min_nu = 0.0_rk
+    endif
+    if (params_acm%penalization .and. params_acm%C_eta > 0.0_rk) then
+        dt_min_vpm = params_acm%CFL_eta*params_acm%C_eta
+    else
+        dt_min_vpm = 0.0_rk
+    endif
+
     ! nice to have this elsewhere in the ACM module:
     params_acm%dx_min = dx_min
 
@@ -442,15 +453,17 @@ subroutine READ_PARAMETERS_ACM( filename, N_mask_components, g )
 
     if (params_acm%mpirank==0) then
       write(*,'(80("<"))')
-      write(*,*) "Some information:"
-      write(*,'("c0=",g12.4," C_eta=",g12.4," CFL=",g12.4)') params_acm%c_0, params_acm%C_eta, params_acm%CFL
-      write(*,'("dx_min=",g12.4," dt(CFL,c0,dx_min)=",g12.4)') dx_min, dt_min
-      write(*,'("if all blocks were at Jmax, the resolution would be nx=",i5)') nx_max
+      write(*,'(A)') "Some information:"
+      write(*,'("   c0=",g12.4," CFL=",g12.4, "CFL_eta=",g12.4, "CFL_nu=",g12.4)') params_acm%c_0, params_acm%CFL, params_acm%CFL_eta, params_acm%CFL_nu
+      write(*,'("   dx_min=",g12.4)') dx_min
+      write(*,'("   dt(CFL,c0,dx_min)=",g12.4)') dt_min_c0
+      write(*,'("   dt(CFL_nu,nu,dx_min)=",g12.4)') dt_min_nu
+      write(*,'("   dt(CFL_vpm,C_eta,dx_min)=",g12.4)') dt_min_vpm
+      write(*,'("   if all blocks were at Jmax, the resolution would be nx=",i5)') nx_max
       if (params_acm%penalization) then
-          write(*,'("C_eta=",es12.4," K_eta=",es12.4)') params_acm%C_eta, sqrt(params_acm%C_eta*params_acm%nu)/dx_min
+          write(*,'("   C_eta=",es12.4," K_eta=",es12.4)') params_acm%C_eta, sqrt(params_acm%C_eta*params_acm%nu)/dx_min
       endif
-      write(*,'("N_mask_components=",i1)') N_mask_components
-      write(*,'("N_scalars=",i2)') params_acm%N_scalars
+      write(*,'("N_mask_components=",i1, " N_scalars=",i1, " N_time_statistics=",i1)') N_mask_components, params_acm%N_scalars, params_acm%N_time_statistics
       write(*,'(80("<"))')
     endif
 

LIB/EQUATION/ACMnew/rhs_ACM.f90

@@ -817,6 +817,7 @@ end subroutine RHS_2D_acm
 
 
 subroutine RHS_3D_acm(g, Bs, dx, x0, phi, order_discretization, time, rhs, mask, n_domain)
+    use module_operators
     implicit none
 
     !> grid parameter
@@ -875,6 +876,9 @@ subroutine RHS_3D_acm(g, Bs, dx, x0, phi, order_discretization, time, rhs, mask,
     real(kind=rk), parameter :: a_FD6(-3:3) = (/-1.0_rk/60.0_rk, 3.0_rk/20.0_rk, -3.0_rk/4.0_rk, 0.0_rk, 3.0_rk/4.0_rk, -3.0_rk/20.0_rk, 1.0_rk/60.0_rk/) ! 1st derivative
     real(kind=rk), parameter :: b_FD6(-3:3) = (/ 1.0_rk/90.0_rk, -3.0_rk/20.0_rk, 3.0_rk/2.0_rk, -49.0_rk/18.0_rk, 3.0_rk/2.0_rk, -3.0_rk/20.0_rk, 1.0_rk/90.0_rk/) ! 2nd derivative
 
+    real(kind=rk), allocatable, dimension(:) :: FD1_l, FD2
+    integer(kind=ik) :: FD1_ls, FD1_le, FD2_s, FD2_e
+
 
     ! set parameters for readability
     c_0         = params_acm%c_0
@@ -1453,7 +1457,131 @@ subroutine RHS_3D_acm(g, Bs, dx, x0, phi, order_discretization, time, rhs, mask,
         endif
 
     case default
-        call abort(441167, "3d Discretization unkown "//order_discretization//", I ll walk into the light now." )
+        ! call abort(44167, "3d Discretization unknown "// trim(order_discretization)//", I'll walk into the light now.")
+
+        ! Flexible operator version for 3D ACM using generalized stencils
+
+        call setup_FD1_left_stencil(order_discretization, FD1_l, FD1_ls, FD1_le)
+        call setup_FD2_stencil(order_discretization, FD2, FD2_s, FD2_e)
+
+        if (params_acm%skew_symmetry) then
+            do iz = g+1, Bs(3)+g
+                do iy = g+1, Bs(2)+g
+                    do ix = g+1, Bs(1)+g
+                        u = phi(ix,iy,iz,1)
+                        v = phi(ix,iy,iz,2)
+                        w = phi(ix,iy,iz,3)
+                        p = phi(ix,iy,iz,4)
+
+                        chi = mask(ix,iy,iz,1) * C_eta_inv
+                        penalx = -chi * (u - mask(ix,iy,iz,2))
+                        penaly = -chi * (v - mask(ix,iy,iz,3))
+                        penalz = -chi * (w - mask(ix,iy,iz,4))
+
+                        ! Generalized first derivatives
+                        u_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,1)) * dx_inv
+                        u_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,1)) * dy_inv
+                        u_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,1)) * dz_inv
+
+                        v_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,2)) * dx_inv
+                        v_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,2)) * dy_inv
+                        v_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,2)) * dz_inv
+
+                        w_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,3)) * dx_inv
+                        w_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,3)) * dy_inv
+                        w_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,3)) * dz_inv
+
+                        p_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,4)) * dx_inv
+                        p_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,4)) * dy_inv
+                        p_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,4)) * dz_inv
+
+                        ! Generalized nonlinear terms
+                        uu_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,1)*phi(ix+FD1_ls:ix+FD1_le,iy,iz,1)) * dx_inv
+                        uv_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,1)*phi(ix,iy+FD1_ls:iy+FD1_le,iz,2)) * dy_inv
+                        uw_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,1)*phi(ix,iy,iz+FD1_ls:iz+FD1_le,3)) * dz_inv
+
+                        vu_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,2)*phi(ix+FD1_ls:ix+FD1_le,iy,iz,1)) * dx_inv
+                        vv_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,2)*phi(ix,iy+FD1_ls:iy+FD1_le,iz,2)) * dy_inv
+                        vw_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,2)*phi(ix,iy,iz+FD1_ls:iz+FD1_le,3)) * dz_inv
+
+                        wu_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,3)*phi(ix+FD1_ls:ix+FD1_le,iy,iz,1)) * dx_inv
+                        wv_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,3)*phi(ix,iy+FD1_ls:iy+FD1_le,iz,2)) * dy_inv
+                        ww_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,3)*phi(ix,iy,iz+FD1_ls:iz+FD1_le,3)) * dz_inv
+
+                        ! Generalized second derivatives
+                        u_dxdx = sum(FD2(FD2_s:FD2_e) * phi(ix+FD2_s:ix+FD2_e,iy,iz,1)) * dx2_inv
+                        u_dydy = sum(FD2(FD2_s:FD2_e) * phi(ix,iy+FD2_s:iy+FD2_e,iz,1)) * dy2_inv
+                        u_dzdz = sum(FD2(FD2_s:FD2_e) * phi(ix,iy,iz+FD2_s:iz+FD2_e,1)) * dz2_inv
+
+                        v_dxdx = sum(FD2(FD2_s:FD2_e) * phi(ix+FD2_s:ix+FD2_e,iy,iz,2)) * dx2_inv
+                        v_dydy = sum(FD2(FD2_s:FD2_e) * phi(ix,iy+FD2_s:iy+FD2_e,iz,2)) * dy2_inv
+                        v_dzdz = sum(FD2(FD2_s:FD2_e) * phi(ix,iy,iz+FD2_s:iz+FD2_e,2)) * dz2_inv
+
+                        w_dxdx = sum(FD2(FD2_s:FD2_e) * phi(ix+FD2_s:ix+FD2_e,iy,iz,3)) * dx2_inv
+                        w_dydy = sum(FD2(FD2_s:FD2_e) * phi(ix,iy+FD2_s:iy+FD2_e,iz,3)) * dy2_inv
+                        w_dzdz = sum(FD2(FD2_s:FD2_e) * phi(ix,iy,iz+FD2_s:iz+FD2_e,3)) * dz2_inv
+
+                        ! Skew-symmetric formulation
+                        rhs(ix,iy,iz,1) = -0.5_rk*(uu_dx + uv_dy + uw_dz   + u*u_dx + v*u_dy + w*u_dz) -p_dx + nu*(u_dxdx + u_dydy + u_dzdz) + penalx
+                        rhs(ix,iy,iz,2) = -0.5_rk*(vu_dx + vv_dy + vw_dz   + u*v_dx + v*v_dy + w*v_dz) -p_dy + nu*(v_dxdx + v_dydy + v_dzdz) + penaly
+                        rhs(ix,iy,iz,3) = -0.5_rk*(wu_dx + wv_dy + ww_dz   + u*w_dx + v*w_dy + w*w_dz) -p_dz + nu*(w_dxdx + w_dydy + w_dzdz) + penalz
+                        rhs(ix,iy,iz,4) = -(c_0**2)*(u_dx + v_dy + w_dz) - gamma*p
+                    end do
+                end do
+            end do
+        else
+            do iz = g+1, Bs(3)+g
+                do iy = g+1, Bs(2)+g
+                    do ix = g+1, Bs(1)+g
+                        u = phi(ix,iy,iz,1)
+                        v = phi(ix,iy,iz,2)
+                        w = phi(ix,iy,iz,3)
+                        p = phi(ix,iy,iz,4)
+
+                        chi = mask(ix,iy,iz,1) * C_eta_inv
+                        penalx = -chi * (u - mask(ix,iy,iz,2))
+                        penaly = -chi * (v - mask(ix,iy,iz,3))
+                        penalz = -chi * (w - mask(ix,iy,iz,4))
+
+                        ! Generalized first derivatives
+                        u_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,1)) * dx_inv
+                        u_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,1)) * dy_inv
+                        u_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,1)) * dz_inv
+
+                        v_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,2)) * dx_inv
+                        v_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,2)) * dy_inv
+                        v_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,2)) * dz_inv
+
+                        w_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,3)) * dx_inv
+                        w_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,3)) * dy_inv
+                        w_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,3)) * dz_inv
+
+                        p_dx = sum(FD1_l(FD1_ls:FD1_le) * phi(ix+FD1_ls:ix+FD1_le,iy,iz,4)) * dx_inv
+                        p_dy = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy+FD1_ls:iy+FD1_le,iz,4)) * dy_inv
+                        p_dz = sum(FD1_l(FD1_ls:FD1_le) * phi(ix,iy,iz+FD1_ls:iz+FD1_le,4)) * dz_inv
+
+                        ! Generalized second derivatives
+                        u_dxdx = sum(FD2(FD2_s:FD2_e) * phi(ix+FD2_s:ix+FD2_e,iy,iz,1)) * dx2_inv
+                        u_dydy = sum(FD2(FD2_s:FD2_e) * phi(ix,iy+FD2_s:iy+FD2_e,iz,1)) * dy2_inv
+                        u_dzdz = sum(FD2(FD2_s:FD2_e) * phi(ix,iy,iz+FD2_s:iz+FD2_e,1)) * dz2_inv
+
+                        v_dxdx = sum(FD2(FD2_s:FD2_e) * phi(ix+FD2_s:ix+FD2_e,iy,iz,2)) * dx2_inv
+                        v_dydy = sum(FD2(FD2_s:FD2_e) * phi(ix,iy+FD2_s:iy+FD2_e,iz,2)) * dy2_inv
+                        v_dzdz = sum(FD2(FD2_s:FD2_e) * phi(ix,iy,iz+FD2_s:iz+FD2_e,2)) * dz2_inv
+
+                        w_dxdx = sum(FD2(FD2_s:FD2_e) * phi(ix+FD2_s:ix+FD2_e,iy,iz,3)) * dx2_inv
+                        w_dydy = sum(FD2(FD2_s:FD2_e) * phi(ix,iy+FD2_s:iy+FD2_e,iz,3)) * dy2_inv
+                        w_dzdz = sum(FD2(FD2_s:FD2_e) * phi(ix,iy,iz+FD2_s:iz+FD2_e,3)) * dz2_inv
+
+                        ! Standard formulation
+                        rhs(ix,iy,iz,1) = (-u*u_dx - v*u_dy - w*u_dz) -p_dx + nu*(u_dxdx + u_dydy + u_dzdz) + penalx
+                        rhs(ix,iy,iz,2) = (-u*v_dx - v*v_dy - w*v_dz) -p_dy + nu*(v_dxdx + v_dydy + v_dzdz) + penaly
+                        rhs(ix,iy,iz,3) = (-u*w_dx - v*w_dy - w*w_dz) -p_dz + nu*(w_dxdx + w_dydy + w_dzdz) + penalz
+                        rhs(ix,iy,iz,4) = -(c_0**2)*(u_dx + v_dy + w_dz) - gamma*p
+                    end do
+                end do
+            end do
+        endif
 
     end select
 
@@ -1483,7 +1611,6 @@ subroutine RHS_3D_acm(g, Bs, dx, x0, phi, order_discretization, time, rhs, mask,
 
     ! --------------------------------------------------------------------------
     ! HIT linear forcing
-    ! ATTENTION! This is at last position because I modify phi to avoid a 3-nested do-loop to subtract mean-flow
     ! --------------------------------------------------------------------------
     if (params_acm%HIT_linear_forcing) then
         G_gain = params_acm%HIT_gain

LIB/EQUATION/ACMnew/sponge.f90

@@ -78,22 +78,30 @@ subroutine sponge_3D(sponge, x0, dx, Bs, g)
 
 
     if (params_acm%sponge_type == "rect") then
-        ! rectangular sponge with 45deg edges
-        do iz = g+1, Bs(3)+g
-            z = dble(iz-(g+1)) * dx(3) + x0(3)
-            do iy = g+1, Bs(2)+g
-                y = dble(iy-(g+1)) * dx(2) + x0(2)
-                do ix = g+1, Bs(1)+g
-                    x = dble(ix-(g+1)) * dx(1) + x0(1)
-
-                    ! distance to borders of domain
-                    tmp = minval( (/x,y,z,-(x-params_acm%domain_size(1)),&
-                         -(y-params_acm%domain_size(2)),-(z-params_acm%domain_size(3))/) )
 
-                    sponge(ix,iy,iz) = smoothstep( tmp, 0.5_rk*params_acm%L_sponge, 0.5_rk*params_acm%L_sponge)
+        ! check if this block is in sponge layer
+        ! find the point in the block closest to the domain boundary
+        tmp = min(minval(x0), minval(params_acm%domain_size - (x0 + dx*dble(Bs))))
+        if (tmp > params_acm%L_sponge) then
+            sponge = 0.0_rk
+        else
+            ! rectangular sponge with 45deg edges
+            do iz = g+1, Bs(3)+g
+                z = dble(iz-(g+1)) * dx(3) + x0(3)
+                do iy = g+1, Bs(2)+g
+                    y = dble(iy-(g+1)) * dx(2) + x0(2)
+                    do ix = g+1, Bs(1)+g
+                        x = dble(ix-(g+1)) * dx(1) + x0(1)
+
+                        ! distance to borders of domain
+                        tmp = minval( (/x,y,z,-(x-params_acm%domain_size(1)),&
+                            -(y-params_acm%domain_size(2)),-(z-params_acm%domain_size(3))/) )
+
+                        sponge(ix,iy,iz) = smoothstep( tmp, 0.5_rk*params_acm%L_sponge, 0.5_rk*params_acm%L_sponge)
+                    end do
                 end do
             end do
-        end do
+        endif
 
         ! sponge for using with symmetry_BC
         ! insect is supposed to be at y=0
@@ -149,6 +157,7 @@ subroutine sponge_3D(sponge, x0, dx, Bs, g)
         ! p-norm sponge. The shape of the sponge is dictated as the p-norm
         ! https://de.wikipedia.org/wiki/P-Norm
         ! which is a nice and simple way to get a rectangle with round corners.
+        ! For 2 it is the eucledian norm (circle), for infinity it is a rectangle and for everything in between a rounded rectangle
 
         if ( maxval(abs(params_acm%domain_size-params_acm%domain_size(1))) > 1.0e-10_rk) then
             call abort(1610184,"ERROR: for the p-norm sponge, the domain has to be same size in all directions.")
@@ -158,21 +167,31 @@ subroutine sponge_3D(sponge, x0, dx, Bs, g)
         pinv = 1.0_rk / p
         offset = 0.5_rk * params_acm%domain_size(1)
 
-        do iz = g+1, Bs(3)+g
-            z = (dble(iz-(g+1)) * dx(3) + x0(3) - offset)**p
-            do iy = g+1, Bs(2)+g
-                y = (dble(iy-(g+1)) * dx(2) + x0(2) - offset)**p
-                do ix = g+1, Bs(1)+g
-                    x = (dble(ix-(g+1)) * dx(1) + x0(1) - offset)**p
-
-                    ! distance to borders of domain
-                    tmp = -( (x + y + z)**pinv - offset)
-
-                    sponge(ix,iy,iz) = smoothstep( tmp, 0.5_rk*params_acm%L_sponge, &
-                    0.5_rk*params_acm%L_sponge)
+        ! check if this block is in sponge layer
+        ! Find the point in the block furthest from the domain center
+        x = max(abs(x0(1)-offset), abs(x0(1)+dx(1)*dble(Bs(1))-offset))**p
+        y = max(abs(x0(2)-offset), abs(x0(2)+dx(2)*dble(Bs(2))-offset))**p
+        z = max(abs(x0(3)-offset), abs(x0(3)+dx(3)*dble(Bs(3))-offset))**p
+        tmp = offset - (x + y + z)**pinv  ! we invert distance to center to distance from boundary
+        if (tmp > params_acm%L_sponge) then
+            sponge = 0.0_rk
+        else
+            do iz = g+1, Bs(3)+g
+                z = (dble(iz-(g+1)) * dx(3) + x0(3) - offset)**p
+                do iy = g+1, Bs(2)+g
+                    y = (dble(iy-(g+1)) * dx(2) + x0(2) - offset)**p
+                    do ix = g+1, Bs(1)+g
+                        x = (dble(ix-(g+1)) * dx(1) + x0(1) - offset)**p
+
+                        ! distance to borders of domain
+                        tmp = -( (x + y + z)**pinv - offset)
+
+                        sponge(ix,iy,iz) = smoothstep( tmp, 0.5_rk*params_acm%L_sponge, &
+                        0.5_rk*params_acm%L_sponge)
+                    end do
                 end do
             end do
-        end do
+        endif
 
     elseif (params_acm%sponge_type == "p-norm-insect-centered") then
         ! p-norm sponge. The shape of the sponge is dictated as the p-norm