From 4da5b7ad6b102e407b3beb35fccd2b08f6b88196 Mon Sep 17 00:00:00 2001 From: David Ketcheson Date: Sat, 15 Nov 2014 16:19:05 +0300 Subject: [PATCH 1/2] Add function space class. The functionSpace class sits between State and Patch. A function space lives on a patch but also has a number of DOFs. It is intended to match the PETSc DMDA abstraction. This will allow us to reduce the number of DA's and Vecs we use, by reusing DA's from the main State object in the _register objects in SharpClaw. It is also anticipated that this abstraction will be useful when implementing DG finite element methods. --- src/petclaw/io/petsc.py | 8 +- src/petclaw/state.py | 178 +++++++++++++++++++-------------- src/pyclaw/sharpclaw/solver.py | 7 +- src/pyclaw/state.py | 40 +++++++- 4 files changed, 149 insertions(+), 84 deletions(-) diff --git a/src/petclaw/io/petsc.py b/src/petclaw/io/petsc.py index 217d0dff3..d2fb6c6e4 100644 --- a/src/petclaw/io/petsc.py +++ b/src/petclaw/io/petsc.py @@ -90,10 +90,10 @@ def write(solution,frame,path='./',file_prefix='claw',write_aux=False, if write_p: state.gpVec.view(q_viewer) else: - state.gqVec.view(q_viewer) + state._q_global_vector.view(q_viewer) if write_aux: - state.gauxVec.view(aux_viewer) + state._aux_global_vector.view(aux_viewer) q_viewer.flush() if write_aux: @@ -188,10 +188,10 @@ def read(solution,frame,path='./',file_prefix='claw',read_aux=False,options={}): # DA View/Load is broken in Petsc-3.1.8, we can load/view the DA if needed in petsc-3.2 # state.q_da.load(q_viewer) - state.gqVec.load(q_viewer) + state._q_global_vector.load(q_viewer) if read_aux: - state.gauxVec.load(aux_viewer) + state._aux_global_vector.load(aux_viewer) solution.states.append(state) patches.append(state.patch) diff --git a/src/petclaw/state.py b/src/petclaw/state.py index 1aa0f671c..76f5d1723 100644 --- a/src/petclaw/state.py +++ b/src/petclaw/state.py @@ -1,5 +1,53 @@ import clawpack.pyclaw +class functionSpace(clawpack.pyclaw.state.functionSpace): + + def __init__(self,patch,num_dof): + self.patch = patch + self.num_dof = num_dof + self.da = self._create_DA() + + def set_num_ghost(self,num_ghost): + self.da = self._create_DA(num_ghost) + + def _create_DA(self,num_ghost=0): + r"""Returns a PETSc DA and associated global Vec. + Note that no local vector is returned. + """ + from petsc4py import PETSc + + #Due to the way PETSc works, we just make the patch always periodic, + #regardless of the boundary conditions actually selected. + #This works because in solver.qbc() we first call globalToLocal() + #and then impose the real boundary conditions (if non-periodic). + + if hasattr(PETSc.DA, 'PeriodicType'): + if self.num_dim == 1: + periodic_type = PETSc.DA.PeriodicType.X + elif self.num_dim == 2: + periodic_type = PETSc.DA.PeriodicType.XY + elif self.num_dim == 3: + periodic_type = PETSc.DA.PeriodicType.XYZ + else: + raise Exception("Invalid number of dimensions") + + DA = PETSc.DA().create(dim=self.patch.num_dim, + dof=self.num_dof, + sizes=self.patch.num_cells_global, + periodic_type = periodic_type, + stencil_width=num_ghost, + comm=PETSc.COMM_WORLD) + else: + DA = PETSc.DA().create(dim=self.patch.num_dim, + dof=self.num_dof, + sizes=self.patch.num_cells_global, + boundary_type = PETSc.DA.BoundaryType.PERIODIC, + stencil_width=num_ghost, + comm=PETSc.COMM_WORLD) + + return DA + + class State(clawpack.pyclaw.State): """Parallel State class""" @@ -8,15 +56,15 @@ class State(clawpack.pyclaw.State): @property def num_eqn(self): r"""(int) - Number of unknowns (components of q)""" - if self.q_da is None: + if self.q_space is None: raise Exception('state.num_eqn has not been set.') - else: return self.q_da.dof + else: return self.q_space.num_dof @property def num_aux(self): r"""(int) - Number of auxiliary fields""" - if self.aux_da is None: return 0 - else: return self.aux_da.dof + if self.aux_space is None: return 0 + else: return self.aux_space.num_dof @property def mp(self): @@ -53,10 +101,10 @@ def q(self): """ shape = self.grid.num_cells shape.insert(0,self.num_eqn) - return self.gqVec.getArray().reshape(shape, order = 'F') + return self._q_global_vector.getArray().reshape(shape, order = 'F') @q.setter def q(self,val): - self.gqVec.setArray(val.reshape([-1], order = 'F')) + self._q_global_vector.setArray(val.reshape([-1], order = 'F')) @property def p(self): @@ -98,19 +146,19 @@ def aux(self): values for the aux array. The global aux vector is used only for outputting the aux values to file; everywhere else we use the local vector. """ - if self.aux_da is None: return None + if self.aux_space is None: return None shape = self.grid.num_cells shape.insert(0,self.num_aux) - aux=self.gauxVec.getArray().reshape(shape, order = 'F') + aux=self._aux_global_vector.getArray().reshape(shape, order = 'F') return aux @aux.setter def aux(self,val): # It would be nice to make this work also for parallel # loading from a file. - if self.aux_da is None: + if self.aux_space is None: num_aux=val.shape[0] - self._init_aux_da(num_aux) - self.gauxVec.setArray(val.reshape([-1], order = 'F')) + self.aux_space = functionSpace(num_aux) + self._aux_global_vector.setArray(val.reshape([-1], order = 'F')) @property def num_dim(self): return self.patch.num_dim @@ -121,6 +169,11 @@ def __init__(self,geom,num_eqn,num_aux=0): Here we don't call super because q and aux must be properties in PetClaw but should not be properties in PyClaw. + The arguments num_eqn and num_aux can be integers, in which case + a new DA with that many DOFs is created. Alternatively, they + can be functionSpaces, in which case the existing function space is + just attached to the state. + :attributes: patch - The patch this state lives on """ @@ -134,8 +187,8 @@ def __init__(self,geom,num_eqn,num_aux=0): raise Exception("""A PetClaw State object must be initialized with a PetClaw Patch or Domain object.""") - self.aux_da = None - self.q_da = None + self.aux_space = None + self.q_space = None self._p_da = None self.gpVec = None @@ -159,66 +212,37 @@ def __init__(self,geom,num_eqn,num_aux=0): stores the values of the corresponding gauge if ``keep_gauges`` is set to ``True``""" - self._init_q_da(num_eqn) - if num_aux>0: self._init_aux_da(num_aux) + if type(num_eqn) is int: + self.q_space = functionSpace(self.patch,num_eqn) + else: + self.function_space = num_eqn + + if (type(num_aux) is int) and num_aux>0: + self.aux_space = functionSpace(self.patch,num_aux) + self._aux_global_vector = self.aux_space.da.createGlobalVector() + elif num_aux != 0: + self.aux_space = num_aux + else: # num_aux==0 + self.aux_space = None + + self._init_global_vecs() - def _init_aux_da(self,num_aux,num_ghost=0): + def _init_global_vecs(self): r""" - Initializes PETSc DA and global & local Vectors for handling the - auxiliary array, aux. - - Initializes aux_da, gauxVec and _aux_local_vector. + Initializes PETSc global Vectors. """ - self.aux_da = self._create_DA(num_aux,num_ghost) - self.gauxVec = self.aux_da.createGlobalVector() - self._aux_local_vector = self.aux_da.createLocalVector() - - def _init_q_da(self,num_eqn,num_ghost=0): + self._q_global_vector = self.q_space.da.createGlobalVector() + if self.num_aux > 0: + self._aux_global_vector = self.aux_space.da.createGlobalVector() + + def _init_local_vecs(self): r""" - Initializes PETSc DA and Vecs for handling the solution, q. - - Initializes q_da, gqVec and _q_local_vector. - """ - self.q_da = self._create_DA(num_eqn,num_ghost) - self.gqVec = self.q_da.createGlobalVector() - self._q_local_vector = self.q_da.createLocalVector() - - def _create_DA(self,dof,num_ghost=0): - r"""Returns a PETSc DA and associated global Vec. - Note that no local vector is returned. + Initializes PETSc local Vectors. """ - from petsc4py import PETSc - - #Due to the way PETSc works, we just make the patch always periodic, - #regardless of the boundary conditions actually selected. - #This works because in solver.qbc() we first call globalToLocal() - #and then impose the real boundary conditions (if non-periodic). - - if hasattr(PETSc.DA, 'PeriodicType'): - if self.num_dim == 1: - periodic_type = PETSc.DA.PeriodicType.X - elif self.num_dim == 2: - periodic_type = PETSc.DA.PeriodicType.XY - elif self.num_dim == 3: - periodic_type = PETSc.DA.PeriodicType.XYZ - else: - raise Exception("Invalid number of dimensions") - - DA = PETSc.DA().create(dim=self.num_dim, - dof=dof, - sizes=self.patch.num_cells_global, - periodic_type = periodic_type, - stencil_width=num_ghost, - comm=PETSc.COMM_WORLD) - else: - DA = PETSc.DA().create(dim=self.num_dim, - dof=dof, - sizes=self.patch.num_cells_global, - boundary_type = PETSc.DA.BoundaryType.PERIODIC, - stencil_width=num_ghost, - comm=PETSc.COMM_WORLD) - - return DA + self._q_local_vector = self.q_space.da.createLocalVector() + if self.num_aux > 0: + self._aux_local_vector = self.aux_space.da.createLocalVector() + def get_qbc_from_q(self,num_ghost,qbc): """ @@ -226,7 +250,7 @@ def get_qbc_from_q(self,num_ghost,qbc): """ shape = [n + 2*num_ghost for n in self.grid.num_cells] - self.q_da.globalToLocal(self.gqVec, self._q_local_vector) + self.q_space.da.globalToLocal(self._q_global_vector, self._q_local_vector) shape.insert(0,self.num_eqn) return self._q_local_vector.getArray().reshape(shape, order = 'F') @@ -236,7 +260,7 @@ def get_auxbc_from_aux(self,num_ghost,auxbc): """ shape = [n + 2*num_ghost for n in self.grid.num_cells] - self.aux_da.globalToLocal(self.gauxVec, self._aux_local_vector) + self.aux_space.da.globalToLocal(self._aux_global_vector, self._aux_local_vector) shape.insert(0,self.num_aux) return self._aux_local_vector.getArray().reshape(shape, order = 'F') @@ -253,12 +277,16 @@ def set_num_ghost(self,num_ghost): but it only happens once so it seems not to be worth it. """ q0 = self.q.copy() - self._init_q_da(self.num_eqn,num_ghost) + self.q_space.set_num_ghost(num_ghost) + if self.num_aux > 0: + self.aux_space.set_num_ghost(num_ghost) + self._init_global_vecs() + self._init_local_vecs() self.q = q0 if self.aux is not None: aux0 = self.aux.copy() - self._init_aux_da(self.num_aux,num_ghost) + self.aux_space.set_num_ghost(num_ghost) self.aux = aux0 def sum_F(self,i): @@ -269,8 +297,8 @@ def get_q_global(self): Returns a copy of the global q array on process 0, otherwise returns None """ from petsc4py import PETSc - q_natural = self.q_da.createNaturalVec() - self.q_da.globalToNatural(self.gqVec, q_natural) + q_natural = self.q_space.da.createNaturalVec() + self.q_space.da.globalToNatural(self._q_global_vector, q_natural) scatter, q0Vec = PETSc.Scatter.toZero(q_natural) scatter.scatter(q_natural, q0Vec, False, PETSc.Scatter.Mode.FORWARD) rank = PETSc.COMM_WORLD.getRank() @@ -292,7 +320,7 @@ def get_aux_global(self): """ from petsc4py import PETSc aux_natural = self.aux_da.createNaturalVec() - self.aux_da.globalToNatural(self.gauxVec, aux_natural) + self.aux_da.globalToNatural(self._aux_global_vector, aux_natural) scatter, aux0Vec = PETSc.Scatter.toZero(aux_natural) scatter.scatter(aux_natural, aux0Vec, False, PETSc.Scatter.Mode.FORWARD) rank = PETSc.COMM_WORLD.getRank() diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index 2f821dc9e..92edd2f6a 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/src/pyclaw/sharpclaw/solver.py @@ -559,11 +559,14 @@ def _allocate_registers(self,solution): self._registers = [] for i in xrange(nregisters): #Maybe should use State.copy() here? - self._registers.append(State(state.patch,state.num_eqn,state.num_aux)) + if state.num_aux > 0: + self._registers.append(State(state.patch,state.q_space,state.aux_space)) + self._registers[-1].aux = state.aux + else: + self._registers.append(State(state.patch,state.q_space)) self._registers[-1].problem_data = state.problem_data self._registers[-1].set_num_ghost(self.num_ghost) self._registers[-1].t = state.t - if state.num_aux > 0: self._registers[-1].aux = state.aux def get_cfl_max(self): diff --git a/src/pyclaw/state.py b/src/pyclaw/state.py index c47dbc4d1..9c1f1e93a 100755 --- a/src/pyclaw/state.py +++ b/src/pyclaw/state.py @@ -7,6 +7,25 @@ David I. Ketcheson -- Initial version (June 2011) """ +class functionSpace(object): + r""" + A PyClaw functionSpace object lives on a patch but also has a concept of + the number of degrees of freedom. Typically, a State object has a functionSpace + for q and another function space for aux, if aux variables are present. + + For now, this class exists mainly to match the DMDA abstraction in PETSc. + It will also be useful in implementing DG finite element methods . + """ + def __init__(self,patch,num_dof): + self.patch = patch + self.num_dof = num_dof + + def set_num_ghost(self,num_ghost): + """ + Virtual routine (does nothing). Overridden in the petclaw.state class. + """ + pass + class State(object): r""" @@ -14,7 +33,8 @@ class State(object): including the unkowns q, the time t, and the auxiliary coefficients aux. The variables num_eqn and num_aux determine the length of the first - dimension of the q and aux arrays. + dimension of the q and aux arrays. The arguments num_eqn and num_aux + to the initialization routine may alternatively be functionSpace objects. :State Data: @@ -151,9 +171,23 @@ def __init__(self,geom,num_eqn,num_aux=0): stores the values of the corresponding gauge if ``keep_gauges`` is set to ``True``""" - - self.q = self.new_array(num_eqn) self.aux = self.new_array(num_aux) + if type(num_eqn) is int: + self.q_space = functionSpace(self.patch,num_eqn) + self.q = self.new_array(num_eqn) + else: + self.q_space = num_eqn + self.q = self.new_array(self.q_space.num_dof) + + if type(num_aux) is int and num_aux>0: + self.aux_space = functionSpace(self.patch,num_eqn) + self.aux = self.new_array(num_aux) + elif num_aux != 0: + self.aux_space = num_aux + self.aux = self.new_array(self.aux_space.num_dof) + else: + self.aux_space = None + def __str__(self): output = "PyClaw State object\n" From 93a8fdfa15a4afd517bc18db81eebb737c7e2e60 Mon Sep 17 00:00:00 2001 From: David Ketcheson Date: Tue, 18 Nov 2014 14:46:40 +0300 Subject: [PATCH 2/2] Move everything to function spaces. --- src/petclaw/io/petsc.py | 2 +- src/petclaw/state.py | 81 ++++++++++++++++------------------ src/pyclaw/classic/solver.py | 2 +- src/pyclaw/sharpclaw/solver.py | 1 - src/pyclaw/state.py | 46 +++++++++---------- 5 files changed, 62 insertions(+), 70 deletions(-) diff --git a/src/petclaw/io/petsc.py b/src/petclaw/io/petsc.py index d2fb6c6e4..5427835ed 100644 --- a/src/petclaw/io/petsc.py +++ b/src/petclaw/io/petsc.py @@ -88,7 +88,7 @@ def write(solution,frame,path='./',file_prefix='claw',write_aux=False, # we will reenable this bad boy when we switch over to petsc-dev # state.q_da.view(q_viewer) if write_p: - state.gpVec.view(q_viewer) + state._p_global_vector.view(q_viewer) else: state._q_global_vector.view(q_viewer) diff --git a/src/petclaw/state.py b/src/petclaw/state.py index 76f5d1723..3cd242cbe 100644 --- a/src/petclaw/state.py +++ b/src/petclaw/state.py @@ -1,14 +1,11 @@ import clawpack.pyclaw -class functionSpace(clawpack.pyclaw.state.functionSpace): +class FunctionSpace(clawpack.pyclaw.state.FunctionSpace): def __init__(self,patch,num_dof): - self.patch = patch - self.num_dof = num_dof + super(FunctionSpace,self).__init__(patch,num_dof) self.da = self._create_DA() - def set_num_ghost(self,num_ghost): - self.da = self._create_DA(num_ghost) def _create_DA(self,num_ghost=0): r"""Returns a PETSc DA and associated global Vec. @@ -69,30 +66,30 @@ def num_aux(self): @property def mp(self): r"""(int) - Number of derived quantities (components of p)""" - if self._p_da is None: + if not hasattr(self,'p_space'): raise Exception('state.mp has not been set.') - else: return self._p_da.dof + else: return self.p_space.num_dof @mp.setter def mp(self,mp): - if self._p_da is not None: + if hasattr(self,'p_space'): raise Exception('You cannot change state.mp after p is initialized.') else: - self._p_da = self._create_DA(mp) - self.gpVec = self._p_da.createGlobalVector() + self.p_space = FunctionSpace(self.patch,mp) + self._p_global_vector = self.p_space.da.createGlobalVector() @property def mF(self): r"""(int) - Number of derived quantities (components of p)""" - if self._F_da is None: + if not hasattr(self,'F_space'): raise Exception('state.mF has not been set.') - else: return self._F_da.dof + else: return self.F_space.num_dof @mF.setter def mF(self,mF): - if self._F_da is not None: - raise Exception('You cannot change state.mp after p is initialized.') + if hasattr(self,'F_space'): + raise Exception('You cannot change state.mF after F is initialized.') else: - self._F_da = self._create_DA(mF) - self.gFVec = self._F_da.createGlobalVector() + self.F_space = FunctionSpace(self.patch,mF) + self._F_global_vector = self.F_space.da.createGlobalVector() @property def q(self): @@ -111,16 +108,16 @@ def p(self): r""" Array containing values of derived quantities for output. """ - if self._p_da is None: return 0 + if not hasattr(self,'p_space'): return 0 shape = self.grid.num_cells shape.insert(0,self.mp) - p=self.gpVec.getArray().reshape(shape, order = 'F') + p=self._p_global_vector.getArray().reshape(shape, order = 'F') return p @p.setter def p(self,val): mp = val.shape[0] - if self.gpVec is None: self.init_p_da(mp) - self.gpVec.setArray(val.reshape([-1], order = 'F')) + if not hasattr(self,'_p_global_vector'): self.mp = mp + self._p_global_vector.setArray(val.reshape([-1], order = 'F')) @property def F(self): @@ -128,16 +125,16 @@ def F(self): Array containing pointwise values (densities) of output functionals. This is just used as temporary workspace before summing. """ - if self._F_da is None: return 0 + if not hasattr(self,'F_space'): return 0 shape = self.grid.num_cells shape.insert(0,self.mF) - F=self.gFVec.getArray().reshape(shape, order = 'F') + F=self._F_global_vector.getArray().reshape(shape, order = 'F') return F @F.setter def fset(self,val): mF = val.shape[0] - if self.gFVec is None: self.init_F_da(mF) - self.gFVec.setArray(val.reshape([-1], order = 'F')) + if not hasattr(self,'_F_global_vector'): self.mF = mF + self._F_global_vector.setArray(val.reshape([-1], order = 'F')) @property def aux(self): @@ -157,7 +154,7 @@ def aux(self,val): # loading from a file. if self.aux_space is None: num_aux=val.shape[0] - self.aux_space = functionSpace(num_aux) + self.aux_space = FunctionSpace(num_aux) self._aux_global_vector.setArray(val.reshape([-1], order = 'F')) @property def num_dim(self): @@ -190,12 +187,6 @@ def __init__(self,geom,num_eqn,num_aux=0): self.aux_space = None self.q_space = None - self._p_da = None - self.gpVec = None - - self._F_da = None - self.gFVec = None - # ========== Attribute Definitions =================================== self.problem_data = {} r"""(dict) - Dictionary of global values for this patch, @@ -213,12 +204,12 @@ def __init__(self,geom,num_eqn,num_aux=0): to ``True``""" if type(num_eqn) is int: - self.q_space = functionSpace(self.patch,num_eqn) + self.q_space = FunctionSpace(self.patch,num_eqn) else: - self.function_space = num_eqn + self.q_space = num_eqn if (type(num_aux) is int) and num_aux>0: - self.aux_space = functionSpace(self.patch,num_aux) + self.aux_space = FunctionSpace(self.patch,num_aux) self._aux_global_vector = self.aux_space.da.createGlobalVector() elif num_aux != 0: self.aux_space = num_aux @@ -226,6 +217,7 @@ def __init__(self,geom,num_eqn,num_aux=0): self.aux_space = None self._init_global_vecs() + self._init_local_vecs() def _init_global_vecs(self): r""" @@ -277,20 +269,23 @@ def set_num_ghost(self,num_ghost): but it only happens once so it seems not to be worth it. """ q0 = self.q.copy() - self.q_space.set_num_ghost(num_ghost) + self.q_space.da = self.q_space._create_DA(num_ghost) + if self.num_aux > 0: - self.aux_space.set_num_ghost(num_ghost) + aux0 = self.aux.copy() + self.aux_space.da = self.aux_space._create_DA(num_ghost) + + # Need new vecs because we have new DAs self._init_global_vecs() self._init_local_vecs() - self.q = q0 - if self.aux is not None: - aux0 = self.aux.copy() - self.aux_space.set_num_ghost(num_ghost) + # Copy old state into new vecs + self.q = q0 + if self.num_aux > 0: self.aux = aux0 def sum_F(self,i): - return self.gFVec.strideNorm(i,0) + return self._F_global_vector.strideNorm(i,0) def get_q_global(self): r""" @@ -319,8 +314,8 @@ def get_aux_global(self): Returns a copy of the global aux array on process 0, otherwise returns None """ from petsc4py import PETSc - aux_natural = self.aux_da.createNaturalVec() - self.aux_da.globalToNatural(self._aux_global_vector, aux_natural) + aux_natural = self.aux_space.da.createNaturalVec() + self.aux_space.da.globalToNatural(self._aux_global_vector, aux_natural) scatter, aux0Vec = PETSc.Scatter.toZero(aux_natural) scatter.scatter(aux_natural, aux0Vec, False, PETSc.Scatter.Mode.FORWARD) rank = PETSc.COMM_WORLD.getRank() diff --git a/src/pyclaw/classic/solver.py b/src/pyclaw/classic/solver.py index 95b6325cb..46e286985 100644 --- a/src/pyclaw/classic/solver.py +++ b/src/pyclaw/classic/solver.py @@ -636,7 +636,7 @@ def _allocate_workspace(self,solution): #with f2py. It involves wastefully allocating three arrays. #f2py seems not able to handle multiple zero-size arrays being passed. # it appears the bug is related to f2py/src/fortranobject.c line 841. - if(aux == None): num_aux=1 + if(aux is None): num_aux=1 grid = state.grid maxmx,maxmy,maxmz = grid.num_cells[0],grid.num_cells[1],grid.num_cells[2] diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index 92edd2f6a..2ecbb4275 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/src/pyclaw/sharpclaw/solver.py @@ -565,7 +565,6 @@ def _allocate_registers(self,solution): else: self._registers.append(State(state.patch,state.q_space)) self._registers[-1].problem_data = state.problem_data - self._registers[-1].set_num_ghost(self.num_ghost) self._registers[-1].t = state.t diff --git a/src/pyclaw/state.py b/src/pyclaw/state.py index 9c1f1e93a..8f64e18b8 100755 --- a/src/pyclaw/state.py +++ b/src/pyclaw/state.py @@ -7,10 +7,10 @@ David I. Ketcheson -- Initial version (June 2011) """ -class functionSpace(object): +class FunctionSpace(object): r""" - A PyClaw functionSpace object lives on a patch but also has a concept of - the number of degrees of freedom. Typically, a State object has a functionSpace + A PyClaw FunctionSpace object lives on a patch but also has a concept of + the number of degrees of freedom. Typically, a State object has a FunctionSpace for q and another function space for aux, if aux variables are present. For now, this class exists mainly to match the DMDA abstraction in PETSc. @@ -20,6 +20,8 @@ def __init__(self,patch,num_dof): self.patch = patch self.num_dof = num_dof + super(FunctionSpace,self).__init__() + def set_num_ghost(self,num_ghost): """ Virtual routine (does nothing). Overridden in the petclaw.state class. @@ -34,7 +36,7 @@ class State(object): The variables num_eqn and num_aux determine the length of the first dimension of the q and aux arrays. The arguments num_eqn and num_aux - to the initialization routine may alternatively be functionSpace objects. + to the initialization routine may alternatively be FunctionSpace objects. :State Data: @@ -118,25 +120,27 @@ def grid(self): @property def mp(self): r"""(int) - Number of derived quantities""" - if self.p is not None: return self.p.shape[0] + if hasattr(self,'p'): return self.p.shape[0] else: return 0 @mp.setter def mp(self,mp): - if self.p is not None: - raise Exception('Cannot change state.mp after aux is initialized.') + if hasattr(self,'p'): + raise Exception('Cannot change state.mp after it is initialized.') else: + self.p_space = FunctionSpace(self.patch,mp) self.p = self.new_array(mp) @property def mF(self): r"""(int) - Number of output functionals""" - if self.F is not None: return self.F.shape[0] + if hasattr(self,'F'): return self.F.shape[0] else: return 0 @mF.setter def mF(self,mF): - if self.F is not None: + if hasattr(self,'F'): raise Exception('Cannot change state.mF after aux is initialized.') else: + self.F_space = FunctionSpace(self.patch,mF) self.F = self.new_array(mF) # ========== Class Methods =============================================== @@ -150,12 +154,6 @@ def __init__(self,geom,num_eqn,num_aux=0): raise Exception("""A PyClaw State object must be initialized with a PyClaw Patch object.""") - # ========== Attribute Definitions =================================== - r"""pyclaw.Patch.patch - The patch this state lives on""" - self.p = None - r"""(ndarray(mp,...)) - Cell averages of derived quantities.""" - self.F = None - r"""(ndarray(mF,...)) - Cell averages of output functional densities.""" self.problem_data = {} r"""(dict) - Dictionary of global values for this patch, ``default = {}``""" @@ -171,22 +169,22 @@ def __init__(self,geom,num_eqn,num_aux=0): stores the values of the corresponding gauge if ``keep_gauges`` is set to ``True``""" - self.aux = self.new_array(num_aux) - if type(num_eqn) is int: - self.q_space = functionSpace(self.patch,num_eqn) - self.q = self.new_array(num_eqn) - else: + if type(num_eqn) is FunctionSpace: self.q_space = num_eqn self.q = self.new_array(self.q_space.num_dof) + else: + self.q_space = FunctionSpace(self.patch,num_eqn) + self.q = self.new_array(num_eqn) - if type(num_aux) is int and num_aux>0: - self.aux_space = functionSpace(self.patch,num_eqn) - self.aux = self.new_array(num_aux) - elif num_aux != 0: + if type(num_aux) is FunctionSpace: self.aux_space = num_aux self.aux = self.new_array(self.aux_space.num_dof) + elif num_aux>0: + self.aux_space = FunctionSpace(self.patch,num_eqn) + self.aux = self.new_array(num_aux) else: self.aux_space = None + self.aux = None def __str__(self):