mwtkrayer
/
suspendtools
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

bugfix: symmetrization chose wrong mirror plane if the lower point coincided with the boundary; also changed default internal names of derivatives

This commit is contained in:
Michael Krayer 2021-07-30 11:13:53 +02:00
parent 352b3ae860
commit e147024374
1 changed files with 21 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
parallel.py
View File

@ -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)
return comm.gather(data,root=0)