diff --git a/src/petclaw/geometry.py b/src/petclaw/geometry.py index e3f29c8c8..0812ca364 100755 --- a/src/petclaw/geometry.py +++ b/src/petclaw/geometry.py @@ -13,11 +13,11 @@ class Patch(pyclaw_geometry.Patch): __doc__ += pyclaw.util.add_parent_doc(pyclaw_geometry.Patch) - def __init__(self,dimensions): + def __init__(self,dimensions,proc_sizes=None): super(Patch,self).__init__(dimensions) - self._da = self._create_DA() + self._da = self._create_DA(proc_sizes) ranges = self._da.getRanges() grid_dimensions = [] for i,nrange in enumerate(ranges): @@ -41,7 +41,7 @@ def __init__(self,dimensions): self.grid = pyclaw_geometry.Grid(grid_dimensions) - def _create_DA(self): + def _create_DA(self,proc_sizes=None): r"""Returns a PETSc DA and associated global Vec. Note that no local vector is returned. """ @@ -62,14 +62,16 @@ def _create_DA(self): sizes=self.num_cells_global, periodic_type = periodic_type, stencil_width=0, - comm=PETSc.COMM_WORLD) + comm=PETSc.COMM_WORLD, + proc_sizes=proc_sizes) else: DA = PETSc.DA().create(dim=self.num_dim, dof=1, sizes=self.num_cells_global, boundary_type = PETSc.DA.BoundaryType.PERIODIC, stencil_width=0, - comm=PETSc.COMM_WORLD) + comm=PETSc.COMM_WORLD, + proc_sizes=proc_sizes) return DA @@ -82,13 +84,13 @@ class Domain(pyclaw_geometry.Domain): __doc__ += pyclaw.util.add_parent_doc(pyclaw.ClawSolver2D) - def __init__(self,geom): + def __init__(self,geom,proc_sizes=None): if not isinstance(geom,list): geom = [geom] if isinstance(geom[0],Patch): self.patches = geom elif isinstance(geom[0],pyclaw_geometry.Dimension): - self.patches = [Patch(geom)] + self.patches = [Patch(geom, proc_sizes)] class Dimension(pyclaw_geometry.Dimension): def __init__(self, lower, upper, num_cells, name='x', diff --git a/src/petclaw/solution.py b/src/petclaw/solution.py index 88a8aad4d..0e3cdd1ba 100644 --- a/src/petclaw/solution.py +++ b/src/petclaw/solution.py @@ -1,5 +1,6 @@ from clawpack import pyclaw from clawpack.pyclaw.solution import Solution +import numpy as np class Solution(Solution): """ Parallel Solution class. @@ -24,4 +25,89 @@ def get_write_func(self, file_format): else: raise ValueError("File format %s not supported." % file_format) + def _init_ds_solution(self): + """ + Initializes a downsampled version of the solution + """ + import clawpack.petclaw as pyclaw + + ds_domain = pyclaw.Domain([pyclaw.Dimension(self.domain.patch.lower_global[i], + self.domain.patch.upper_global[i], + self.domain.patch.num_cells_global[i]/self.downsampling_factors[i], + self.domain.grid.dimensions[i].name) + for i in range(self.domain.num_dim)], proc_sizes=self.state.q_da.getProcSizes()) + ds_state = self._init_ds_state(self.state) + self._ds_solution = pyclaw.Solution(ds_state, ds_domain) + + def _init_ds_state(self, state): + """ + Returns a downsampled version of the given state object + """ + import clawpack.petclaw as pyclaw + + ds_domain = pyclaw.Domain([pyclaw.Dimension(self.domain.patch.lower_global[i], + self.domain.patch.upper_global[i], + self.domain.patch.num_cells_global[i]/self.downsampling_factors[i], + self.domain.grid.dimensions[i].name) + for i in range(self.domain.num_dim)], proc_sizes=state.q_da.getProcSizes()) + ds_state = pyclaw.State(ds_domain,state.num_eqn,state.num_aux) + + return ds_state + + def downsample(self, write_aux,write_p): + """ + Returns a downsampled version of the solution by local averaging over the downsampling factors + """ + for i in range(len(self.states)): + state = self.states[i] + if i > 0: + self.ds_solution.states.append(self._init_ds_state(self.downsample_factors, state)) + ds_state = self.ds_solution.states[i] + + # Downsample q + if write_p: + ds_state.p = self._downsample_global_to_local(state.get_q_da_for_ds(np.max(self.downsampling_factors)), + state.gqVec, state.num_eqn, ds_state.patch._da.getRanges()) + else: + ds_state.q = self._downsample_global_to_local(state.get_q_da_for_ds(np.max(self.downsampling_factors)), + state.gqVec, state.num_eqn, ds_state.patch._da.getRanges()) + # Dowsample aux + if write_aux: + ds_state.aux = self._downsample_global_to_local(state.get_aux_da_for_ds(np.max(self.downsampling_factors)), + state.gauxVec, state.num_aux, ds_state.patch._da.getRanges()) + + return self.ds_solution + + def _downsample_global_to_local(self, da_for_ds, gVec, num_eqn, ds_domain_ranges): + """ + Returns a locally averaged array and handles ranges mapping between the original & downsampled solution objects + + Input: + - da_for_ds: the DA object of the original solution (q or aux) with stencil width adjusted as appropriate to + handle the boundaries of the array to be averaged because after domain decomposition, there are no guarantees that + the downsampling factors will evenly divide the the sub-domain in each direction + - gVec: The global vector associated with q or aux + - num_eqn: number of components of q or aux + - ds_domain_ranges: global ranges of the downsampled solution + """ + + from skimage.transform import downscale_local_mean + + # Create local array with ghost cells + q_local_vec = da_for_ds.createLocalVec() + da_for_ds.globalToLocal(gVec, q_local_vec) + q_local_array = q_local_vec.array + shape = [end-start for start,end in da_for_ds.getGhostRanges()] + shape.insert(0, num_eqn) + q_local_array = q_local_array.reshape(shape, order='f') + + # Map global ranges in the original DA object to local ranges for local averaging + q_da_ghost_ranges = da_for_ds.getGhostRanges() + new_global_ranges = tuple([(s * self.downsampling_factors[i], e * self.downsampling_factors[i]) for i, (s,e) in enumerate(ds_domain_ranges)]) + new_local_slices = (slice(0, q_local_array.shape[0]),) + new_local_slices = new_local_slices + tuple([slice(s - gs,q_local_array.shape[i+1] - (ge - e)) for i, ((s,e), (gs, ge)) in enumerate(zip(new_global_ranges, q_da_ghost_ranges))]) + q_local_array = q_local_array[new_local_slices] + + # Compute local mean + return downscale_local_mean(q_local_array, (1,) + self.downsampling_factors) diff --git a/src/petclaw/state.py b/src/petclaw/state.py index 1809b97b5..3f1bb125c 100644 --- a/src/petclaw/state.py +++ b/src/petclaw/state.py @@ -107,7 +107,7 @@ def aux(self): 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_da is None: num_aux=val.shape[0] self._init_aux_da(num_aux) self.gauxVec.setArray(val.reshape([-1], order = 'F')) @@ -136,6 +136,8 @@ def __init__(self,geom,num_eqn,num_aux=0): self.aux_da = None self.q_da = None + self._q_da_for_ds = None + self._aux_da_for_ds = None self._p_da = None self.gpVec = None @@ -145,14 +147,14 @@ def __init__(self,geom,num_eqn,num_aux=0): # ========== Attribute Definitions =================================== self.problem_data = {} - r"""(dict) - Dictionary of global values for this patch, + r"""(dict) - Dictionary of global values for this patch, ``default = {}``""" self.t=0. - r"""(float) - Current time represented on this patch, + r"""(float) - Current time represented on this patch, ``default = 0.0``""" self.index_capa = -1 self.keep_gauges = False - r"""(bool) - Keep gauge values in memory for every time step, + r"""(bool) - Keep gauge values in memory for every time step, ``default = False``""" self.gauge_data = [] r"""(list) - List of numpy.ndarray objects. Each element of the list @@ -165,18 +167,18 @@ def __init__(self,geom,num_eqn,num_aux=0): def _init_aux_da(self,num_aux,num_ghost=0): r""" Initializes PETSc DA and global & local Vectors for handling the - auxiliary array, aux. - + auxiliary array, aux. + Initializes aux_da, gauxVec and _aux_local_vector. """ 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): r""" - Initializes PETSc DA and Vecs for handling the solution, q. - + 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) @@ -225,7 +227,7 @@ def get_qbc_from_q(self,num_ghost,qbc): Returns q with ghost cells attached, by accessing the local vector. """ shape = [n + 2*num_ghost for n in self.grid.num_cells] - + self.q_da.globalToLocal(self.gqVec, self._q_local_vector) shape.insert(0,self.num_eqn) return self._q_local_vector.getArray().reshape(shape, order = 'F') @@ -235,7 +237,7 @@ def get_auxbc_from_aux(self,num_ghost,auxbc): Returns aux with ghost cells attached, by accessing the local vector. """ shape = [n + 2*num_ghost for n in self.grid.num_cells] - + self.aux_da.globalToLocal(self.gauxVec, self._aux_local_vector) shape.insert(0,self.num_aux) return self._aux_local_vector.getArray().reshape(shape, order = 'F') @@ -316,3 +318,27 @@ def __deepcopy__(self,memo={}): result.set_num_ghost(self.q_da.stencil_width) return result + + def get_q_da_for_ds(self, ds_stencil_width): + """ + Returns a q_da object with an adjusted width if needed in order to compute local averages for downsampling + """ + if self._q_da_for_ds is None: + if self.q_da.getStencilWidth() < ds_stencil_width: + self._q_da_for_ds = self.q_da.duplicate(stencil_width=ds_stencil_width) + else: + self._q_da_for_ds = self.q_da + + return self._q_da_for_ds + + def get_aux_da_for_ds(self, ds_stencil_width): + """ + Returns an aux_da object with an adjusted width if needed in order to compute local averages for downsampling + """ + if self._aux_da_for_ds is None: + if self.aux_da.getStencilWidth() < ds_stencil_width: + self._aux_da_for_ds = self.aux_da.duplicate(stencil_width=ds_stencil_width) + else: + self._aux_da_for_ds = self.aux_da + + return self._aux_da_for_ds \ No newline at end of file diff --git a/src/pyclaw/controller.py b/src/pyclaw/controller.py index 1ffdd9c86..603bfe1e7 100644 --- a/src/pyclaw/controller.py +++ b/src/pyclaw/controller.py @@ -193,6 +193,7 @@ def __init__(self): r"""(function) - Function that computes density of functional F""" self.F_file_name = 'F' r"""(string) - Name of text file containing functionals""" + self.downsampling_factors = None # ========== Access methods =============================================== def __str__(self): @@ -326,11 +327,20 @@ def run(self): if self.t == self.tfinal: print "Simulation has already reached tfinal." return None - + + self.solution.downsampling_factors = self.downsampling_factors + # Output and save initial frame if self.keep_copy: self.frames.append(copy.deepcopy(self.solution)) if self.output_format is not None: + if not (self.downsampling_factors is None): + if self.solution.domain.num_dim != len(self.downsampling_factors): + raise ValueError("Invalid number of downsampling factors. The number of downsampling factors should match the number of dimensions.") + for i, factor in enumerate(self.downsampling_factors): + if self.solution.domain.patch.dimensions[i].num_cells % factor != 0: + raise ValueError("Invalid downsampling factors. The downsampling factors should evenly divide the grid in each dimension.") + if os.path.exists(self.outdir) and self.overwrite==False: raise Exception("Refusing to overwrite existing output data. \ \nEither delete/move the directory or set controller.overwrite=True.") @@ -341,7 +351,7 @@ def run(self): self.file_prefix_p, write_aux = False, options = self.output_options, - write_p = True) + write_p = True) write_aux = (self.write_aux_always or self.write_aux_init) self.solution.write(frame,self.outdir, @@ -374,7 +384,7 @@ def run(self): self.file_prefix_p, write_aux = False, options = self.output_options, - write_p = True) + write_p = True) self.solution.write(frame,self.outdir, self.output_format, diff --git a/src/pyclaw/solution.py b/src/pyclaw/solution.py index 3e404c09e..161fcac87 100755 --- a/src/pyclaw/solution.py +++ b/src/pyclaw/solution.py @@ -8,6 +8,8 @@ import logging from .geometry import Patch, Dimension, Domain +from .state import State +import numpy as np # ============================================================================ # Solution Class @@ -82,7 +84,7 @@ def __setattr__(self, key, value): self.__dict__[key] = value # ========== Attributes ================================================== - + # ========== Properties ================================================== @property def state(self): @@ -100,17 +102,25 @@ def start_frame(self): return self._start_frame _start_frame = 0 + @property + def ds_solution(self): + """ + Lazily instantiates and returns a downsampled version of the solution + """ + if self._ds_solution is None and not (self.downsampling_factors is None): + self._init_ds_solution() + return self._ds_solution # ========== Class Methods =============================================== def __init__(self,*arg,**kargs): r"""Solution Initiatlization Routine - + See :class:`Solution` for more info. """ # select package to build solver objects from, by default this will be # the package that contains the module implementing the derived class - # for example, if Solution is implemented in 'clawpack.petclaw.solution', then + # for example, if Solution is implemented in 'clawpack.petclaw.solution', then # the computed claw_package will be 'clawpack.petclaw' import sys @@ -179,16 +189,17 @@ def get_clawpack_dot_xxx(modname): return modname.rpartition('.')[0] self.states.append(State(self.domain,arg[0])) if self.states == [] or self.domain is None: raise Exception("Invalid argument list") - - + + self.downsampling_factors = None + self._ds_solution = None def is_valid(self): r""" Checks to see if this solution is valid - + The Solution checks to make sure it is valid by checking each of its states. If an invalid state is found, a message is logged what specifically made this solution invalid. - + :Output: - (bool) - True if valid, false otherwise """ @@ -201,51 +212,51 @@ def __str__(self): for state in self.states: output = output + str(state) return str(output) - - + + def set_all_states(self,attr,value,overwrite=True): r""" Sets all member states attribute 'attr' to value - + :Input: - *attr* - (string) Attribute name to be set - *value* - (id) Value for attribute - - *overwrite* - (bool) Whether to overwrite the attribute if it + - *overwrite* - (bool) Whether to overwrite the attribute if it already exists. ``default = True`` """ for state in self.states: if getattr(state,attr) is None or overwrite: - setattr(state,attr,value) - - + setattr(state,attr,value) + + def _get_base_state_attribute(self, name): r""" Return base state attribute - + :Output: - (id) - Value of attribute from ``states[0]`` """ return getattr(self.states[0],name) - - + + def __copy__(self): return self.__class__(self) - - + + def __deepcopy__(self,memo={}): import copy # Create basic container result = self.__class__() result.__init__() - + # Populate the states for state in self.states: result.states.append(copy.deepcopy(state)) result.domain = copy.deepcopy(self.domain) - + return result - - + + # ========== IO Functions ================================================ def write(self,frame,path='./',file_format='ascii',file_prefix=None, write_aux=False,options={},write_p=False): @@ -258,15 +269,15 @@ def write(self,frame,path='./',file_format='ascii',file_prefix=None, :Input: - *frame* - (int) Frame number to append to the file output - - *path* - (string) Root path, will try and create the path if it + - *path* - (string) Root path, will try and create the path if it does not already exist. ``default = './'`` - - *format* - (string or list of strings) a string or list of strings + - *format* - (string or list of strings) a string or list of strings containing the desired output formats. ``default = 'ascii'`` - *file_prefix* - (string) Prefix for the file name. Defaults to the particular io modules default. - - *write_aux* - (book) Write the auxillary array out as well if + - *write_aux* - (book) Write the auxillary array out as well if present. ``default = False`` - - *options* - (dict) Dictionary of optional arguments dependent on + - *options* - (dict) Dictionary of optional arguments dependent on which format is being used. ``default = {}`` """ # Determine if we need to create the path @@ -275,7 +286,7 @@ def write(self,frame,path='./',file_format='ascii',file_prefix=None, try: os.makedirs(path) except OSError: - print "directory already exists, ignoring" + print "directory already exists, ignoring" # Call the correct write function based on the output format if isinstance(file_format,str): @@ -283,54 +294,59 @@ def write(self,frame,path='./',file_format='ascii',file_prefix=None, elif isinstance(file_format,list): format_list = file_format + # Downsample the solution if needed + if self.downsampling_factors is None: + solution = self + else: + solution = self.downsample(write_aux,write_p) # Loop over list of formats requested for form in format_list: write_func = self.get_write_func(form) if file_prefix is None: - write_func(self,frame,path,write_aux=write_aux, + write_func(solution,frame,path,write_aux=write_aux, options=options,write_p=write_p) else: - write_func(self,frame,path,file_prefix=file_prefix, + write_func(solution,frame,path,file_prefix=file_prefix, write_aux=write_aux,options=options, write_p=write_p) msg = "Wrote out solution in format %s for time t=%s" % (form,self.t) logging.getLogger('pyclaw.io').info(msg) - - def read(self, frame, path='./_output', file_format='ascii', + + def read(self, frame, path='./_output', file_format='ascii', file_prefix=None, read_aux=True, options={}, **kargs): r""" Reads in a Solution object from a file - - Reads in and initializes this Solution with the data specified. This - function will raise an IOError if it was unsuccessful. + + Reads in and initializes this Solution with the data specified. This + function will raise an IOError if it was unsuccessful. Any format must conform to the following call signiture and return True if the file has been successfully read into the given solution or False otherwise. Options is a dictionary of parameters that each format can specify. See the ascii module for an example.:: - + read_(solution,path,frame,file_prefix,options={}) - + ```` is the name of the format in question. - + :Input: - *frame* - (int) Frame number to be read in - - *path* - (string) Base path to the files to be read. + - *path* - (string) Base path to the files to be read. ``default = './_output'`` - - *file_format* - (string) Format of the file, should match on of the + - *file_format* - (string) Format of the file, should match on of the modules inside of the io package. ``default = 'ascii'`` - - *file_prefix* - (string) Name prefix in front of all the files, + - *file_prefix* - (string) Name prefix in front of all the files, defaults to whatever the format defaults to, e.g. fort for ascii - - *options* - (dict) Dictionary of optional arguments dependent on + - *options* - (dict) Dictionary of optional arguments dependent on the format being read in. ``default = {}`` - + :Output: - (bool) - True if read was successful, False otherwise """ - + read_func = self.get_read_func(file_format) path = os.path.expandvars(os.path.expanduser(path)) @@ -356,7 +372,7 @@ def get_write_func(self, file_format): from clawpack.pyclaw import io return getattr(getattr(io,file_format),'write') - + def plot(self): r""" Plot the solution @@ -365,6 +381,47 @@ def plot(self): "implemented as of yet, please refer to the plotting module for" + " how to plot solutions.") + def _init_ds_solution(self): + """ + Initializes a downsampled version of the solution + """ + domain = Domain(self.domain.grid.lower,self.domain.grid.upper,self.domain.grid.num_cells/np.array(self.downsampling_factors)) + state = State(domain,self.state.num_eqn,self.state.num_aux) + self._ds_solution = Solution(state, domain) + + def _init_ds_state(self, state): + """ + Returns a downsampled version of the given state object + """ + ds_domain = Domain(self.domain.grid.lower,self.domain.grid.upper,self.domain.grid.num_cells/np.array(self.downsampling_factors)) + ds_state = State(ds_domain,self.state.num_eqn,self.state.num_aux) + + return ds_state + + def downsample(self, write_aux,write_p): + """ + Returns a downsampled version of the solution by local averaging over the downsampling factors + """ + from skimage.transform import downscale_local_mean + + for i in range(len(self.states)): + state = self.states[i] + if i > 0: + self.ds_solution.states.append(self._init_ds_state(self.downsample_factors, state)) + ds_state = self.ds_solution.states[i] + + # Downsample q + if write_p: + ds_state.p = downscale_local_mean(state.p, (1,) + self.downsampling_factors) + else: + ds_state.q = downscale_local_mean(state.q, (1,) + self.downsampling_factors) + + # Dowsample aux + if write_aux: + ds_state.aux = downscale_local_mean(state.aux, (1,) + self.downsampling_factors) + + return self.ds_solution + if __name__ == "__main__": import doctest doctest.testmod()