Merge pull request #17 from rjleveque/q_out_vars

WIP: Allow specifying which components of q to output on each fgout grid

examples/radial_slide/pyvista_plotting.py

@@ -0,0 +1,205 @@
+from pylab import *
+import glob, os
+import pyvista as pv
+from clawpack.geoclaw import topotools
+from clawpack.visclaw import animation_tools
+from clawpack.geoclaw import fgout_tools
+from clawpack.visclaw import colormaps
+
+warpfactor = 5  # vertical amplification of elevations in 3D plots
+
+outdir  = '_output'
+fgno = 1 # fgout grid number
+
+if 1:
+    # determine how many fgout frames are in outdir:
+    fgout_frames = glob.glob(os.path.join(outdir, \
+                              'fgout%s.t*' % str(fgno).zfill(4)))
+    nfgout = len(fgout_frames)
+    fgframenos = array(range(1, nfgout+1))
+    print('Found %i fgout frames' % nfgout)
+else:
+    # specify which frames to use:
+    fgframenos = array(range(1,50,2))
+
+
+# where to save a png file for each frame, for constructing an animation:
+framedir = '_frames'
+os.system('mkdir -p %s' % framedir)
+os.system('rm %s/*' % framedir)  # remove frames from previous version
+
+window_size = (1200,1200)
+
+tan = [0.8,0.5,0.2]
+light_blue = [0.0,1.0,1.0]
+cmap_massfrac = colormaps.make_colormap({0:light_blue, 0.1:light_blue,
+                                         1:tan})
+
+def plot_topo():
+    """
+    plot topo from topofile
+    """
+    topo = topotools.Topography('basal_topo.tt3')
+    z = array([0.])
+    x = topo.x
+    y = -topo.y
+    X,Y,Z = meshgrid(x, y, z, indexing='ij')
+    topoxyz = pv.StructuredGrid(X,Y,Z)
+
+    B = flipud(topo.Z)
+    B = fliplr(B)
+    topoxyz.point_data['B'] = B.flatten(order='C')
+
+    topowarp = topoxyz.warp_by_scalar('B', factor=warpfactor)
+
+    p = pv.Plotter()
+    p.add_mesh(topowarp,cmap='gist_earth',clim=(-20,400))
+    p.add_title('Topography')
+    p.show(window_size=window_size)
+
+def plot_topo_fgout():
+    """
+    plot topo from an fgout file
+    """
+    fgno = 1  # which fgout grid
+
+    fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format, qmap='dclaw')
+
+    fgout = fgout_grid.read_frame(1)
+
+    x = fgout.x
+    y = fgout.y
+    z = array([0.])
+    X,Y,Z = meshgrid(x, y, z, indexing='ij')
+    topoxyz = pv.StructuredGrid(X,Y,Z)
+
+    B = fgout.B
+    topoxyz.point_data['B'] = B.flatten(order='C')
+    warpfactor = 5  # amplification of elevations
+    topowarp = topoxyz.warp_by_scalar('B', factor=warpfactor)
+
+    h = fgout.h
+    topoxyz.point_data['h'] = h.flatten(order='C')
+    warpfactor = 5  # amplification of elevations
+    hwarp = topoxyz.warp_by_scalar('h', factor=warpfactor)
+
+    eta = fgout.eta
+    topoxyz.point_data['eta'] = B.flatten(order='C')
+    warpfactor = 5  # amplification of elevations
+    etawarp = topoxyz.warp_by_scalar('eta', factor=warpfactor)
+
+    p = pv.Plotter()
+    p.window_size = window_size
+    p.add_mesh(topowarp,cmap='gist_earth',clim=(-20,400))
+    #p.add_mesh(topowarp,color='g')
+    p.add_title('Topography fgout.B')
+    p.show()
+
+
+def plot_fgout(make_animation=False):
+    """
+    If make_animation == False, plot one fgout frame with a slider bar to
+    change frame number. Also prints the camera position after changing frame,
+    which you can copy and paste into this function in order to set the
+    camera position for making an animation.
+    
+    If make_animation == True, make an mp4 animation of all frames
+    specified by fgframenos set above.
+    """
+
+    global etamesh
+
+    fgno = 1  # which fgout grid
+    warpfactor = 5  # amplification of elevations
+
+    # Instantiate object for reading fgout frames:
+    fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format, qmap='dclaw')
+
+    fgout = fgout_grid.read_frame(1)
+
+    x = fgout.x
+    y = fgout.y
+    z = array([0.])
+    X,Y,Z = meshgrid(x, y, z, indexing='ij')
+    topoxyz = pv.StructuredGrid(X,Y,Z)
+
+    B = fgout.B
+    topoxyz.point_data['B'] = B.flatten(order='C')
+    topowarp = topoxyz.warp_by_scalar('B', factor=warpfactor)
+
+    p = pv.Plotter(off_screen=make_animation, lighting='three lights')
+    p.window_size = window_size
+    #p.add_mesh(topowarp,cmap='gist_earth',clim=(-20,400))
+    p.add_mesh(topowarp,color='lightgreen')
+    etamesh = None
+
+
+    def set_frameno(fgframeno):
+        global etamesh
+        fgframeno = int(round(fgframeno))
+        fgout = fgout_grid.read_frame(fgframeno)
+        tsec = fgout.t
+        print('Frame %i, t = %.1f seconds' % (fgframeno, fgout.t))
+
+        # replace land surface by nan in eta so it only shows water:
+        # (big tolerance for landslide, would normally be smaller)
+        eta = where(fgout.h>1, fgout.eta, nan)
+
+        topoxyz.point_data['eta'] = eta.flatten(order='C')
+        etawarp = topoxyz.warp_by_scalar('eta', factor=warpfactor)
+        if etamesh:
+            p.remove_actor(etamesh)
+
+        # color the mesh based on mass fraction, using cmap_massfrac
+        h = fgout.h
+        massfrac = divide(fgout.hm, h, where=h>0, out=nan*ones(h.shape))
+        etamesh = p.add_mesh(etawarp,scalars=massfrac,
+                             colormap=cmap_massfrac, clim=(0,0.6))
+
+        p.add_title('Frame %i at time %.1f seconds' % (fgframeno,tsec))
+
+        if not make_animation:
+            # print camera position so that this can be copied and pasted
+            # into this script after adjusting (and then sliding frameno)
+            print('p.camera_position = ', p.camera_position)
+
+
+    # initial camera position:
+    p.camera_position =  [
+         (9034.796206317897, -2484.6171987620633, 3703.6128928929047),
+         (1500.0, 1500.0, 950.89111328125),
+         (-0.29984440650740857, 0.08916816605502331, 0.949811755059182)]
+
+
+    if not make_animation:
+        fgfr1 = fgframenos[0]
+        fgfr2 = fgframenos[-1]
+        p.add_slider_widget(set_frameno, [fgfr1,fgfr2],
+                        value=fgfr1, title='Frame',
+                        pointa=(0.4,0.85), pointb=(0.9,0.85), color='blue',
+                        slider_width=0.02, tube_width=0.005)
+
+        p.show()
+
+    else:
+
+        # make a png file for each frame:
+        for fgframeno in fgframenos:
+            set_frameno(fgframeno)
+            fname_png = '%s/PyVistaFrame%s.png' \
+                        % (framedir, str(fgframeno).zfill(4))
+            p.screenshot(fname_png)
+            print('Created ',fname_png)
+
+        p.close()
+
+        # combine png files into mp4 animation:
+        anim = animation_tools.make_anim(framedir,
+               fname_pattern='PyVistaFrame*.png')
+        fname_mp4 = 'radial_slide_animation.mp4'
+        animation_tools.make_mp4(anim, fname_mp4)
+        print('Created ',fname_mp4)
+
+if __name__=='__main__':
+
+    plot_fgout(make_animation=False)

examples/radial_slide/setrun.py

@@ -16,6 +16,8 @@
     raise Exception("*** Must first set CLAW environment variable")
 
 from clawpack.amrclaw.data import FlagRegion
+from clawpack.geoclaw import fgout_tools
+
 
 
 #------------------------------
@@ -124,8 +126,8 @@ def setrun(claw_pkg='dclaw'):
 
     if clawdata.output_style==1:
         # Output nout frames at equally spaced times up to tfinal:
-        clawdata.num_output_times = 30 #240
-        clawdata.tfinal = 120. #240.
+        clawdata.num_output_times = 10 #240
+        clawdata.tfinal = 100. #240.
         clawdata.output_t0 = True  # output at initial (or restart) time?
 
     elif clawdata.output_style == 2:
@@ -280,7 +282,7 @@ def setrun(claw_pkg='dclaw'):
     amrdata = rundata.amrdata
 
     # max number of refinement levels:
-    amrdata.amr_levels_max = 1
+    amrdata.amr_levels_max = 2
 
     # List of refinement ratios at each level (length at least mxnest-1)
     # dx = dy = 2', 10", 2", 1/3":
@@ -394,8 +396,8 @@ def setrun(claw_pkg='dclaw'):
     etafile = 'surface_topo.tt3'
     qinitdclaw_data.qinitfiles.append([3, 8, 1, 2, etafile])
 
-    #mfile = 'mass_frac.tt3'
-    mfile = 'mass_frac.tt3' # with m0 = 0 below
+    mfile = 'mass_frac.tt3' 
+    #mfile = 'mass_frac0.tt3' # with m0 = 0 below
     qinitdclaw_data.qinitfiles.append([3, 4, 1, 2, mfile])
 
     #hfile = 'landslide_depth.tt3'
@@ -460,7 +462,32 @@ def setrun(claw_pkg='dclaw'):
     #flowgrades_data.flowgrades.append([1.0e-6, 2, 1, 1])
     #flowgrades_data.flowgrades.append([1.0e-6, 1, 1, 1])
 
-    #  ----- For developers -----
+
+    # == fgout_grids.data values ==
+    # NEW IN v5.9.0
+    # Set rundata.fgout_data.fgout_grids to be a list of
+    # objects of class clawpack.geoclaw.fgout_tools.FGoutGrid:
+    fgout_grids = rundata.fgout_data.fgout_grids  # empty list initially
+
+    fgout = fgout_tools.FGoutGrid()
+    fgout.fgno = 1
+    fgout.point_style = 2       # will specify a 2d grid of points
+    #fgout.output_format = 'binary32'  # 4-byte, float32
+    fgout.output_format = 'ascii'  # 4-byte, float32
+    fgout.nx = 300
+    fgout.ny = 300
+    fgout.x1 = 0.  # specify edges (fgout pts will be cell centers)
+    fgout.x2 = 3e3
+    fgout.y1 = 0.
+    fgout.y2 = 3e3
+    fgout.tstart = 0.
+    fgout.tend = 100.
+    fgout.nout = 101
+    fgout.q_out_vars = [1,4,8]
+    fgout_grids.append(fgout)    # written to fgout_grids.data
+
+
+    #  ----- For developers ----- 
     # Toggle debugging print statements:
     amrdata.dprint = False      # print domain flags
     amrdata.eprint = False      # print err est flags

src/2d/dig/Makefile.dclaw

@@ -23,14 +23,14 @@ COMMON_MODULES += \
  $(GEOLIB)/multilayer/multilayer_module.f90 \
  $(GEOLIB)/friction_module.f90 \
  $(DIGLIB)/refinement_module.f90 \
- $(GEOLIB)/adjointsup_module.f90
+ $(GEOLIB)/adjointsup_module.f90 \
+ $(GEOLIB)/fgout_module.f90 \
 
 #list of common modules needed for dclaw codes
 COMMON_MODULES += \
  $(DIGLIB)/digclaw_module.f90 \
  $(DIGLIB)/qinit_module.f90 \
  $(DIGLIB)/auxinit_module.f90 \
- $(DIGLIB)/fgout_module.f90 \
 
 
 # list of source files from AMR library.

src/2d/dig/amr2.f90

@@ -491,7 +491,7 @@ program amr2
         call read_dtopo_settings()        ! specifies file with dtopo from earthquake
         call read_topo_settings(rest)     ! specifies topography (bathymetry) files
         call set_auxinit()            ! specifies file(s) for auxiliary variables
-        call set_fgout(rest)            ! Fixed grid settings
+        call set_fgout(rest,nvar)            ! Fixed grid settings
         call set_regions()                ! Set refinement regions
         call set_dig(naux)
         call set_pinit()
@@ -526,7 +526,7 @@ program amr2
         call read_dtopo_settings()        ! specifies file with dtopo from earthquake
         call read_topo_settings(rest)     ! specifies topography (bathymetry) files
         call set_auxinit()            ! specifies file(s) for auxiliary variables
-        call set_fgout(rest)            ! Fixed grid settings
+        call set_fgout(rest,nvar)            ! Fixed grid settings
         call set_regions()                ! Set refinement regions
         call set_gauges(rest, nvar, naux) ! Set gauge output
         call set_fgmax()