Merge pull request #212 from clawpack/batchGaugeOutput

Batch gauge output

examples/multi-layer/plane_wave/Makefile

@@ -81,7 +81,6 @@ MODULES = \
  ./qinit_module.f90
 
 SOURCES = \
-  $(GEOLIB)/multilayer/setprob.f90 \
   ./qinit.f90 \
   $(GEOLIB)/multilayer/setaux.f90 \
   $(GEOLIB)/multilayer/b4step2.f90 \

examples/multi-layer/plane_wave/qinit_module.f90

@@ -2,6 +2,8 @@ module qinit_module
 
     implicit none
     save
+
+    logical, private :: module_setup = .false.
 
     ! Type of q initialization
     integer, public :: qinit_type
@@ -33,6 +35,70 @@ module qinit_module
 
 contains
 
+    subroutine set_qinit(fname)
+
+        use geoclaw_module, only: GEO_PARM_UNIT
+
+        implicit none
+
+        ! Subroutine arguments
+        character(len=*), optional, intent(in) :: fname
+
+        ! File handling
+        integer, parameter :: unit = 7
+        character(len=150) :: qinit_fname
+
+        if (.not.module_setup) then
+
+            write(GEO_PARM_UNIT,*) ' '
+            write(GEO_PARM_UNIT,*) '--------------------------------------------'
+            write(GEO_PARM_UNIT,*) 'SETQINIT:'
+            write(GEO_PARM_UNIT,*) '-------------'
+
+            ! Open the data file
+            if (present(fname)) then
+                call opendatafile(unit,fname)
+            else
+                call opendatafile(unit,"qinit.data")
+            endif
+
+            read(unit,"(i1)") qinit_type
+            if (qinit_type == 0) then
+                ! No perturbation specified
+                write(GEO_PARM_UNIT,*)  '  qinit_type = 0, no perturbation'
+                print *,'  qinit_type = 0, no perturbation'
+                return
+            else if (qinit_type > 0 .and. qinit_type < 5) then
+                read(unit,*) qinit_fname
+                read(unit,"(2i2)") min_level_qinit, max_level_qinit
+
+                write(GEO_PARM_UNIT,*) '   min_level, max_level, qinit_fname:'
+                write(GEO_PARM_UNIT,*)  min_level_qinit, max_level_qinit, qinit_fname
+
+                call read_qinit(qinit_fname)
+            else if (qinit_type >= 5) then
+                read(unit,*) epsilon
+                read(unit,*) init_location
+                read(unit,*) wave_family
+                read(unit,*) angle
+                read(unit,*) sigma
+
+                write(GEO_PARM_UNIT,*) " epsilon = ",  epsilon
+                write(GEO_PARM_UNIT,*) " init_location = ",  init_location
+                write(GEO_PARM_UNIT,*) " wave_family = ",  wave_family
+                write(GEO_PARM_UNIT,*) " angle = ",  angle
+                write(GEO_PARM_UNIT,*) " sigma = ",  sigma
+            endif
+
+            close(unit)
+
+            module_setup = .true.
+        end if
+
+    end subroutine set_qinit
+
+
+
     subroutine add_perturbation(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux)
 
         use geoclaw_module, only: sea_level, pi, g => grav
@@ -164,64 +230,6 @@ subroutine add_perturbation(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux)
 
     end subroutine add_perturbation
 
-
-    subroutine set_qinit(fname)
-
-        use geoclaw_module, only: GEO_PARM_UNIT
-
-        implicit none
-
-        ! Subroutine arguments
-        character(len=*), optional, intent(in) :: fname
-
-        ! File handling
-        integer, parameter :: unit = 7
-        character(len=150) :: qinit_fname
-
-        write(GEO_PARM_UNIT,*) ' '
-        write(GEO_PARM_UNIT,*) '--------------------------------------------'
-        write(GEO_PARM_UNIT,*) 'SETQINIT:'
-        write(GEO_PARM_UNIT,*) '-------------'
-
-        ! Open the data file
-        if (present(fname)) then
-            call opendatafile(unit,fname)
-        else
-            call opendatafile(unit,"qinit.data")
-        endif
-
-        read(unit,"(i1)") qinit_type
-        if (qinit_type == 0) then
-            ! No perturbation specified
-            write(GEO_PARM_UNIT,*)  '  qinit_type = 0, no perturbation'
-            print *,'  qinit_type = 0, no perturbation'
-            return
-        else if (qinit_type > 0 .and. qinit_type < 5) then
-            read(unit,*) qinit_fname
-            read(unit,"(2i2)") min_level_qinit, max_level_qinit
-
-            write(GEO_PARM_UNIT,*) '   min_level, max_level, qinit_fname:'
-            write(GEO_PARM_UNIT,*)  min_level_qinit, max_level_qinit, qinit_fname
-
-            call read_qinit(qinit_fname)
-        else if (qinit_type >= 5) then
-            read(unit,*) epsilon
-            read(unit,*) init_location
-            read(unit,*) wave_family
-            read(unit,*) angle
-            read(unit,*) sigma
-
-            write(GEO_PARM_UNIT,*) " epsilon = ",  epsilon
-            write(GEO_PARM_UNIT,*) " init_location = ",  init_location
-            write(GEO_PARM_UNIT,*) " wave_family = ",  wave_family
-            write(GEO_PARM_UNIT,*) " angle = ",  angle
-            write(GEO_PARM_UNIT,*) " sigma = ",  sigma
-        endif
-
-        close(unit)
-
-    end subroutine set_qinit
-
 
     ! currently only supports one file type:
     ! x,y,z values, one per line in standard order from NW corner to SE

examples/multi-layer/plane_wave/setplot.py

@@ -27,21 +27,18 @@ def setplot(plotdata,  bathy_location=0.15,  bathy_angle=0.0,
 
     import clawpack.clawutil.data as clawutil
     import clawpack.amrclaw.data as amrclaw
-    import clawpack.geoclaw.data as geodata
+    import clawpack.geoclaw.data
 
-    import clawpack.geoclaw.multilayer.data as ml_data
     import clawpack.geoclaw.multilayer.plot as ml_plot
 
     # Load data from output
     clawdata = clawutil.ClawInputData(2)
     clawdata.read(os.path.join(plotdata.outdir,'claw.data'))
     amrdata = amrclaw.AmrclawInputData(clawdata)
     amrdata.read(os.path.join(plotdata.outdir,'amr.data'))
-    geodata = geodata.GeoClawData()
+    geodata = clawpack.geoclaw.data.GeoClawData()
     geodata.read(os.path.join(plotdata.outdir,'geoclaw.data'))
-    # surge_data = surge.data.SurgeData()
-    # surge_data.read(os.path.join(plotdata.outdir,'surge.data'))
-    multilayer_data = ml_data.MultilayerData()
+    multilayer_data = clawpack.geoclaw.data.MultilayerData()
     multilayer_data.read(os.path.join(plotdata.outdir,'multilayer.data'))
 
     def transform_c2p(x,y,x0,y0,theta):
@@ -567,7 +564,7 @@ def bottom_speed(current_data):
 
     # Plot surface as blue curve:
     plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
-    plotitem.plot_var = 6
+    plotitem.plot_var = 3
     plotitem.plotstyle = 'b-'
 
     # Bottom
@@ -579,7 +576,7 @@ def bottom_speed(current_data):
     # Set up for axes in this figure:
     plotaxes = plotfigure.new_plotaxes()
     plotaxes.xlimits = [0.0,1.0]
-    plotaxes.ylimits = internal_surface_limits
+    # plotaxes.ylimits = internal_surface_limits
     plotaxes.title = 'Bottom Surface'
 
     # Plot surface as blue curve:

examples/multi-layer/plane_wave/setrun.py

@@ -8,9 +8,10 @@
 
 import numpy as numpy
 
-import clawpack.geoclaw.multilayer.data as multilayer
+import clawpack.geoclaw.data
 import clawpack.geoclaw.topotools as tt
 
+
 # Rotation transformations
 def transform_c2p(x,y,x0,y0,theta):
     return ((x+x0)*numpy.cos(theta) - (y+y0)*numpy.sin(theta),
@@ -21,6 +22,54 @@ def transform_p2c(x,y,x0,y0,theta):
             -x*numpy.sin(theta) + y*numpy.cos(theta) - y0)
 
 
+# Class containing some setup for the qinit especially for multilayer tests 
+class QinitMultilayerData(clawpack.geoclaw.data.QinitData):
+    r"""
+    Modified Qinit data object for multiple layers
+
+    """
+
+    def __init__(self):
+        super(QinitMultilayerData, self).__init__()
+
+        # Test qinit data > 4
+        self.add_attribute("init_location", [0.0, 0.0])
+        self.add_attribute("wave_family", 1)
+        self.add_attribute("epsilon", 0.02)
+        self.add_attribute("angle", 0.0)
+        self.add_attribute("sigma", 0.02)
+
+    def write(self, data_source='setrun.py'):
+
+        # Initial perturbation
+        self.open_data_file('qinit.data',data_source)
+        self.data_write('qinit_type')
+
+        # Perturbation requested
+        if self.qinit_type == 0:
+            pass
+        elif 0 < self.qinit_type < 5:
+            # Check to see if each qinit file is present and then write the data
+            for tfile in self.qinitfiles:
+                try:
+                    fname = "'%s'" % os.path.abspath(tfile[-1])
+                except:
+                    raise Warning("File %s was not found." % fname)
+                    # raise MissingFile("file not found")
+                self._out_file.write("\n%s  \n" % fname)
+                self._out_file.write("%3i %3i \n" % tuple(tfile[:-1]))
+        elif self.qinit_type >= 5 and self.qinit_type <= 9:
+            self.data_write('epsilon')
+            self.data_write("init_location")
+            self.data_write("wave_family")
+            self.data_write("angle")
+            self.data_write("sigma")
+        else:
+            raise ValueError("Invalid qinit_type parameter %s." % self.qinit_type)
+        self.close_data_file()
+
+
+
 #------------------------------
 def setrun(claw_pkg='geoclaw'):
 #------------------------------
@@ -47,9 +96,7 @@ def setrun(claw_pkg='geoclaw'):
     # GeoClaw specific parameters:
     #------------------------------------------------------------------
     rundata = setgeo(rundata)
-
-    rundata.add_data(multilayer.MultilayerData(), 'multilayer_data')
-    set_multilayer(rundata)
+    rundata = set_multilayer(rundata)
 
     #------------------------------------------------------------------
     # Standard Clawpack parameters to be written to claw.data:
@@ -432,15 +479,16 @@ def set_multilayer(rundata):
     # data.wave_tolerance = [0.1,0.1]
     # data.dry_limit = True
 
-    # Set special initial conditions for qinit
-    rundata.replace_data('qinit_data', multilayer.QinitMultilayerData())
+    rundata.replace_data('qinit_data', QinitMultilayerData())
     rundata.qinit_data.qinit_type = 6
     rundata.qinit_data.epsilon = 0.02
     rundata.qinit_data.angle = 0.0
     rundata.qinit_data.sigma = 0.02
     rundata.qinit_data.wave_family = 4
     rundata.qinit_data.init_location = [-0.1,0.0]
 
+    return rundata
+
 
 def bathy_step(x, y, location=0.15, angle=0.0, left=-1.0, right=-0.2):
     x_c,y_c = transform_p2c(x, y, location, 0.0, angle)

src/2d/shallow/Makefile.geoclaw

@@ -22,6 +22,7 @@ COMMON_MODULES += \
  $(GEOLIB)/surge/stommel_storm_module.f90 \
  $(GEOLIB)/surge/constant_storm_module.f90 \
  $(GEOLIB)/surge/storm_module.f90 \
+ $(GEOLIB)/multilayer/multilayer_module.f90 \
  $(GEOLIB)/friction_module.f90
 
 # list of source files from AMR library. 

src/2d/shallow/advanc.f

@@ -166,7 +166,7 @@ subroutine prepgrids(listgrids,num, level)
       subroutine par_advanc (mptr,mitot,mjtot,nvar,naux,dtnew)
 c
       use amr_module
-      use gauges_module, only: print_gauges, num_gauges
+      use gauges_module, only: update_gauges, num_gauges
       implicit double precision (a-h,o-z)
 
 
@@ -240,9 +240,11 @@ subroutine par_advanc (mptr,mitot,mjtot,nvar,naux,dtnew)
 c        # now has boundary conditions filled in.
 
 c     should change the way print_gauges does io - right now is critical section
-
+c     NOW changed, mjb 2/6/2015.
+c     NOTE that gauge subr called before stepgrid, so never get
+c     the very last gauge time at end of run.
       if (num_gauges > 0) then
-           call print_gauges(alloc(locnew:locnew+nvar*mitot*mjtot), 
+           call update_gauges(alloc(locnew:locnew+nvar*mitot*mjtot), 
      .                       alloc(locaux:locnew+nvar*mitot*mjtot),
      .                       xlow,ylow,nvar,mitot,mjtot,naux,mptr)
            endif