From 0b595de01b554af20f52c18a0f3873435c33f42e Mon Sep 17 00:00:00 2001 From: Wes Lowrie ARA/SED Date: Wed, 22 Nov 2017 16:30:03 -0800 Subject: [PATCH] Various changes to pyclaw. Updates to properly handle restarts using hdf5 --- examples/euler_1d/shocksine.py | 7 ++++++- src/pyclaw/classic/setup.py | 2 +- src/pyclaw/classic/solver.py | 25 +++++++++++++++++++++++-- src/pyclaw/io/hdf5.py | 4 ++-- src/pyclaw/sharpclaw/solver.py | 21 +++++++++++++++++---- src/pyclaw/solver.py | 21 ++++++++++++++++++++- src/pyclaw/state.py | 2 +- 7 files changed, 70 insertions(+), 12 deletions(-) diff --git a/examples/euler_1d/shocksine.py b/examples/euler_1d/shocksine.py index 57d9229c8..b770f3ec9 100755 --- a/examples/euler_1d/shocksine.py +++ b/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 diff --git a/src/pyclaw/classic/setup.py b/src/pyclaw/classic/setup.py index a50de681b..d5ff10de7 100644 --- a/src/pyclaw/classic/setup.py +++ b/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 diff --git a/src/pyclaw/classic/solver.py b/src/pyclaw/classic/solver.py index 2b91c234e..57fdf6474 100644 --- a/src/pyclaw/classic/solver.py +++ b/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) diff --git a/src/pyclaw/io/hdf5.py b/src/pyclaw/io/hdf5.py index 390e5e4ee..8dbde118a 100644 --- a/src/pyclaw/io/hdf5.py +++ b/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) diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index 3a42282b7..4158d5861 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/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) diff --git a/src/pyclaw/solver.py b/src/pyclaw/solver.py index 5bc708c6a..e73209f3b 100644 --- a/src/pyclaw/solver.py +++ b/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 diff --git a/src/pyclaw/state.py b/src/pyclaw/state.py index c47dbc4d1..9851fd657 100755 --- a/src/pyclaw/state.py +++ b/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