diff --git a/examples/acoustics_2d_radial/Makefile_abl b/examples/acoustics_2d_radial/Makefile_abl new file mode 100644 index 000000000..a37d340fd --- /dev/null +++ b/examples/acoustics_2d_radial/Makefile_abl @@ -0,0 +1,61 @@ + +# 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 = amrclaw # Clawpack package to use +EXE = xamr # 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: +# --------------------------------- + +AMRLIB = $(CLAW)/amrclaw/src/2d +include $(AMRLIB)/Makefile.amr_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.f \ + $(CLAW)/riemann/src/rpt2_acoustics.f + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + diff --git a/examples/acoustics_2d_radial/README.rst b/examples/acoustics_2d_radial/README.rst index 12208367f..49c027044 100644 --- a/examples/acoustics_2d_radial/README.rst +++ b/examples/acoustics_2d_radial/README.rst @@ -13,3 +13,17 @@ with the 1d solution. Three gauges are specified, Gauge 0 is at the origin and Gauges 1 and 2 are are both at radius 0.7, along the x-axis and along the diagonal. +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 + diff --git a/examples/acoustics_2d_radial/setplot_abl.py b/examples/acoustics_2d_radial/setplot_abl.py new file mode 100644 index 000000000..6b1cbc10c --- /dev/null +++ b/examples/acoustics_2d_radial/setplot_abl.py @@ -0,0 +1,175 @@ + +""" +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. + + """ + + + from clawpack.visclaw import colormaps + from clawpack.clawutil.data import ClawData + + if plotdata is None: + from clawpack.visclaw.data import ClawPlotData + plotdata = ClawPlotData() + + # determine original computational domain + abldata = ClawData() + abldata.read(plotdata.outdir + '/abl.data', force=True) + clawdata = ClawData() + clawdata.read(plotdata.outdir + '/claw.data', force=True) + 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 + + + # 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.add_colorbar = True + plotitem.show = True # show on plot? + plotitem.pcolor_cmin = -2.0 + plotitem.pcolor_cmax = 2.0 + plotitem.amr_patchedges_show = [1,1,1] + plotitem.amr_celledges_show = [1,0,0] + + def plot_original_domain_and_add_gauges(current_data): + 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') + addgauges(current_data) + + plotaxes.afteraxes = plot_original_domain_and_add_gauges + + # Figure for scatter plot + # ----------------------- + + plotfigure = plotdata.new_plotfigure(name='scatter', figno=3) + plotfigure.show = (qref_dir is not None) + + # 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 < 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 + + plotitem.map_2d_to_1d = p_vs_r + plotitem.plot_var = 0 + plotitem.plotstyle = 'o' + plotitem.color = 'b' + plotitem.show = True # 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? + plotaxes.afteraxes = "import pylab; pylab.legend(('2d data (interior only)', '1d reference solution'))" + + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='q', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Pressure' + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.plotstyle = 'b-' + + + # 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? + plotdata.parallel = True # make multiple frame png's at once + + return plotdata + + + +# To plot gauge locations on pcolor or contour plot, use this as +# an afteraxis function: + +def addgauges(current_data): + from clawpack.visclaw import gaugetools + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos='all', format_string='ko', add_labels=True) diff --git a/examples/acoustics_2d_radial/setrun_abl.py b/examples/acoustics_2d_radial/setrun_abl.py new file mode 100644 index 000000000..0753e447a --- /dev/null +++ b/examples/acoustics_2d_radial/setrun_abl.py @@ -0,0 +1,366 @@ +""" +Module to set up run time parameters for Clawpack. + +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='amrclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "amrclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + + assert claw_pkg.lower() == 'amrclaw', "Expected claw_pkg = 'amrclaw'" + + 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' + + clawdata.output_q_components = 'all' # could be list such as [True,True] + clawdata.output_aux_components = 'none' # could be list + clawdata.output_aux_onlyonce = True # output aux arrays only at t0 + + + # --------------------------------------------------- + # 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? + 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 + + + # --------------- + # Gauges: + # --------------- + rundata.gaugedata.gauges = [] + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + rundata.gaugedata.gauges.append([0, 0.0, 0.0, 0., 10.]) + rundata.gaugedata.gauges.append([1, 0.7, 0.0, 0., 10.]) + rundata.gaugedata.gauges.append([2, 0.7/np.sqrt(2.), 0.7/np.sqrt(2.), 0., 10.]) + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif clawdata.checkpt_style == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = [0.1,0.15] + + elif clawdata.checkpt_style == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + + # --------------- + # AMR parameters: + # --------------- + + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 3 + + # List of refinement ratios at each level (length at least amr_level_max-1) + amrdata.refinement_ratios_x = [2, 2] + amrdata.refinement_ratios_y = [2, 2] + amrdata.refinement_ratios_t = [2, 2] + + + # Specify type of each aux variable in clawdata.auxtype. + # This must be a list of length num_aux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + amrdata.aux_type = [] + + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = False # use Richardson? + amrdata.flag_richardson_tol = 0.001000e+00 # Richardson tolerance + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = True # use this? + amrdata.flag2refine_tol = 0.2 # tolerance used in this routine + # User can modify flag2refine to change the criterion for flagging. + # Default: check maximum absolute difference of first component of q + # between a cell and each of its neighbors. + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 2 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.7 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 0 + + + # --------------- + # Regions: + # --------------- + rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + 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.amr_2d b/src/2d/Makefile.amr_2d index 75b015816..45f001a0f 100644 --- a/src/2d/Makefile.amr_2d +++ b/src/2d/Makefile.amr_2d @@ -6,7 +6,8 @@ COMMON_MODULES += \ $(CLAW)/classic/src/utility_module.f90 \ $(AMRLIB)/amr_module.f90 \ $(AMRLIB)/gauges_module.f90 \ - $(AMRLIB)/regions_module.f90 + $(AMRLIB)/regions_module.f90 \ + $(AMRLIB)/abl_module.f90 COMMON_SOURCES += \ $(AMRLIB)/qinit.f \ diff --git a/src/2d/abl_module.f90 b/src/2d/abl_module.f90 new file mode 100644 index 000000000..4349ead36 --- /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.f b/src/2d/flux2.f index a8e873c2c..5d421ddea 100644 --- a/src/2d/flux2.f +++ b/src/2d/flux2.f @@ -76,6 +76,7 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, c c use amr_module + use abl_module, only: abl_type implicit double precision (a-h,o-z) external rpn2, rpt2 dimension q1d(meqn,1-mbc:maxm+mbc) @@ -124,6 +125,20 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, c call rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,q1d,q1d, & aux2,aux2,wave,s,amdq,apdq) + +c # apply ABL scaling to propagation speed and update fluctuations if needed + if (abl_type .ne. 0) then + do i = 1-mbc,mx+mbc + do m = 1,mwaves + s(m,i) = s(m,i)*aux2(maux-2+ixy,i) + end do + do m = 1,meqn + amdq(m,i) = amdq(m,i)*aux2(maux-2+ixy,i) + apdq(m,i) = apdq(m,i)*aux2(maux-2+ixy,i) + end do + end do + end if + c c # Set fadd for the donor-cell upwind method (Godunov) do 40 i=1,mx+1 @@ -202,6 +217,18 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, call rpt2(ixy,1,maxm,meqn,mwaves,maux,mbc,mx, & q1d,q1d,aux1,aux2,aux3, & amdq,bmasdq,bpasdq) + +c +c # Adjust for ABL if needed + if (abl_type .ne. 0) then + do i = 1,mx+mbc + do m = 1,meqn + bmasdq(m,i) = bmasdq(m,i)*aux2(maux+1-ixy,i-1) + bpasdq(m,i) = bpasdq(m,i)*aux3(maux+1-ixy,i-1) + end do + end do + end if + c c # modify flux below and above by B^- A^- Delta q and B^+ A^- Delta q: do 160 i = 1, mx+1 @@ -219,6 +246,16 @@ subroutine flux2(ixy,maxm,meqn,maux,mbc,mx, call rpt2(ixy,2,maxm,meqn,mwaves,maux,mbc,mx, & q1d,q1d,aux1,aux2,aux3, & apdq,bmasdq,bpasdq) + +c # Adjust for ABL if needed + if (abl_type .ne. 0) then + do i = 1,mx+mbc + do m = 1,meqn + bmasdq(m,i) = bmasdq(m,i)*aux2(maux+1-ixy,i) + bpasdq(m,i) = bpasdq(m,i)*aux3(maux+1-ixy,i) + end do + end do + end if c c # modify flux below and above by B^- A^+ Delta q and B^+ A^+ Delta q: do 180 i = 1, mx+1 diff --git a/src/2d/setaux.f90 b/src/2d/setaux.f90 index 9dfa290f6..73ab08470 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(inout) :: 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 diff --git a/src/python/amrclaw/data.py b/src/python/amrclaw/data.py index c1439def9..175015ae3 100755 --- a/src/python/amrclaw/data.py +++ b/src/python/amrclaw/data.py @@ -22,7 +22,7 @@ def __init__(self, clawdata): # Need to have a pointer to this so we can get num_dim and num_aux self._clawdata = clawdata - + # Refinement control self.add_attribute('amr_levels_max',1) self.add_attribute('refinement_ratios_x',[1]) @@ -45,7 +45,7 @@ def __init__(self, clawdata): self.add_attribute('clustering_cutoff',0.7) self.add_attribute('verbosity_regrid',3) - + # debugging flags: self.add_attribute('dprint',False) self.add_attribute('eprint',False) @@ -62,7 +62,7 @@ def __init__(self, clawdata): def write(self, out_file='amr.data', data_source='setrun.py'): self.open_data_file(out_file, data_source) - + self.data_write('amr_levels_max') num_ratios = max(abs(self.amr_levels_max)-1, 1) @@ -89,14 +89,24 @@ def write(self, out_file='amr.data', data_source='setrun.py'): if len(self.aux_type) < self._clawdata.num_aux: raise ValueError("*** length of aux_type array should be " +\ "same as num_aux = %i" % self._clawdata.num_aux) - if self._clawdata.num_aux > 0: + + # If necessary, adjust data for absorbing boundary layer + if self._clawdata.abldata.abltype != 0: + aux_type = self.aux_type + aux_type.append('xleft') + if self._clawdata.num_dim > 1: + aux_type.append('yleft') + if self._clawdata.num_dim > 2: + aux_type.append('zleft') + self.data_write('', value=aux_type, alt_name='aux_type') + elif self._clawdata.num_aux > 0: self.data_write('aux_type') self.data_write() self.data_write('flag_richardson') if self.flag_richardson == False: # Still need to add flag_richardson to fortran, for now this works: - self.flag_richardson_tol = -1. + self.flag_richardson_tol = -1. self.data_write('flag_richardson_tol') self.data_write('flag2refine') self.data_write('flag2refine_tol') @@ -149,7 +159,7 @@ def write(self,out_file='regions.data',data_source='setrun.py'): # write minlevel,maxlevel as integers, remaining floats in e format: self._out_file.write("%i %i " % (regions[0], regions[1])) self._out_file.write((len(regions)-2)*"%20.14e " % tuple(regions[2:]) +"\n") - + self.close_data_file() @@ -169,7 +179,7 @@ def read(self, path, force=False): else: break - # Read in number of regions + # Read in number of regions num_regions, tail = line.split("=:") num_regions = int(num_regions) varname = tail.split()[0] @@ -186,7 +196,7 @@ def read(self, path, force=False): # ============================================================================== - + # ============================================================================== # Gauge data object class GaugeData(clawpack.clawutil.data.ClawData): @@ -225,12 +235,12 @@ def __str__(self): output = "\t".join((output, "display: %s" % self.display_format)) output = "\n\t".join((output, "q fields: %s" % self.q_out_fields)) output = "\n\t".join((output, "aux fields: %s" % self.aux_out_fields)) - output = "\n\t".join((output, "min. time increment: %s" + output = "\n\t".join((output, "min. time increment: %s" % self.min_time_increment)) return output - def write(self, num_eqn, num_aux, out_file='gauges.data', + def write(self, num_eqn, num_aux, out_file='gauges.data', data_source='setrun.py'): r"""Write out gauge information data file.""" @@ -246,7 +256,7 @@ def write(self, num_eqn, num_aux, out_file='gauges.data', format = "%4i" + (len(gauge)-3) * " %17.10e" + 2 * " %13.6e" + "\n" self._out_file.write(format % tuple(gauge)) self.data_write() - + # Expand all gauge format option dictionaries for key in self.defaults.keys(): self.expand_gauge_format_option(key) @@ -282,16 +292,16 @@ def write(self, num_eqn, num_aux, out_file='gauges.data', if isinstance(self.q_out_fields[gauge_num], six.string_types): if self.q_out_fields[gauge_num].lower() == 'all': self._out_file.write("%s\n" % " ".join(['True'] * num_eqn)) - elif self.q_out_fields[gauge_num].lower() == 'none': + elif self.q_out_fields[gauge_num].lower() == 'none': self._out_file.write("%s\n" % " ".join(['False'] * num_eqn)) else: - raise ValueError("Unknown q field string specified, '%s'" + raise ValueError("Unknown q field string specified, '%s'" % self.q_out_fields[gauge_num]) else: # Specify by field number if not isinstance(self.q_out_fields[gauge_num], list): self.q_out_fields[gauge_num] = [self.q_out_fields[gauge_num]] - bool_list = [n in self.q_out_fields[gauge_num] + bool_list = [n in self.q_out_fields[gauge_num] for n in range(num_eqn)] bool_list = [str(value) for value in bool_list] self._out_file.write("%s\n" % (" ".join(bool_list))) @@ -305,16 +315,16 @@ def write(self, num_eqn, num_aux, out_file='gauges.data', if isinstance(self.aux_out_fields[gauge_num], six.string_types): if self.aux_out_fields[gauge_num].lower() == 'all': self._out_file.write("%s\n" % " ".join(['True'] * num_aux)) - elif self.aux_out_fields[gauge_num].lower() == 'none': + elif self.aux_out_fields[gauge_num].lower() == 'none': self._out_file.write("%s\n" % " ".join(['False'] * num_aux)) else: - raise ValueError("Unknown q field string specified, '%s'" + raise ValueError("Unknown q field string specified, '%s'" % self.aux_out_fields[gauge_num]) else: # Specify by field number if not isinstance(self.aux_out_fields[gauge_num], list): self.aux_out_fields[gauge_num] = [self.aux_out_fields[gauge_num]] - bool_list = [n in self.aux_out_fields[gauge_num] + bool_list = [n in self.aux_out_fields[gauge_num] for n in range(num_aux)] bool_list = [str(value) for value in bool_list] self._out_file.write("%s\n" % (" ".join(bool_list))) @@ -385,4 +395,3 @@ def read(self, data_path="./", file_name='gauges.data'): if __name__ == '__main__': raise Exception("Not unit tests have been defined for this module.") -