From 3bcb6541017c399c80f3f530e57ba7c40fa6adf4 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 9 Jul 2018 10:47:53 -0700 Subject: [PATCH 01/10] added absorbing layer source code and updated acoustics_2d_radial example to use ABL --- .gitignore | 2 + examples/acoustics_2d_radial/Makefile | 24 +- examples/acoustics_2d_radial/setrun_abl.py | 260 +++++++++++++++++++++ src/2d/Makefile.classic_2d | 5 + src/2d/abl_module.f90 | 155 ++++++++++++ src/2d/flux2.f90 | 38 ++- src/2d/setaux.f90 | 27 ++- 7 files changed, 493 insertions(+), 18 deletions(-) create mode 100644 examples/acoustics_2d_radial/setrun_abl.py create mode 100644 src/2d/abl_module.f90 diff --git a/.gitignore b/.gitignore index ff041f9..b967f32 100644 --- a/.gitignore +++ b/.gitignore @@ -3,7 +3,9 @@ xclaw *.mod *~ .data +.data_abl .output +.output_abl _output/* _output _plots/* diff --git a/examples/acoustics_2d_radial/Makefile b/examples/acoustics_2d_radial/Makefile index 4ea5c6c..a77d3a0 100644 --- a/examples/acoustics_2d_radial/Makefile +++ b/examples/acoustics_2d_radial/Makefile @@ -24,7 +24,7 @@ OVERWRITE ?= True # False ==> make a copy of OUTDIR first # Environment variable FC should be set to fortran compiler, e.g. gfortran # Compiler flags can be specified here or set as an environment variable -FFLAGS ?= +FFLAGS ?= # --------------------------------- # package sources for this program: @@ -34,7 +34,7 @@ include $(CLAW)/classic/src/2d/Makefile.classic_2d # --------------------------------------- # package sources specifically to exclude -# (i.e. if a custom replacement source +# (i.e. if a custom replacement source # under a different name is provided) # --------------------------------------- @@ -53,9 +53,27 @@ SOURCES = \ qinit.f \ setprob.f \ $(CLAW)/riemann/src/rpn2_acoustics.f90 \ - $(CLAW)/riemann/src/rpt2_acoustics.f90 + $(CLAW)/riemann/src/rpt2_acoustics.f90 #------------------------------------------------------------------- # Include Makefile containing standard definitions and make options: include $(CLAWMAKE) +#------------------------------------------------------------------- +# Add make options that use absorbing boundary layer setrun +.data_abl: setrun_abl.py $(MAKEFILE_LIST) ; + $(MAKE) data_abl -f $(MAKEFILE_LIST) + +data_abl: $(MAKEFILE_LIST); + -rm -f .data_abl + $(CLAW_PYTHON) setrun_abl.py $(CLAW_PKG) + touch .data_abl + +.output_abl: $(EXE) .data_abl $(MAKEFILE_LIST); + $(MAKE) output_abl -f $(MAKEFILE_LIST) + +output_abl: $(MAKEFILE_LIST); + -rm -f .output_abl + $(CLAW_PYTHON) $(CLAW)/clawutil/src/python/clawutil/runclaw.py $(EXE) $(OUTDIR) \ + $(OVERWRITE) $(RESTART) . $(GIT_STATUS) $(NOHUP) $(NICE) + @echo $(OUTDIR) > .output_abl diff --git a/examples/acoustics_2d_radial/setrun_abl.py b/examples/acoustics_2d_radial/setrun_abl.py new file mode 100644 index 0000000..4ee25d9 --- /dev/null +++ b/examples/acoustics_2d_radial/setrun_abl.py @@ -0,0 +1,260 @@ +""" +Module to set up run time parameters for Clawpack -- classic code. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +from __future__ import absolute_import +import os +import numpy as np + +#------------------------------ +def setrun(claw_pkg='classic'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "classic" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + + assert claw_pkg.lower() == 'classic', "Expected claw_pkg = 'classic'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + probdata.add_param('rho', 1., 'density of medium') + probdata.add_param('bulk', 4., 'bulk modulus') + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amrclaw.data for AMR) + #------------------------------------------------------------------ + + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + clawdata.lower[0] = -1.000000e+00 # xlower + clawdata.upper[0] = 1.000000e+00 # xupper + clawdata.lower[1] = -1.000000e+00 # ylower + clawdata.upper[1] = 1.000000e+00 # yupper + + # Number of grid cells: + clawdata.num_cells[0] = 50 # mx + clawdata.num_cells[1] = 50 # my + + # ------------------------- + # Absorbing boundary layer: + # ------------------------- + abldata = clawdata.abldata + # Select grid mapping for layer + # 'trigonometric', 'appelo_colonius', 'petersson_sjogreen', user defined # + abldata.abltype = 'trigonometric' + # Set depth of layer at each computation boundary + abldata.depth_lower[0] = 0.2 + abldata.depth_upper[0] = 0.2 + abldata.depth_lower[1] = 0.2 + abldata.depth_upper[1] = 0.2 + # Set additional parameters if needed + # 1) trigonometric - no additional parameters needed + # 2) appelo_colonius - epsilon, p, q (e.g. {'epsilon': 1e-4, 'p': 4, 'q': 4}) + # 3) petterson_sjogreen - epsilon (e.g. {'epsilon': 1e-4}) + abldata.parameters = {} + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 0 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 0 + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.000000 + + + # Restart from checkpoint file of a previous run? + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = 'fort.chk00006' # File to use for restart data + + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output ntimes frames at equally spaced times up to tfinal: + # Can specify num_output_times = 0 for no output + clawdata.num_output_times = 20 + clawdata.tfinal = 1.0 + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list or numpy array of output times: + # Include t0 if you want output at the initial time. + clawdata.output_times = [0., 0.1] + + elif clawdata.output_style == 3: + # Output every step_interval timesteps over total_steps timesteps: + clawdata.output_step_interval = 2 + clawdata.total_steps = 4 + clawdata.output_t0 = True # output at initial (or restart) time? + + + clawdata.output_format = 'ascii' # 'ascii', 'binary', 'netcdf' + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 0 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==True: variable time steps used based on cfl_desired, + # if dt_variable==False: fixed time steps dt = dt_initial always used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # (If dt_variable==0 then dt=dt_initial for all steps) + clawdata.dt_initial = 1.00000e-02 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1.000000e+99 + + # Desired Courant number if variable dt used + clawdata.cfl_desired = 0.900000 + # max Courant number to allow without retaking step with a smaller dt: + clawdata.cfl_max = 1.000000 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 50000 + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + + # Number of waves in the Riemann solution: + clawdata.num_waves = 2 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'vanleer' ==> van Leer + # 4 or 'mc' ==> MC limiter + clawdata.limiter = ['vanleer','vanleer'] + + clawdata.use_fwaves = False # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never + # called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, + # not recommended. + clawdata.source_split = 0 + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 or 'user' => user specified (must modify bcNamr.f to use this + # option) + # 1 or 'extrap' => extrapolation (non-reflecting outflow) + # 2 or 'periodic' => periodic (must specify this at both boundaries) + # 3 or 'wall' => solid wall for systems where q(2) is normal + # velocity + + clawdata.bc_lower[0] = 'extrap' # at xlower + clawdata.bc_upper[0] = 'extrap' # at xupper + + clawdata.bc_lower[1] = 'extrap' # at ylower + clawdata.bc_upper[1] = 'extrap' # at yupper + + return rundata + + # end of function setrun + # ---------------------- + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() diff --git a/src/2d/Makefile.classic_2d b/src/2d/Makefile.classic_2d index 37bc02b..c0be700 100644 --- a/src/2d/Makefile.classic_2d +++ b/src/2d/Makefile.classic_2d @@ -1,6 +1,11 @@ #get the directory of this makefile LIB:=$(CLAW)/classic/src/2d +#list of common modules files for classic 2d codes +COMMON_MODULES = \ + $(LIB)/abl_module.f90 + + #list of common source files for classic 2d codes COMMON_SOURCES = \ $(LIB)/qinit.f \ diff --git a/src/2d/abl_module.f90 b/src/2d/abl_module.f90 new file mode 100644 index 0000000..4349ead --- /dev/null +++ b/src/2d/abl_module.f90 @@ -0,0 +1,155 @@ +module abl_module + + implicit none + integer, parameter :: NO_ABL = 0 + integer, parameter :: TRIG = 1 + integer, parameter :: AC = 2 ! Appelo & Colonius + integer, parameter :: PS = 3 ! Petersson & Sjogreen + + integer, parameter :: MAX_PARAM = 3 + + real(kind=8), parameter :: PIH = 2.d0*datan(1.d0) + + + integer :: abl_type = 0 + integer :: p_AC, q_AC + real(kind=8) :: eps_AC, eps_PS + real(kind=8) :: xdepth(2), ydepth(2), xpos(2), ypos(2) + logical :: initialized = .false. + +contains + + subroutine initialize() + + implicit none + + real(kind=8) :: xlower, xupper, ylower, yupper + + ! Read in ABL data + ! open the unit with new routine from Clawpack 4.4 to skip over + ! comment lines starting with #: + call opendatafile(7, 'abl.data') + + ! Read in parameters: + read(7,*) abl_type + if (abl_type == TRIG) then + print *, "ABL is trigonometric" + else if (abl_type == AC) then + read(7,*) eps_AC + read(7,*) p_AC + read(7,*) q_AC + print *, "ABL is Appelo/Colonius with (eps,p,q)=", eps_AC, p_AC, q_AC + else if (abl_type == PS) then + read(7,*) eps_PS + print *, "ABL is Petersson/Sjogreen with eps=", eps_PS + end if + read(7,*) xdepth(1), ydepth(1) + read(7,*) xdepth(2), ydepth(2) + close(7) + + ! Read in grid parameters + ! open the unit with new routine from Clawpack 4.4 to skip over + ! comment lines starting with #: + call opendatafile(7, 'claw.data') + + read(7,*) xlower ! num_dim + read(7,*) xlower, ylower ! lower + read(7,*) xupper, yupper ! upper + close(7) + + ! Compute ABL positions + xpos(1) = xlower + xdepth(1) + xpos(2) = xupper - xdepth(2) + ypos(1) = ylower + ydepth(1) + ypos(2) = yupper - ydepth(2) + + initialized = .true. + + end subroutine initialize + + subroutine set_scaling_factor(mbc,mx,my,xlower,ylower,dx,dy,maux,aux) + + implicit none + + integer, intent(in) :: mbc, mx, my, maux + real(kind=8), intent(in) :: xlower, ylower, dx, dy + real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) + + integer :: i, j + real(kind=8) :: xcenter, ycenter, lower(2) + + ! initialize module if needed + if (.not. initialized) then + call initialize() + end if + + if (abl_type == 0) then + return + end if + + ! Loop over all cell edges + do j=1-mbc,my+mbc + + ycenter = ylower + (j-0.5d0)*dy + + do i=1-mbc,mx+mbc + + xcenter = xlower + (i-0.5d0)*dx + + ! determine lower cell edge locations + ! constrained to the computational domain without ghost cells + lower(1) = max(xcenter-0.5d0*dx,xpos(1)-xdepth(1)) + lower(1) = min(lower(1),xpos(2)+xdepth(2)-dx) + lower(2) = max(ycenter-0.5d0*dy,ypos(1)-ydepth(1)) + lower(2) = min(lower(2),ypos(2)+ydepth(2)-dy) + + ! Calculate inverse of Jacobian 1/(1 + g'(x_k-1/2)) & 1/(1 + g'(y_k-1/2)) + aux(maux-1,i,j) = 1.d0/(1.d0 + gp(lower(1),1)) + aux(maux,i,j) = 1.d0/(1.d0 + gp(lower(2),2)) + + end do + end do + + end subroutine set_scaling_factor + + function gp(z,direction) + + implicit none + + real(kind=8) :: gp, z + integer :: direction + + real(kind=8) :: pos(2), depth(2), sig, val + + gp = 0.d0 + sig = 0.d0 + if (direction == 1) then + pos = xpos + depth = xdepth + else if (direction == 2) then + pos = ypos + depth = ydepth + end if + + if (z < pos(1)) then + val = (pos(1)-z)/depth(1) + else if (z > pos(2)) then + val = (z-pos(2))/depth(2) + else + return + end if + + if (abl_type == TRIG) then + gp = dtan(PIH*val)**2 + else if (abl_type == AC) then + sig = (1.0 - eps_AC)*(1.0 - (1.0 - val)**p_AC)**q_AC + gp = sig/(1.0-sig) + else if (abl_type == PS) then + sig = (1.0 - eps_PS)* & + val**6*(462.0 - 1980.0*val + 3465.0*val**2 - 3080.0*val**3 + 1386.0*val**4 - 252.0*val**5) + gp = sig/(1.0-sig) + end if + + end function gp + +end module abl_module diff --git a/src/2d/flux2.f90 b/src/2d/flux2.f90 index a805e41..02cae5f 100644 --- a/src/2d/flux2.f90 +++ b/src/2d/flux2.f90 @@ -57,6 +57,8 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & !f2py external rpn2, rpt2 !f2py intent(callback) rpn2,rpt2 + use abl_module, only: abl_type + implicit double precision (a-h,o-z) integer num_aux external rpn2,rpt2 @@ -71,6 +73,7 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & dimension gadd(num_eqn, 2, 1-num_ghost:maxm+num_ghost) dimension dtdx1d(1-num_ghost:maxm+num_ghost) + dimension aux1(num_aux,1-num_ghost:maxm+num_ghost) dimension aux2(num_aux,1-num_ghost:maxm+num_ghost) dimension aux3(num_aux,1-num_ghost:maxm+num_ghost) @@ -84,7 +87,7 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & limit = .false. do mw=1,num_waves - if (mthlim(mw) > 0) limit = .TRUE. + if (mthlim(mw) > 0) limit = .TRUE. end do ! # initialize flux increments: @@ -110,6 +113,17 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & call rpn2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, & aux2,aux2,wave,s,amdq,apdq) + ! apply ABL scaling to propagation speed and update fluctuations if needed + if (abl_type .ne. 0) then + do mw=1,num_waves + s(mw,:) = s(mw,:)*aux2(num_aux-2+ixy,:) + end do + do m=1,num_eqn + amdq(m,:) = amdq(m,:)*aux2(num_aux-2+ixy,:) + apdq(m,:) = apdq(m,:)*aux2(num_aux-2+ixy,:) + end do + end if + ! # Set qadd for the donor-cell upwind method (Godunov) do i = 1, mx+1 do m = 1, num_eqn ! qadd(:,i-1) is still in cache from last cycle of outer loop @@ -146,7 +160,6 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & ! # Second order corrections: if (use_fwave.eqv. .FALSE. ) then do i = 2-num_ghost, mx+num_ghost - ! # modified in Version 4.3 to use average only in cqxx, not transverse dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) do m = 1, num_eqn cqxx(m,i) = 0.d0 @@ -204,6 +217,18 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & call rpt2(ixy,1,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, & aux1,aux2,aux3,amdq,bmasdq,bpasdq) +! # Adjust for ABL if needed + if (abl_type .ne. 0) then + do m=1,num_eqn + bmasdq(m,2-num_ghost:mx+num_ghost) = & + bmasdq(m,2-num_ghost:mx+num_ghost) * & + aux2(num_aux+1-ixy,1-num_ghost:mx+num_ghost-1) + bpasdq(m,2-num_ghost:mx+num_ghost) = & + bpasdq(m,2-num_ghost:mx+num_ghost) * & + aux3(num_aux+1-ixy,1-num_ghost:mx+num_ghost-1) + end do + end if + ! # modify flux below and above by B^- A^- Delta q and B^+ A^- Delta q: do i = 1, mx+1 ! Having two inner loops here allows traversal of gadd in memory-contiguous order @@ -214,12 +239,19 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & gadd(m,2,i-1) = -0.5d0*dtdx1d(i-1) * bpasdq(m,i) end do end do - ! # split the right-going flux difference into down-going and up-going: call rpt2(ixy,2,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, & aux1,aux2,aux3,apdq,bmasdq,bpasdq) +! # Adjust for ABL if needed + if (abl_type .ne. 0) then + do m=1,num_eqn + bmasdq(m,:) = bmasdq(m,:) * aux2(num_aux+1-ixy,:) + bpasdq(m,:) = bpasdq(m,:) * aux3(num_aux+1-ixy,:) + end do + end if + ! # modify flux below and above by B^- A^+ Delta q and B^+ A^+ Delta q: do i = 1, mx+1 do m = 1, num_eqn diff --git a/src/2d/setaux.f90 b/src/2d/setaux.f90 index ded13c0..73ab084 100644 --- a/src/2d/setaux.f90 +++ b/src/2d/setaux.f90 @@ -1,16 +1,19 @@ subroutine setaux(mbc,mx,my,xlower,ylower,dx,dy,maux,aux) - ! Called at start of computation before calling qinit, and - ! when AMR is used, also called every time a new grid patch is created. - ! Use to set auxiliary arrays aux(1:maux, 1-mbc:mx+mbc, 1-mbc:my+mbc). - ! Note that ghost cell values may need to be set if the aux arrays - ! are used by the Riemann solver(s). - ! - ! This default version does nothing. - - implicit none - integer, intent(in) :: mbc,mx,my,maux - real(kind=8), intent(in) :: xlower,ylower,dx,dy - real(kind=8), intent(out) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) + ! Called at start of computation before calling qinit, and + ! when AMR is used, also called every time a new grid patch is created. + ! Use to set auxiliary arrays aux(1:maux, 1-mbc:mx+mbc, 1-mbc:my+mbc). + ! Note that ghost cell values may need to be set if the aux arrays + ! are used by the Riemann solver(s). + + use abl_module, only: set_scaling_factor + + implicit none + integer, intent(in) :: mbc,mx,my,maux + real(kind=8), intent(in) :: xlower,ylower,dx,dy + real(kind=8), intent(out) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) + + aux(:,:,:) = 0.d0 + call set_scaling_factor(mbc,mx,my,xlower,ylower,dx,dy,maux,aux) end subroutine setaux From 32681eb58ad6e61ebf8f93d660e41a2f01bcdcd2 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Tue, 10 Jul 2018 15:38:01 -0700 Subject: [PATCH 02/10] reverted Makefile and created Makefile_abl and setplot_abl instead --- examples/acoustics_2d_radial/Makefile | 24 +--- examples/acoustics_2d_radial/Makefile_abl | 60 ++++++++ examples/acoustics_2d_radial/setplot_abl.py | 150 ++++++++++++++++++++ 3 files changed, 213 insertions(+), 21 deletions(-) create mode 100644 examples/acoustics_2d_radial/Makefile_abl create mode 100644 examples/acoustics_2d_radial/setplot_abl.py diff --git a/examples/acoustics_2d_radial/Makefile b/examples/acoustics_2d_radial/Makefile index a77d3a0..4ea5c6c 100644 --- a/examples/acoustics_2d_radial/Makefile +++ b/examples/acoustics_2d_radial/Makefile @@ -24,7 +24,7 @@ OVERWRITE ?= True # False ==> make a copy of OUTDIR first # Environment variable FC should be set to fortran compiler, e.g. gfortran # Compiler flags can be specified here or set as an environment variable -FFLAGS ?= +FFLAGS ?= # --------------------------------- # package sources for this program: @@ -34,7 +34,7 @@ include $(CLAW)/classic/src/2d/Makefile.classic_2d # --------------------------------------- # package sources specifically to exclude -# (i.e. if a custom replacement source +# (i.e. if a custom replacement source # under a different name is provided) # --------------------------------------- @@ -53,27 +53,9 @@ SOURCES = \ qinit.f \ setprob.f \ $(CLAW)/riemann/src/rpn2_acoustics.f90 \ - $(CLAW)/riemann/src/rpt2_acoustics.f90 + $(CLAW)/riemann/src/rpt2_acoustics.f90 #------------------------------------------------------------------- # Include Makefile containing standard definitions and make options: include $(CLAWMAKE) -#------------------------------------------------------------------- -# Add make options that use absorbing boundary layer setrun -.data_abl: setrun_abl.py $(MAKEFILE_LIST) ; - $(MAKE) data_abl -f $(MAKEFILE_LIST) - -data_abl: $(MAKEFILE_LIST); - -rm -f .data_abl - $(CLAW_PYTHON) setrun_abl.py $(CLAW_PKG) - touch .data_abl - -.output_abl: $(EXE) .data_abl $(MAKEFILE_LIST); - $(MAKE) output_abl -f $(MAKEFILE_LIST) - -output_abl: $(MAKEFILE_LIST); - -rm -f .output_abl - $(CLAW_PYTHON) $(CLAW)/clawutil/src/python/clawutil/runclaw.py $(EXE) $(OUTDIR) \ - $(OVERWRITE) $(RESTART) . $(GIT_STATUS) $(NOHUP) $(NICE) - @echo $(OUTDIR) > .output_abl diff --git a/examples/acoustics_2d_radial/Makefile_abl b/examples/acoustics_2d_radial/Makefile_abl new file mode 100644 index 0000000..e821b7c --- /dev/null +++ b/examples/acoustics_2d_radial/Makefile_abl @@ -0,0 +1,60 @@ + +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = classic # Clawpack package to use +EXE = xclaw # Executable to create +SETRUN_FILE = setrun_abl.py # File containing function to make data +OUTDIR = _output_abl # Directory for output +SETPLOT_FILE = setplot_abl.py # File containing function to set plots +PLOTDIR = _plots_abl # Directory for plots + +OVERWRITE ?= True # False ==> make a copy of OUTDIR first + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +include $(CLAW)/classic/src/2d/Makefile.classic_2d + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + + +# --------------------------------- +# List of custom sources for this program: +# --------------------------------- + +MODULES = \ + +SOURCES = \ + qinit.f \ + setprob.f \ + $(CLAW)/riemann/src/rpn2_acoustics.f90 \ + $(CLAW)/riemann/src/rpt2_acoustics.f90 + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) diff --git a/examples/acoustics_2d_radial/setplot_abl.py b/examples/acoustics_2d_radial/setplot_abl.py new file mode 100644 index 0000000..486b063 --- /dev/null +++ b/examples/acoustics_2d_radial/setplot_abl.py @@ -0,0 +1,150 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import os +if os.path.exists('./1drad/_output'): + qref_dir = os.path.abspath('./1drad/_output') +else: + qref_dir = None + print("Directory ./1drad/_output not found") + + +#-------------------------- +def setplot(plotdata=None): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of clawpack.visclaw.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + if plotdata is None: + from clawpack.visclaw.data import ClawPlotData + plotdata = ClawPlotData() + + + from clawpack.visclaw import colormaps + from clawpack.clawutil.data import ClawData + + # determine original computational domain + abldata = ClawData() + abldata.read(plotdata.outdir + '/abl.data', force=True) + clawdata = ClawData() + clawdata.read(plotdata.outdir + '/claw.data', force=True) + clawdata.lower[0] += abldata.depth_lower[0] + clawdata.upper[0] -= abldata.depth_upper[0] + clawdata.lower[1] += abldata.depth_lower[1] + clawdata.upper[1] -= abldata.depth_upper[1] + + plotdata.clearfigures() # clear any old figures,axes,items data + + + # Figure for pressure + # ------------------- + + plotfigure = plotdata.new_plotfigure(name='Pressure', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Pressure' + plotaxes.scaled = True # so aspect ratio is 1 + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = 0 + plotitem.pcolor_cmap = colormaps.blue_yellow_red + plotitem.pcolor_cmin = -2.0 + plotitem.pcolor_cmax = 2.0 + plotitem.add_colorbar = True + + def plot_original_domain(current_data): + from matplotlib.pyplot import gca + ax = gca() + x = [clawdata.lower[0], clawdata.upper[0], clawdata.upper[0], \ + clawdata.lower[0], clawdata.lower[0]] + y = [clawdata.lower[1], clawdata.lower[1], clawdata.upper[1], \ + clawdata.upper[1], clawdata.lower[1]] + ax.plot(x, y, '--k') + + plotaxes.afteraxes = plot_original_domain + + # Figure for scatter plot + # ----------------------- + + plotfigure = plotdata.new_plotfigure(name='scatter', figno=3) + plotfigure.show = (qref_dir is not None) # don't plot if 1d solution is missing + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,1.5] + plotaxes.ylimits = [-2.,4.] + plotaxes.title = 'Scatter plot' + + # Set up for item on these axes: scatter of 2d data + plotitem = plotaxes.new_plotitem(plot_type='1d_from_2d_data') + + def p_vs_r(current_data): + # Return radius of each grid cell and p value in the cell + from numpy.ma import MaskedArray, masked_where + from pylab import sqrt + x = current_data.x + y = current_data.y + r = sqrt(x**2 + y**2) + r = masked_where(x < clawdata.lower[0], r) + r = masked_where(x > clawdata.upper[0], r) + r = masked_where(y < clawdata.lower[1], r) + r = masked_where(y > clawdata.upper[1], r) + q = current_data.q + p = MaskedArray(q[0,:,:], mask=r.mask) + return r,p + + plotitem.map_2d_to_1d = p_vs_r + plotitem.plot_var = 0 + plotitem.plotstyle = 'o' + plotitem.color = 'b' + plotitem.show = (qref_dir is not None) # show on plot? + + # Set up for item on these axes: 1d reference solution + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.outdir = qref_dir + plotitem.plot_var = 0 + plotitem.plotstyle = '-' + plotitem.color = 'r' + plotitem.kwargs = {'linewidth': 2} + plotitem.show = True # show on plot? + + def make_legend(current_data): + import matplotlib.pyplot as plt + plt.legend(('2d data', '1d reference solution')) + + plotaxes.afteraxes = make_legend + + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via clawpack.visclaw.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.html_movie = 'JSAnimation' # new style, or "4.x" for old style + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + + return plotdata From 0a52c78b62197b8dbed47ea6166388f7e3f7e71c Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Tue, 24 Jul 2018 17:47:44 -0700 Subject: [PATCH 03/10] for clarity don't redefine clawdata attributes in setplot_abl.py --- examples/acoustics_2d_radial/setplot_abl.py | 22 ++++++++++----------- 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/examples/acoustics_2d_radial/setplot_abl.py b/examples/acoustics_2d_radial/setplot_abl.py index 486b063..2ae5cd6 100644 --- a/examples/acoustics_2d_radial/setplot_abl.py +++ b/examples/acoustics_2d_radial/setplot_abl.py @@ -41,10 +41,10 @@ def setplot(plotdata=None): abldata.read(plotdata.outdir + '/abl.data', force=True) clawdata = ClawData() clawdata.read(plotdata.outdir + '/claw.data', force=True) - clawdata.lower[0] += abldata.depth_lower[0] - clawdata.upper[0] -= abldata.depth_upper[0] - clawdata.lower[1] += abldata.depth_lower[1] - clawdata.upper[1] -= abldata.depth_upper[1] + x1 = clawdata.lower[0] + abldata.depth_lower[0] + x2 = clawdata.upper[0] - abldata.depth_upper[0] + y1 = clawdata.lower[1] + abldata.depth_lower[1] + y2 = clawdata.upper[1] - abldata.depth_upper[1] plotdata.clearfigures() # clear any old figures,axes,items data @@ -72,10 +72,8 @@ def setplot(plotdata=None): def plot_original_domain(current_data): from matplotlib.pyplot import gca ax = gca() - x = [clawdata.lower[0], clawdata.upper[0], clawdata.upper[0], \ - clawdata.lower[0], clawdata.lower[0]] - y = [clawdata.lower[1], clawdata.lower[1], clawdata.upper[1], \ - clawdata.upper[1], clawdata.lower[1]] + x = [x1,x2,x2,x1,x1] + y = [y1,y1,y2,y2,y1] ax.plot(x, y, '--k') plotaxes.afteraxes = plot_original_domain @@ -102,10 +100,10 @@ def p_vs_r(current_data): x = current_data.x y = current_data.y r = sqrt(x**2 + y**2) - r = masked_where(x < clawdata.lower[0], r) - r = masked_where(x > clawdata.upper[0], r) - r = masked_where(y < clawdata.lower[1], r) - r = masked_where(y > clawdata.upper[1], r) + r = masked_where(x < x1, r) + r = masked_where(x > x2, r) + r = masked_where(y < y1, r) + r = masked_where(y > y2, r) q = current_data.q p = MaskedArray(q[0,:,:], mask=r.mask) return r,p From 5f7a38f761fcb149f42b36542b547f6267f068f2 Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Tue, 24 Jul 2018 17:49:18 -0700 Subject: [PATCH 04/10] added text to plot of ABL --- examples/acoustics_2d_radial/setplot_abl.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/acoustics_2d_radial/setplot_abl.py b/examples/acoustics_2d_radial/setplot_abl.py index 2ae5cd6..8aeb65c 100644 --- a/examples/acoustics_2d_radial/setplot_abl.py +++ b/examples/acoustics_2d_radial/setplot_abl.py @@ -70,11 +70,12 @@ def setplot(plotdata=None): plotitem.add_colorbar = True def plot_original_domain(current_data): - from matplotlib.pyplot import gca + from matplotlib.pyplot import gca, text ax = gca() x = [x1,x2,x2,x1,x1] y = [y1,y1,y2,y2,y1] ax.plot(x, y, '--k') + text(-0.6,1.05,'Absorbing Boundary Layer') plotaxes.afteraxes = plot_original_domain From 87f7745be277574b9a1270f293b9eeeefad1ad34 Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Tue, 24 Jul 2018 17:53:39 -0700 Subject: [PATCH 05/10] modify legend for scatter plot --- examples/acoustics_2d_radial/setplot_abl.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/acoustics_2d_radial/setplot_abl.py b/examples/acoustics_2d_radial/setplot_abl.py index 8aeb65c..18ceea0 100644 --- a/examples/acoustics_2d_radial/setplot_abl.py +++ b/examples/acoustics_2d_radial/setplot_abl.py @@ -126,7 +126,7 @@ def p_vs_r(current_data): def make_legend(current_data): import matplotlib.pyplot as plt - plt.legend(('2d data', '1d reference solution')) + plt.legend(('2d data (interior only)', '1d reference solution')) plotaxes.afteraxes = make_legend From 74df71c25a8dd68537e1fafa3e6fd9baf7f04be3 Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Tue, 24 Jul 2018 17:58:19 -0700 Subject: [PATCH 06/10] started adding info on ABL to README --- examples/acoustics_2d_radial/README.rst | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/examples/acoustics_2d_radial/README.rst b/examples/acoustics_2d_radial/README.rst index ff8878a..c5cc8b2 100644 --- a/examples/acoustics_2d_radial/README.rst +++ b/examples/acoustics_2d_radial/README.rst @@ -10,3 +10,21 @@ compute the "true solution" and then setplot.py contains code to produce a scatter plot of the computed 2d pressure vs. distance from origin and compare with the 1d solution. +Extrapolation BCs +------------------ + +The code is set up to use extrapolation boundary conditions at all +boundaries. This does a reasonably good job of providing non-reflecting +boundaries, but there are some artifacts visible at later times. + +Absorbing boundary layer +------------------------ + +New cababilities have been added to Clawpack 5.5.0 to allow extending the +computational domain with an aborbing boundary layer that does a better job +of eliminating artificial reflections. [Add more discussion.] + +To try this version:: + + make all -f Makefile_abl + From 0209fbc73446c9ec2a683054c482571ac9645628 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Wed, 8 Aug 2018 09:05:57 -0700 Subject: [PATCH 07/10] no longer need .output_abl or .data_abl in .gitignore as they will not be generated --- .gitignore | 2 -- 1 file changed, 2 deletions(-) diff --git a/.gitignore b/.gitignore index b967f32..ff041f9 100644 --- a/.gitignore +++ b/.gitignore @@ -3,9 +3,7 @@ xclaw *.mod *~ .data -.data_abl .output -.output_abl _output/* _output _plots/* From 33d6bca320b555a11def2a80d28df7c055100ae7 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Tue, 5 Mar 2019 19:37:57 -0800 Subject: [PATCH 08/10] implemented aux-less absorbing boundary layer --- examples/acoustics_2d_radial/setrun_abl.py | 2 +- src/2d/abl_module.f90 | 71 +++++++--- src/2d/claw2.f | 2 +- src/2d/claw2ez.f | 5 + src/2d/flux2.f90 | 80 +++++++----- src/2d/setaux.f90 | 27 ++-- src/2d/step2.f90 | 145 ++++++++++++++------- 7 files changed, 218 insertions(+), 114 deletions(-) diff --git a/examples/acoustics_2d_radial/setrun_abl.py b/examples/acoustics_2d_radial/setrun_abl.py index 4ee25d9..02568e4 100644 --- a/examples/acoustics_2d_radial/setrun_abl.py +++ b/examples/acoustics_2d_radial/setrun_abl.py @@ -97,7 +97,7 @@ def setrun(claw_pkg='classic'): clawdata.num_eqn = 3 # Number of auxiliary variables in the aux array (initialized in setaux) - clawdata.num_aux = 0 + clawdata.num_aux = 2 # Index of aux array corresponding to capacity function, if there is one: clawdata.capa_index = 0 diff --git a/src/2d/abl_module.f90 b/src/2d/abl_module.f90 index 4349ead..49b436d 100644 --- a/src/2d/abl_module.f90 +++ b/src/2d/abl_module.f90 @@ -67,51 +67,82 @@ subroutine initialize() end subroutine initialize - subroutine set_scaling_factor(mbc,mx,my,xlower,ylower,dx,dy,maux,aux) + subroutine set_scaling_factor(xlower, ylower, dx, dy, mx, my, num_ghost, & + abl_center, abl_edge, & + i_abl_lower, i_abl_upper, j_abl_lower, j_abl_upper) implicit none - integer, intent(in) :: mbc, mx, my, maux + integer, intent(in) :: mx, my, num_ghost real(kind=8), intent(in) :: xlower, ylower, dx, dy - real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) + integer, intent(out) :: i_abl_lower, i_abl_upper, j_abl_lower, j_abl_upper + real(kind=8), intent(out) :: abl_center(1-num_ghost:mx+num_ghost, & + 1-num_ghost:my+num_ghost, 2) + real(kind=8), intent(out) :: abl_edge(1-num_ghost:mx+num_ghost, & + 1-num_ghost:my+num_ghost, 2) integer :: i, j - real(kind=8) :: xcenter, ycenter, lower(2) + real(kind=8) :: xcenter, ycenter, xedge, yedge ! initialize module if needed if (.not. initialized) then call initialize() end if - if (abl_type == 0) then - return - end if - - ! Loop over all cell edges - do j=1-mbc,my+mbc + ! determine which indices are in the layer + i_abl_lower = int((1.d0 + 1.d-14)*(xpos(1)-xlower)/dx) + i_abl_upper = int((1.d0 + 1.d-14)*(xpos(2)-xlower)/dx) + j_abl_lower = int((1.d0 + 1.d-14)*(ypos(1)-ylower)/dy) + j_abl_upper = int((1.d0 + 1.d-14)*(ypos(2)-ylower)/dy) + ! Loop over all cells + do j=1-num_ghost,my+num_ghost ycenter = ylower + (j-0.5d0)*dy - do i=1-mbc,mx+mbc - + do i=1-num_ghost,mx+num_ghost xcenter = xlower + (i-0.5d0)*dx ! determine lower cell edge locations ! constrained to the computational domain without ghost cells - lower(1) = max(xcenter-0.5d0*dx,xpos(1)-xdepth(1)) - lower(1) = min(lower(1),xpos(2)+xdepth(2)-dx) - lower(2) = max(ycenter-0.5d0*dy,ypos(1)-ydepth(1)) - lower(2) = min(lower(2),ypos(2)+ydepth(2)-dy) + xedge = max(xcenter-0.5d0*dx,xpos(1)-xdepth(1)) + xedge = min(xedge,xpos(2)+xdepth(2)-dx) + yedge = max(ycenter-0.5d0*dy,ypos(1)-ydepth(1)) + yedge = min(yedge,ypos(2)+ydepth(2)-dy) - ! Calculate inverse of Jacobian 1/(1 + g'(x_k-1/2)) & 1/(1 + g'(y_k-1/2)) - aux(maux-1,i,j) = 1.d0/(1.d0 + gp(lower(1),1)) - aux(maux,i,j) = 1.d0/(1.d0 + gp(lower(2),2)) + ! Calculate inverse of Jacobian at cell edges + abl_edge(i,j,1) = 1.d0/(1.d0 + gp(xedge,1)) + abl_edge(i,j,2) = 1.d0/(1.d0 + gp(yedge,2)) + + ! Calculate cell width ratio at cell centers + ! (use inverse of Jacobian for methods that dont have g) + abl_center(i,j,1) = 1.d0/(1.d0 + gp(xcenter,1)) + abl_center(i,j,2) = 1.d0/(1.d0 + gp(ycenter,2)) end do + end do end subroutine set_scaling_factor + subroutine scale_for_abl(abl_factor,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & + scaled_value,offset) + + implicit none + integer, intent(in) :: num_ghost, maxm, mx, i_abl_lower, i_abl_upper, offset + double precision, intent(in) :: abl_factor(1-num_ghost:maxm+num_ghost) + double precision, intent(inout) :: scaled_value(1-num_ghost:maxm+num_ghost) + + integer :: i + + do i = 1-num_ghost, i_abl_lower + scaled_value(i) = abl_factor(i-offset)*scaled_value(i) + end do + do i = i_abl_upper, mx+num_ghost + scaled_value(i) = abl_factor(i-offset)*scaled_value(i) + end do + + end subroutine scale_for_abl + function gp(z,direction) implicit none @@ -152,4 +183,6 @@ function gp(z,direction) end function gp + + end module abl_module diff --git a/src/2d/claw2.f b/src/2d/claw2.f index 3de21fa..32debca 100644 --- a/src/2d/claw2.f +++ b/src/2d/claw2.f @@ -661,7 +661,7 @@ subroutine claw2(meqn,mwaves,maux,mbc,mx,my, c call step2(maxm,meqn,mwaves,maux,mbc,mx,my, & work(i0qwrk1),q,aux, - & dx,dy,dt,method,mthlim,cfl, + & xlower,ylower,dx,dy,dt,method,mthlim,cfl, & work(i0qadd),work(i0fadd),work(i0gadd), & work(i0q1d),work(i0dtdx1),work(i0dtdy1), & work(i0aux1),work(i0aux2),work(i0aux3), diff --git a/src/2d/claw2ez.f b/src/2d/claw2ez.f index 78441bf..9afc68a 100644 --- a/src/2d/claw2ez.f +++ b/src/2d/claw2ez.f @@ -11,6 +11,8 @@ subroutine claw2ez ! No arguments c c Authors: Randall J. LeVeque, Grady I. Lemoine c + use abl_module, only: initialize_abl=>initialize + implicit double precision (a-h,o-z) external bc2,rpn2,rpt2,src2,b4step2 @@ -233,6 +235,9 @@ subroutine claw2ez ! No arguments go to 900 end if +c # initialize abl module + call initialize_abl() + c c # set initial conditions: c diff --git a/src/2d/flux2.f90 b/src/2d/flux2.f90 index 02cae5f..6d83345 100644 --- a/src/2d/flux2.f90 +++ b/src/2d/flux2.f90 @@ -5,7 +5,8 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & q1d,dtdx1d,aux1,aux2,aux3,method,mthlim, & qadd,fadd,gadd,cfl1d,wave,s, & - amdq,apdq,cqxx,bmasdq,bpasdq,rpn2,rpt2,use_fwave) + amdq,apdq,cqxx,bmasdq,bpasdq,rpn2,rpt2,use_fwave, & + abl_center,abl_edge1,abl_edge2,abl_edge3,i_abl_lower,i_abl_upper) ! ===================================================== ! # Compute the modification to fluxes f and g that are generated by @@ -57,7 +58,7 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & !f2py external rpn2, rpt2 !f2py intent(callback) rpn2,rpt2 - use abl_module, only: abl_type + use abl_module, only: abl_type, scale_for_abl implicit double precision (a-h,o-z) integer num_aux @@ -73,7 +74,8 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & dimension gadd(num_eqn, 2, 1-num_ghost:maxm+num_ghost) dimension dtdx1d(1-num_ghost:maxm+num_ghost) - + dimension dtdx1d_abl(1-num_ghost:maxm+num_ghost) + dimension dtdxave_abl(1-num_ghost:maxm+num_ghost) dimension aux1(num_aux,1-num_ghost:maxm+num_ghost) dimension aux2(num_aux,1-num_ghost:maxm+num_ghost) dimension aux3(num_aux,1-num_ghost:maxm+num_ghost) @@ -81,6 +83,11 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & dimension s(num_waves,1-num_ghost:maxm+num_ghost) dimension wave(num_eqn, num_waves, 1-num_ghost:maxm+num_ghost) + dimension abl_center(1-num_ghost:maxm+num_ghost) + dimension abl_edge1(1-num_ghost:maxm+num_ghost) + dimension abl_edge2(1-num_ghost:maxm+num_ghost) + dimension abl_edge3(1-num_ghost:maxm+num_ghost) + dimension method(7),mthlim(num_waves) logical :: limit, use_fwave common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom @@ -113,24 +120,23 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & call rpn2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, & aux2,aux2,wave,s,amdq,apdq) - ! apply ABL scaling to propagation speed and update fluctuations if needed - if (abl_type .ne. 0) then - do mw=1,num_waves - s(mw,:) = s(mw,:)*aux2(num_aux-2+ixy,:) - end do - do m=1,num_eqn - amdq(m,:) = amdq(m,:)*aux2(num_aux-2+ixy,:) - apdq(m,:) = apdq(m,:)*aux2(num_aux-2+ixy,:) - end do +! adjust dtdx1d for abl if needed + do i = 0, mx+1 + dtdx1d_abl(i) = dtdx1d(i) + end do + if (abl_type /= 0) then + call scale_for_abl(abl_center,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & + dtdx1d_abl,0) end if + ! # Set qadd for the donor-cell upwind method (Godunov) do i = 1, mx+1 do m = 1, num_eqn ! qadd(:,i-1) is still in cache from last cycle of outer loop - qadd(m,i-1) = qadd(m,i-1) - dtdx1d(i-1)*amdq(m,i) + qadd(m,i-1) = qadd(m,i-1) - dtdx1d_abl(i-1)*amdq(m,i) end do do m = 1, num_eqn - qadd(m,i) = qadd(m,i) - dtdx1d(i)*apdq(m,i) + qadd(m,i) = qadd(m,i) - dtdx1d_abl(i)*apdq(m,i) end do end do @@ -159,15 +165,23 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & ! # Second order corrections: if (use_fwave.eqv. .FALSE. ) then + ! # modified in Version 4.3 to use average only in cqxx, not transverse + ! adjust dtdxave for abl (if ndeed) + do i = 2-num_ghost, mx+num_ghost + dtdxave_abl(i) = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) + end do + if (abl_type /= 0) then + call scale_for_abl(abl_edge2,num_ghost,maxm,mx,i_abl_lower, & + i_abl_upper,dtdxave_abl,0) + end if do i = 2-num_ghost, mx+num_ghost - dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) do m = 1, num_eqn cqxx(m,i) = 0.d0 end do do mw=1,num_waves ! Traverse the wave array in memory-contiguous fashion do m=1,num_eqn cqxx(m,i) = cqxx(m,i) + dabs(s(mw,i)) & - * (1.d0 - dabs(s(mw,i))*dtdxave) * wave(m,mw,i) + * (1.d0 - dabs(s(mw,i))*dtdxave_abl(i)) * wave(m,mw,i) enddo enddo do m = 1, num_eqn @@ -175,6 +189,10 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & end do enddo else + if (abl_type /= 0) then + print *, "ABL not yet implemented for fwave splitting" + stop + end if do i = 2-num_ghost, mx+num_ghost dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) do m = 1, num_eqn @@ -217,16 +235,12 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & call rpt2(ixy,1,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, & aux1,aux2,aux3,amdq,bmasdq,bpasdq) -! # Adjust for ABL if needed - if (abl_type .ne. 0) then - do m=1,num_eqn - bmasdq(m,2-num_ghost:mx+num_ghost) = & - bmasdq(m,2-num_ghost:mx+num_ghost) * & - aux2(num_aux+1-ixy,1-num_ghost:mx+num_ghost-1) - bpasdq(m,2-num_ghost:mx+num_ghost) = & - bpasdq(m,2-num_ghost:mx+num_ghost) * & - aux3(num_aux+1-ixy,1-num_ghost:mx+num_ghost-1) - end do +! adjust transverse fluxes for abl if needed + if (abl_type /= 0) then + call scale_for_abl(abl_edge1,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & + bmasdq,1) + call scale_for_abl(abl_edge3,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & + bpasdq,1) end if ! # modify flux below and above by B^- A^- Delta q and B^+ A^- Delta q: @@ -240,18 +254,20 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & end do end do + ! # split the right-going flux difference into down-going and up-going: call rpt2(ixy,2,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, & aux1,aux2,aux3,apdq,bmasdq,bpasdq) -! # Adjust for ABL if needed - if (abl_type .ne. 0) then - do m=1,num_eqn - bmasdq(m,:) = bmasdq(m,:) * aux2(num_aux+1-ixy,:) - bpasdq(m,:) = bpasdq(m,:) * aux3(num_aux+1-ixy,:) - end do +! adjust transverse fluxes for abl if needed + if (abl_type /= 0) then + call scale_for_abl(abl_edge1,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & + bmasdq,0) + call scale_for_abl(abl_edge3,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & + bpasdq,0) end if + ! # modify flux below and above by B^- A^+ Delta q and B^+ A^+ Delta q: do i = 1, mx+1 do m = 1, num_eqn diff --git a/src/2d/setaux.f90 b/src/2d/setaux.f90 index 73ab084..ded13c0 100644 --- a/src/2d/setaux.f90 +++ b/src/2d/setaux.f90 @@ -1,19 +1,16 @@ subroutine setaux(mbc,mx,my,xlower,ylower,dx,dy,maux,aux) - ! Called at start of computation before calling qinit, and - ! when AMR is used, also called every time a new grid patch is created. - ! Use to set auxiliary arrays aux(1:maux, 1-mbc:mx+mbc, 1-mbc:my+mbc). - ! Note that ghost cell values may need to be set if the aux arrays - ! are used by the Riemann solver(s). - - use abl_module, only: set_scaling_factor - - implicit none - integer, intent(in) :: mbc,mx,my,maux - real(kind=8), intent(in) :: xlower,ylower,dx,dy - real(kind=8), intent(out) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) - - aux(:,:,:) = 0.d0 - call set_scaling_factor(mbc,mx,my,xlower,ylower,dx,dy,maux,aux) + ! Called at start of computation before calling qinit, and + ! when AMR is used, also called every time a new grid patch is created. + ! Use to set auxiliary arrays aux(1:maux, 1-mbc:mx+mbc, 1-mbc:my+mbc). + ! Note that ghost cell values may need to be set if the aux arrays + ! are used by the Riemann solver(s). + ! + ! This default version does nothing. + + implicit none + integer, intent(in) :: mbc,mx,my,maux + real(kind=8), intent(in) :: xlower,ylower,dx,dy + real(kind=8), intent(out) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) end subroutine setaux diff --git a/src/2d/step2.f90 b/src/2d/step2.f90 index c4cd3e1..20ceceb 100644 --- a/src/2d/step2.f90 +++ b/src/2d/step2.f90 @@ -1,6 +1,6 @@ ! ========================================================== subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & - qold,qnew,aux,dx,dy,dt,method,mthlim,cfl, & + qold,qnew,aux,xlower,ylower,dx,dy,dt,method,mthlim,cfl, & qadd,fadd,gadd,q1d,dtdx1d,dtdy1d, & aux1,aux2,aux3,work,mwork,use_fwave,rpn2,rpt2) ! ========================================================== @@ -15,6 +15,7 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & ! # fadd and gadd are used to return flux increments from flux2. ! # See the flux2 documentation for more information. + use abl_module, only: abl_type, set_scaling_factor implicit double precision (a-h,o-z) external :: rpn2,rpt2 @@ -38,6 +39,19 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & logical :: use_fwave double precision :: work(mwork) + ! added for abl implementation + integer :: i_abl_lower, i_abl_upper, j_abl_lower, j_abl_upper + double precision :: xlower, ylower + double precision :: abl_center(1-num_ghost:mx+num_ghost, & + 1-num_ghost:my+num_ghost, 2) + double precision :: abl_edge(1-num_ghost:mx+num_ghost, & + 1-num_ghost:my+num_ghost, 2) + double precision :: abl_center2(1-num_ghost:maxm+num_ghost) + double precision :: abl_edge1(1-num_ghost:maxm+num_ghost) + double precision :: abl_edge2(1-num_ghost:maxm+num_ghost) + double precision :: abl_edge3(1-num_ghost:maxm+num_ghost) + + common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom !f2py intent(out) cfl @@ -85,25 +99,38 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & end do endif +! compute abl scaling factors if needed + if (abl_type /= 0) then + call set_scaling_factor(xlower,ylower,dx,dy,mx,my,num_ghost, & + abl_center,abl_edge, & + i_abl_lower,i_abl_upper,j_abl_lower,j_abl_upper) + else + i_abl_lower = -num_ghost + j_abl_lower = -num_ghost + i_abl_upper = maxm+num_ghost + j_abl_upper = maxm+num_ghost + abl_center(:,:,:) = 1.d0 + end if + ! # perform x-sweeps ! ================== do j = 0,my+1 - + ! # copy data along a slice into 1d arrays: do i = 1-num_ghost, mx+num_ghost do m=1,num_eqn q1d(m,i) = qold(m,i,j) end do end do - + if (index_capa > 0) then do i = 1-num_ghost, mx+num_ghost dtdx1d(i) = dtdx/aux(index_capa,i,j) end do endif - + if (num_aux > 0) then do i = 1-num_ghost, mx+num_ghost do ma=1,num_aux @@ -113,52 +140,65 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & end do end do endif - + + if (abl_type /= 0) then + do i = 1-num_ghost, mx+num_ghost + abl_center2(i) = abl_center(i,j,1) + abl_edge1(i) = abl_edge(i,j-1,2) + abl_edge2(i) = abl_edge(i,j,1) + abl_edge3(i) = abl_edge(i,j+1,2) + end do + end if + ! # Store the value of j along this slice in the common block ! # comxyt in case it is needed in the Riemann solver (for ! # variable coefficient problems) jcom = j - + ! # compute modifications fadd and gadd to fluxes along this slice: call flux2(1,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & q1d,dtdx1d,aux1,aux2,aux3,method,mthlim, & qadd,fadd,gadd,cfl1d, & work(i0wave),work(i0s),work(i0amdq),work(i0apdq), & work(i0cqxx),work(i0bmadq),work(i0bpadq),rpn2,rpt2, & - use_fwave) + use_fwave,abl_center2,abl_edge1,abl_edge2,abl_edge3, & + i_abl_lower,i_abl_upper) cfl = dmax1(cfl,cfl1d) - + ! # update qnew by flux differencing. ! # (rather than maintaining arrays f and g for the total fluxes, ! # the modifications are used immediately to update qnew ! # in order to save storage.) - + if (index_capa == 0) then - + ! # no capa array. Standard flux differencing: do i=1,mx do m=1,num_eqn qnew(m,i,j) = qnew(m,i,j) + qadd(m,i) & - - dtdx * (fadd(m,i+1) - fadd(m,i)) & - - dtdy * (gadd(m,2,i) - gadd(m,1,i)) - qnew(m,i,j-1) = qnew(m,i,j-1) - dtdy * gadd(m,1,i) - qnew(m,i,j+1) = qnew(m,i,j+1) + dtdy * gadd(m,2,i) + - dtdx * abl_center(i,j,1) * (fadd(m,i+1) - fadd(m,i)) & + - dtdy * abl_center(i,j,2) * (gadd(m,2,i) - gadd(m,1,i)) + qnew(m,i,j-1) = qnew(m,i,j-1) & + - dtdy * abl_center(i,j-1,2) * gadd(m,1,i) + qnew(m,i,j+1) = qnew(m,i,j+1) & + + dtdy * abl_center(i,j+1,2) * gadd(m,2,i) end do end do - else - + ! # with capa array. do i=1,mx do m=1,num_eqn qnew(m,i,j) = qnew(m,i,j) + qadd(m,i) & - - (dtdx * (fadd(m,i+1) - fadd(m,i)) & - + dtdy * (gadd(m,2,i) - gadd(m,1,i))) & + - (dtdx * abl_center(i,j,1) * (fadd(m,i+1) - fadd(m,i)) & + + dtdy * abl_center(i,j,2) * (gadd(m,2,i) - gadd(m,1,i))) & / aux(index_capa,i,j) - qnew(m,i,j-1) = qnew(m,i,j-1) - dtdy * gadd(m,1,i) & - / aux(index_capa,i,j-1) - qnew(m,i,j+1) = qnew(m,i,j+1) + dtdy * gadd(m,2,i) & - / aux(index_capa,i,j+1) + qnew(m,i,j-1) = qnew(m,i,j-1) & + - dtdy * abl_center(i,j-1,2) * gadd(m,1,i) & + / aux(index_capa,i,j-1) + qnew(m,i,j+1) = qnew(m,i,j+1) & + + dtdy * abl_center(i,j+1,2) * gadd(m,2,i) & + / aux(index_capa,i,j+1) end do end do endif @@ -171,20 +211,20 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & do i = 0, mx+1 - + ! # copy data along a slice into 1d arrays: do j = 1-num_ghost, my+num_ghost do m=1,num_eqn q1d(m,j) = qold(m,i,j) end do end do - + if (index_capa > 0) then do j = 1-num_ghost, my+num_ghost dtdy1d(j) = dtdy / aux(index_capa,i,j) end do endif - + if (num_aux > 0) then do j = 1-num_ghost, my+num_ghost do ma=1,num_aux @@ -194,54 +234,67 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & end do end do endif - - + + if (abl_type /= 0) then + do j = 1-num_ghost, my+num_ghost + abl_center2(j) = abl_center(i,j,2) + abl_edge1(j) = abl_edge(i-1,j,1) + abl_edge2(j) = abl_edge(i,j,2) + abl_edge3(j) = abl_edge(i+1,j,1) + end do + end if + ! # Store the value of i along this slice in the common block ! # comxyt in case it is needed in the Riemann solver (for ! # variable coefficient problems) icom = i - + ! # compute modifications fadd and gadd to fluxes along this slice: call flux2(2,maxm,num_eqn,num_waves,num_aux,num_ghost,my, & q1d,dtdy1d,aux1,aux2,aux3,method,mthlim, & qadd,fadd,gadd,cfl1d, & work(i0wave),work(i0s),work(i0amdq),work(i0apdq), & work(i0cqxx),work(i0bmadq),work(i0bpadq),rpn2,rpt2, & - use_fwave) - + use_fwave,abl_center2,abl_edge1,abl_edge2,abl_edge3, & + j_abl_lower,j_abl_upper) + cfl = dmax1(cfl,cfl1d) - + ! # update qnew by flux differencing. ! # Note that the roles of fadd and gadd are reversed for ! # the y-sweeps -- fadd is the modification to g-fluxes and ! # gadd is the modification to f-fluxes to the left and right. - + if (index_capa == 0) then - + ! # no capa array. Standard flux differencing: do j=1,my do m=1,num_eqn - qnew(m,i,j) = qnew(m,i,j) + (qadd(m,j) & - - dtdy * (fadd(m,j+1) - fadd(m,j)) & - - dtdx * (gadd(m,2,j) - gadd(m,1,j))) - qnew(m,i-1,j) = qnew(m,i-1,j) - dtdx * gadd(m,1,j) - qnew(m,i+1,j) = qnew(m,i+1,j) + dtdx * gadd(m,2,j) + qnew(m,i,j) = qnew(m,i,j) + qadd(m,j) & + - dtdy * abl_center(i,j,2) * (fadd(m,j+1) - fadd(m,j)) & + - dtdx * abl_center(i,j,1) * (gadd(m,2,j) - gadd(m,1,j)) + qnew(m,i-1,j) = qnew(m,i-1,j) & + - dtdx * abl_center(i-1,j,1) * gadd(m,1,j) + qnew(m,i+1,j) = qnew(m,i+1,j) & + + dtdx * abl_center(i+1,j,1) * gadd(m,2,j) end do end do - + else - + ! # with capa array. do j=1,my do m=1,num_eqn qnew(m,i,j) = qnew(m,i,j) + qadd(m,j) & - - (dtdy * (fadd(m,j+1) - fadd(m,j)) & - + dtdx * (gadd(m,2,j) - gadd(m,1,j))) & + - (dtdy * abl_center(i,j,2) * (fadd(m,j+1) - fadd(m,j)) & + + dtdx * abl_center(i,j,1) * (gadd(m,2,j) - gadd(m,1,j))) & / aux(index_capa,i,j) - qnew(m,i-1,j) = qnew(m,i-1,j) - dtdx * gadd(m,1,j) & - / aux(index_capa,i-1,j) - qnew(m,i+1,j) = qnew(m,i+1,j) + dtdx * gadd(m,2,j) & - / aux(index_capa,i+1,j) + qnew(m,i-1,j) = qnew(m,i-1,j) & + - dtdx * abl_center(i-1,j,1) * gadd(m,1,j) & + / aux(index_capa,i-1,j) + qnew(m,i+1,j) = qnew(m,i+1,j) & + + dtdx * abl_Center(i+1,j,1) * gadd(m,2,j) & + / aux(index_capa,i+1,j) end do end do endif From 1d764719e112dfa7fc27f98b14783a318ec60408 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Tue, 5 Mar 2019 19:51:16 -0800 Subject: [PATCH 09/10] minor cleanup of step2 for abl --- src/2d/step2.f90 | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/2d/step2.f90 b/src/2d/step2.f90 index 20ceceb..a784ff5 100644 --- a/src/2d/step2.f90 +++ b/src/2d/step2.f90 @@ -101,14 +101,10 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & ! compute abl scaling factors if needed if (abl_type /= 0) then - call set_scaling_factor(xlower,ylower,dx,dy,mx,my,num_ghost, & - abl_center,abl_edge, & - i_abl_lower,i_abl_upper,j_abl_lower,j_abl_upper) + call set_scaling_factor(xlower,ylower,dx,dy,mx,my,num_ghost, & + abl_center,abl_edge, & + i_abl_lower,i_abl_upper,j_abl_lower,j_abl_upper) else - i_abl_lower = -num_ghost - j_abl_lower = -num_ghost - i_abl_upper = maxm+num_ghost - j_abl_upper = maxm+num_ghost abl_center(:,:,:) = 1.d0 end if From 95dc031f8ee96e62f7e96da873a99b0500292830 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Tue, 23 Apr 2019 00:58:03 -0700 Subject: [PATCH 10/10] many improvements to ABL implementation, but need to look in abl_R or remove it --- src/2d/abl_module.f90 | 172 +++++++++++++++++++++++++--------- src/2d/flux2.f90 | 77 ++++++---------- src/2d/step2.f90 | 209 ++++++++++++++++++++++++++++-------------- 3 files changed, 295 insertions(+), 163 deletions(-) diff --git a/src/2d/abl_module.f90 b/src/2d/abl_module.f90 index 49b436d..304de6f 100644 --- a/src/2d/abl_module.f90 +++ b/src/2d/abl_module.f90 @@ -68,7 +68,7 @@ subroutine initialize() end subroutine initialize subroutine set_scaling_factor(xlower, ylower, dx, dy, mx, my, num_ghost, & - abl_center, abl_edge, & + ablX_M, ablX_J, ablX_R, ablY_M, ablY_J, ablY_R, & i_abl_lower, i_abl_upper, j_abl_lower, j_abl_upper) implicit none @@ -76,10 +76,12 @@ subroutine set_scaling_factor(xlower, ylower, dx, dy, mx, my, num_ghost, & integer, intent(in) :: mx, my, num_ghost real(kind=8), intent(in) :: xlower, ylower, dx, dy integer, intent(out) :: i_abl_lower, i_abl_upper, j_abl_lower, j_abl_upper - real(kind=8), intent(out) :: abl_center(1-num_ghost:mx+num_ghost, & - 1-num_ghost:my+num_ghost, 2) - real(kind=8), intent(out) :: abl_edge(1-num_ghost:mx+num_ghost, & - 1-num_ghost:my+num_ghost, 2) + real(kind=8), intent(out) :: ablX_M(1-num_ghost:mx+num_ghost) + real(kind=8), intent(out) :: ablX_J(1-num_ghost:mx+num_ghost) + real(kind=8), intent(out) :: ablX_R(1-num_ghost:mx+num_ghost) + real(kind=8), intent(out) :: ablY_M(1-num_ghost:my+num_ghost) + real(kind=8), intent(out) :: ablY_J(1-num_ghost:my+num_ghost) + real(kind=8), intent(out) :: ablY_R(1-num_ghost:my+num_ghost) integer :: i, j real(kind=8) :: xcenter, ycenter, xedge, yedge @@ -89,59 +91,113 @@ subroutine set_scaling_factor(xlower, ylower, dx, dy, mx, my, num_ghost, & call initialize() end if - ! determine which indices are in the layer + ! x direction + + ! determine which indices are in the layer and loop over those i_abl_lower = int((1.d0 + 1.d-14)*(xpos(1)-xlower)/dx) - i_abl_upper = int((1.d0 + 1.d-14)*(xpos(2)-xlower)/dx) - j_abl_lower = int((1.d0 + 1.d-14)*(ypos(1)-ylower)/dy) - j_abl_upper = int((1.d0 + 1.d-14)*(ypos(2)-ylower)/dy) + i_abl_upper = 1 + int((1.d0 + 1.d-14)*(xpos(2)-xlower)/dx) - ! Loop over all cells - do j=1-num_ghost,my+num_ghost - ycenter = ylower + (j-0.5d0)*dy + do i=1-num_ghost,i_abl_lower + + xcenter = xlower + (i-0.5d0)*dx + + ! Calculate inverse of Jacobian at cell edges + xedge = max(xcenter-0.5d0*dx, xpos(1)-xdepth(1)) + xedge = min(xedge, xpos(2)+xdepth(2)) + ablX_J(i) = 1.d0/(1.d0 + gp(xedge,1)) + + ! Calculate grid map scaling and ratio at cell centers + xcenter = max(xcenter, xpos(1)-xdepth(1)+0.5d0*dx) + xcenter = min(xcenter, xpos(2)+xdepth(2)-0.5d0*dx) + if (abl_type == TRIG) then + xedge = min(xedge, xpos(2)+xdepth(2)-dx) + ablX_M(i) = 1.d0/(1.d0 + (g(xedge+dx,1)-g(xedge,1))/dx) + ablX_R(i) = ablX_M(i)*(1.d0 + gp(xcenter,1)) + else + ablX_M(i) = 1.d0/(1.d0 + gp(xcenter,1)) + ablX_R(i) = 1.d0 + end if - do i=1-num_ghost,mx+num_ghost - xcenter = xlower + (i-0.5d0)*dx + end do - ! determine lower cell edge locations - ! constrained to the computational domain without ghost cells - xedge = max(xcenter-0.5d0*dx,xpos(1)-xdepth(1)) - xedge = min(xedge,xpos(2)+xdepth(2)-dx) - yedge = max(ycenter-0.5d0*dy,ypos(1)-ydepth(1)) - yedge = min(yedge,ypos(2)+ydepth(2)-dy) + do i=i_abl_upper,mx+num_ghost - ! Calculate inverse of Jacobian at cell edges - abl_edge(i,j,1) = 1.d0/(1.d0 + gp(xedge,1)) - abl_edge(i,j,2) = 1.d0/(1.d0 + gp(yedge,2)) + xcenter = xlower + (i-0.5d0)*dx - ! Calculate cell width ratio at cell centers - ! (use inverse of Jacobian for methods that dont have g) - abl_center(i,j,1) = 1.d0/(1.d0 + gp(xcenter,1)) - abl_center(i,j,2) = 1.d0/(1.d0 + gp(ycenter,2)) + ! Calculate inverse of Jacobian at cell edges + xedge = max(xcenter-0.5d0*dx, xpos(1)-xdepth(1)) + xedge = min(xedge, xpos(2)+xdepth(2)) + ablX_J(i) = 1.d0/(1.d0 + gp(xedge,1)) - end do + ! Calculate grid map scaling and ratio at cell centers + xcenter = max(xcenter, xpos(1)-xdepth(1)+0.5d0*dx) + xcenter = min(xcenter, xpos(2)+xdepth(2)-0.5d0*dx) + if (abl_type == TRIG) then + xedge = min(xedge, xpos(2)+xdepth(2)-dx) + ablX_M(i) = 1.d0/(1.d0 + (g(xedge+dx,1)-g(xedge,1))/dx) + ablX_R(i) = ablX_M(i)*(1.d0 + gp(xcenter,1)) + else + ablX_M(i) = 1.d0/(1.d0 + gp(xcenter,1)) + ablX_R(i) = 1.d0 + end if end do - end subroutine set_scaling_factor - subroutine scale_for_abl(abl_factor,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & - scaled_value,offset) + ! y direction - implicit none - integer, intent(in) :: num_ghost, maxm, mx, i_abl_lower, i_abl_upper, offset - double precision, intent(in) :: abl_factor(1-num_ghost:maxm+num_ghost) - double precision, intent(inout) :: scaled_value(1-num_ghost:maxm+num_ghost) + ! determine which indices are in the layer and loop over those + j_abl_lower = int((1.d0 + 1.d-14)*(ypos(1)-ylower)/dy) + j_abl_upper = 1 + int((1.d0 + 1.d-14)*(ypos(2)-ylower)/dy) + + do j=1-num_ghost,j_abl_lower - integer :: i + ycenter = ylower + (j-0.5d0)*dy + + ! Calculate inverse of Jacobian at cell edges + yedge = max(ycenter-0.5d0*dy, ypos(1)-ydepth(1)) + yedge = min(yedge, ypos(2)+ydepth(2)) + ablY_J(j) = 1.d0/(1.d0 + gp(yedge,2)) + + ! Calculate grid map scaling and ratio at cell centers + ycenter = max(ycenter, ypos(1)-ydepth(1)+0.5d0*dy) + ycenter = min(ycenter, ypos(2)+ydepth(2)-0.5d0*dy) + if (abl_type == TRIG) then + yedge = min(yedge, ypos(2)+ydepth(2)-dy) + ablY_M(j) = 1.d0/(1.d0 + (g(yedge+dy,2)-g(yedge,2))/dy) + ablY_R(j) = ablY_M(j)*(1.d0 + gp(ycenter,2)) + else + ablY_M(j) = 1.d0/(1.d0 + gp(ycenter,2)) + ablY_R(j) = 1.d0 + end if - do i = 1-num_ghost, i_abl_lower - scaled_value(i) = abl_factor(i-offset)*scaled_value(i) end do - do i = i_abl_upper, mx+num_ghost - scaled_value(i) = abl_factor(i-offset)*scaled_value(i) + + do j=j_abl_upper,my+num_ghost + + ycenter = ylower + (j-0.5d0)*dy + + ! Calculate inverse of Jacobian at cell edges + yedge = max(ycenter-0.5d0*dy, ypos(1)-ydepth(1)) + yedge = min(yedge, ypos(2)+ydepth(2)) + ablY_J(j) = 1.d0/(1.d0 + gp(yedge,2)) + + ! Calculate grid map scaling and ratio at cell centers + ycenter = max(ycenter, ypos(1)-ydepth(1)+0.5d0*dy) + ycenter = min(ycenter, ypos(2)+ydepth(2)-0.5d0*dy) + if (abl_type == TRIG) then + yedge = min(yedge, ypos(2)+ydepth(2)-dy) + ablY_M(j) = 1.d0/(1.d0 + (g(yedge+dy,2)-g(yedge,2))/dy) + ablY_R(j) = ablY_M(j)*(1.d0 + gp(ycenter,2)) + else + ablY_M(j) = 1.d0/(1.d0 + gp(ycenter,2)) + ablY_R(j) = 1.d0 + end if + end do - end subroutine scale_for_abl + + end subroutine set_scaling_factor function gp(z,direction) @@ -183,6 +239,40 @@ function gp(z,direction) end function gp + function g(z,direction) + + implicit none + + real(kind=8) :: g, z + integer :: direction + + real(kind=8) :: pos(2), depth(2), val, depth_scalar + + g = 0.d0 + if (direction == 1) then + pos = xpos + depth = xdepth + else if (direction == 2) then + pos = ypos + depth = ydepth + end if + + if (z < pos(1)) then + val = (z-pos(1))/depth(1) + depth_scalar = depth(1) + else if (z > pos(2)) then + val = (z-pos(2))/depth(2) + depth_scalar = depth(2) + else + return + end if + + if (abl_type == TRIG) then + g = depth_scalar*(dtan(PIH*val)/PIH - val) + end if + + end function g + end module abl_module diff --git a/src/2d/flux2.f90 b/src/2d/flux2.f90 index 6d83345..10c6a6a 100644 --- a/src/2d/flux2.f90 +++ b/src/2d/flux2.f90 @@ -5,8 +5,7 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & q1d,dtdx1d,aux1,aux2,aux3,method,mthlim, & qadd,fadd,gadd,cfl1d,wave,s, & - amdq,apdq,cqxx,bmasdq,bpasdq,rpn2,rpt2,use_fwave, & - abl_center,abl_edge1,abl_edge2,abl_edge3,i_abl_lower,i_abl_upper) + amdq,apdq,cqxx,bmasdq,bpasdq,rpn2,rpt2,use_fwave,abl,i_abl_lower,i_abl_upper) ! ===================================================== ! # Compute the modification to fluxes f and g that are generated by @@ -58,7 +57,7 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & !f2py external rpn2, rpt2 !f2py intent(callback) rpn2,rpt2 - use abl_module, only: abl_type, scale_for_abl + use abl_module, only: abl_type implicit double precision (a-h,o-z) integer num_aux @@ -74,8 +73,6 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & dimension gadd(num_eqn, 2, 1-num_ghost:maxm+num_ghost) dimension dtdx1d(1-num_ghost:maxm+num_ghost) - dimension dtdx1d_abl(1-num_ghost:maxm+num_ghost) - dimension dtdxave_abl(1-num_ghost:maxm+num_ghost) dimension aux1(num_aux,1-num_ghost:maxm+num_ghost) dimension aux2(num_aux,1-num_ghost:maxm+num_ghost) dimension aux3(num_aux,1-num_ghost:maxm+num_ghost) @@ -83,10 +80,8 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & dimension s(num_waves,1-num_ghost:maxm+num_ghost) dimension wave(num_eqn, num_waves, 1-num_ghost:maxm+num_ghost) - dimension abl_center(1-num_ghost:maxm+num_ghost) - dimension abl_edge1(1-num_ghost:maxm+num_ghost) - dimension abl_edge2(1-num_ghost:maxm+num_ghost) - dimension abl_edge3(1-num_ghost:maxm+num_ghost) + dimension abl(1-num_ghost:maxm+num_ghost) + dimension dtdxave(1-num_ghost:mx+num_ghost) dimension method(7),mthlim(num_waves) logical :: limit, use_fwave @@ -108,6 +103,22 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & end do end do + ! Store dtdxave for second order corrections so that it may be scaled for + ! ABL if needed + do i = 2-num_ghost, mx+num_ghost + dtdxave(i) = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) + end do + + ! apply ABL if needed + if (abl_type /= 0) then + do i = 2-num_ghost, i_abl_lower + dtdxave(i) = dtdxave(i) * abl(i) + end do + do i = i_abl_upper, mx+num_ghost + dtdxave(i) = dtdxave(i) * abl(i) + end do + end if + ! Unnecessary to init fadd to zero here because it only gets ! modified by the second-order corrections @@ -120,23 +131,13 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & call rpn2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, & aux2,aux2,wave,s,amdq,apdq) -! adjust dtdx1d for abl if needed - do i = 0, mx+1 - dtdx1d_abl(i) = dtdx1d(i) - end do - if (abl_type /= 0) then - call scale_for_abl(abl_center,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & - dtdx1d_abl,0) - end if - - ! # Set qadd for the donor-cell upwind method (Godunov) do i = 1, mx+1 do m = 1, num_eqn ! qadd(:,i-1) is still in cache from last cycle of outer loop - qadd(m,i-1) = qadd(m,i-1) - dtdx1d_abl(i-1)*amdq(m,i) + qadd(m,i-1) = qadd(m,i-1) - amdq(m,i) end do do m = 1, num_eqn - qadd(m,i) = qadd(m,i) - dtdx1d_abl(i)*apdq(m,i) + qadd(m,i) = qadd(m,i) - apdq(m,i) end do end do @@ -165,23 +166,15 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & ! # Second order corrections: if (use_fwave.eqv. .FALSE. ) then - ! # modified in Version 4.3 to use average only in cqxx, not transverse - ! adjust dtdxave for abl (if ndeed) - do i = 2-num_ghost, mx+num_ghost - dtdxave_abl(i) = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) - end do - if (abl_type /= 0) then - call scale_for_abl(abl_edge2,num_ghost,maxm,mx,i_abl_lower, & - i_abl_upper,dtdxave_abl,0) - end if do i = 2-num_ghost, mx+num_ghost + ! # modified in Version 4.3 to use average only in cqxx, not transverse do m = 1, num_eqn cqxx(m,i) = 0.d0 end do do mw=1,num_waves ! Traverse the wave array in memory-contiguous fashion do m=1,num_eqn cqxx(m,i) = cqxx(m,i) + dabs(s(mw,i)) & - * (1.d0 - dabs(s(mw,i))*dtdxave_abl(i)) * wave(m,mw,i) + * (1.d0 - dabs(s(mw,i))*dtdxave(i)) * wave(m,mw,i) enddo enddo do m = 1, num_eqn @@ -190,11 +183,10 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & enddo else if (abl_type /= 0) then - print *, "ABL not yet implemented for fwave splitting" + print *, "ABL not implemented yet for F-wave splitting" stop end if do i = 2-num_ghost, mx+num_ghost - dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) do m = 1, num_eqn cqxx(m,i) = 0.d0 end do @@ -202,7 +194,7 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & do m=1,num_eqn ! # second order corrections: cqxx(m,i) = cqxx(m,i) + dsign(1.d0,s(mw,i)) & - * (1.d0 - dabs(s(mw,i))*dtdxave) * wave(m,mw,i) + * (1.d0 - dabs(s(mw,i))*dtdxave(i)) * wave(m,mw,i) enddo enddo do m = 1, num_eqn @@ -235,14 +227,6 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & call rpt2(ixy,1,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, & aux1,aux2,aux3,amdq,bmasdq,bpasdq) -! adjust transverse fluxes for abl if needed - if (abl_type /= 0) then - call scale_for_abl(abl_edge1,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & - bmasdq,1) - call scale_for_abl(abl_edge3,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & - bpasdq,1) - end if - ! # modify flux below and above by B^- A^- Delta q and B^+ A^- Delta q: do i = 1, mx+1 ! Having two inner loops here allows traversal of gadd in memory-contiguous order @@ -259,15 +243,6 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & call rpt2(ixy,2,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, & aux1,aux2,aux3,apdq,bmasdq,bpasdq) -! adjust transverse fluxes for abl if needed - if (abl_type /= 0) then - call scale_for_abl(abl_edge1,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & - bmasdq,0) - call scale_for_abl(abl_edge3,num_ghost,maxm,mx,i_abl_lower,i_abl_upper, & - bpasdq,0) - end if - - ! # modify flux below and above by B^- A^+ Delta q and B^+ A^+ Delta q: do i = 1, mx+1 do m = 1, num_eqn diff --git a/src/2d/step2.f90 b/src/2d/step2.f90 index a784ff5..ce44ca9 100644 --- a/src/2d/step2.f90 +++ b/src/2d/step2.f90 @@ -41,16 +41,13 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & ! added for abl implementation integer :: i_abl_lower, i_abl_upper, j_abl_lower, j_abl_upper - double precision :: xlower, ylower - double precision :: abl_center(1-num_ghost:mx+num_ghost, & - 1-num_ghost:my+num_ghost, 2) - double precision :: abl_edge(1-num_ghost:mx+num_ghost, & - 1-num_ghost:my+num_ghost, 2) - double precision :: abl_center2(1-num_ghost:maxm+num_ghost) - double precision :: abl_edge1(1-num_ghost:maxm+num_ghost) - double precision :: abl_edge2(1-num_ghost:maxm+num_ghost) - double precision :: abl_edge3(1-num_ghost:maxm+num_ghost) - + double precision :: xlower, ylower, dtdxm, dtdxc, dtdxp, dtdyp, dtdyc, dtdym + double precision :: ablX_M(1-num_ghost:mx+num_ghost) + double precision :: ablX_J(1-num_ghost:mx+num_ghost) + double precision :: ablX_R(1-num_ghost:mx+num_ghost) + double precision :: ablY_M(1-num_ghost:my+num_ghost) + double precision :: ablY_J(1-num_ghost:my+num_ghost) + double precision :: ablY_R(1-num_ghost:my+num_ghost) common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom @@ -90,6 +87,14 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & cfl = 0.d0 dtdx = dt/dx dtdy = dt/dy + dtdxm = dtdx + dtdxc = dtdx + dtdxp = dtdx + dtdym = dtdy + dtdyc = dtdy + dtdyp = dtdy + + if (index_capa == 0) then ! # no capa array: @@ -102,13 +107,15 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & ! compute abl scaling factors if needed if (abl_type /= 0) then call set_scaling_factor(xlower,ylower,dx,dy,mx,my,num_ghost, & - abl_center,abl_edge, & + ablX_M,ablX_J,ablX_R,ablY_M,ablY_J,ablY_R, & i_abl_lower,i_abl_upper,j_abl_lower,j_abl_upper) else - abl_center(:,:,:) = 1.d0 + i_abl_lower = 0 + j_abl_lower = 0 + i_abl_upper = mx+1 + j_abl_upper = my+1 end if - ! # perform x-sweeps ! ================== @@ -137,28 +144,39 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & end do endif - if (abl_type /= 0) then - do i = 1-num_ghost, mx+num_ghost - abl_center2(i) = abl_center(i,j,1) - abl_edge1(i) = abl_edge(i,j-1,2) - abl_edge2(i) = abl_edge(i,j,1) - abl_edge3(i) = abl_edge(i,j+1,2) - end do - end if - ! # Store the value of j along this slice in the common block ! # comxyt in case it is needed in the Riemann solver (for ! # variable coefficient problems) jcom = j + ! adjust dtdy* if ABL is present + if (abl_type /= 0) then + if (j-1 <= j_abl_lower .or. j-1 >= j_abl_upper) then + dtdym = dtdy * ablY_R(j-1) * ablY_M(j-1) + else + dtdym = dtdy + end if + if (j <= j_abl_lower .or. j >= j_abl_upper) then + dtdyc = dtdy * ablY_R(j) * ablY_M(j) + dtdxc = dtdx * ablY_R(j) + else + dtdyc = dtdy + dtdxc = dtdx + end if + if (j+1 <= j_abl_lower .or. j+1 >= j_abl_upper) then + dtdyp = dtdy * ablY_R(j+1) * ablY_M(j+1) + else + dtdyp = dtdy + end if + end if + ! # compute modifications fadd and gadd to fluxes along this slice: call flux2(1,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, & q1d,dtdx1d,aux1,aux2,aux3,method,mthlim, & qadd,fadd,gadd,cfl1d, & work(i0wave),work(i0s),work(i0amdq),work(i0apdq), & work(i0cqxx),work(i0bmadq),work(i0bpadq),rpn2,rpt2, & - use_fwave,abl_center2,abl_edge1,abl_edge2,abl_edge3, & - i_abl_lower,i_abl_upper) + use_fwave,ablX_J,i_abl_lower,i_abl_upper) cfl = dmax1(cfl,cfl1d) ! # update qnew by flux differencing. @@ -169,32 +187,51 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & if (index_capa == 0) then ! # no capa array. Standard flux differencing: - do i=1,mx + do i=i_abl_lower+1,i_abl_upper-1 do m=1,num_eqn - qnew(m,i,j) = qnew(m,i,j) + qadd(m,i) & - - dtdx * abl_center(i,j,1) * (fadd(m,i+1) - fadd(m,i)) & - - dtdy * abl_center(i,j,2) * (gadd(m,2,i) - gadd(m,1,i)) - qnew(m,i,j-1) = qnew(m,i,j-1) & - - dtdy * abl_center(i,j-1,2) * gadd(m,1,i) - qnew(m,i,j+1) = qnew(m,i,j+1) & - + dtdy * abl_center(i,j+1,2) * gadd(m,2,i) + qnew(m,i,j) = qnew(m,i,j) + dtdxc * qadd(m,i) & + - dtdxc * (fadd(m,i+1) - fadd(m,i)) & + - dtdyc * (gadd(m,2,i) - gadd(m,1,i)) + qnew(m,i,j-1) = qnew(m,i,j-1) - dtdym * gadd(m,1,i) + qnew(m,i,j+1) = qnew(m,i,j+1) + dtdyp * gadd(m,2,i) + end do + end do + ! these loops are skipped if no ABL present + do i=1,i_abl_lower + do m=1,num_eqn + qnew(m,i,j) = qnew(m,i,j) + ablX_M(i)*( dtdxc * qadd(m,i) & + - dtdxc * (fadd(m,i+1) - fadd(m,i)) & + - dtdyc * (gadd(m,2,i) - gadd(m,1,i)) ) + qnew(m,i,j-1) = qnew(m,i,j-1) - dtdym*ablX_M(i) * gadd(m,1,i) + qnew(m,i,j+1) = qnew(m,i,j+1) + dtdyp*ablX_M(i) * gadd(m,2,i) + end do + end do + do i=i_abl_upper,mx + do m=1,num_eqn + qnew(m,i,j) = qnew(m,i,j) + ablX_M(i)*( dtdxc * qadd(m,i) & + - dtdxc * (fadd(m,i+1) - fadd(m,i)) & + - dtdyc * (gadd(m,2,i) - gadd(m,1,i)) ) + qnew(m,i,j-1) = qnew(m,i,j-1) - dtdym*ablX_M(i) * gadd(m,1,i) + qnew(m,i,j+1) = qnew(m,i,j+1) + dtdyp*ablX_M(i) * gadd(m,2,i) end do end do - else + else ! # with capa array. + if (abl_type /= 0) then + print *, "ABL not implemented yet for capa array" + stop + end if do i=1,mx do m=1,num_eqn - qnew(m,i,j) = qnew(m,i,j) + qadd(m,i) & - - (dtdx * abl_center(i,j,1) * (fadd(m,i+1) - fadd(m,i)) & - + dtdy * abl_center(i,j,2) * (gadd(m,2,i) - gadd(m,1,i))) & + qnew(m,i,j) = qnew(m,i,j) + dtdx * qadd(m,i) & + - (dtdx * (fadd(m,i+1) - fadd(m,i)) & + + dtdy * (gadd(m,2,i) - gadd(m,1,i))) & / aux(index_capa,i,j) - qnew(m,i,j-1) = qnew(m,i,j-1) & - - dtdy * abl_center(i,j-1,2) * gadd(m,1,i) & - / aux(index_capa,i,j-1) - qnew(m,i,j+1) = qnew(m,i,j+1) & - + dtdy * abl_center(i,j+1,2) * gadd(m,2,i) & - / aux(index_capa,i,j+1) + qnew(m,i,j-1) = qnew(m,i,j-1) - dtdy * gadd(m,1,i) & + / aux(index_capa,i,j-1) + qnew(m,i,j+1) = qnew(m,i,j+1) + dtdy * gadd(m,2,i) & + / aux(index_capa,i,j+1) end do end do endif @@ -205,7 +242,6 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & ! # perform y sweeps ! ================== - do i = 0, mx+1 ! # copy data along a slice into 1d arrays: @@ -231,28 +267,40 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & end do endif - if (abl_type /= 0) then - do j = 1-num_ghost, my+num_ghost - abl_center2(j) = abl_center(i,j,2) - abl_edge1(j) = abl_edge(i-1,j,1) - abl_edge2(j) = abl_edge(i,j,2) - abl_edge3(j) = abl_edge(i+1,j,1) - end do - end if - ! # Store the value of i along this slice in the common block ! # comxyt in case it is needed in the Riemann solver (for ! # variable coefficient problems) icom = i + ! adjust dtdx* if ABL is present + if (abl_type /= 0) then + if (i-1 <= i_abl_lower .or. i-1 >= i_abl_upper) then + dtdxm = dtdx * ablX_R(i-1) * ablX_M(i-1) + else + dtdxm = dtdx + end if + if (i <= i_abl_lower .or. i >= i_abl_upper) then + dtdxc = dtdx * ablX_R(i) * ablX_M(i) + dtdyc = dtdy * ablX_R(i) + else + dtdxc = dtdx + dtdyc = dtdy + end if + if (i+1 <= i_abl_lower .or. i+1 >= i_abl_upper) then + dtdxp = dtdx * ablX_R(i+1) * ablX_M(i+1) + else + dtdxp = dtdx + end if + end if + + ! # compute modifications fadd and gadd to fluxes along this slice: call flux2(2,maxm,num_eqn,num_waves,num_aux,num_ghost,my, & q1d,dtdy1d,aux1,aux2,aux3,method,mthlim, & qadd,fadd,gadd,cfl1d, & work(i0wave),work(i0s),work(i0amdq),work(i0apdq), & work(i0cqxx),work(i0bmadq),work(i0bpadq),rpn2,rpt2, & - use_fwave,abl_center2,abl_edge1,abl_edge2,abl_edge3, & - j_abl_lower,j_abl_upper) + use_fwave,ablY_J,j_abl_lower,j_abl_upper) cfl = dmax1(cfl,cfl1d) @@ -264,33 +312,52 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, & if (index_capa == 0) then ! # no capa array. Standard flux differencing: - do j=1,my + do j=j_abl_lower+1,j_abl_upper-1 + do m=1,num_eqn + qnew(m,i,j) = qnew(m,i,j) + dtdyc * qadd(m,j) & + - dtdyc * (fadd(m,j+1) - fadd(m,j)) & + - dtdxc * (gadd(m,2,j) - gadd(m,1,j)) + qnew(m,i-1,j) = qnew(m,i-1,j) - dtdxm * gadd(m,1,j) + qnew(m,i+1,j) = qnew(m,i+1,j) + dtdxp * gadd(m,2,j) + end do + end do + ! these loops are skipped if no ABL present + do j=1,j_abl_lower + do m=1,num_eqn + qnew(m,i,j) = qnew(m,i,j) + ablY_M(j)*( dtdyc * qadd(m,j) & + - dtdyc * (fadd(m,j+1) - fadd(m,j)) & + - dtdxc * (gadd(m,2,j) - gadd(m,1,j)) ) + qnew(m,i-1,j) = qnew(m,i-1,j) - dtdxm*ablY_M(j) * gadd(m,1,j) + qnew(m,i+1,j) = qnew(m,i+1,j) + dtdxp*ablY_M(j) * gadd(m,2,j) + end do + end do + do j=j_abl_upper,my do m=1,num_eqn - qnew(m,i,j) = qnew(m,i,j) + qadd(m,j) & - - dtdy * abl_center(i,j,2) * (fadd(m,j+1) - fadd(m,j)) & - - dtdx * abl_center(i,j,1) * (gadd(m,2,j) - gadd(m,1,j)) - qnew(m,i-1,j) = qnew(m,i-1,j) & - - dtdx * abl_center(i-1,j,1) * gadd(m,1,j) - qnew(m,i+1,j) = qnew(m,i+1,j) & - + dtdx * abl_center(i+1,j,1) * gadd(m,2,j) + qnew(m,i,j) = qnew(m,i,j) + ablY_M(j)*( dtdyc * qadd(m,j) & + - dtdyc * (fadd(m,j+1) - fadd(m,j)) & + - dtdxc * (gadd(m,2,j) - gadd(m,1,j)) ) + qnew(m,i-1,j) = qnew(m,i-1,j) - dtdxm*ablY_M(j) * gadd(m,1,j) + qnew(m,i+1,j) = qnew(m,i+1,j) + dtdxp*ablY_M(j) * gadd(m,2,j) end do end do else ! # with capa array. + if (abl_type /= 0) then + print *, "ABL not implemented yet for capa array" + stop + end if do j=1,my do m=1,num_eqn - qnew(m,i,j) = qnew(m,i,j) + qadd(m,j) & - - (dtdy * abl_center(i,j,2) * (fadd(m,j+1) - fadd(m,j)) & - + dtdx * abl_center(i,j,1) * (gadd(m,2,j) - gadd(m,1,j))) & + qnew(m,i,j) = qnew(m,i,j) + dtdy * qadd(m,j) & + - (dtdy * (fadd(m,j+1) - fadd(m,j)) & + + dtdx * (gadd(m,2,j) - gadd(m,1,j))) & / aux(index_capa,i,j) - qnew(m,i-1,j) = qnew(m,i-1,j) & - - dtdx * abl_center(i-1,j,1) * gadd(m,1,j) & - / aux(index_capa,i-1,j) - qnew(m,i+1,j) = qnew(m,i+1,j) & - + dtdx * abl_Center(i+1,j,1) * gadd(m,2,j) & - / aux(index_capa,i+1,j) + qnew(m,i-1,j) = qnew(m,i-1,j) - dtdx * gadd(m,1,j) & + / aux(index_capa,i-1,j) + qnew(m,i+1,j) = qnew(m,i+1,j) + dtdx * gadd(m,2,j) & + / aux(index_capa,i+1,j) end do end do endif