Various changes to pyclaw. Updates to properly handle restarts using hdf5

Author: weslowrie

examples/euler_1d/shocksine.py

@@ -52,10 +52,15 @@ def setup(use_petsc=False,iplot=False,htmlplot=False,outdir='./_output',solver_t
 
     if solver_type=='sharpclaw':
         solver = pyclaw.SharpClawSolver1D(rs)
-        solver.time_integrator = 'RK'
+        solver.time_integrator = 'Euler'
         solver.a, solver.b, solver.c = a, b, c
         solver.cfl_desired = 0.6
         solver.cfl_max = 0.7
+
+        solver.weno_order = 2
+        solver.num_ghost = 2
+        solver.lim_type = 1
+        
         if use_char_decomp:
             try:
                 import sharpclaw1               # Import custom Fortran code

src/pyclaw/classic/setup.py

@@ -13,7 +13,7 @@ def configuration(parent_package='',top_path=None):
                          ['limiter.f90','philim.f90','flux2.f90','step2ds.f90','step2.f90'],f2py_options=['--quiet'])
 
     config.add_extension('classic3',
-                         ['limiter.f90','philim.f90','flux3.f90','step3ds.f90','step3.f90'],f2py_options=['--quiet'],extra_f90_compile_args=['-fopenmp'],libraries=['gomp'])
+                         ['limiter.f90','philim.f90','flux3.f90','step3ds.f90','step3.f90'],f2py_options=['--quiet'],extra_f90_compile_args=['-fopenmp -Warray-bounds -fcheck=all'],libraries=['gomp'])
 
     return config
 

src/pyclaw/classic/solver.py

@@ -141,7 +141,7 @@ def step(self,solution,take_one_step,tstart,tend):
             # Godunov Splitting
             if self.source_split == 1:
                 self.step_source(self,solution.states[0],self.dt)
-                
+
         return True
 
     def _check_cfl_settings(self):
@@ -684,7 +684,7 @@ def step_hyperbolic(self,solution):
             if self.dimensional_split:
                 #Right now only Godunov-dimensional-splitting is implemented.
                 #Strang-dimensional-splitting could be added following dimsp3.f in Clawpack.
-
+                #"""
                 q, cfl_x = self.fmod.step3ds(maxm,self.nthreads,self.num_ghost,\
                     mx,my,mz,qold,qnew,self.auxbc,dx,dy,dz,self.dt,self._method,\
                     self._mthlim,self.aux1,self.aux2,self.aux3,self.work,1,\
@@ -699,6 +699,27 @@ def step_hyperbolic(self,solution):
                     mx,my,mz,q,q,self.auxbc,dx,dy,dz,self.dt,self._method,\
                     self._mthlim,self.aux1,self.aux2,self.aux3,self.work,3,\
                     self.fwave,rpn3,rpt3,rptt3)
+                #"""
+
+                """
+                #print "No Dim Splitting..."
+                # No Dimensional Splitting, Similar to SharpClaw
+                dq = qnew.copy('F')
+                dq[...] = 0.
+                dq, cfl_x = self.fmod.step3ds(maxm,self.nthreads,self.num_ghost,\
+                    mx,my,mz,qold,dq,self.auxbc,dx,dy,dz,self.dt,self._method,\
+                    self._mthlim,self.aux1,self.aux2,self.aux3,self.work,1,\
+                    self.fwave,rpn3,rpt3,rptt3)
+                dq, cfl_y = self.fmod.step3ds(maxm,self.nthreads,self.num_ghost,\
+                    mx,my,mz,qold,dq,self.auxbc,dx,dy,dz,self.dt,self._method,\
+                    self._mthlim,self.aux1,self.aux2,self.aux3,self.work,2,\
+                    self.fwave,rpn3,rpt3,rptt3)
+                dq, cfl_z = self.fmod.step3ds(maxm,self.nthreads,self.num_ghost,\
+                    mx,my,mz,qold,dq,self.auxbc,dx,dy,dz,self.dt,self._method,\
+                    self._mthlim,self.aux1,self.aux2,self.aux3,self.work,3,\
+                    self.fwave,rpn3,rpt3,rptt3)
+                self.qbc += dq
+                """
 
                 cfl = max(cfl_x,cfl_y,cfl_z)
 

src/pyclaw/io/hdf5.py

@@ -193,11 +193,11 @@ def read(solution,frame,path='./',file_prefix='claw',read_aux=True,
                 state = pyclaw.state.State(pyclaw_patch, \
                          patch.attrs['num_eqn'],patch.attrs['num_aux'])
                 state.t = patch.attrs['t']
-                state.q = patch['q'][:].reshape(state.q.shape,order='F')
+                state.q = patch['q'][:].ravel(order='F').reshape(state.q.shape,order='F')
 
                 # Read in aux if applicable
                 if read_aux and patch.get('aux',None) is not None:
-                    state.aux = patch['aux'][:].reshape(state.aux.shape,order='F')
+                    state.aux = patch['aux'][:].ravel(order='F').reshape(state.aux.shape,order='F')
 
                 solution.states.append(state)
                 patches.append(pyclaw_patch)

src/pyclaw/sharpclaw/solver.py

@@ -284,6 +284,7 @@ def step(self,solution,take_one_step,tstart,tend):
         Take one step with a Runge-Kutta or multistep method as specified by
         `solver.time_integrator`.
         """
+        
         state = solution.states[0]
         step_index = self.status['numsteps'] + 1
         if self.accept_step == True:
@@ -515,9 +516,19 @@ def check_3rd_ord_cond(self,state,step_index,dtFE):
     def _set_mthlim(self):
         self._mthlim = self.limiters
         if not isinstance(self.limiters,list): self._mthlim=[self._mthlim]
-        if len(self._mthlim)==1: self._mthlim = self._mthlim * self.num_waves
-        if len(self._mthlim)!=self.num_waves:
-            raise Exception('Length of solver.limiters is not equal to 1 or to solver.num_waves')
+        if self.lim_type==1:
+            if self.char_decomp == 0 or self.char_decomp == 3:
+                if len(self._mthlim)==1: self._mthlim = self._mthlim * self.num_eqn
+                if len(self._mthlim)!=self.num_eqn:
+                    raise Exception('Length of solver.limiters is not equal to 1 or to solver.num_eqn')
+            elif self.char_decomp == 1 or self.char_decomp == 4:
+                if len(self._mthlim)==1: self._mthlim = self._mthlim * self.num_waves
+                if len(self._mthlim)!=self.num_waves:
+                    raise Exception('Length of solver.limiters is not equal to 1 or to solver.num_waves')
+        else:
+            if len(self._mthlim)==1: self._mthlim = self._mthlim * self.num_waves
+            if len(self._mthlim)!=self.num_waves:
+                raise Exception('Length of solver.limiters is not equal to 1 or to solver.num_waves')
 
 
     def dq(self,state):
@@ -526,7 +537,9 @@ def dq(self,state):
         """
 
         deltaq = self.dq_hyperbolic(state)
-
+        
+        #state.q += deltaq
+        #deltaq[...] = 0.
         if self.dq_src is not None:
             deltaq+=self.dq_src(self,state,self.dt)
 

src/pyclaw/solver.py

@@ -32,6 +32,14 @@ def before_step(solver,solution):
     """
     pass
 
+def after_step(solver,solution):
+    r"""
+    Dummy routine called after each step
+    
+    Replace this routine if you want to do something after each time step.
+    """
+    pass
+
 class Solver(object):
     r"""
     Pyclaw solver superclass.
@@ -85,6 +93,13 @@ class Solver(object):
         
         def before_step(solver,solution)
 
+    .. attribute:: after_step
+    
+        Function called after each time step is taken.
+        The required signature for this function is:
+        
+        def after_step(solver,solution)
+
     .. attribute:: dt_variable
     
         Whether to allow the time step to vary, ``default = True``.
@@ -179,6 +194,7 @@ def __init__(self,riemann_solver=None,claw_package=None):
         self._use_old_bc_sig = False
         self.accept_step = True
         self.before_step = None
+        self.after_step = None
 
         # select package to build solver objects from, by default this will be
         # the package that contains the module implementing the derived class
@@ -610,6 +626,9 @@ def evolve_to_time(self,solution,tend=None):
             # Note that the solver may alter dt during the step() routine
             self.step(solution,take_one_step,tstart,tend)
 
+            if self.after_step is not None:
+                self.after_step(self,solution.states[0])
+
             # Check to make sure that the Courant number was not too large
             cfl = self.cfl.get_cached_max()
             self.accept_step = self.accept_reject_step(state)
@@ -642,7 +661,7 @@ def evolve_to_time(self,solution,tend=None):
                     self.status['cflmax'] = \
                         max(cfl, self.status['cflmax'])
                     raise Exception('CFL too large, giving up!')
- 
+                
             # See if we are finished yet
             if solution.t >= tend or take_one_step:
                 break

src/pyclaw/state.py

@@ -189,7 +189,7 @@ def is_valid(self):
             valid = False
         if self.aux is not None:
             if not self.aux.flags['F_CONTIGUOUS']:
-                logger.debug('q array is not Fortran contiguous.')
+                logger.debug('aux array is not Fortran contiguous.')
                 valid = False
         return valid