From e14702437498ce246a241a3d0fdddde83178acd4 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Fri, 30 Jul 2021 11:13:53 +0200 Subject: [PATCH] bugfix: symmetrization chose wrong mirror plane if the lower point coincided with the boundary; also changed default internal names of derivatives --- parallel.py | 32 +++++++++++++++++++++----------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/parallel.py b/parallel.py index b4da8ba..5b565bc 100644 --- a/parallel.py +++ b/parallel.py @@ -198,7 +198,7 @@ class PPP: def differentiate(self,key,axis,key_out=None,on_pressure_grid=False): assert axis<3, "'axis' must be one of 0,1,2." if key_out is None: - key_out = 'd'+key+('/dx','/dy','/dz')[axis] + key_out = 'd'+key+('dx','dy','dz')[axis] key_grid = self._key_grid(key) origin = list(self.origin[key_grid]) shifting_state = self.shifting_state(key,axis=axis) @@ -270,7 +270,7 @@ class PPP: keyx = ('x','y','z') keyu = ('u','v','w') for ii in range(3): - key_tmp = 'd'+keyu[ii]+'/d'+keyx[ii] + key_tmp = 'd'+keyu[ii]+'d'+keyx[ii] parprint(key_tmp) self.differentiate(keyu[ii],axis=ii,key_out=key_tmp,on_pressure_grid=True) self.field[key_out] += self.field[key_tmp]**2 @@ -278,8 +278,8 @@ class PPP: self.delete(key_tmp) # iip1 = (ii+1)%3 - key_tmp1 = 'd'+keyu[ii]+'/d'+keyx[iip1] - key_tmp2 = 'd'+keyu[iip1]+'/d'+keyx[ii] + key_tmp1 = 'd'+keyu[ii]+'d'+keyx[iip1] + key_tmp2 = 'd'+keyu[iip1]+'d'+keyx[ii] parprint(key_tmp1,key_tmp2) self.differentiate(keyu[ii],axis=iip1,key_out=key_tmp1,on_pressure_grid=True) self.differentiate(keyu[iip1],axis=ii,key_out=key_tmp2,on_pressure_grid=True) @@ -983,10 +983,15 @@ class PPP: for axis in range(3): if positionBd[axis]==-1: # lower boundary # index of first point within the domain - idx = self.field[key].nearest_gridpoint(self.bounds[2*axis],axis=axis,lower=True)+1 + idx = self.field[key].nearest_gridpoint(self.bounds[2*axis],axis=axis) + dist = self.field[key].coordinate(idx,axis=axis)-self.bounds[2*axis] + if np.abs(dist)>self.field[key].eps_collapse and np.sign(dist)<0.0: + idx += 1 + dist = self.field[key].coordinate(idx,axis=axis)-self.bounds[2*axis] # distance first point to wall: should be either 0 or dx/2 - dist = np.abs(self.field[key].coordinate(idx,axis=axis)-self.bounds[2*axis]) - if dist>=0.25*self.spacing[axis]: # no point on boundary + assert (dist<=self.field[key].eps_collapse or + (dist-0.5*self.spacing[axis])<=self.field[key].eps_collapse) + if np.abs(dist)>=0.25*self.spacing[axis]: # no point on boundary sl_dst[axis] = slice(0,idx,1) sl_src[axis] = slice(2*idx-1,idx-1,-1) else: # point on boundary @@ -995,10 +1000,15 @@ class PPP: break elif positionBd[axis]==1: # upper boundary # index of last point within the domain - idx = self.field[key].nearest_gridpoint(self.bounds[2*axis+1],axis=axis,lower=True) + idx = self.field[key].nearest_gridpoint(self.bounds[2*axis+1],axis=axis) + dist = self.field[key].coordinate(idx,axis=axis)-self.bounds[2*axis+1] + if np.abs(dist)>self.field[key].eps_collapse and np.sign(dist)>0.0: + idx -= 1 + dist = self.field[key].coordinate(idx,axis=axis)-self.bounds[2*axis+1] # distance last point to wall: should be either 0 or -dx/2 - dist = np.abs(self.field[key].coordinate(idx,axis=axis)-self.bounds[2*axis+1]) - if dist>=0.25*self.spacing[axis]: # no point on boundary + assert (dist<=self.field[key].eps_collapse or + (dist+0.5*self.spacing[axis])<=self.field[key].eps_collapse) + if np.abs(dist)>=0.25*self.spacing[axis]: # no point on boundary sl_dst[axis] = slice(idx+1,self.field[key].dim(axis=axis),1) sl_src[axis] = slice(idx,2*idx+1-self.field[key].dim(axis=axis),-1) else: # point on boundary @@ -1218,4 +1228,4 @@ def is_root(comm=None,root=0): def gather(data,comm=None): from mpi4py import MPI comm = MPI.COMM_WORLD if comm is None else comm - return comm.gather(data,root=0) \ No newline at end of file + return comm.gather(data,root=0)