diff --git a/paper2_examples/README.md b/paper2_examples/README.md index d62f7c8..57e0e72 100644 --- a/paper2_examples/README.md +++ b/paper2_examples/README.md @@ -1,10 +1,18 @@ -The code in this repository supplements a paper that is currently in production. +# Computational Experiments to accompany the paper -## Computational Experiments for the paper Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack +### "Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack" -The code in this repository supplements the paper (currently in production): "Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack." +by Brisa N. Davis and Randall J. LeVeque +This version goes with version 2 of the paper, currently submitted for publication. + +Check out the `paper2v1` branch for the version of the code that accompanies the original preprint https://arxiv.org/abs/1810.00927. + +This version is updated in the following ways: + - It uses Clawpack v5.6.1, which has the adjoint flagging now built in. + - Example 1 from the original paper has been dropped from v2 of the paper, so the subdirectories have been renamed. + The code is designed to work with the Clawpack 5 software. For more information visit [the Clawpack webpage](http://www.clawpack.org/ ). The Clawpack software is open source, and the code is openly available at @@ -18,10 +26,10 @@ the use of the adjoint method for guiding adaptive mesh refinement in Clawpack i * gfortran 5.0.0 * python 2.7.10 -* clawpack 5.3.0 : [Clawpack](http://www.clawpack.org/ ) +* clawpack 5.6.1 : [Clawpack](http://www.clawpack.org/ ) * git: to clone this repository -Note: Clawpack 5.0 or higher is required, other versions of gfortran and python might work as well. +Note: Later versions of Clawpack should work, other versions of gfortran and python might work as well. **Operating system.** @@ -30,24 +38,24 @@ Note: Clawpack 5.0 or higher is required, other versions of gfortran and python ### Examples Included and Folder Organization -* **acoustics_1d_ex1:** contains the code to run the first one-dimensional variable coefficients -linear acoustics equations example. +* **acoustics_1d_ex1_old:** contains the code to run the first one-dimensional variable coefficients +linear acoustics example from v1 of the paper, not used in v2 of the paper. The internal folder: **adjoint**, contains the code for the adjoint problem. -* **acoustics_1d_ex2:** contains the code to run the second one-dimensional variable coefficients -linear acoustics equations example. +* **acoustics_1d_example1:** contains the code to run the one-dimensional variable coefficients +linear acoustics example that is Example 1 in v2 of the paper. The internal folder: **adjoint**, contains the code for the adjoint problem. -* **acoustics_2d_ex1:** contains the code to run the first two-dimensional variable coefficients -linear acoustics equations example. +* **acoustics_2d_example2:** contains the code to run the first two-dimensional variable coefficients +linear acoustics problem that is Example 2 in v2 of the paper. The internal folder: **adjoint**, contains the code for the adjoint problem. -* **acoustics_2d_ex2:** contains the code to run the second two-dimensional variable coefficients -linear acoustics equations example. +* **acoustics_2d_example3:** contains the code to run the second two-dimensional variable coefficients +linear acoustics problem that is Example 3 in v2 of the paper. The internal folder: **adjoint**, contains the code for the adjoint problem. ### Running the Code -* **Install Clawpack 5.3.0:** +* **Install Clawpack 5.6.1:** - Go to: http://www.clawpack.org/installing.html#installation-instructions - Follow the download and installation intructions in the section: Install all Clawpack packages - Follow the setup instructions in the section: Set environment variables. @@ -60,9 +68,11 @@ The internal folder: **adjoint**, contains the code for the adjoint problem. - Create a copy of the adjoint repository by running the code: ``` -git://github.com/BrisaDavis/adjoint.git + git clone https://github.com/clawpack/adjoint.git + cd adjoint + git checkout paper2v2 ``` * **Run the code** -* Go any of the folders and follow the instructions in the README.md in that folder. +* Go any of the folders in `paper2_examples` and follow the instructions in the README.md in that folder. diff --git a/paper2_examples/acoustics_1d_ex1/setplot.py b/paper2_examples/acoustics_1d_ex1/setplot.py deleted file mode 100644 index e0fbcba..0000000 --- a/paper2_examples/acoustics_1d_ex1/setplot.py +++ /dev/null @@ -1,105 +0,0 @@ - -""" -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 - -#-------------------------- -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() - - plotdata.clearfigures() # clear any old figures,axes,items data - - def fix_plot(current_data): - from pylab import plot - from pylab import xticks,yticks,xlabel,ylabel,savefig,ylim,title - t = current_data.t - plot([0., 0.], [-1000., 1000.], 'k--') - title('Pressure at t = %5.3f seconds' % t, fontsize=26) - yticks(fontsize=23) - xticks(fontsize=23) - - def fix_plot_innerprod(current_data): - from pylab import plot - from pylab import xticks,yticks,xlabel,ylabel,savefig,ylim,title - t = current_data.t - plot([0., 0.], [-1000., 1000.], 'k--') - title('Pressure at t = %5.3f seconds' % t, fontsize=26) - yticks(fontsize=23) - xticks(fontsize=23) - - - # Figure for q[0] - plotfigure = plotdata.new_plotfigure(name='Pressure', figno=1) - plotfigure.kwargs = {'figsize': (10,3.5)} - # Set up for axes in this figure: - plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = [-12,12] - plotaxes.ylimits = [-1.1,1.1] - plotaxes.title = 'Pressure' - plotaxes.afteraxes = fix_plot - - # Set up for item on these axes: - plotitem = plotaxes.new_plotitem(plot_type='1d_plot') - plotitem.plot_var = 0 - plotitem.amr_color = 'b' - plotitem.amr_plotstyle = 'o' - plotitem.amr_kwargs = [{'linewidth':2}] - plotitem.amr_kwargs = [{'markersize':4}] - - # Figure for inner product, q[2] - - plotfigure = plotdata.new_plotfigure(name='Inner Product', figno=10) - - # Set up for axes in this figure: - plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = 'auto' - #plotaxes.ylimits = [-.5,1.1] # use when taking inner product with forward solution - plotaxes.ylimits = [0,0.02] # use when taking inner product with Richardson error - plotaxes.title = 'Inner Product' - plotaxes.afteraxes = fix_plot_innerprod - - # Set up for item on these axes: - plotitem = plotaxes.new_plotitem(plot_type='1d') - plotitem.plot_var = 2 - plotitem.amr_color = 'b' - plotitem.amr_plotstyle = 'o' - plotitem.amr_kwargs = [{'linewidth':2}] - plotitem.amr_data_show = [1,1,1,0] - plotitem.show = True # show on plot? - - # 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' - 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 - - diff --git a/paper2_examples/acoustics_1d_ex2/Makefile b/paper2_examples/acoustics_1d_ex1_old/Makefile similarity index 95% rename from paper2_examples/acoustics_1d_ex2/Makefile rename to paper2_examples/acoustics_1d_ex1_old/Makefile index c270adc..04c1e79 100644 --- a/paper2_examples/acoustics_1d_ex2/Makefile +++ b/paper2_examples/acoustics_1d_ex1_old/Makefile @@ -31,8 +31,8 @@ FFLAGS ?= # package sources for this program: # --------------------------------- -AMRLIB = $(CLAW)/amrclaw/src/1d_adjoint -include $(AMRLIB)/Makefile.adjoint_amr_1d +AMRLIB = $(CLAW)/amrclaw/src/1d +include $(AMRLIB)/Makefile.amr_1d # --------------------------------------- # package sources specifically to exclude diff --git a/paper2_examples/acoustics_1d_ex1/README.md b/paper2_examples/acoustics_1d_ex1_old/README.md similarity index 79% rename from paper2_examples/acoustics_1d_ex1/README.md rename to paper2_examples/acoustics_1d_ex1_old/README.md index 4dd3a50..564f5ac 100644 --- a/paper2_examples/acoustics_1d_ex1/README.md +++ b/paper2_examples/acoustics_1d_ex1_old/README.md @@ -1,4 +1,6 @@ -This code produces Example 1, the first one-dimensional variable coefficient acoustics example, presented in ''Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack''. +This code was for the original Example 1, the first one-dimensional variable +coefficient acoustics example, presented in the preprint version of +''Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack''. ### Folder Organization * **adjoint:** diff --git a/paper2_examples/acoustics_1d_ex1/adjoint/Makefile b/paper2_examples/acoustics_1d_ex1_old/adjoint/Makefile similarity index 100% rename from paper2_examples/acoustics_1d_ex1/adjoint/Makefile rename to paper2_examples/acoustics_1d_ex1_old/adjoint/Makefile diff --git a/paper2_examples/acoustics_1d_ex1/adjoint/qinit.f b/paper2_examples/acoustics_1d_ex1_old/adjoint/qinit.f similarity index 100% rename from paper2_examples/acoustics_1d_ex1/adjoint/qinit.f rename to paper2_examples/acoustics_1d_ex1_old/adjoint/qinit.f diff --git a/paper2_examples/acoustics_1d_ex1/adjoint/setaux.f b/paper2_examples/acoustics_1d_ex1_old/adjoint/setaux.f similarity index 100% rename from paper2_examples/acoustics_1d_ex1/adjoint/setaux.f rename to paper2_examples/acoustics_1d_ex1_old/adjoint/setaux.f diff --git a/paper2_examples/acoustics_1d_ex1/adjoint/setplot.py b/paper2_examples/acoustics_1d_ex1_old/adjoint/setplot.py similarity index 100% rename from paper2_examples/acoustics_1d_ex1/adjoint/setplot.py rename to paper2_examples/acoustics_1d_ex1_old/adjoint/setplot.py diff --git a/paper2_examples/acoustics_1d_ex1/adjoint/setprob.f b/paper2_examples/acoustics_1d_ex1_old/adjoint/setprob.f similarity index 100% rename from paper2_examples/acoustics_1d_ex1/adjoint/setprob.f rename to paper2_examples/acoustics_1d_ex1_old/adjoint/setprob.f diff --git a/paper2_examples/acoustics_1d_ex1/adjoint/setrun.py b/paper2_examples/acoustics_1d_ex1_old/adjoint/setrun.py similarity index 100% rename from paper2_examples/acoustics_1d_ex1/adjoint/setrun.py rename to paper2_examples/acoustics_1d_ex1_old/adjoint/setrun.py diff --git a/paper2_examples/acoustics_1d_ex1/qinit.f b/paper2_examples/acoustics_1d_ex1_old/qinit.f similarity index 100% rename from paper2_examples/acoustics_1d_ex1/qinit.f rename to paper2_examples/acoustics_1d_ex1_old/qinit.f diff --git a/paper2_examples/acoustics_1d_ex2/setaux.f b/paper2_examples/acoustics_1d_ex1_old/setaux.f similarity index 88% rename from paper2_examples/acoustics_1d_ex2/setaux.f rename to paper2_examples/acoustics_1d_ex1_old/setaux.f index 35af59c..5301706 100644 --- a/paper2_examples/acoustics_1d_ex2/setaux.f +++ b/paper2_examples/acoustics_1d_ex1_old/setaux.f @@ -22,10 +22,15 @@ subroutine setaux(mbc,mx,xlower,dx,maux,aux) common /comaux/ Zl, cl, Zr, cr double precision Zl, cl, Zr, cr - integer i,ii + integer i double precision xcell do i=1-mbc,mx+mbc + + if (aux(1,i) .eq. NEEDS_TO_BE_SET) then + aux(innerprod_index,i) = 0.d0 + endif + xcell = xlower + (i-0.5d0)*dx if (xcell .lt. 0.0d0) then aux(1,i) = Zl @@ -36,10 +41,5 @@ subroutine setaux(mbc,mx,xlower,dx,maux,aux) endif enddo - ! Initialize innerproduct to zero. - do ii=1-mbc,mx+mbc - aux(innerprod_index,ii) = 0.d0 - enddo - return end diff --git a/paper2_examples/acoustics_1d_ex1_old/setplot.py b/paper2_examples/acoustics_1d_ex1_old/setplot.py new file mode 100644 index 0000000..9f16b04 --- /dev/null +++ b/paper2_examples/acoustics_1d_ex1_old/setplot.py @@ -0,0 +1,165 @@ + +""" +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 + +amr_color = ['g','b','r','c','k'] +amr_marker = ['x','s','o','^','.'] +amr_linestyle = ['-','-','-','-','-'] +amr_plotstyle = [ma+li for ma,li in zip(amr_marker,amr_linestyle)] +print('amr_plotstyle = ',amr_plotstyle) +#amr_plotstyle = ['x-','s-','o-','^-','.-'] + +#-------------------------- +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() + + plotdata.clearfigures() # clear any old figures,axes,items data + + def draw_interface_add_legend(current_data): + from pylab import plot + plot([0., 0.], [-1000., 1000.], 'k--') + try: + from clawpack.visclaw import legend_tools + labels = ['Level 1','Level 2', 'Level 3', 'Level 4', 'Level 5'] + legend_tools.add_legend(labels, colors=amr_color, + markers=amr_marker, linestyles=amr_linestyle, + loc='upper left') + except: + pass + + def draw_interface_add_legend_innerprod(current_data): + from pylab import plot + plot([0., 0.], [-1000., 1000.], 'k--') + try: + from clawpack.visclaw import legend_tools + labels = ['Level 3','Level 4'] + legend_tools.add_legend(labels, colors=['r','c'], + markers=['o','^'], linestyles=['',''], + loc='upper left') + except: + pass + + + # Figure for q[0] + plotfigure = plotdata.new_plotfigure(name='Pressure and Velocity', figno=1) + plotfigure.kwargs = {'figsize': (8,8)} + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.axescmd = 'subplot(2,1,1)' # top figure + plotaxes.xlimits = [-12,12] + plotaxes.ylimits = [-.5,1.1] + plotaxes.title = 'Pressure' + plotaxes.afteraxes = draw_interface_add_legend + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.amr_color = amr_color + plotitem.amr_plotstyle = amr_plotstyle + plotitem.amr_data_show = [1,1,1] + plotitem.amr_kwargs = [{'markersize':5},{'markersize':4},{'markersize':3}] + + # Figure for q[1] + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.axescmd = 'subplot(2,1,2)' # bottom figure + plotaxes.xlimits = [-12,12] + plotaxes.ylimits = [-.5,1.1] + plotaxes.title = 'Velocity' + plotaxes.afteraxes = draw_interface_add_legend + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 1 + plotitem.amr_color = amr_color + plotitem.amr_plotstyle = amr_plotstyle + plotitem.amr_data_show = [1,1,1] + plotitem.amr_kwargs = [{'markersize':5},{'markersize':4},{'markersize':3}] + + # Figure for inner product, q[2] + + plotfigure = plotdata.new_plotfigure(name='Inner Product', figno=10) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [-12,12] + plotaxes.ylimits = [-0.1,1] # use when taking inner product with forward solution + #plotaxes.ylimits = [-0.01,0.02] # use when taking inner product with Richardson error + plotaxes.title = 'Inner Product' + plotaxes.afteraxes = draw_interface_add_legend_innerprod + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d') + plotitem.plot_var = plot_innerprod + plotitem.amr_color = amr_color + plotitem.amr_plotstyle = amr_plotstyle + plotitem.amr_data_show = [0,0,1,1,0] + plotitem.show = True # show on plot? + + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='q', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + plotfigure.kwargs = {'figsize': (10,10)} + + plotaxes = plotfigure.new_plotaxes() + plotaxes.axescmd = 'subplot(211)' + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Pressure' + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.plotstyle = 'b-' + + plotaxes = plotfigure.new_plotaxes() + plotaxes.axescmd = 'subplot(212)' + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Velocity' + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 1 + 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' + 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 + +def plot_innerprod(current_data): + return current_data.aux[2,:] + diff --git a/paper2_examples/acoustics_1d_ex1/setprob.f b/paper2_examples/acoustics_1d_ex1_old/setprob.f similarity index 92% rename from paper2_examples/acoustics_1d_ex1/setprob.f rename to paper2_examples/acoustics_1d_ex1_old/setprob.f index 4f3a86d..1d0a44a 100644 --- a/paper2_examples/acoustics_1d_ex1/setprob.f +++ b/paper2_examples/acoustics_1d_ex1_old/setprob.f @@ -35,7 +35,5 @@ subroutine setprob read(7,*) cr Zr = rhor*cr - call read_adjoint_data() !# Read adjoint info and solution - return end diff --git a/paper2_examples/acoustics_1d_ex1/setrun.py b/paper2_examples/acoustics_1d_ex1_old/setrun.py similarity index 78% rename from paper2_examples/acoustics_1d_ex1/setrun.py rename to paper2_examples/acoustics_1d_ex1_old/setrun.py index 1c51e09..ae947f1 100644 --- a/paper2_examples/acoustics_1d_ex1/setrun.py +++ b/paper2_examples/acoustics_1d_ex1_old/setrun.py @@ -10,26 +10,7 @@ import numpy as np -#----------------------------------------------- -# Set these parameters for adjoint flagging.... - -# location of output from computing adjoint: -adjoint_output = os.path.abspath('adjoint/_output') -print('Will flag using adjoint solution from %s' % adjoint_output) - -# Time period of interest: -t1 = 33.5 -t2 = 34.0 - -# Determining type of adjoint flagging: - -# taking inner product with forward solution or Richardson error: -flag_forward_adjoint = True -flag_richardson_adjoint = False - -# tolerance for adjoint flagging: -adjoint_flag_tolerance = 1e-3 # suggested if using forward solution -#adjoint_flag_tolerance = 1e-5 # suggested if using Richardson error +# This version is set up for adjoint flagging #----------------------------------------------- @@ -50,7 +31,8 @@ def setrun(claw_pkg='amrclaw'): """ from clawpack.clawutil import data - + + assert claw_pkg.lower() == 'amrclaw', "Expected claw_pkg = 'amrclaw'" num_dim = 1 @@ -100,7 +82,8 @@ def setrun(claw_pkg='amrclaw'): clawdata.num_eqn = 2 # Number of auxiliary variables in the aux array (initialized in setaux) - # see setadjoint + rundata.clawdata.num_aux = 2 + # Note: as required for original problem - modified below for adjoint # Index of aux array corresponding to capacity function, if there is one: clawdata.capa_index = 0 @@ -155,8 +138,8 @@ def setrun(claw_pkg='amrclaw'): 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 + clawdata.output_aux_components = 'all' # could be list + clawdata.output_aux_onlyonce = False # output aux arrays only at t0 # --------------------------------------------------- @@ -292,16 +275,18 @@ def setrun(claw_pkg='amrclaw'): # This must be a list of length num_aux, each element of which is one of: # 'center', 'capacity', or 'xleft' (see documentation). - # need 3 values, set in setadjoint + rundata.amrdata.aux_type = ['center','center'] + # Note: as required for original problem - modified below for adjoint + # set tolerances appropriate for adjoint flagging: # Flag for refinement based on Richardson error estimater: amrdata.flag_richardson = False + amrdata.flag_richardson_tol = 1e-5 # suggested if using adjoint-error # Flag for refinement using routine flag2refine: - amrdata.flag2refine = False - - # see setadjoint to set tolerance and flag for adjoint flagging + amrdata.flag2refine = True + rundata.amrdata.flag2refine_tol = 1e-3 # suggested if using adjoint-mag # steps to take on each level L between regriddings of level L+1: amrdata.regrid_interval = 2 @@ -328,7 +313,23 @@ def setrun(claw_pkg='amrclaw'): #------------------------------------------------------------------ # Adjoint specific data: #------------------------------------------------------------------ - rundata = setadjoint(rundata) + # Also need to set flagging method and appropriate tolerances above + + adjointdata = rundata.adjointdata + adjointdata.use_adjoint = True + + # location of adjoint solution, must first be created: + adjointdata.adjoint_outdir = os.path.abspath('adjoint/_output') + + # time period of interest: + adjointdata.t1 = 33.5 + adjointdata.t2 = 34.0 + + if adjointdata.use_adjoint: + # need an additional aux variable for inner product: + rundata.amrdata.aux_type.append('center') + rundata.clawdata.num_aux = len(rundata.amrdata.aux_type) + adjointdata.innerprod_index = len(rundata.amrdata.aux_type) # ----- For developers ----- @@ -350,61 +351,6 @@ def setrun(claw_pkg='amrclaw'): # end of function setrun # ---------------------- -#------------------- -def setadjoint(rundata): - #------------------- - - """ - Setting up adjoint variables and - reading in all of the checkpointed Adjoint files - """ - - import glob - import os - - # Set these parameters at top of this file: - # adjoint_flag_tolerance, t1, t2, adjoint_output - # Then you don't need to modify this function... - - # flag and tolerance for adjoint flagging: - if flag_forward_adjoint == True: - # setting up taking inner product with forward solution - rundata.amrdata.flag2refine = True - rundata.amrdata.flag2refine_tol = adjoint_flag_tolerance - elif flag_richardson_adjoint == True: - # setting up taking inner product with Richardson error - rundata.amrdata.flag_richardson = True - rundata.amrdata.flag_richardson_tol = adjoint_flag_tolerance - else: - print("No refinement flag set!") - - rundata.clawdata.num_aux = 3 # 3 required for adjoint flagging - rundata.amrdata.aux_type = ['center','center','center'] - - adjointdata = rundata.new_UserData(name='adjointdata',fname='adjoint.data') - adjointdata.add_param('adjoint_output',adjoint_output,'adjoint_output') - adjointdata.add_param('t1',t1,'t1, start time of interest') - adjointdata.add_param('t2',t2,'t2, final time of interest') - - files = glob.glob(os.path.join(adjoint_output,"fort.b*")) - files.sort() - - if (len(files) == 0): - print("No binary files found for adjoint output!") - - adjointdata.add_param('numadjoints', len(files), 'Number of adjoint checkpoint files.') - adjointdata.add_param('innerprod_index', 3, 'Index for innerproduct data in aux array.') - - counter = 1 - for fname in files: - f = open(fname) - adjointdata.add_param('file' + str(counter), fname, 'Binary file' + str(counter)) - counter = counter + 1 - - return rundata -# end of function setadjoint -# ---------------------- - if __name__ == '__main__': # Set up run-time parameters and write all data files. diff --git a/paper2_examples/acoustics_1d_ex2/setplot.py b/paper2_examples/acoustics_1d_ex2/setplot.py deleted file mode 100644 index 3524b41..0000000 --- a/paper2_examples/acoustics_1d_ex2/setplot.py +++ /dev/null @@ -1,105 +0,0 @@ - -""" -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 - -#-------------------------- -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() - - plotdata.clearfigures() # clear any old figures,axes,items data - - def fix_plot(current_data): - from pylab import plot - from pylab import xticks,yticks,xlabel,ylabel,savefig,ylim,title - t = current_data.t - plot([0., 0.], [-1000., 1000.], 'k--') - title('Pressure at t = %5.3f seconds' % t, fontsize=26) - yticks(fontsize=23) - xticks(fontsize=23) - - def fix_plot_innerprod(current_data): - from pylab import plot - from pylab import xticks,yticks,xlabel,ylabel,savefig,ylim,title - t = current_data.t - plot([0., 0.], [-1000., 1000.], 'k--') - title('Inner Product at t = %5.3f seconds' % t, fontsize=26) - yticks(fontsize=23) - xticks(fontsize=23) - - - # Figure for q[0] - plotfigure = plotdata.new_plotfigure(name='Pressure', figno=1) - plotfigure.kwargs = {'figsize': (10,3.5)} - # Set up for axes in this figure: - plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = [-12,12] - plotaxes.ylimits = [-1.1,1.1] - plotaxes.title = 'Pressure' - plotaxes.afteraxes = fix_plot - - # Set up for item on these axes: - plotitem = plotaxes.new_plotitem(plot_type='1d_plot') - plotitem.plot_var = 0 - plotitem.amr_color = 'b' - plotitem.amr_plotstyle = 'o' - plotitem.amr_kwargs = [{'linewidth':2}] - plotitem.amr_kwargs = [{'markersize':4}] - - # Figure for inner product, q[2] - - plotfigure = plotdata.new_plotfigure(name='Inner Product', figno=10) - - # Set up for axes in this figure: - plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = 'auto' - #plotaxes.ylimits = [-.5,1.1] # use when taking inner product with forward solution - plotaxes.ylimits = [-0.01,0.02] # use when taking inner product with Richardson error - plotaxes.title = 'Inner Product' - plotaxes.afteraxes = fix_plot_innerprod - - # Set up for item on these axes: - plotitem = plotaxes.new_plotitem(plot_type='1d') - plotitem.plot_var = 2 - plotitem.amr_color = 'b' - plotitem.amr_plotstyle = 'o' - plotitem.amr_kwargs = [{'linewidth':2}] - plotitem.amr_kwargs = [{'markersize':4}] - plotitem.show = True # show on plot? - - # 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' - 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 - - diff --git a/paper2_examples/acoustics_1d_ex1/Makefile b/paper2_examples/acoustics_1d_example1/Makefile similarity index 94% rename from paper2_examples/acoustics_1d_ex1/Makefile rename to paper2_examples/acoustics_1d_example1/Makefile index ce57bdc..04c1e79 100644 --- a/paper2_examples/acoustics_1d_ex1/Makefile +++ b/paper2_examples/acoustics_1d_example1/Makefile @@ -25,14 +25,14 @@ RESTART ?= False # Should = clawdata.restart in setrun # 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 ?= -O2 -fopenmp +FFLAGS ?= # --------------------------------- # package sources for this program: # --------------------------------- -AMRLIB = $(CLAW)/amrclaw/src/1d_adjoint -include $(AMRLIB)/Makefile.adjoint_amr_1d +AMRLIB = $(CLAW)/amrclaw/src/1d +include $(AMRLIB)/Makefile.amr_1d # --------------------------------------- # package sources specifically to exclude diff --git a/paper2_examples/acoustics_2d_ex4/README.md b/paper2_examples/acoustics_1d_example1/README.md similarity index 71% rename from paper2_examples/acoustics_2d_ex4/README.md rename to paper2_examples/acoustics_1d_example1/README.md index f439cb8..0e7ef6e 100644 --- a/paper2_examples/acoustics_2d_ex4/README.md +++ b/paper2_examples/acoustics_1d_example1/README.md @@ -1,4 +1,4 @@ -This code produces Example 4, the second two-dimensional variable coefficient acoustics example, presented in ''Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack''. +This code produces Example 1, the one-dimensional variable coefficient acoustics example, presented in ''Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack''. ### Folder Organization * **adjoint:** @@ -20,7 +20,7 @@ The code will produce two new folders: _output and _plots. The first one contains all the output files, while the latter one contains the plots and interactive visualization apps. -* Go to the main folder **acoustics_2d_ex3** and run in the terminal: +* Go to the main folder **acoustics_1d_example2** and run in the terminal: ``` make new diff --git a/paper2_examples/acoustics_1d_ex2/adjoint/Makefile b/paper2_examples/acoustics_1d_example1/adjoint/Makefile similarity index 100% rename from paper2_examples/acoustics_1d_ex2/adjoint/Makefile rename to paper2_examples/acoustics_1d_example1/adjoint/Makefile diff --git a/paper2_examples/acoustics_1d_ex2/adjoint/qinit.f b/paper2_examples/acoustics_1d_example1/adjoint/qinit.f similarity index 100% rename from paper2_examples/acoustics_1d_ex2/adjoint/qinit.f rename to paper2_examples/acoustics_1d_example1/adjoint/qinit.f diff --git a/paper2_examples/acoustics_1d_ex2/adjoint/setaux.f b/paper2_examples/acoustics_1d_example1/adjoint/setaux.f similarity index 100% rename from paper2_examples/acoustics_1d_ex2/adjoint/setaux.f rename to paper2_examples/acoustics_1d_example1/adjoint/setaux.f diff --git a/paper2_examples/acoustics_1d_ex2/adjoint/setplot.py b/paper2_examples/acoustics_1d_example1/adjoint/setplot.py similarity index 100% rename from paper2_examples/acoustics_1d_ex2/adjoint/setplot.py rename to paper2_examples/acoustics_1d_example1/adjoint/setplot.py diff --git a/paper2_examples/acoustics_1d_ex2/adjoint/setprob.f b/paper2_examples/acoustics_1d_example1/adjoint/setprob.f similarity index 100% rename from paper2_examples/acoustics_1d_ex2/adjoint/setprob.f rename to paper2_examples/acoustics_1d_example1/adjoint/setprob.f diff --git a/paper2_examples/acoustics_1d_ex2/adjoint/setrun.py b/paper2_examples/acoustics_1d_example1/adjoint/setrun.py similarity index 100% rename from paper2_examples/acoustics_1d_ex2/adjoint/setrun.py rename to paper2_examples/acoustics_1d_example1/adjoint/setrun.py diff --git a/paper2_examples/acoustics_1d_ex2/qinit.f b/paper2_examples/acoustics_1d_example1/qinit.f similarity index 100% rename from paper2_examples/acoustics_1d_ex2/qinit.f rename to paper2_examples/acoustics_1d_example1/qinit.f diff --git a/paper2_examples/acoustics_1d_ex1/setaux.f b/paper2_examples/acoustics_1d_example1/setaux.f similarity index 88% rename from paper2_examples/acoustics_1d_ex1/setaux.f rename to paper2_examples/acoustics_1d_example1/setaux.f index 35af59c..5301706 100644 --- a/paper2_examples/acoustics_1d_ex1/setaux.f +++ b/paper2_examples/acoustics_1d_example1/setaux.f @@ -22,10 +22,15 @@ subroutine setaux(mbc,mx,xlower,dx,maux,aux) common /comaux/ Zl, cl, Zr, cr double precision Zl, cl, Zr, cr - integer i,ii + integer i double precision xcell do i=1-mbc,mx+mbc + + if (aux(1,i) .eq. NEEDS_TO_BE_SET) then + aux(innerprod_index,i) = 0.d0 + endif + xcell = xlower + (i-0.5d0)*dx if (xcell .lt. 0.0d0) then aux(1,i) = Zl @@ -36,10 +41,5 @@ subroutine setaux(mbc,mx,xlower,dx,maux,aux) endif enddo - ! Initialize innerproduct to zero. - do ii=1-mbc,mx+mbc - aux(innerprod_index,ii) = 0.d0 - enddo - return end diff --git a/paper2_examples/acoustics_1d_example1/setplot.py b/paper2_examples/acoustics_1d_example1/setplot.py new file mode 100644 index 0000000..9f16b04 --- /dev/null +++ b/paper2_examples/acoustics_1d_example1/setplot.py @@ -0,0 +1,165 @@ + +""" +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 + +amr_color = ['g','b','r','c','k'] +amr_marker = ['x','s','o','^','.'] +amr_linestyle = ['-','-','-','-','-'] +amr_plotstyle = [ma+li for ma,li in zip(amr_marker,amr_linestyle)] +print('amr_plotstyle = ',amr_plotstyle) +#amr_plotstyle = ['x-','s-','o-','^-','.-'] + +#-------------------------- +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() + + plotdata.clearfigures() # clear any old figures,axes,items data + + def draw_interface_add_legend(current_data): + from pylab import plot + plot([0., 0.], [-1000., 1000.], 'k--') + try: + from clawpack.visclaw import legend_tools + labels = ['Level 1','Level 2', 'Level 3', 'Level 4', 'Level 5'] + legend_tools.add_legend(labels, colors=amr_color, + markers=amr_marker, linestyles=amr_linestyle, + loc='upper left') + except: + pass + + def draw_interface_add_legend_innerprod(current_data): + from pylab import plot + plot([0., 0.], [-1000., 1000.], 'k--') + try: + from clawpack.visclaw import legend_tools + labels = ['Level 3','Level 4'] + legend_tools.add_legend(labels, colors=['r','c'], + markers=['o','^'], linestyles=['',''], + loc='upper left') + except: + pass + + + # Figure for q[0] + plotfigure = plotdata.new_plotfigure(name='Pressure and Velocity', figno=1) + plotfigure.kwargs = {'figsize': (8,8)} + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.axescmd = 'subplot(2,1,1)' # top figure + plotaxes.xlimits = [-12,12] + plotaxes.ylimits = [-.5,1.1] + plotaxes.title = 'Pressure' + plotaxes.afteraxes = draw_interface_add_legend + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.amr_color = amr_color + plotitem.amr_plotstyle = amr_plotstyle + plotitem.amr_data_show = [1,1,1] + plotitem.amr_kwargs = [{'markersize':5},{'markersize':4},{'markersize':3}] + + # Figure for q[1] + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.axescmd = 'subplot(2,1,2)' # bottom figure + plotaxes.xlimits = [-12,12] + plotaxes.ylimits = [-.5,1.1] + plotaxes.title = 'Velocity' + plotaxes.afteraxes = draw_interface_add_legend + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 1 + plotitem.amr_color = amr_color + plotitem.amr_plotstyle = amr_plotstyle + plotitem.amr_data_show = [1,1,1] + plotitem.amr_kwargs = [{'markersize':5},{'markersize':4},{'markersize':3}] + + # Figure for inner product, q[2] + + plotfigure = plotdata.new_plotfigure(name='Inner Product', figno=10) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [-12,12] + plotaxes.ylimits = [-0.1,1] # use when taking inner product with forward solution + #plotaxes.ylimits = [-0.01,0.02] # use when taking inner product with Richardson error + plotaxes.title = 'Inner Product' + plotaxes.afteraxes = draw_interface_add_legend_innerprod + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d') + plotitem.plot_var = plot_innerprod + plotitem.amr_color = amr_color + plotitem.amr_plotstyle = amr_plotstyle + plotitem.amr_data_show = [0,0,1,1,0] + plotitem.show = True # show on plot? + + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='q', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + plotfigure.kwargs = {'figsize': (10,10)} + + plotaxes = plotfigure.new_plotaxes() + plotaxes.axescmd = 'subplot(211)' + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Pressure' + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 0 + plotitem.plotstyle = 'b-' + + plotaxes = plotfigure.new_plotaxes() + plotaxes.axescmd = 'subplot(212)' + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Velocity' + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 1 + 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' + 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 + +def plot_innerprod(current_data): + return current_data.aux[2,:] + diff --git a/paper2_examples/acoustics_1d_ex2/setprob.f b/paper2_examples/acoustics_1d_example1/setprob.f similarity index 92% rename from paper2_examples/acoustics_1d_ex2/setprob.f rename to paper2_examples/acoustics_1d_example1/setprob.f index 4f3a86d..1d0a44a 100644 --- a/paper2_examples/acoustics_1d_ex2/setprob.f +++ b/paper2_examples/acoustics_1d_example1/setprob.f @@ -35,7 +35,5 @@ subroutine setprob read(7,*) cr Zr = rhor*cr - call read_adjoint_data() !# Read adjoint info and solution - return end diff --git a/paper2_examples/acoustics_1d_ex2/setrun.py b/paper2_examples/acoustics_1d_example1/setrun.py similarity index 78% rename from paper2_examples/acoustics_1d_ex2/setrun.py rename to paper2_examples/acoustics_1d_example1/setrun.py index 385ffa7..df4a1ad 100644 --- a/paper2_examples/acoustics_1d_ex2/setrun.py +++ b/paper2_examples/acoustics_1d_example1/setrun.py @@ -10,26 +10,7 @@ import numpy as np -#----------------------------------------------- -# Set these parameters for adjoint flagging.... - -# location of output from computing adjoint: -adjoint_output = os.path.abspath('adjoint/_output') -print('Will flag using adjoint solution from %s' % adjoint_output) - -# Time period of interest: -t1 = 51.5 -t2 = 52.0 - -# Determining type of adjoint flagging: - -# taking inner product with forward solution or Richardson error: -flag_forward_adjoint = True -flag_richardson_adjoint = False - -# tolerance for adjoint flagging: -adjoint_flag_tolerance = 1e-3 # suggested if using forward solution -#adjoint_flag_tolerance = 1e-5 # suggested if using Richardson error +# This version is set up for adjoint flagging #----------------------------------------------- @@ -50,7 +31,8 @@ def setrun(claw_pkg='amrclaw'): """ from clawpack.clawutil import data - + + assert claw_pkg.lower() == 'amrclaw', "Expected claw_pkg = 'amrclaw'" num_dim = 1 @@ -100,7 +82,8 @@ def setrun(claw_pkg='amrclaw'): clawdata.num_eqn = 2 # Number of auxiliary variables in the aux array (initialized in setaux) - # see setadjoint + rundata.clawdata.num_aux = 2 + # Note: as required for original problem - modified below for adjoint # Index of aux array corresponding to capacity function, if there is one: clawdata.capa_index = 0 @@ -155,8 +138,8 @@ def setrun(claw_pkg='amrclaw'): 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 + clawdata.output_aux_components = 'all' # could be list + clawdata.output_aux_onlyonce = False # output aux arrays only at t0 # --------------------------------------------------- @@ -292,16 +275,18 @@ def setrun(claw_pkg='amrclaw'): # This must be a list of length num_aux, each element of which is one of: # 'center', 'capacity', or 'xleft' (see documentation). - # need 3 values, set in setadjoint + rundata.amrdata.aux_type = ['center','center'] + # Note: as required for original problem - modified below for adjoint + # set tolerances appropriate for adjoint flagging: # Flag for refinement based on Richardson error estimater: amrdata.flag_richardson = False + amrdata.flag_richardson_tol = 1e-5 # suggested if using adjoint-error # Flag for refinement using routine flag2refine: - amrdata.flag2refine = False - - # see setadjoint to set tolerance and flag for adjoint flagging + amrdata.flag2refine = True + rundata.amrdata.flag2refine_tol = 1e-3 # suggested if using adjoint-mag # steps to take on each level L between regriddings of level L+1: amrdata.regrid_interval = 2 @@ -328,7 +313,23 @@ def setrun(claw_pkg='amrclaw'): #------------------------------------------------------------------ # Adjoint specific data: #------------------------------------------------------------------ - rundata = setadjoint(rundata) + # Also need to set flagging method and appropriate tolerances above + + adjointdata = rundata.adjointdata + adjointdata.use_adjoint = True + + # location of adjoint solution, must first be created: + adjointdata.adjoint_outdir = os.path.abspath('adjoint/_output') + + # time period of interest: + adjointdata.t1 = 51.5 + adjointdata.t2 = 52.0 + + if adjointdata.use_adjoint: + # need an additional aux variable for inner product: + rundata.amrdata.aux_type.append('center') + rundata.clawdata.num_aux = len(rundata.amrdata.aux_type) + adjointdata.innerprod_index = len(rundata.amrdata.aux_type) # ----- For developers ----- @@ -350,61 +351,6 @@ def setrun(claw_pkg='amrclaw'): # end of function setrun # ---------------------- -#------------------- -def setadjoint(rundata): - #------------------- - - """ - Setting up adjoint variables and - reading in all of the checkpointed Adjoint files - """ - - import glob - import os - - # Set these parameters at top of this file: - # adjoint_flag_tolerance, t1, t2, adjoint_output - # Then you don't need to modify this function... - - # flag and tolerance for adjoint flagging: - if flag_forward_adjoint == True: - # setting up taking inner product with forward solution - rundata.amrdata.flag2refine = True - rundata.amrdata.flag2refine_tol = adjoint_flag_tolerance - elif flag_richardson_adjoint == True: - # setting up taking inner product with Richardson error - rundata.amrdata.flag_richardson = True - rundata.amrdata.flag_richardson_tol = adjoint_flag_tolerance - else: - print("No refinement flag set!") - - rundata.clawdata.num_aux = 3 # 3 required for adjoint flagging - rundata.amrdata.aux_type = ['center','center','center'] - - adjointdata = rundata.new_UserData(name='adjointdata',fname='adjoint.data') - adjointdata.add_param('adjoint_output',adjoint_output,'adjoint_output') - adjointdata.add_param('t1',t1,'t1, start time of interest') - adjointdata.add_param('t2',t2,'t2, final time of interest') - - files = glob.glob(os.path.join(adjoint_output,"fort.b*")) - files.sort() - - if (len(files) == 0): - print("No binary files found for adjoint output!") - - adjointdata.add_param('numadjoints', len(files), 'Number of adjoint checkpoint files.') - adjointdata.add_param('innerprod_index', 3, 'Index for innerproduct data in aux array.') - - counter = 1 - for fname in files: - f = open(fname) - adjointdata.add_param('file' + str(counter), fname, 'Binary file' + str(counter)) - counter = counter + 1 - - return rundata -# end of function setadjoint -# ---------------------- - if __name__ == '__main__': # Set up run-time parameters and write all data files. diff --git a/paper2_examples/acoustics_2d_ex3/README.md b/paper2_examples/acoustics_2d_ex3/README.md deleted file mode 100644 index f049781..0000000 --- a/paper2_examples/acoustics_2d_ex3/README.md +++ /dev/null @@ -1,39 +0,0 @@ -This code produces Example 3, the first two-dimensional variable coefficient acoustics example, presented in ''Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack''. - -### Folder Organization -* **adjoint:** - -Contains code to solve the adjoint problem. - -The start and end times specified in this directory should agree with those for the forward code. - -### Running the Code - -* Go to the folder **adjoint** and run in a terminal: - -``` -make new -make .plots -``` - -The code will produce two new folders: _output and _plots. -The first one contains all the output files, while the latter one contains the plots and interactive -visualization apps. - -* Go to the main folder **acoustics_2d_ex3** and run in the terminal: - -``` -make new -make .plots -``` - -The code will produce two new folders: _output and _plots. -The first one contains all the output files, while the latter one contains the plots and interactive -visualization apps. - -### Extra Python Scripts - -If you run the script generate_tolplots.py in the terminal you will generate -the tolerance, error, and CPU timing plots that appear in the paper. - - diff --git a/paper2_examples/acoustics_2d_ex3/generate_tolplots.py b/paper2_examples/acoustics_2d_ex3/generate_tolplots.py deleted file mode 100644 index 53abea0..0000000 --- a/paper2_examples/acoustics_2d_ex3/generate_tolplots.py +++ /dev/null @@ -1,151 +0,0 @@ -from numpy import * -from matplotlib.pyplot import * -from pylab import * - -# Setting up local variables -tols = ['1e-0','6e-1','3e-1', - '1e-1','6e-2','3e-2', - '1e-2','6e-3','3e-3', - '1e-3','6e-4','3e-4', - '1e-4','6e-5','3e-5', - '1e-5'] - -## --------------------------------- -## Setting up vectors of amount of work done: -## --------------------------------- - -# Timing and cells from adjoint-magnitude flagging -amag_regrid = [ 168.3828 , 160.60822222, 150.9964 , 138.6031 , - 82.6206 , 47.2141 , 35.9217 , 14.6139 , - 2.756 , 0.796 , 0.3159 , 0.2619 , - 0.2592 , 0.2268 , 0.2229 , 0.2215 ] -amag_times = [ 1.58622970e+03, 1.49281189e+03, 1.34246280e+03, - 1.16435800e+03, 6.55311900e+02, 3.61060100e+02, - 2.66386500e+02, 1.03079700e+02, 1.77373000e+01, - 5.11060000e+00, 1.97270000e+00, 1.55870000e+00, - 1.49700000e+00, 1.25520000e+00, 1.24660000e+00, - 1.23730000e+00] -amag_regrid = amag_regrid[::-1] -amag_times = amag_times[::-1] - -# Timing and cells from adjoint-error flagging -aerr_regrid = [ 222.4333, 221.5709, 219.9372, 219.7524, 215.9752, 213.4921, - 213.1097, 202.8152, 104.2194, 0.3299, 0.3157, 0.2973, - 0.322 , 0.3236, 0.3348, 0.3322] -aerr_times = [ 1.65396950e+03, 1.64503170e+03, 1.63447290e+03, - 1.63348660e+03, 1.58072580e+03, 1.55252160e+03, - 1.54025240e+03, 1.39645700e+03, 6.43393800e+02, - 1.42510000e+00, 1.40220000e+00, 1.29830000e+00, - 1.38490000e+00, 1.37970000e+00, 1.42630000e+00, - 1.41970000e+00] -aerr_regrid = aerr_regrid[::-1] -aerr_times = aerr_times[::-1] - -# Timing and cells from difference-flagging -diff_regrid = [ 101.3334, 100.4365, 99.8504, 101.3604, 99.7821, 96.9741, - 94.8886, 88.6307, 80.675 , 76.9986, 49.4936, 29.1913, - 14.6271, 1.3023, 0.2764, 0.1175] -diff_times = [ 1.71767820e+03, 1.71282980e+03, 1.70584980e+03, - 1.70690280e+03, 1.71527140e+03, 1.67923300e+03, - 1.63604370e+03, 1.54464870e+03, 1.39813830e+03, - 1.29744650e+03, 7.76082100e+02, 4.37399700e+02, - 2.10789500e+02, 1.55218000e+01, 3.33560000e+00, - 1.47840000e+00] -diff_regrid = diff_regrid[::-1] -diff_times = diff_times[::-1] - -# Timing and cells from error-flagging -err_regrid = [ 207.3158, 194.2371, 182.5131, 164.8369, 103.2788, 59.2029, - 22.7793, 2.6414, 0.7793, 0.446 , 0.2799, 0.2688, - 0.2807, 0.2664, 0.2709, 0.2642] -err_times = [ 1.70307200e+03, 1.48634670e+03, 1.28684730e+03, - 1.04123930e+03, 5.06128100e+02, 2.31970000e+02, - 7.67219000e+01, 8.73310000e+00, 2.72430000e+00, - 1.68660000e+00, 1.16740000e+00, 1.12400000e+00, - 1.18160000e+00, 1.16350000e+00, 1.19180000e+00, - 1.16100000e+00] -err_regrid = err_regrid[::-1] -err_times = err_times[::-1] - -## --------------------------------- -## Setting up vectors with error in solutions: -## --------------------------------- -diff_errors_fine = [ 1.88803323e-02, 1.80756249e-02, 1.31571836e-02, - 2.44708545e-03, 6.67765406e-04, 1.13649307e-04, - 4.30387844e-05, 2.34284826e-06, 8.14584907e-06, - 5.12725708e-07, 5.49981515e-06, 3.11343017e-07, - 3.09194313e-07, 3.08818669e-07, 3.08488960e-07, - 3.08261548e-07] - -amag_errors_fine = [ 1.66309793e-02, 1.66309793e-02, 1.66309793e-02, - 1.63869141e-02, 1.63830704e-02, 1.63455552e-02, - 1.76997206e-02, 1.46060597e-02, 1.13565015e-03, - 5.12797344e-04, 2.74694425e-04, 3.40168469e-04, - 3.38413701e-05, 1.26387074e-07, 8.35988978e-06, - 0.00000000e+00] - -aerr_errors_fine = [ 1.91503471e-02, 1.91503471e-02, 1.91503471e-02, - 1.91503471e-02, 1.91503471e-02, 1.91503471e-02, - 1.91180102e-02, 8.06758744e-04, 4.41458490e-05, - 7.03947247e-06, 7.14602929e-06, 1.95601111e-06, - 1.14262929e-06, 2.76824135e-06, 2.12155203e-06, - 1.96862724e-06] - -err_errors_fine = [ 1.66309793e-02, 1.66309793e-02, 1.66309793e-02, - 1.66309793e-02, 1.66309163e-02, 1.66234431e-02, - 1.64147606e-02, 1.50676298e-02, 1.10137446e-02, - 7.27707216e-03, 2.06864462e-03, 6.11217733e-04, - 8.21253224e-05, 5.67292759e-05, 8.46845563e-06, - 4.05137217e-06] - -size = 20 -fig = figure(1,(10,8)) -axes([0.11,0.23,0.85,0.72]) -loglog(tols,err_errors_fine,'r',label='Error-Flagging',linewidth=2) -loglog(tols,diff_errors_fine,'b',label='Difference-Flagging',linewidth=2) -loglog(tols,aerr_errors_fine,'k',label='Adjoint-Error Flagging',linewidth=2) -loglog(tols,amag_errors_fine,'g',label='Adjoint-Magnitude Flagging',linewidth=2) -legend(bbox_to_anchor=(0.95,0), loc="lower right",bbox_transform=fig.transFigure, ncol=2,fontsize=size) -title("Tolerance vs. Error in J",fontsize=size) -xlabel("Tolerance",fontsize=size) -ylabel("Error in J",fontsize=size) -tick_params(axis='y',labelsize=size) -tick_params(axis='x',labelsize=size) -plt.axis([10**(-4), 10**(0), 5*10**(-4), 2*10**(-2)]) -savefig('tolvserror_2d_ex3.png') -clf() - -fig2 = figure(1,(10,8)) -axes([0.11,0.23,0.85,0.72]) -plt.semilogy(err_times[:],err_errors_fine[:],'r',label='Error-Flagging',linewidth=2) -plt.semilogy(diff_times[:],diff_errors_fine[:],'b',label='Difference-Flagging',linewidth=2) -plt.semilogy(aerr_times[:],aerr_errors_fine[:],'k',label='Adjoint-Error Flagging',linewidth=2) -plt.semilogy(amag_times,amag_errors_fine,'g',label='Adjoint-Magnitude Flagging',linewidth=2) -plt.title('Error vs. Total CPU Time', fontsize=size) -plt.legend(loc=3,fontsize=size) -xlabel("CPU Time",fontsize=size) -ylabel("Error",fontsize=size) -plt.axis([0, 800, 5*10**(-4), 2*10**(-2)]) -tick_params(axis='y',labelsize=size) -tick_params(axis='x',labelsize=size) -legend(bbox_to_anchor=(0.95,0), loc="lower right",bbox_transform=fig2.transFigure, ncol=2,fontsize=size) -plt.savefig('errvstime_2d_ex3.png') -clf() - -fig3 = figure(1,(10,8)) -axes([0.11,0.23,0.85,0.72]) -plt.semilogy(err_regrid,err_errors_fine,'r',label='Error-Flagging',linewidth=2) -plt.semilogy(diff_regrid,diff_errors_fine,'b',label='Difference-Flagging',linewidth=2) -plt.semilogy(aerr_regrid,aerr_errors_fine,'k',label='Adjoint-Error Flagging',linewidth=2) -plt.semilogy(amag_regrid,amag_errors_fine,'g',label='Adjoint-Magnitude Flagging',linewidth=2) -plt.title('Error vs. Regridding CPU Time', fontsize=size) -plt.legend(loc=3,fontsize=size) -xlabel("CPU Time",fontsize=size) -ylabel("Error",fontsize=size) -plt.axis([0, 140, 5*10**(-4), 2*10**(-2)]) -tick_params(axis='y',labelsize=size) -tick_params(axis='x',labelsize=size) -legend(bbox_to_anchor=(0.95,0), loc="lower right",bbox_transform=fig3.transFigure, ncol=2,fontsize=size) -plt.savefig('errvsregridtime_2d_ex3.png') -clf() - diff --git a/paper2_examples/acoustics_2d_ex4/Makefile b/paper2_examples/acoustics_2d_ex4_old/Makefile similarity index 90% rename from paper2_examples/acoustics_2d_ex4/Makefile rename to paper2_examples/acoustics_2d_ex4_old/Makefile index 0cd4a4a..972b578 100644 --- a/paper2_examples/acoustics_2d_ex4/Makefile +++ b/paper2_examples/acoustics_2d_ex4_old/Makefile @@ -24,14 +24,14 @@ 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 ?= -O2 -fopenmp +FFLAGS ?= # --------------------------------- # package sources for this program: # --------------------------------- -ADJLIB = $(CLAW)/amrclaw/src/2d/adjoint -include $(ADJLIB)/Makefile.adjoint_amr_2d +AMRLIB = $(CLAW)/amrclaw/src/2d +include $(AMRLIB)/Makefile.amr_2d # --------------------------------------- # package sources specifically to exclude @@ -53,8 +53,8 @@ SOURCES += \ qinit.f \ setprob.f \ setaux.f \ - $(CLAW)/riemann/src/rpn2_vc_acoustics.f \ - $(CLAW)/riemann/src/rpt2_vc_acoustics.f \ + rpn2_vc_acoustics.f90 \ + $(CLAW)/riemann/src/rpt2_vc_acoustics.f90 \ #------------------------------------------------------------------- # Include Makefile containing standard definitions and make options: diff --git a/paper2_examples/acoustics_1d_ex2/README.md b/paper2_examples/acoustics_2d_ex4_old/README.md similarity index 69% rename from paper2_examples/acoustics_1d_ex2/README.md rename to paper2_examples/acoustics_2d_ex4_old/README.md index 62d6afd..ac1dda0 100644 --- a/paper2_examples/acoustics_1d_ex2/README.md +++ b/paper2_examples/acoustics_2d_ex4_old/README.md @@ -1,4 +1,6 @@ -This code produces Example 2, the second one-dimensional variable coefficient acoustics example, presented in ''Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack''. +This code was for the original Example 4, the second two-dimensional +variable coefficient acoustics example, presented in the preprint version of +''Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack''. ### Folder Organization * **adjoint:** @@ -20,7 +22,7 @@ The code will produce two new folders: _output and _plots. The first one contains all the output files, while the latter one contains the plots and interactive visualization apps. -* Go to the main folder **acoustics_1d_ex2** and run in the terminal: +* Go to the main folder **acoustics_2d_ex4** and run in the terminal: ``` make new diff --git a/paper2_examples/acoustics_2d_ex3/adjoint/Makefile b/paper2_examples/acoustics_2d_ex4_old/adjoint/Makefile similarity index 97% rename from paper2_examples/acoustics_2d_ex3/adjoint/Makefile rename to paper2_examples/acoustics_2d_ex4_old/adjoint/Makefile index 87852de..cfa9ed9 100644 --- a/paper2_examples/acoustics_2d_ex3/adjoint/Makefile +++ b/paper2_examples/acoustics_2d_ex4_old/adjoint/Makefile @@ -53,7 +53,7 @@ SOURCES = \ qinit.f \ setprob.f \ setaux.f \ - $(CLAW)/riemann/src/rpn2_vc_acoustics_adjoint.f90 \ + rpn2_vc_acoustics_adjoint.f90 \ $(CLAW)/riemann/src/rpt2_vc_acoustics_adjoint.f90 \ #------------------------------------------------------------------- diff --git a/paper2_examples/acoustics_2d_ex4/adjoint/qinit.f b/paper2_examples/acoustics_2d_ex4_old/adjoint/qinit.f similarity index 87% rename from paper2_examples/acoustics_2d_ex4/adjoint/qinit.f rename to paper2_examples/acoustics_2d_ex4_old/adjoint/qinit.f index 490a7fe..6f2f0bb 100644 --- a/paper2_examples/acoustics_2d_ex4/adjoint/qinit.f +++ b/paper2_examples/acoustics_2d_ex4_old/adjoint/qinit.f @@ -9,6 +9,9 @@ subroutine qinit(meqn,mbc,mx,my,xlower,ylower, implicit double precision (a-h,o-z) dimension q(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc) + area = dx*dy + write(6,*) 'dx,dy,phi: ',dx,dy,1.d0/area + do 20 i=1,mx xcell = xlower + (i-0.5d0)*dx @@ -17,7 +20,7 @@ subroutine qinit(meqn,mbc,mx,my,xlower,ylower, if (dabs(xcell-3.5d0) .le. dx & .and. dabs(ycell-0.5d0) .le. dy) then - pressure = 1.d0 + pressure = 1.d0 / area else pressure = 0.d0 endif diff --git a/paper2_examples/acoustics_2d_ex4_old/adjoint/rpn2_vc_acoustics_adjoint.f90 b/paper2_examples/acoustics_2d_ex4_old/adjoint/rpn2_vc_acoustics_adjoint.f90 new file mode 100644 index 0000000..643f420 --- /dev/null +++ b/paper2_examples/acoustics_2d_ex4_old/adjoint/rpn2_vc_acoustics_adjoint.f90 @@ -0,0 +1,135 @@ +! ===================================================== +subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq) +! ===================================================== +! # Adjoint Riemann solver for the acoustics equations in 2d, with varying +! # material properties rho and kappa + +! waves: 2 +! equations: 3 +! aux fields: 2 + +! Conserved quantities: +! 1 pressure +! 2 x_momentum +! 3 y_momentum + +! Auxiliary variables: +! 1 density +! 2 sound_speed + +! # Note that although there are 3 eigenvectors, the second eigenvalue +! # is always zero and so we only need to compute 2 waves. +! # +! # solve Riemann problems along one slice of data. + +! # On input, ql contains the state vector at the left edge of each cell +! # qr contains the state vector at the right edge of each cell + +! # Here it is assumed that auxl=auxr gives the cell values. + +! # This data is along a slice in the x-direction if ixy=1 +! # or the y-direction if ixy=2. +! # On output, wave contains the waves, +! # s the speeds, +! # amdq the left-going flux difference A^- \Delta q +! # apdq the right-going flux difference A^+ \Delta q + + +! # Note that the i'th Riemann problem has left state qr(i-1,:) +! # and right state ql(i,:) +! # From the basic clawpack routines, this routine is called with ql = qr + + implicit double precision (a-h,o-z) + + dimension fwave(meqn, mwaves, 1-mbc:maxm+mbc) + dimension s(mwaves, 1-mbc:maxm+mbc) + dimension ql(meqn, 1-mbc:maxm+mbc) + dimension qr(meqn, 1-mbc:maxm+mbc) + dimension apdq(meqn, 1-mbc:maxm+mbc) + dimension amdq(meqn, 1-mbc:maxm+mbc) + dimension auxl(maux, 1-mbc:maxm+mbc) + dimension auxr(maux, 1-mbc:maxm+mbc) + +! local arrays +! ------------ + dimension delta(3) + +! # set mu to point to the component of the system that corresponds +! # to velocity in the direction of this slice, mv to the orthogonal +! # velocity: + + if (ixy == 1) then + mu = 2 + mv = 3 + else + mu = 3 + mv = 2 + endif + +! # note that notation for u and v reflects assumption that the +! # Riemann problems are in the x-direction with u in the normal +! # direciton and v in the orthogonal direcion, but with the above +! # definitions of mu and mv the routine also works with ixy=2 +! # in which case waves come from the +! # Riemann problems u_t + g(u)_y = 0 in the y-direction. + + +! # split the jump in f(q) at each interface into waves +! # The jump is split into a leftgoing wave traveling at speed -c +! # relative to the material properties to the left of the interface, +! # and a rightgoing wave traveling at speed +c +! # relative to the material properties to the right of the interface, + +! # find b1 and b2, the coefficients of the 2 eigenvectors: + do 20 i = 2-mbc, mx+mbc +! # material properties + +! # impedances: + zi = auxl(1,i)*auxl(2,i) + zim = auxr(1,i-1)*auxr(2,i-1) + +! # sound speed + ci = auxl(2,i) + cim = auxr(2,i-1) + +! # density + rhoi = auxl(1,i) + rhoim = auxr(1,i-1) + +! # bulk modulus + bulki = rhoi*ci**2 + bulkim = rhoim*cim**2 + +! # f-wave splitting + delta(1) = -(ql(mu,i)/rhoi - qr(mu,i-1)/rhoim) + delta(2) = -(bulki*ql(1,i) - bulkim*qr(1,i-1)) + + beta1 = (zi*delta(1) + delta(2)) / (zi+zim) + beta2 = (zim*delta(1) - delta(2)) / (zi+zim) + +! # Compute the waves. Eigenvectors switched for adjoint: + + fwave(1,1,i) = beta1 + fwave(mu,1,i) = beta1*zim + fwave(mv,1,i) = 0.d0 + s(1,i) = -cim + + fwave(1,2,i) = beta2 + fwave(mu,2,i) = -beta2*zi + fwave(mv,2,i) = 0.d0 + s(2,i) = ci + + 20 END DO + + +! # compute the leftgoing and rightgoing flux differences: +! # Note s(i,1) < 0 and s(i,2) > 0. + +! # f-wave splitting, so do not multiply by wave speeds: + forall (m=1:meqn, i=2-mbc: mx+mbc) + amdq(m,i) = fwave(m,1,i) + apdq(m,i) = fwave(m,2,i) + end forall + + return + end subroutine rpn2 diff --git a/paper2_examples/acoustics_2d_ex3/adjoint/setaux.f b/paper2_examples/acoustics_2d_ex4_old/adjoint/setaux.f similarity index 100% rename from paper2_examples/acoustics_2d_ex3/adjoint/setaux.f rename to paper2_examples/acoustics_2d_ex4_old/adjoint/setaux.f diff --git a/paper2_examples/acoustics_2d_ex4/adjoint/setplot.py b/paper2_examples/acoustics_2d_ex4_old/adjoint/setplot.py similarity index 97% rename from paper2_examples/acoustics_2d_ex4/adjoint/setplot.py rename to paper2_examples/acoustics_2d_ex4_old/adjoint/setplot.py index 1e3792d..90dace2 100644 --- a/paper2_examples/acoustics_2d_ex4/adjoint/setplot.py +++ b/paper2_examples/acoustics_2d_ex4_old/adjoint/setplot.py @@ -45,8 +45,8 @@ def setplot(plotdata): plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') plotitem.plot_var = 0 plotitem.pcolor_cmap = colormaps.blue_white_red - plotitem.pcolor_cmin = -0.005 - plotitem.pcolor_cmax = 0.005 + plotitem.pcolor_cmin = -0.5 + plotitem.pcolor_cmax = 0.5 plotitem.add_colorbar = False diff --git a/paper2_examples/acoustics_2d_ex3/adjoint/setprob.f b/paper2_examples/acoustics_2d_ex4_old/adjoint/setprob.f similarity index 100% rename from paper2_examples/acoustics_2d_ex3/adjoint/setprob.f rename to paper2_examples/acoustics_2d_ex4_old/adjoint/setprob.f diff --git a/paper2_examples/acoustics_2d_ex4/adjoint/setrun.py b/paper2_examples/acoustics_2d_ex4_old/adjoint/setrun.py similarity index 98% rename from paper2_examples/acoustics_2d_ex4/adjoint/setrun.py rename to paper2_examples/acoustics_2d_ex4_old/adjoint/setrun.py index 075b08f..802222f 100644 --- a/paper2_examples/acoustics_2d_ex4/adjoint/setrun.py +++ b/paper2_examples/acoustics_2d_ex4_old/adjoint/setrun.py @@ -68,8 +68,8 @@ def setrun(claw_pkg='amrclaw'): clawdata.upper[1] = 11.000000e+00 # yupper # Number of grid cells: - clawdata.num_cells[0] = 200 # mx - clawdata.num_cells[1] = 200 # my + clawdata.num_cells[0] = 160 # mx + clawdata.num_cells[1] = 120 # my # --------------- @@ -114,8 +114,8 @@ def setrun(claw_pkg='amrclaw'): 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 = 21*8 - clawdata.tfinal = 21.0 + clawdata.num_output_times = 22*8 + clawdata.tfinal = 22.0 clawdata.output_t0 = True # output at initial (or restart) time? elif clawdata.output_style == 2: diff --git a/paper2_examples/acoustics_2d_ex3/qinit.f b/paper2_examples/acoustics_2d_ex4_old/qinit.f similarity index 100% rename from paper2_examples/acoustics_2d_ex3/qinit.f rename to paper2_examples/acoustics_2d_ex4_old/qinit.f diff --git a/paper2_examples/acoustics_2d_ex4_old/rpn2_vc_acoustics.f90 b/paper2_examples/acoustics_2d_ex4_old/rpn2_vc_acoustics.f90 new file mode 100644 index 0000000..5886898 --- /dev/null +++ b/paper2_examples/acoustics_2d_ex4_old/rpn2_vc_acoustics.f90 @@ -0,0 +1,122 @@ +! ===================================================== +subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq) +! ===================================================== + +! Riemann solver for the acoustics equations in 2d, with varying +! material properties rho and kappa + +! waves: 2 +! equations: 3 +! aux fields: 2 + +! Conserved quantities: +! 1 pressure +! 2 x_momentum +! 3 y_momentum + +! Auxiliary variables: +! 1 density +! 2 sound_speed + +! Note that although there are 3 eigenvectors, the second eigenvalue +! is always zero and so we only need to compute 2 waves. +! +! solve Riemann problems along one slice of data. + +! On input, ql contains the state vector at the left edge of each cell +! qr contains the state vector at the right edge of each cell + +! Here it is assumed that auxl=auxr gives the cell values. + +! On output, wave contains the waves, +! s the speeds, +! amdq the left-going flux difference A^- \Delta q +! apdq the right-going flux difference A^+ \Delta q + +! This data is along a slice in the x-direction if ixy=1 +! or the y-direction if ixy=2. + +! Note that the i'th Riemann problem has left state qr(:,i-1) +! and right state ql(:,i) +! From the basic clawpack routines, this routine is called with ql = qr + + + implicit double precision (a-h,o-z) + + dimension wave(meqn, mwaves, 1-mbc:maxm+mbc) + dimension s(mwaves, 1-mbc:maxm+mbc) + dimension ql(meqn, 1-mbc:maxm+mbc) + dimension qr(meqn, 1-mbc:maxm+mbc) + dimension apdq(meqn, 1-mbc:maxm+mbc) + dimension amdq(meqn, 1-mbc:maxm+mbc) + dimension auxl(maux, 1-mbc:maxm+mbc) + dimension auxr(maux, 1-mbc:maxm+mbc) + +! local arrays +! ------------ + dimension delta(3) + +! # set mu to point to the component of the system that corresponds +! # to velocity in the direction of this slice, mv to the orthogonal +! # velocity. + + + if (ixy == 1) then + mu = 2 + mv = 3 + else + mu = 3 + mv = 2 + endif + +! # note that notation for u and v reflects assumption that the +! # Riemann problems are in the x-direction with u in the normal +! # direciton and v in the orthogonal direcion, but with the above +! # definitions of mu and mv the routine also works with ixy=2 + + +! # split the jump in q at each interface into waves +! # The jump is split into a leftgoing wave traveling at speed -c +! # relative to the material properties to the left of the interface, +! # and a rightgoing wave traveling at speed +c +! # relative to the material properties to the right of the interface, + +! # find a1 and a2, the coefficients of the 2 eigenvectors: + do 20 i = 2-mbc, mx+mbc + delta(1) = ql(1,i) - qr(1,i-1) + delta(2) = ql(mu,i) - qr(mu,i-1) + ! # impedances: + zi = auxl(1,i)*auxl(2,i) + zim = auxr(1,i-1)*auxr(2,i-1) + + a1 = (-delta(1) + zi*delta(2)) / (zim + zi) + a2 = (delta(1) + zim*delta(2)) / (zim + zi) + + + ! # Compute the waves. + + wave(1,1,i) = -a1*zim + wave(mu,1,i) = a1 + wave(mv,1,i) = 0.d0 + s(1,i) = -auxr(2,i-1) + + wave(1,2,i) = a2*zi + wave(mu,2,i) = a2 + wave(mv,2,i) = 0.d0 + s(2,i) = auxl(2,i) + + 20 END DO + + + +! # compute the leftgoing and rightgoing flux differences: +! # Note s(1,i) < 0 and s(2,i) > 0. + + do 220 m=1,meqn + do 220 i = 2-mbc, mx+mbc + amdq(m,i) = s(1,i)*wave(m,1,i) + apdq(m,i) = s(2,i)*wave(m,2,i) + 220 END DO + + return + end subroutine rpn2 diff --git a/paper2_examples/acoustics_2d_ex3/setaux.f b/paper2_examples/acoustics_2d_ex4_old/setaux.f similarity index 85% rename from paper2_examples/acoustics_2d_ex3/setaux.f rename to paper2_examples/acoustics_2d_ex4_old/setaux.f index b70900d..f10c236 100644 --- a/paper2_examples/acoustics_2d_ex3/setaux.f +++ b/paper2_examples/acoustics_2d_ex4_old/setaux.f @@ -11,7 +11,7 @@ subroutine setaux(mbc,mx,my,xlower,ylower, c c use amr_module, only : NEEDS_TO_BE_SET - use adjoint_module, only: innerprod_index + use adjoint_module, only: innerprod_index, adjoint_flagging implicit double precision (a-h,o-z) double precision aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) @@ -22,7 +22,8 @@ subroutine setaux(mbc,mx,my,xlower,ylower, do i = 1-mbc,mx+mbc xcell = xlower + (i-0.5d0)*dx - if (aux(1,i,j) .eq. NEEDS_TO_BE_SET) then + if (adjoint_flagging .and. + & (aux(1,i,j) .eq. NEEDS_TO_BE_SET)) then aux(innerprod_index,i,j) = 0.d0 endif diff --git a/paper2_examples/acoustics_2d_ex4/setplot.py b/paper2_examples/acoustics_2d_ex4_old/setplot.py similarity index 72% rename from paper2_examples/acoustics_2d_ex4/setplot.py rename to paper2_examples/acoustics_2d_ex4_old/setplot.py index 0639c7a..2ac322e 100644 --- a/paper2_examples/acoustics_2d_ex4/setplot.py +++ b/paper2_examples/acoustics_2d_ex4_old/setplot.py @@ -7,6 +7,9 @@ """ +outdir_fine = '_output_1600x1200' + + #-------------------------- def setplot(plotdata): #-------------------------- @@ -20,8 +23,13 @@ def setplot(plotdata): from clawpack.visclaw import colormaps + from clawpack.clawutil.data import ClawData + adjoint_data = ClawData() + adjoint_data.read('adjoint.data', force=True) + print('use_adjoint = ', adjoint_data.use_adjoint) + plotdata.clearfigures() # clear any old figures,axes,items data - plotdata.format = 'binary' # 'ascii', 'binary', 'netcdf' + plotdata.format = 'binary' # 'ascii', 'binary' # Figure for pressure @@ -42,7 +50,7 @@ def setplot(plotdata): plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') plotitem.plot_var = 0 plotitem.pcolor_cmap = colormaps.blue_white_red - plotitem.add_colorbar = False + plotitem.add_colorbar = True plotitem.show = True # show on plot? plotitem.pcolor_cmin = -0.3 plotitem.pcolor_cmax = 0.3 @@ -54,6 +62,7 @@ def setplot(plotdata): #----------------------------------------- plotfigure = plotdata.new_plotfigure(name='Inner Product', figno=1) plotfigure.kwargs = {'figsize': (5.5,4)} + plotfigure.show = adjoint_data.use_adjoint # Set up for axes in this figure: plotaxes = plotfigure.new_plotaxes() @@ -70,14 +79,35 @@ def setplot(plotdata): plotitem.pcolor_cmap = colormaps.white_red plotitem.add_colorbar = False plotitem.show = True # show on plot? - plotitem.pcolor_cmin = 0.01 # use when plotting inner product with q - #plotitem.pcolor_cmin = 0.0 # use when plotting inner product with error - plotitem.pcolor_cmax = 0.12 # use when plotting inner product with q - #plotitem.pcolor_cmax = 0.005 # use when plotting inner product with error + plotitem.pcolor_cmin = 0.0 # use when plotting inner product with q + #plotitem.pcolor_cmin = 0.0 # use when plotting inner product with error + #plotitem.pcolor_cmax = 0.12 # use when plotting inner product with q + plotitem.pcolor_cmax = 0.001 # for adjoint-error + #plotitem.pcolor_cmax = 0.0001 # for adjoint-mag plotitem.amr_patchedges_show = [0,0,0] plotitem.amr_celledges_show = [0,0,0] plotitem.amr_data_show = [1,1,1,1,0] + + + #----------------------------------------- + # Figure for innerproduct vs x + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Inner Product slice', figno=20) + #plotfigure.kwargs = {'figsize': (5.5,4)} + plotfigure.show = adjoint_data.use_adjoint + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.title = 'Inner Product vs x' + plotaxes.xlimits = [-8,8] + #plotaxes.ylimits = [-1,11] + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='1d_from_2d_data') + plotitem.map_2d_to_1d = plot_innerprod_x + plotitem.plotstyle = 'kx' + #----------------------------------------- # Figure for levels #----------------------------------------- @@ -94,8 +124,8 @@ def setplot(plotdata): # Water plotitem = plotaxes.new_plotitem(plot_type='2d_patch') - plotitem.amr_patch_bgcolor = [[1,1,1], [0.8,0.8,0.8], [0.8,1,0.8], [1,.7,.7],[0.6,0.6,1]] - plotitem.amr_patchedges_color = ['k','k','g','r','b'] + plotitem.amr_patch_bgcolor = [[1,1,1], [0.8,0.8,0.8], [0.8,1,0.8], [1,.7,.7],[0.6,0.6,1],[0.9,0.9,0]] + plotitem.amr_patchedges_color = ['k','k','g','r','b','yellow'] plotitem.amr_celledges_show = [0] plotitem.amr_patchedges_show = [0,1,1,1,1] @@ -109,14 +139,24 @@ def setplot(plotdata): plotaxes = plotfigure.new_plotaxes() plotaxes.xlimits = 'auto' plotaxes.ylimits = 'auto' - plotaxes.xlimits = [1,3] - plotaxes.ylimits = [-0.4,0.5] + plotaxes.xlimits = [0,21] + plotaxes.ylimits = [-0.4,0.6] plotaxes.title = 'Pressure' + plotaxes.afteraxes = fixup_gauge + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') plotitem.plot_var = 0 plotitem.plotstyle = 'b-' plotitem.kwargs = {'linewidth': 3} + + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.outdir = outdir_fine + plotitem.plot_var = 0 + plotitem.plotstyle = 'k-' + plotitem.kwargs = {'linewidth': 1} plotaxes.afteraxes = fixup_gauge + + # Parameters used only when creating html and/or latex hardcopy # e.g., via clawpack.visclaw.frametools.printframes: @@ -132,12 +172,21 @@ def setplot(plotdata): 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 frames in parallel with omp? return plotdata + def plot_innerprod(current_data): return current_data.aux[2,:,:] +def plot_innerprod_x(current_data): + x = current_data.x[:,0] + ipx = current_data.aux[2,:,0] + print('+++ x.shape, ipx.shape: ',x.shape, ipx.shape) + print('+++ abs(ipx).max() = ',abs(ipx).max()) + return x,ipx + # Afteraxis functions: def addgauges(current_data): @@ -165,6 +214,7 @@ def fixup_innerprod(current_data): pylab.yticks([0, 5, 10], fontsize=size) plot([0., 0.], [-1000., 1000.], 'k--') + def plot_rectangle(current_data): from clawpack.visclaw.plottools import plotbox x1 = 3.18; x2 = 3.82; y1 = 0.26; y2 = 0.74 @@ -173,7 +223,8 @@ def plot_rectangle(current_data): def fixup_gauge(current_data): import pylab - size = 34 + size = 15 + pylab.grid(True) pylab.title('Pressure at Gauge 0', fontsize=size) - pylab.xticks([1.0, 1.5, 2.0, 2.5, 3.0], fontsize=size) - pylab.yticks([-0.3, -0.1, 0.1, 0.3, 0.5], fontsize=size) + #pylab.xticks([1.0, 1.5, 2.0, 2.5, 3.0], fontsize=size) + #pylab.yticks([-0.3, -0.1, 0.1, 0.3, 0.5], fontsize=size) diff --git a/paper2_examples/acoustics_2d_ex3/setprob.f b/paper2_examples/acoustics_2d_ex4_old/setprob.f similarity index 100% rename from paper2_examples/acoustics_2d_ex3/setprob.f rename to paper2_examples/acoustics_2d_ex4_old/setprob.f diff --git a/paper2_examples/acoustics_2d_ex4/setrun.py b/paper2_examples/acoustics_2d_ex4_old/setrun.py similarity index 76% rename from paper2_examples/acoustics_2d_ex4/setrun.py rename to paper2_examples/acoustics_2d_ex4_old/setrun.py index b709d92..8d5e2e4 100644 --- a/paper2_examples/acoustics_2d_ex4/setrun.py +++ b/paper2_examples/acoustics_2d_ex4_old/setrun.py @@ -9,32 +9,16 @@ import os import numpy as np -#----------------------------------------------- -# Set these parameters for adjoint flagging.... - -# location of output from computing adjoint: -adjoint_output = os.path.abspath('adjoint/_output') -print('Will flag using adjoint solution from %s' % adjoint_output) - -# Time period of interest: -t1 = 20.5 -t2 = 21. -# Determining type of adjoint flagging: +# This version is set up for adjoint flagging +#----------------------------------------------- -# taking inner product with forward solution or Richardson error: -flag_forward_adjoint = False -flag_richardson_adjoint = True -# tolerance for adjoint flagging: -#adjoint_flag_tolerance = 3e-4 # suggested if using forward solution -adjoint_flag_tolerance = 3e-3 # suggested if using Richardson error -#----------------------------------------------- #------------------------------ def setrun(claw_pkg='amrclaw'): #------------------------------ - + """ Define the parameters used for running Clawpack. @@ -90,8 +74,8 @@ def setrun(claw_pkg='amrclaw'): clawdata.upper[1] = 11.000000e+00 # yupper # Number of grid cells: - clawdata.num_cells[0] = 50 # mx - clawdata.num_cells[1] = 50 # my + clawdata.num_cells[0] = 50 #1600 #50 # mx + clawdata.num_cells[1] = 40 #1200 #50 # my # --------------- @@ -102,7 +86,8 @@ def setrun(claw_pkg='amrclaw'): clawdata.num_eqn = 3 # Number of auxiliary variables in the aux array (initialized in setaux) - # see setadjoint + rundata.clawdata.num_aux = 2 + # Note: as required for original problem - modified below for adjoint # Index of aux array corresponding to capacity function, if there is one: clawdata.capa_index = 0 @@ -131,24 +116,24 @@ def setrun(claw_pkg='amrclaw'): # 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 + clawdata.output_style = 2 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 = 42 - clawdata.tfinal = 21.0 + clawdata.num_output_times = 44 + clawdata.tfinal = 22.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] + clawdata.output_times = [0., 2.5, 6., 15., 20.4, 22.] 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_step_interval = 1 + clawdata.total_steps = 10 clawdata.output_t0 = True # output at initial (or restart) time? @@ -258,8 +243,7 @@ def setrun(claw_pkg='amrclaw'): # --------------- rundata.gaugedata.gauges = [] # for gauges append lines of the form [gaugeno, x, y, t1, t2] - rundata.gaugedata.gauges.append([0, 3.5, 0.5, 18., 21.]) - #rundata.gaugedata.gauges.append([1, 3.6, 0.5, 2.7, 2.85]) + rundata.gaugedata.gauges.append([0, 3.5, 0.5, 0., 1e9]) # -------------- # Checkpointing: @@ -268,7 +252,7 @@ def setrun(claw_pkg='amrclaw'): # Specify when checkpoint files should be created that can be # used to restart a computation. - clawdata.checkpt_style = 0 + clawdata.checkpt_style = 1 if clawdata.checkpt_style == 0: # Do not checkpoint at all @@ -308,17 +292,12 @@ def setrun(claw_pkg='amrclaw'): # This must be a list of length num_aux, each element of which is one of: # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). - # need 1 value, set in setadjoint - - - # Flag for refinement based on Richardson error estimater: - amrdata.flag_richardson = False + rundata.amrdata.aux_type = ['center','center'] + # Note: as required for original problem - modified below for adjoint + + # set tolerances for flagging in adjoint section below - # Flag for refinement using routine flag2refine: - amrdata.flag2refine = False - # see setadjoint to set tolerance for adjoint flagging - # steps to take on each level L between regriddings of level L+1: amrdata.regrid_interval = 2 @@ -340,13 +319,43 @@ def setrun(claw_pkg='amrclaw'): rundata.regiondata.regions = [] # to specify regions of refinement append lines of the form # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] - + #------------------------------------------------------------------ # Adjoint specific data: #------------------------------------------------------------------ - rundata = setadjoint(rundata) + # Also need to set flagging method and appropriate tolerances above + + flag_method = 'adjoint-error' + + assert flag_method in ['forward-diff','forward-error',\ + 'adjoint-mag','adjoint-error'], '*** Unsupported flag_method' + + adjointdata = rundata.adjointdata + adjointdata.use_adjoint = (flag_method in ['adjoint-mag','adjoint-error']) + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = (flag_method in ['adjoint-error','forward-error']) + amrdata.flag_richardson_tol = 5. # suggested if using adjoint-error flag + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = (flag_method in ['adjoint-mag','forward-diff']) + #amrdata.flag2refine_tol = 3e-4 # suggested if using adjoint-mag flag + amrdata.flag2refine_tol = 0.05 # suggested if using forward-diff + # location of adjoint solution, must first be created: + adjointdata.adjoint_outdir = os.path.abspath('adjoint/_output') + # time period of interest: + adjointdata.t1 = 20. + adjointdata.t2 = 21. + + if adjointdata.use_adjoint: + # need an additional aux variable for inner product: + rundata.amrdata.aux_type.append('center') + rundata.clawdata.num_aux = len(rundata.amrdata.aux_type) + adjointdata.innerprod_index = len(rundata.amrdata.aux_type) + + # ----- For developers ----- # Toggle debugging print statements: amrdata.dprint = False # print domain flags @@ -365,60 +374,6 @@ def setrun(claw_pkg='amrclaw'): # end of function setrun # ---------------------- -#------------------- -def setadjoint(rundata): - #------------------- - - """ - Setting up adjoint variables and - reading in all of the checkpointed Adjoint files - """ - - import glob - - - # Set these parameters at top of this file: - # adjoint_flag_tolerance, t1, t2, adjoint_output - # Then you don't need to modify this function... - - # flag and tolerance for adjoint flagging: - if flag_forward_adjoint == True: - # setting up taking inner product with forward solution - rundata.amrdata.flag2refine = True - rundata.amrdata.flag2refine_tol = adjoint_flag_tolerance - elif flag_richardson_adjoint == True: - # setting up taking inner product with Richardson error - rundata.amrdata.flag_richardson = True - rundata.amrdata.flag_richardson_tol = adjoint_flag_tolerance - else: - print("No refinement flag set!") - - rundata.clawdata.num_aux = 3 # 3 required for adjoint flagging - rundata.amrdata.aux_type = ['center','center','center'] - - adjointdata = rundata.new_UserData(name='adjointdata',fname='adjoint.data') - adjointdata.add_param('adjoint_output',adjoint_output,'adjoint_output') - adjointdata.add_param('t1',t1,'t1, start time of interest') - adjointdata.add_param('t2',t2,'t2, final time of interest') - - files = glob.glob(os.path.join(adjoint_output,"fort.b*")) - files.sort() - - if (len(files) == 0): - print("No binary files found for adjoint output!") - - adjointdata.add_param('numadjoints', len(files), 'Number of adjoint checkpoint files.') - adjointdata.add_param('innerprod_index', 3, 'Index for innerproduct data in aux array.') - - counter = 1 - for fname in files: - f = open(fname) - adjointdata.add_param('file' + str(counter), fname, 'Binary file' + str(counter)) - counter = counter + 1 - - return rundata -# end of function setadjoint -# ---------------------- if __name__ == '__main__': # Set up run-time parameters and write all data files. diff --git a/paper2_examples/acoustics_2d_ex3/Makefile b/paper2_examples/acoustics_2d_example2/Makefile similarity index 90% rename from paper2_examples/acoustics_2d_ex3/Makefile rename to paper2_examples/acoustics_2d_example2/Makefile index 2f92075..972b578 100644 --- a/paper2_examples/acoustics_2d_ex3/Makefile +++ b/paper2_examples/acoustics_2d_example2/Makefile @@ -24,14 +24,14 @@ 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 ?= -O2 -fopenmp +FFLAGS ?= # --------------------------------- # package sources for this program: # --------------------------------- -ADJLIB = $(CLAW)/amrclaw/src/2d/adjoint -include $(ADJLIB)/Makefile.adjoint_amr_2d +AMRLIB = $(CLAW)/amrclaw/src/2d +include $(AMRLIB)/Makefile.amr_2d # --------------------------------------- # package sources specifically to exclude @@ -53,8 +53,8 @@ SOURCES += \ qinit.f \ setprob.f \ setaux.f \ - $(CLAW)/riemann/src/rpn2_vc_acoustics.f \ - $(CLAW)/riemann/src/rpt2_vc_acoustics.f \ + rpn2_vc_acoustics.f90 \ + $(CLAW)/riemann/src/rpt2_vc_acoustics.f90 \ #------------------------------------------------------------------- # Include Makefile containing standard definitions and make options: diff --git a/paper2_examples/acoustics_2d_example2/Makefile_fine b/paper2_examples/acoustics_2d_example2/Makefile_fine new file mode 100644 index 0000000..86d7e6d --- /dev/null +++ b/paper2_examples/acoustics_2d_example2/Makefile_fine @@ -0,0 +1,62 @@ + +# 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_fine.py # File containing function to make data +OUTDIR = _output_fine # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots_fine # 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 \ + setaux.f \ + rpn2_vc_acoustics.f90 \ + $(CLAW)/riemann/src/rpt2_vc_acoustics.f90 \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + diff --git a/paper2_examples/acoustics_2d_example2/README.md b/paper2_examples/acoustics_2d_example2/README.md new file mode 100644 index 0000000..ea71d98 --- /dev/null +++ b/paper2_examples/acoustics_2d_example2/README.md @@ -0,0 +1,55 @@ +This code produces Example 2, the two-dimensional variable coefficient acoustics example, presented in ''Analysis and Performance Evaluation of Adjoint-Guided Adaptive Mesh Refinement Using Clawpack''. + +### Folder Organization +* **adjoint:** + +Contains code to solve the adjoint problem. + +The start and end times specified in this directory should agree with those for the forward code. + +### Running the Code + +* Go to the folder **adjoint** and run in a terminal: + +``` +make new +make .plots +``` + +The code will produce two new folders: _output and _plots. +The first one contains all the output files, while the latter one contains the plots and interactive +visualization apps. + +* Go to the main folder **acoustics_2d_example3** and run in the terminal: + +``` +make new +make .plots +``` + +The code will produce two new folders: _output and _plots. +The first one contains all the output files, while the latter one contains the plots and interactive +visualization apps. + +### Extra Python Scripts + +To run all the tests reported in the paper: + + # first run the code in adjoint subdirectory if you haven't already + + make .output -f Makefile_fine # create reference solution + python run_tests.py # runs all tests + +The results will be summarized in timings_errors.txt, listing both CPU +times and error estimates over [t1, t2]. + +The script make_cpu_vs_error_plot.py creates the performance plot in the +paper. + +To make the gauge plots, the script plotgauge.py can be used. This requires +running the fine grid code and also the adjoint-error code with at least 3 +tolerances and saving the results into directories named as indicated in +the script. + + + diff --git a/paper2_examples/acoustics_2d_ex4/adjoint/Makefile b/paper2_examples/acoustics_2d_example2/adjoint/Makefile similarity index 97% rename from paper2_examples/acoustics_2d_ex4/adjoint/Makefile rename to paper2_examples/acoustics_2d_example2/adjoint/Makefile index 87852de..cfa9ed9 100644 --- a/paper2_examples/acoustics_2d_ex4/adjoint/Makefile +++ b/paper2_examples/acoustics_2d_example2/adjoint/Makefile @@ -53,7 +53,7 @@ SOURCES = \ qinit.f \ setprob.f \ setaux.f \ - $(CLAW)/riemann/src/rpn2_vc_acoustics_adjoint.f90 \ + rpn2_vc_acoustics_adjoint.f90 \ $(CLAW)/riemann/src/rpt2_vc_acoustics_adjoint.f90 \ #------------------------------------------------------------------- diff --git a/paper2_examples/acoustics_2d_ex3/adjoint/qinit.f b/paper2_examples/acoustics_2d_example2/adjoint/qinit.f similarity index 100% rename from paper2_examples/acoustics_2d_ex3/adjoint/qinit.f rename to paper2_examples/acoustics_2d_example2/adjoint/qinit.f diff --git a/paper2_examples/acoustics_2d_example2/adjoint/rpn2_vc_acoustics_adjoint.f90 b/paper2_examples/acoustics_2d_example2/adjoint/rpn2_vc_acoustics_adjoint.f90 new file mode 100644 index 0000000..643f420 --- /dev/null +++ b/paper2_examples/acoustics_2d_example2/adjoint/rpn2_vc_acoustics_adjoint.f90 @@ -0,0 +1,135 @@ +! ===================================================== +subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq) +! ===================================================== +! # Adjoint Riemann solver for the acoustics equations in 2d, with varying +! # material properties rho and kappa + +! waves: 2 +! equations: 3 +! aux fields: 2 + +! Conserved quantities: +! 1 pressure +! 2 x_momentum +! 3 y_momentum + +! Auxiliary variables: +! 1 density +! 2 sound_speed + +! # Note that although there are 3 eigenvectors, the second eigenvalue +! # is always zero and so we only need to compute 2 waves. +! # +! # solve Riemann problems along one slice of data. + +! # On input, ql contains the state vector at the left edge of each cell +! # qr contains the state vector at the right edge of each cell + +! # Here it is assumed that auxl=auxr gives the cell values. + +! # This data is along a slice in the x-direction if ixy=1 +! # or the y-direction if ixy=2. +! # On output, wave contains the waves, +! # s the speeds, +! # amdq the left-going flux difference A^- \Delta q +! # apdq the right-going flux difference A^+ \Delta q + + +! # Note that the i'th Riemann problem has left state qr(i-1,:) +! # and right state ql(i,:) +! # From the basic clawpack routines, this routine is called with ql = qr + + implicit double precision (a-h,o-z) + + dimension fwave(meqn, mwaves, 1-mbc:maxm+mbc) + dimension s(mwaves, 1-mbc:maxm+mbc) + dimension ql(meqn, 1-mbc:maxm+mbc) + dimension qr(meqn, 1-mbc:maxm+mbc) + dimension apdq(meqn, 1-mbc:maxm+mbc) + dimension amdq(meqn, 1-mbc:maxm+mbc) + dimension auxl(maux, 1-mbc:maxm+mbc) + dimension auxr(maux, 1-mbc:maxm+mbc) + +! local arrays +! ------------ + dimension delta(3) + +! # set mu to point to the component of the system that corresponds +! # to velocity in the direction of this slice, mv to the orthogonal +! # velocity: + + if (ixy == 1) then + mu = 2 + mv = 3 + else + mu = 3 + mv = 2 + endif + +! # note that notation for u and v reflects assumption that the +! # Riemann problems are in the x-direction with u in the normal +! # direciton and v in the orthogonal direcion, but with the above +! # definitions of mu and mv the routine also works with ixy=2 +! # in which case waves come from the +! # Riemann problems u_t + g(u)_y = 0 in the y-direction. + + +! # split the jump in f(q) at each interface into waves +! # The jump is split into a leftgoing wave traveling at speed -c +! # relative to the material properties to the left of the interface, +! # and a rightgoing wave traveling at speed +c +! # relative to the material properties to the right of the interface, + +! # find b1 and b2, the coefficients of the 2 eigenvectors: + do 20 i = 2-mbc, mx+mbc +! # material properties + +! # impedances: + zi = auxl(1,i)*auxl(2,i) + zim = auxr(1,i-1)*auxr(2,i-1) + +! # sound speed + ci = auxl(2,i) + cim = auxr(2,i-1) + +! # density + rhoi = auxl(1,i) + rhoim = auxr(1,i-1) + +! # bulk modulus + bulki = rhoi*ci**2 + bulkim = rhoim*cim**2 + +! # f-wave splitting + delta(1) = -(ql(mu,i)/rhoi - qr(mu,i-1)/rhoim) + delta(2) = -(bulki*ql(1,i) - bulkim*qr(1,i-1)) + + beta1 = (zi*delta(1) + delta(2)) / (zi+zim) + beta2 = (zim*delta(1) - delta(2)) / (zi+zim) + +! # Compute the waves. Eigenvectors switched for adjoint: + + fwave(1,1,i) = beta1 + fwave(mu,1,i) = beta1*zim + fwave(mv,1,i) = 0.d0 + s(1,i) = -cim + + fwave(1,2,i) = beta2 + fwave(mu,2,i) = -beta2*zi + fwave(mv,2,i) = 0.d0 + s(2,i) = ci + + 20 END DO + + +! # compute the leftgoing and rightgoing flux differences: +! # Note s(i,1) < 0 and s(i,2) > 0. + +! # f-wave splitting, so do not multiply by wave speeds: + forall (m=1:meqn, i=2-mbc: mx+mbc) + amdq(m,i) = fwave(m,1,i) + apdq(m,i) = fwave(m,2,i) + end forall + + return + end subroutine rpn2 diff --git a/paper2_examples/acoustics_2d_ex4/adjoint/setaux.f b/paper2_examples/acoustics_2d_example2/adjoint/setaux.f similarity index 100% rename from paper2_examples/acoustics_2d_ex4/adjoint/setaux.f rename to paper2_examples/acoustics_2d_example2/adjoint/setaux.f diff --git a/paper2_examples/acoustics_2d_ex3/adjoint/setplot.py b/paper2_examples/acoustics_2d_example2/adjoint/setplot.py similarity index 93% rename from paper2_examples/acoustics_2d_ex3/adjoint/setplot.py rename to paper2_examples/acoustics_2d_example2/adjoint/setplot.py index 76bf86c..039357f 100644 --- a/paper2_examples/acoustics_2d_ex3/adjoint/setplot.py +++ b/paper2_examples/acoustics_2d_example2/adjoint/setplot.py @@ -27,10 +27,10 @@ def setplot(plotdata): plotdata.format = 'binary' # 'ascii', 'binary', 'netcdf' - # Figure for pressure + # Figure for adjoint # ------------------- - plotfigure = plotdata.new_plotfigure(name='Pressure', figno=0) + plotfigure = plotdata.new_plotfigure(name='Adjoint', figno=0) plotfigure.kwargs = {'figsize': (5.5,4)} # Set up for axes in this figure: @@ -79,6 +79,7 @@ def fixup_adjoint(current_data): def plot_rectangle(current_data): from clawpack.visclaw.plottools import plotbox - x1 = 0.68; x2 = 1.32; y1 = 5.26; y2 = 5.74 + #x1 = 0.68; x2 = 1.32; y1 = 5.26; y2 = 5.74 + x1 = 0.9; x2 = 1.1; y1 = 5.4; y2 = 5.6 xy = [x1, x2, y1, y2] plotbox(xy, kwargs={'color': 'k', 'linewidth': 2}) diff --git a/paper2_examples/acoustics_2d_ex4/adjoint/setprob.f b/paper2_examples/acoustics_2d_example2/adjoint/setprob.f similarity index 100% rename from paper2_examples/acoustics_2d_ex4/adjoint/setprob.f rename to paper2_examples/acoustics_2d_example2/adjoint/setprob.f diff --git a/paper2_examples/acoustics_2d_ex3/adjoint/setrun.py b/paper2_examples/acoustics_2d_example2/adjoint/setrun.py similarity index 98% rename from paper2_examples/acoustics_2d_ex3/adjoint/setrun.py rename to paper2_examples/acoustics_2d_example2/adjoint/setrun.py index 075b08f..802222f 100644 --- a/paper2_examples/acoustics_2d_ex3/adjoint/setrun.py +++ b/paper2_examples/acoustics_2d_example2/adjoint/setrun.py @@ -68,8 +68,8 @@ def setrun(claw_pkg='amrclaw'): clawdata.upper[1] = 11.000000e+00 # yupper # Number of grid cells: - clawdata.num_cells[0] = 200 # mx - clawdata.num_cells[1] = 200 # my + clawdata.num_cells[0] = 160 # mx + clawdata.num_cells[1] = 120 # my # --------------- @@ -114,8 +114,8 @@ def setrun(claw_pkg='amrclaw'): 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 = 21*8 - clawdata.tfinal = 21.0 + clawdata.num_output_times = 22*8 + clawdata.tfinal = 22.0 clawdata.output_t0 = True # output at initial (or restart) time? elif clawdata.output_style == 2: diff --git a/paper2_examples/acoustics_2d_example2/collect_results.py b/paper2_examples/acoustics_2d_example2/collect_results.py new file mode 100644 index 0000000..2238975 --- /dev/null +++ b/paper2_examples/acoustics_2d_example2/collect_results.py @@ -0,0 +1,131 @@ + +""" +Initial tests of various ways to plot errors. +For version used in paper, see make_cpu_vs_error_plot.py +""" + +from pylab import * +import clawpack.pyclaw.gauges as gauges +from scipy.interpolate import interp1d + +fs = 15 # fontsize + +gaugeno=0 +t1 = 20.5 +t2 = 21.5 + +truns = ['1600x1200','0025','005','01','02', + 'lte0025','lte005','lte01','lte02', + 'fb0025','fb005','fb01','fb02', + '800x600','400x300','1200x900','600x450','300x224','150x112'] +outdirs = ['_output_%s' % trun for trun in truns] + +tunif = linspace(t1,t2,10000) + +punif = {} +perr = {} +cpu = {} + +for k in range(len(truns)): + outdir = outdirs[k] + trun = truns[k] + + timings = open('%s/timing.txt' % outdir).readlines() + for line in timings: + if 'Total time:' in line: + tokens = line.split() + cpu[k] = float(tokens[-1]) + print('%s: CPU = %s' % (trun,tokens[-1])) + gauge = gauges.GaugeSolution(gaugeno, outdir) + t = gauge.t + q = gauge.q + p = q[0,:] + j = where(logical_and(t>t1-0.1,t0: + pfine = punif[truns[0]] + perr[trun] = abs(punif[trun] - pfine).max() + print('%s: diff = %.4f' % (trun,perr[trun])) + +if 1: + tol = array([.0025,.005,.01,.02]) + cpu_lte = [1329,581,225,63] + cpu_a = [395+5, 123+5,17+5,2+5] #means #[407+5,125+5,17+5,2+5] + err_a = [0.0027,.0083,.0354,.0837] + err_lte = [.0031,.0067,.0122,.0347] + cpu_fine = 2823.5 + + cpu_fa = [1031,559,283,140] + err_fa = [0.0017,0.0054,0.0116,0.0171] + cpu_fb = [702,357,193,82] # means #[756,384,189,84] + err_fb = [0.0034,0.0087,0.0124,0.0294] + + cpu_unif = [1097,345,43] #means #[1212,339,37] + err_unif = [0.0037,0.0126,0.0445] + #cpu_unif = [1212,152,15,2] + #err_unif = [0.0037,0.0232,0.0622,0.0867] + + + figure(21) + clf() + loglog(tol,err_lte,'bo-',label='forward error flagging') + loglog(tol/1.5,err_fb,'go-',label='forward 1-step flagging') + loglog(tol,err_a,'rs-',label='adjoint error flagging') + axis('scaled'); grid(True) + axis([1e-3,1e-1,1e-3,1e-1]) + loglog([1e-3,0.1],[1e-3,0.1],'k--') + tick_params(axis='both', which='major', labelsize=fs) + xlabel('tolerance specified',fontsize=fs) + ylabel('max error in p over [t1,t2]',fontsize=fs) + legend(loc='lower right',framealpha=1,fontsize=12) + if 0: + fname = 'error_plot.png' + savefig(fname, bbox_inches='tight') + print('Created ',fname) + + figure(22, figsize=(8,5)) + clf() + semilogx(tol,cpu_fine*ones(len(tol)),'g--',lw=3,label='uniform fine grid') + semilogx(tol,cpu_lte,'bo-',label='forward error flagging') + semilogx(tol,cpu_fb,'go-',label='forward 1-step flagging') + semilogx(tol,cpu_a,'rs-',label='adjoint error flagging') + + axis([1e-3,0.1,0,3000]) + tick_params(axis='both', which='major', labelsize=fs) + #text(tol[0],cpu_fine+20, 'uniform fine grid',fontsize=fs) + grid(True) + xlabel('tolerance specified',fontsize=fs) + ylabel('CPU time (seconds)',fontsize=fs) + legend(loc='upper right',framealpha=1,fontsize=12) + if 0: + fname = 'cpu_plot.png' + savefig(fname, bbox_inches='tight') + print('Created ',fname) + + figure(23, figsize=(8,5)) + clf() + #semilogx(tol,cpu_fine*ones(len(tol)),'g--',lw=3,label='uniform fine grid') + loglog(err_unif,cpu_unif,'g^-',label='uniform grids') + #loglog(err_lte,cpu_lte,'bo-',label='forward error flagging') + loglog(err_fb,cpu_fb,'bo-',label='forward error flagging') + loglog(err_a,cpu_a,'rs-',label='adjoint error flagging') + loglog([1e-3,1e-1],[0.1*(1e-3)**(-1.5), 0.1* (1e-1)**(-1.5)],'k--', + label='reference line') + + xlim(1e-3,0.1) + #ylim(0,3000) + tick_params(axis='both', which='major', labelsize=fs) + #text(tol[0],cpu_fine+20, 'uniform fine grid',fontsize=fs) + grid(True) + xlabel('error achieved',fontsize=fs) + ylabel('CPU time (seconds)',fontsize=fs) + legend(loc='upper right',framealpha=1,fontsize=12) + + if 0: + fname = 'cpu_vs_error.png' + savefig(fname, bbox_inches='tight') + print('Created ',fname) \ No newline at end of file diff --git a/paper2_examples/acoustics_2d_example2/make_cpu_vs_error_plot.py b/paper2_examples/acoustics_2d_example2/make_cpu_vs_error_plot.py new file mode 100644 index 0000000..a724338 --- /dev/null +++ b/paper2_examples/acoustics_2d_example2/make_cpu_vs_error_plot.py @@ -0,0 +1,37 @@ + +from pylab import * + +fs = 15 # fontsize + +# adjoint-error: +cpu_a = [395+5, 123+5,17+5,2+5] #means + adjoint time +err_a = [0.0027,.0083,.0354,.0837] + +# forward-error: +cpu_fb = [702,357,193,82] # means +err_fb = [0.0034,0.0087,0.0124,0.0294] + +# uniform: +cpu_unif = [1097,345,43] # means +err_unif = [0.0037,0.0126,0.0445] + +figure(23, figsize=(8,5)) +clf() +loglog(err_unif,cpu_unif,'g^-',label='uniform grids') +loglog(err_fb,cpu_fb,'bo-',label='forward error flagging') +loglog(err_a,cpu_a,'rs-',label='adjoint error flagging') +loglog([1e-3,1e-1],[0.1*(1e-3)**(-1.5), 0.1* (1e-1)**(-1.5)],'k--', + label='reference line') + +xlim(1e-3,0.1) +#ylim(0,3000) +tick_params(axis='both', which='major', labelsize=fs) +grid(True) +xlabel('error achieved',fontsize=fs) +ylabel('CPU time (seconds)',fontsize=fs) +legend(loc='upper right',framealpha=1,fontsize=12) + +if 0: + fname = 'cpu_vs_error.png' + savefig(fname, bbox_inches='tight') + print('Created ',fname) diff --git a/paper2_examples/acoustics_2d_example2/plotgauge.py b/paper2_examples/acoustics_2d_example2/plotgauge.py new file mode 100644 index 0000000..bccd3d6 --- /dev/null +++ b/paper2_examples/acoustics_2d_example2/plotgauge.py @@ -0,0 +1,54 @@ + +from pylab import * +import clawpack.pyclaw.gauges as gauges + +gaugeno=0 +t1 = 20.5 +t2 = 21.5 + +labels = ['Fine grid','tol = 0.005','tol = 0.01','tol = 0.02','tol = 0.0025'] +outdirs = ['_output_fine', '_output_005', '_output_01', + '_output_02', '_output_0025'] +colors = ['r','k','m','b','m'] +lws = [2.5,1.5,1.5,1.5,1.5] + +figure(400, figsize=(11,8)) +clf() +ax1 = subplot(211) +ax2 = subplot(212) + +for k in range(4): + outdir = outdirs[k] + gauge = gauges.GaugeSolution(gaugeno, outdir) + t = gauge.t + q = gauge.q + p = q[0,:] + + ax1.plot(t, p, color=colors[k], lw=lws[k], label=labels[k]) + ax2.plot(t, p, color=colors[k], lw=lws[k], label=labels[k]) + + pphi = where(logical_and(t>t1,t 0. + + do 220 m=1,meqn + do 220 i = 2-mbc, mx+mbc + amdq(m,i) = s(1,i)*wave(m,1,i) + apdq(m,i) = s(2,i)*wave(m,2,i) + 220 END DO + + return + end subroutine rpn2 diff --git a/paper2_examples/acoustics_2d_example2/run_tests.py b/paper2_examples/acoustics_2d_example2/run_tests.py new file mode 100644 index 0000000..713ea3f --- /dev/null +++ b/paper2_examples/acoustics_2d_example2/run_tests.py @@ -0,0 +1,177 @@ + +from pylab import * +import clawpack.pyclaw.gauges as gauges +from scipy.interpolate import interp1d +from clawpack.clawutil import runclaw +import setrun_cases + +# number of reps to run for each case for timing: +num_reps = 5 + +# gauge and time interval for error estimates: +gaugeno=0 +t1 = 20.5 +t2 = 21.5 + +# direct output from all runs to the same temp output directory: +outdir = '_output_temp' + + +# Read in fine grid reference results for estimating errors +# Assumes code was already run with uniform fine grid, via +# make .output -f Makefile_fine + +gauge = gauges.GaugeSolution(gaugeno, '_output_fine') +t = gauge.t +q = gauge.q +p = q[0,:] +j = where(logical_and(t>t1-0.1,tt1-0.1,t 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] = 'wall' # at xlower + clawdata.bc_upper[0] = 'wall' # at xupper + + clawdata.bc_lower[1] = 'extrap' # at ylower + clawdata.bc_upper[1] = 'wall' # at yupper + + + # --------------- + # Gauges: + # --------------- + rundata.gaugedata.gauges = [] + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + rundata.gaugedata.gauges.append([0, 1.0, 5.5, 0., 1e9]) + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 1 + + 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 = 6 + + # List of refinement ratios at each level (length at least amr_level_max-1) + amrdata.refinement_ratios_x = [2,2,2,2,2,2,2,2,2] + amrdata.refinement_ratios_y = [2,2,2,2,2,2,2,2,2] + amrdata.refinement_ratios_t = [2,2,2,2,2,2,2,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). + + rundata.amrdata.aux_type = ['center','center'] + # Note: as required for original problem - modified below for adjoint + + # set tolerances for flagging in adjoint section below + + + # 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] + + #------------------------------------------------------------------ + # Adjoint specific data: + #------------------------------------------------------------------ + # Also need to set flagging method and appropriate tolerances above + + flag_method = 'adjoint-error' + + assert flag_method in ['forward-diff','forward-error',\ + 'adjoint-mag','adjoint-error'], '*** Unsupported flag_method' + + adjointdata = rundata.adjointdata + adjointdata.use_adjoint = (flag_method in ['adjoint-mag','adjoint-error']) + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = (flag_method in ['adjoint-error','forward-error']) + amrdata.flag_richardson_tol = 0.005 # suggested if using adjoint-error flag + amrdata.flag_richardson_tol = 0.005*0.015 # suggested if using forward-error flag + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = (flag_method in ['adjoint-mag','forward-diff']) + amrdata.flag2refine_tol = 0.001 # suggested if using adjoint-mag flag + #amrdata.flag2refine_tol = 0.05 # suggested if using forward-diff + + # location of adjoint solution, must first be created: + adjointdata.adjoint_outdir = os.path.abspath('adjoint/_output') + + # time period of interest: + adjointdata.t1 = 20.5 + adjointdata.t2 = 21.5 + + if adjointdata.use_adjoint: + # need an additional aux variable for inner product: + rundata.amrdata.aux_type.append('center') + rundata.clawdata.num_aux = len(rundata.amrdata.aux_type) + adjointdata.innerprod_index = len(rundata.amrdata.aux_type) + + + # ----- 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/paper2_examples/acoustics_2d_ex3/setrun.py b/paper2_examples/acoustics_2d_example2/setrun_cases.py similarity index 76% rename from paper2_examples/acoustics_2d_ex3/setrun.py rename to paper2_examples/acoustics_2d_example2/setrun_cases.py index 83712f2..bc57869 100644 --- a/paper2_examples/acoustics_2d_ex3/setrun.py +++ b/paper2_examples/acoustics_2d_example2/setrun_cases.py @@ -9,32 +9,16 @@ import os import numpy as np -#----------------------------------------------- -# Set these parameters for adjoint flagging.... - -# location of output from computing adjoint: -adjoint_output = os.path.abspath('adjoint/_output') -print('Will flag using adjoint solution from %s' % adjoint_output) - -# Time period of interest: -t1 = 20.5 -t2 = 21. -# Determining type of adjoint flagging: +# This version is set up for adjoint flagging +#----------------------------------------------- -# taking inner product with forward solution or Richardson error: -flag_forward_adjoint = True -flag_richardson_adjoint = False -# tolerance for adjoint flagging: -adjoint_flag_tolerance = 3e-4 # suggested if using forward solution -#adjoint_flag_tolerance = 3e-3 # suggested if using Richardson error -#----------------------------------------------- #------------------------------ def setrun(claw_pkg='amrclaw'): #------------------------------ - + """ Define the parameters used for running Clawpack. @@ -90,8 +74,8 @@ def setrun(claw_pkg='amrclaw'): clawdata.upper[1] = 11.000000e+00 # yupper # Number of grid cells: - clawdata.num_cells[0] = 50 # mx - clawdata.num_cells[1] = 50 # my + clawdata.num_cells[0] = 50 # mx + clawdata.num_cells[1] = 40 # my # --------------- @@ -102,7 +86,8 @@ def setrun(claw_pkg='amrclaw'): clawdata.num_eqn = 3 # Number of auxiliary variables in the aux array (initialized in setaux) - # see setadjoint + rundata.clawdata.num_aux = 2 + # Note: as required for original problem - modified below for adjoint # Index of aux array corresponding to capacity function, if there is one: clawdata.capa_index = 0 @@ -136,19 +121,19 @@ def setrun(claw_pkg='amrclaw'): 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 = 42 - clawdata.tfinal = 21.0 + clawdata.num_output_times = 1 + clawdata.tfinal = 22.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] + clawdata.output_times = [0., 2.5, 6., 15., 21., 22.] 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_step_interval = 1 + clawdata.total_steps = 10 clawdata.output_t0 = True # output at initial (or restart) time? @@ -180,7 +165,7 @@ def setrun(claw_pkg='amrclaw'): # Initial time step for variable dt. # (If dt_variable==0 then dt=dt_initial for all steps) - clawdata.dt_initial = 1.00000e-03 + clawdata.dt_initial = 0.13 # Max time step to be allowed if variable dt used: clawdata.dt_max = 1.000000e+99 @@ -258,8 +243,7 @@ def setrun(claw_pkg='amrclaw'): # --------------- rundata.gaugedata.gauges = [] # for gauges append lines of the form [gaugeno, x, y, t1, t2] - rundata.gaugedata.gauges.append([0, 1.0, 5.5, 18., 21.]) - #rundata.gaugedata.gauges.append([1, 3.6, 0.5, 2.7, 2.85]) + rundata.gaugedata.gauges.append([0, 1.0, 5.5, 0., 1e9]) # -------------- # Checkpointing: @@ -268,7 +252,7 @@ def setrun(claw_pkg='amrclaw'): # Specify when checkpoint files should be created that can be # used to restart a computation. - clawdata.checkpt_style = 0 + clawdata.checkpt_style = 1 if clawdata.checkpt_style == 0: # Do not checkpoint at all @@ -296,7 +280,7 @@ def setrun(claw_pkg='amrclaw'): amrdata = rundata.amrdata # max number of refinement levels: - amrdata.amr_levels_max = 5 + amrdata.amr_levels_max = 6 # List of refinement ratios at each level (length at least amr_level_max-1) amrdata.refinement_ratios_x = [2,2,2,2,2,2,2,2,2] @@ -308,17 +292,12 @@ def setrun(claw_pkg='amrclaw'): # This must be a list of length num_aux, each element of which is one of: # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). - # need 1 value, set in setadjoint - - - # Flag for refinement based on Richardson error estimater: - amrdata.flag_richardson = False + rundata.amrdata.aux_type = ['center','center'] + # Note: as required for original problem - modified below for adjoint + + # set tolerances for flagging in adjoint section below - # Flag for refinement using routine flag2refine: - amrdata.flag2refine = False - # see setadjoint to set tolerance for adjoint flagging - # steps to take on each level L between regriddings of level L+1: amrdata.regrid_interval = 2 @@ -340,13 +319,43 @@ def setrun(claw_pkg='amrclaw'): rundata.regiondata.regions = [] # to specify regions of refinement append lines of the form # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] - + #------------------------------------------------------------------ # Adjoint specific data: #------------------------------------------------------------------ - rundata = setadjoint(rundata) + # Also need to set flagging method and appropriate tolerances above + + flag_method = 'forward-error' + + assert flag_method in ['forward-diff','forward-error',\ + 'adjoint-mag','adjoint-error'], '*** Unsupported flag_method' + + adjointdata = rundata.adjointdata + adjointdata.use_adjoint = (flag_method in ['adjoint-mag','adjoint-error']) + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = (flag_method in ['adjoint-error','forward-error']) + amrdata.flag_richardson_tol = 0.02*0.015 # suggested if using adjoint-error flag + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = (flag_method in ['adjoint-mag','forward-diff']) + amrdata.flag2refine_tol = 0.001 # suggested if using adjoint-mag flag + #amrdata.flag2refine_tol = 0.05 # suggested if using forward-diff + # location of adjoint solution, must first be created: + adjointdata.adjoint_outdir = os.path.abspath('adjoint/_output') + # time period of interest: + adjointdata.t1 = 20.5 + adjointdata.t2 = 21.5 + + if adjointdata.use_adjoint: + # need an additional aux variable for inner product: + rundata.amrdata.aux_type.append('center') + rundata.clawdata.num_aux = len(rundata.amrdata.aux_type) + adjointdata.innerprod_index = len(rundata.amrdata.aux_type) + + # ----- For developers ----- # Toggle debugging print statements: amrdata.dprint = False # print domain flags @@ -365,60 +374,6 @@ def setrun(claw_pkg='amrclaw'): # end of function setrun # ---------------------- -#------------------- -def setadjoint(rundata): - #------------------- - - """ - Setting up adjoint variables and - reading in all of the checkpointed Adjoint files - """ - - import glob - - - # Set these parameters at top of this file: - # adjoint_flag_tolerance, t1, t2, adjoint_output - # Then you don't need to modify this function... - - # flag and tolerance for adjoint flagging: - if flag_forward_adjoint == True: - # setting up taking inner product with forward solution - rundata.amrdata.flag2refine = True - rundata.amrdata.flag2refine_tol = adjoint_flag_tolerance - elif flag_richardson_adjoint == True: - # setting up taking inner product with Richardson error - rundata.amrdata.flag_richardson = True - rundata.amrdata.flag_richardson_tol = adjoint_flag_tolerance - else: - print("No refinement flag set!") - - rundata.clawdata.num_aux = 3 # 3 required for adjoint flagging - rundata.amrdata.aux_type = ['center','center','center'] - - adjointdata = rundata.new_UserData(name='adjointdata',fname='adjoint.data') - adjointdata.add_param('adjoint_output',adjoint_output,'adjoint_output') - adjointdata.add_param('t1',t1,'t1, start time of interest') - adjointdata.add_param('t2',t2,'t2, final time of interest') - - files = glob.glob(os.path.join(adjoint_output,"fort.b*")) - files.sort() - - if (len(files) == 0): - print("No binary files found for adjoint output!") - - adjointdata.add_param('numadjoints', len(files), 'Number of adjoint checkpoint files.') - adjointdata.add_param('innerprod_index', 3, 'Index for innerproduct data in aux array.') - - counter = 1 - for fname in files: - f = open(fname) - adjointdata.add_param('file' + str(counter), fname, 'Binary file' + str(counter)) - counter = counter + 1 - - return rundata -# end of function setadjoint -# ---------------------- if __name__ == '__main__': # Set up run-time parameters and write all data files. diff --git a/paper2_examples/acoustics_2d_example2/setrun_fine.py b/paper2_examples/acoustics_2d_example2/setrun_fine.py new file mode 100644 index 0000000..595a4c0 --- /dev/null +++ b/paper2_examples/acoustics_2d_example2/setrun_fine.py @@ -0,0 +1,383 @@ +""" +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. + +""" + +import os +import numpy as np + + +# This version is set up for adjoint flagging +#----------------------------------------------- + + + +#------------------------------ +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('p1', 2., 'density on left') + probdata.add_param('c1', 2., 'sound speed on left') + probdata.add_param('p2', 2., 'density on right') + probdata.add_param('c2', 1., 'sound speed on right') + + #------------------------------------------------------------------ + # 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] = -8.000000e+00 # xlower + clawdata.upper[0] = 8.000000e+00 # xupper + clawdata.lower[1] = -1.000000e+00 # ylower + clawdata.upper[1] = 11.000000e+00 # yupper + + # Number of grid cells: + clawdata.num_cells[0] = 1600 # mx + clawdata.num_cells[1] = 1200 # my + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + rundata.clawdata.num_aux = 2 + # Note: as required for original problem - modified below for adjoint + + # 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 = 2 + + 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 = 1 + clawdata.tfinal = 22.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., 2.5, 6., 15., 21., 22.] + + elif clawdata.output_style == 3: + # Output every step_interval timesteps over total_steps timesteps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 10 + clawdata.output_t0 = True # output at initial (or restart) time? + + + clawdata.output_format = 'binary' # 'ascii', 'binary', 'netcdf' + + clawdata.output_q_components = 'all' # could be list such as [True,True] + clawdata.output_aux_components = 'all' # could be list + clawdata.output_aux_onlyonce = False # 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 = 0.13 / 32. + + # 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] = 'wall' # at xlower + clawdata.bc_upper[0] = 'wall' # at xupper + + clawdata.bc_lower[1] = 'extrap' # at ylower + clawdata.bc_upper[1] = 'wall' # at yupper + + + # --------------- + # Gauges: + # --------------- + rundata.gaugedata.gauges = [] + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + rundata.gaugedata.gauges.append([0, 1.0, 5.5, 0., 1e9]) + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 1 + + 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 = 1 + + # List of refinement ratios at each level (length at least amr_level_max-1) + amrdata.refinement_ratios_x = [2,2,2,2,2,2,2,2,2] + amrdata.refinement_ratios_y = [2,2,2,2,2,2,2,2,2] + amrdata.refinement_ratios_t = [2,2,2,2,2,2,2,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). + + rundata.amrdata.aux_type = ['center','center'] + # Note: as required for original problem - modified below for adjoint + + # set tolerances for flagging in adjoint section below + + + # 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] + + #------------------------------------------------------------------ + # Adjoint specific data: + #------------------------------------------------------------------ + # Also need to set flagging method and appropriate tolerances above + + flag_method = 'forward-error' + + assert flag_method in ['forward-diff','forward-error',\ + 'adjoint-mag','adjoint-error'], '*** Unsupported flag_method' + + adjointdata = rundata.adjointdata + adjointdata.use_adjoint = (flag_method in ['adjoint-mag','adjoint-error']) + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = (flag_method in ['adjoint-error','forward-error']) + amrdata.flag_richardson_tol = 0.02*0.015 # suggested if using adjoint-error flag + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = (flag_method in ['adjoint-mag','forward-diff']) + amrdata.flag2refine_tol = 0.001 # suggested if using adjoint-mag flag + #amrdata.flag2refine_tol = 0.05 # suggested if using forward-diff + + # location of adjoint solution, must first be created: + adjointdata.adjoint_outdir = os.path.abspath('adjoint/_output') + + # time period of interest: + adjointdata.t1 = 20.5 + adjointdata.t2 = 21.5 + + if adjointdata.use_adjoint: + # need an additional aux variable for inner product: + rundata.amrdata.aux_type.append('center') + rundata.clawdata.num_aux = len(rundata.amrdata.aux_type) + adjointdata.innerprod_index = len(rundata.amrdata.aux_type) + + + # ----- 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() +