verbose mode for some functions

field.py

@@ -150,7 +150,7 @@ class Field3d:
 
     def insert_subfield(self,subfield):
         assert all([abs(subfield.spacing[ii]-self.spacing[ii])<self.eps_collapse 
-                        for ii in range(3)]), "spacing differs."
+                        for ii in range(3)]), "spacing differs. Got {}, have {}".format(subfield.spacing,self.spacing)
         assert all([self.distance_to_nearest_gridpoint(subfield.origin[ii],axis=ii)<self.eps_collapse 
                         for ii in range(3)]), "subfield has shifted origin."
         assert all(self.is_within_bounds(subfield.origin,axis=None)), "subfield origin is out-of-bounds."

parallel.py

@@ -128,13 +128,13 @@ class PPP:
             proc_grid_ext,nxp_ext,nyp_ext,nzp_ext,
             nghbr,field,symmetries)
 
-    def load_field(self,key,io_limit=None,dtype=np.float64,verbose=False):
+    def load_field(self,key,io_limit=None,dtype=np.float64,verbose=True):
         '''Loads the required chunks from file'''
         from .field import Field3d
         import numpy as np
         # Verbose output
         if verbose and self.rank==0:
-            print('[load_field] loading {} (io_limit={})'.format(key,io_limit))
+            print('[load_field] key={}, io_limit={}'.format(key,io_limit))
         # Block execution of some processors if IO is limited
         sb = SequentialBlock(io_limit,comm=self.comm,per_node=True)
         # Determine which chunks are to be loaded by the current processor
@@ -305,7 +305,7 @@ class PPP:
             self.field[key_out] *= self.field[key2]
         return
 
-    def gaussian_filter(self,key,sigma,truncate=4.0,key_out=None,iterate=False):
+    def gaussian_filter(self,key,sigma,truncate=4.0,key_out=None,iterate=False,verbose=True):
         '''Applies a gaussian filter to a field as in-place operation. Sigma is the std of the filter in terms of grid width.'''
         import numpy as np
         from mpi4py import MPI
@@ -315,6 +315,8 @@ class PPP:
             self.copy(key,key_out,skip_data=True)
         # Compute radius of Gaussian filter
         radius = self.field[key].gaussian_filter_radius(sigma,truncate=truncate)
+        if verbose and self.rank==0:
+            print('[gaussian_filter] key={}, stencil radius={}'.format(key,radius))
         if not iterate:
             # Assert that we have sufficient amount of ghost cells
             assert all([self.num_ghost[ii]>=radius[ii] for ii in range(3)]),\
@@ -331,9 +333,11 @@ class PPP:
             assert all([radius[ii]>0 if sigma[ii]>0.0 else True for ii in range(3)]),\
                 "Iterative procedure leads to invalid stencil radius: "\
                 "increase number of ghost cells. {}".format(radius)
-            parprint('Gaussian filter: iterations={}, stencil radius={}'.format(niter,radius))
+            if verbose and self.rank==0:
+                print('[gaussian_filter] iterations={}'.format(niter))
         for iiter in range(niter):
-            parprint('Iter #{:d}: '.format(iiter),end='')
+            if verbose and self.rank==0:
+                print('[gaussian_filter] iter #{:d}: '.format(iiter),end='')
             tbeg = MPI.Wtime()
             # Filter field: if key_out is None, perform operation inplace
             self.field[key_out] = self.field[key].gaussian_filter(sigma,
@@ -343,7 +347,8 @@ class PPP:
             self.impose_boundary_conditions(key_out)
             # Iterate inplace from now on
             key = key_out
-            parprint('{:g} sec'.format(MPI.Wtime()-tbeg))
+            if verbose and self.rank==0:
+                print('{:g} sec'.format(MPI.Wtime()-tbeg))
 
     def broadcast(self,key,operation,arg,key_out=None):
         '''Broadcasts an operation involving a scalar or matrix on 
@@ -541,9 +546,11 @@ class PPP:
         gc.collect()
         return
 
-    def save_state(self,file,parallel=False):
+    def save_state(self,file,parallel=False,verbose=True):
         import h5py
         from mpi4py import MPI
+        if verbose and self.rank==0:
+            print('[save_state] file={}'.format(file))
         tbeg = MPI.Wtime()
         ascii_type = h5py.string_dtype('ascii',32)
         # Only use parallel IO if flag is set and h5py has MPIIO support
@@ -624,11 +631,11 @@ class PPP:
         if not parallel: sb.proceed()
         self.comm.Barrier()
         tend = MPI.Wtime()
-        if self.rank==0:
-            print("[save_state] Elapsed time: {:f}".format(tend-tbeg))
+        if verbose and self.rank==0:
+            print("[save_state] elapsed time: {:f}".format(tend-tbeg))
         return
 
-    def save_for_vtk(self,file,key,stride=(1,1,1),truncate=True,merge_at_root=True,on_pressure_grid=True):
+    def save_for_vtk(self,file,key,stride=(1,1,1),truncate=True,merge_at_root=True,on_pressure_grid=True,verbose=True):
         '''Saves a field for visualization purposes. This means it will only have a single
         lower ghost cell if there is an upper neighbor, and both a single and an upper 
         ghost cell if there is no upper neighbor (or merged).'''
@@ -643,6 +650,8 @@ class PPP:
             return
         # Since the data is usually much smaller than a full 'save_state', I only
         # implement sequential IO for now.
+        if verbose and self.rank==0:
+            print('[save_for_vtk] key={}, file={}'.format(key,file))
         tbeg = MPI.Wtime()
         # If flag is set, shift data onto pressure grid first. Use a temporary field for this.
         name = key
@@ -666,8 +675,8 @@ class PPP:
         self.comm.Barrier()
         # Print timing
         tend = MPI.Wtime()
-        if self.rank==0:
-            print("[save_for_vtk] Elapsed time: {:f}".format(tend-tbeg))
+        if verbose and self.rank==0:
+            print("[save_for_vtk] elapsed time: {:f}".format(tend-tbeg))
         return
 
     def to_vtk(self,key,stride=(1,1,1),merge_at_root=False):

visu.py

@@ -1,67 +1,35 @@
-def add_domain_bounds(plotter,bounds,color='black',line_width=2):
+def add_domain_bounds(plotter,bounds,color='black',line_width=2.):
     import pyvista
     domain = pyvista.Box(bounds=bounds)
     plotter.add_mesh(domain,color=color,style='wireframe',line_width=line_width)
-    return 
-
-def add_particles(plotter,part,
-                        name=None,
-                        color='white',
-                        scalars=None,cmap=None,clim=None,
-                        opacity=1.0,
-                        theta_resolution=10,phi_resolution=10):
-    '''
-    name (str, optional) – The name for the added mesh/actor so that it can be easily updated. If an actor of this name already exists in the rendering window, it will be replaced by the new actor.
-
-    color (string or 3 item list, optional, defaults to white) – 
-    Use to make the entire mesh have a single solid color. Either a string, RGB list, or hex color string. For example: color='white', color='w', color=[1, 1, 1], or color='#FFFFFF'. Color will be overridden if scalars are specified.
-
-    scalars (str, optional) – 
-    Scalars used to “color” the mesh. Accepts a key of the col dictionary. If both color and scalars are None, then the active scalars are used.
-
-    cmap (str, list, optional) –
-    Name of the Matplotlib colormap to use when mapping the scalars. See available Matplotlib colormaps. Only applicable for when displaying scalars. Requires Matplotlib to be installed. colormap is also an accepted alias for this. If colorcet or cmocean are installed, their colormaps can be specified by name.
-    You can also specify a list of colors to override an existing colormap with a custom one. For example, to create a three color colormap you might specify ['green', 'red', 'blue']
-
-    clim (2 item list, optional) – 
-    Color bar range for scalars. Defaults to minimum and maximum of scalars array. Example: [-1, 2]. rng is also an accepted alias for this.
-
-    opacity (float, str, array-like) – 
-    Opacity of the mesh. If a single float value is given, it will be the global opacity of the mesh and uniformly applied everywhere - should be between 0 and 1. A string can also be specified to map the scalars range to a predefined opacity transfer function (options include: ‘linear’, ‘linear_r’, ‘geom’, ‘geom_r’). A string could also be used to map a scalars array from the mesh to the opacity (must have same number of elements as the scalars argument). Or you can pass a custom made transfer function that is an array either n_colors in length or shorter.
-    '''
-    import pyvista
-    assert part.has_attribute('r')
-    pts = pyvista.PolyData(part.get_position().transpose())
-    pts['r'] = part.attr['r']
-
-    if scalars:
-        assert part.has_attribute(scalars)
-        pts[scalars] = part.attr[scalars]
-        color=None
-
-    geom = pyvista.Sphere(radius=1,theta_resolution=theta_resolution,phi_resolution=phi_resolution)
-    glyphs = pts.glyph(scale='r',factor=1,geom=geom)
-    plotter.add_mesh(glyphs,name=name,
-        color=color,
-        scalars=scalars,cmap=cmap,clim=clim,
-        opacity=opacity)
     return
 
-def add_trajectories(plotter,traj,
-                        name=None,
-                        color='black',
-                        scalars=None,cmap=None,clim=None,
-                        opacity=1.0):
-    import pyvista
-    segments = traj.get_trajectories_segmented()
-    lines = pyvista.PolyData()
-    for part in segments:
-        for seg in part:
-            lines += pyvista.helpers.lines_from_points(seg.transpose())
-    plotter.add_mesh(lines,name=name,
-            color=color,
-            scalars=scalars,cmap=cmap,clim=clim,
-            opacity=opacity)
+def setup_camera(pl,bounds,viewvec,dist=None,parallel_projection=False,
+                rel_focus=(0.5,0.5,0.5),abs_focus=None,viewup=(0,1,0)):
+    import numpy as np
+    if not abs_focus:
+        focal_point = np.array(
+            (rel_focus[0]*(bounds[1]+bounds[0]),
+             rel_focus[1]*(bounds[3]+bounds[2]),
+             rel_focus[2]*(bounds[5]+bounds[4])))
+    else:
+        focal_point = np.array(abs_focus)
+    viewvec  = np.array(viewvec)
+    viewup   = np.array(viewup)
+    if dist is None:
+        pl.set_focus(tuple(focal_point))
+        pl.view_vector(tuple(-viewvec),tuple(viewup))
+    else:
+        position = focal_point-viewvec/np.linalg.norm(viewvec)*dist
+        # https://github.com/pyvista/pyvista-support/issues/40
+        pl.camera_position = (
+            tuple(position),
+            tuple(focal_point),
+            tuple(viewup))
+    if parallel_projection:
+        pl.enable_parallel_projection()
+    else:
+        pl.disable_parallel_projection()
     return
 
 def chunk_to_pvmesh(chunk,gridg):
@@ -82,15 +50,6 @@ def chunk_to_pvmesh(chunk,gridg):
     mesh.point_arrays['values'] = chunk['data'].flatten(order='F')  # Flatten the array!
     return mesh
 
-def field_to_pvmesh(field):
-    import pyvista
-    mesh = pyvista.UniformGrid()
-    mesh.dimensions = field.dim(axis=None)
-    mesh.origin     = field.origin  # The bottom left corner of the data set
-    mesh.spacing    = field.spacing # These are the cell sizes along each axis
-    mesh.point_arrays['values'] = field.data.flatten(order='F') # TBD: check order
-    return mesh
-
 def clip_out_of_bounds(mesh,bounds,axis=0,inplace=True,clip_lower=True,clip_upper=True):
     assert axis<3, "axis must be in [0,1,2]."
     new_bounds = bounds[2*axis:2*axis+2]
@@ -122,39 +81,6 @@ def clip_out_of_bounds(mesh,bounds,axis=0,inplace=True,clip_lower=True,clip_uppe
     else:
         return mesh
 
-# def common_mesh(gridg)
-#     import pyvista
-#     mesh = pyvista.UniformGrid()
-#     mesh.dimensions = (
-#         chunk['nxl']+2*chunk['ighost'],
-#         chunk['nyl']+2*chunk['ighost'],
-#         chunk['nzl']+2*chunk['ighost']
-#     )
-#     xg,yg,zg = gridg
-#     dx,dy,dz = xg[2]-xg[1],yg[2]-yg[1],zg[2]-zg[1]
-#     x0 = xg[chunk['ibeg']]-chunk['ighost']*dx
-#     y0 = yg[chunk['jbeg']]-chunk['ighost']*dy
-#     z0 = zg[chunk['kbeg']]-chunk['ighost']*dz
-#     mesh.origin  = (x0,y0,z0)  # The bottom left corner of the data set
-#     mesh.spacing = (dx,dy,dz)  # These are the cell sizes along each axis
-#     return mesh
-
-# def mesh_qcriterion(chunks_uvw,gridg_uvwp):
-#     for chunk in chunks_uvw:
-#         assert(chunk['ighost']==1)
-#     import pyvista
-#     dx = gridg_uvwp[0][0][1]-gridg_uvwp[0][0][1]
-#     dy = gridg_uvwp[0][1][1]-gridg_uvwp[0][1][1]
-#     dz = gridg_uvwp[0][2][1]-gridg_uvwp[0][2][1]
-#     xq0 = gridg_uvwp[3][0][chunks]
-
-
-#     xg,yg,zg = gridg
-#     dx,dy,dz = xg[2]-xg[1],yg[2]-yg[1],zg[2]-zg[1]
-#     x0 = xg[chunk['ibeg']]-chunk['ighost']*dx
-#     y0 = yg[chunk['jbeg']]-chunk['ighost']*dy
-#     z0 = zg[chunk['kbeg']]-chunk['ighost']*dz
-
 def translate_circular(pd,translation,bounds,axis=0):
     '''Translates pyvista PolyData objects while taking into account
     the bounding box'''