Commit 75fdc106 authored by ylan's avatar ylan
Browse files

[350k] add missing files and update chk_txt.py and README

- fix chk_txt.py
parent 61451fe2
......@@ -31,10 +31,14 @@ ann352625 N, Version 3.3, E98,782,067
Rmin 24.1010
Rmax 82.3009
dtouch 1.997404 << min distance between any two spheres' centers
dsep 1.999376 << "average" length to seperate spheres, = kissing number at 2.5
id: s1 s2 1874160 11416 11410 << find min sphere pair, for fun?
s1 loc: [-76.8514634 -25.9385264 -48.2381829]
s2 loc: [-77.9093191 -26.3022988 -49.8929437]
chk: 1.9974039016285345 1.9974039016285345
dsep 2.000000 << "average" length to seperate spheres, = kissing number at 2.5
```
(`modulde load python` on Summit)
(It takes some time to compute dtouch)
(It takes some time to compute dtouch. I believe it's a O(N log N) algo.)
dsep is about to 2, which makes sense for the radius=0.91945
......
......@@ -44,10 +44,11 @@ print('Rmax %10.4f'%max(rad_cyl))
# cheap pdist
def comp_dist(P):
tri = Delaunay(P, qhull_options="Qbb C-0")
arr1 = P[tri.simplices[:,1],:] - P[tri.simplices[:,0],:]
arr2 = P[tri.simplices[:,2],:] - P[tri.simplices[:,1],:]
arr3 = P[tri.simplices[:,0],:] - P[tri.simplices[:,2],:]
arr4 = P[tri.simplices[:,3],:] - P[tri.simplices[:,0],:]
arr2 = P[tri.simplices[:,2],:] - P[tri.simplices[:,0],:]
arr3 = P[tri.simplices[:,3],:] - P[tri.simplices[:,0],:]
arr4 = P[tri.simplices[:,2],:] - P[tri.simplices[:,1],:]
arr5 = P[tri.simplices[:,3],:] - P[tri.simplices[:,1],:]
arr6 = P[tri.simplices[:,3],:] - P[tri.simplices[:,2],:]
......@@ -57,11 +58,33 @@ def comp_dist(P):
arr=np.append(arr,arr4,axis=0)
arr=np.append(arr,arr5,axis=0)
arr=np.append(arr,arr6,axis=0)
arr,uid=np.unique(arr,axis=0,return_index=True)
dist=np.sqrt(arr[:,0]**2+arr[:,1]**2+arr[:,2]**2)
return dist
dist = comp_dist(centers)
pr1 = np.array([tri.simplices[:,1],tri.simplices[:,0]])
pr2 = np.array([tri.simplices[:,2],tri.simplices[:,0]])
pr3 = np.array([tri.simplices[:,3],tri.simplices[:,0]])
pr4 = np.array([tri.simplices[:,2],tri.simplices[:,1]])
pr5 = np.array([tri.simplices[:,3],tri.simplices[:,1]])
pr6 = np.array([tri.simplices[:,3],tri.simplices[:,2]])
pr = pr1
pr = np.append(pr,pr2,axis=1)
pr = np.append(pr,pr3,axis=1)
pr = np.append(pr,pr4,axis=1)
pr = np.append(pr,pr5,axis=1)
pr = np.append(pr,pr6,axis=1)
pr = pr[:,uid]
return dist,pr
dist,pr = comp_dist(centers)
print('dtouch %f'%min(dist))
im = np.argmin(dist)
print(' id: s1 s2 ',im,pr[0,im],pr[1,im])
print(' s1 loc: ',centers[pr[0,im],:])
print(' s2 loc: ',centers[pr[1,im],:])
arr=centers[pr[0,im],:]-centers[pr[1,im],:]
print(' chk: ',np.sqrt(arr[0]**2+arr[1]**2+arr[2]**2),dist[im])
dist=np.sort(dist); Nplt=dist.shape[0];
tx = np.linspace(1,Nplt,num=Nplt,endpoint=True); tx = tx/Nsph*2
......
c dependencies
c my_outpost.f
c my_outpost
c-----------------------------------------------------------------------
subroutine my_mesh_metric
call xm1toxc
call mesh_metrics ! read jacm, xc
return
end
c-----------------------------------------------------------------------
subroutine my_mshchk(s6,ifdump)
implicit none
include 'SIZE'
include 'TOTAL'
character*6 s6
integer ierr,ierr1,ierr2,ierr3
integer e,nel,nxyz,iwk(lelt),elist(lelt)
integer inei,nnei,nplt0,nplt1,iglsum
real tmp(lx1,ly1,lz1,lelt),dummy(lx1,ly1,lz1,lelt), vlmax
logical ifdump
nel=nelv
nxyz=lx1*ly1*lz1
nnei=3
ierr=0
call izero(elist,lelt)
! rhs
call izero(iwk,lelt)
call my_chkrhs(ierr2,iwk)
if(ierr2.gt.0) then
ierr=1
call my_iadd2(elist,iwk,nel)
endif
! jac
call izero(iwk,lelt)
call my_chkjac(ierr3,iwk)
if(ierr3.gt.0) then
ierr=1
call my_iadd2(elist,iwk,nel)
endif
! summary
if(ierr.eq.0) then
if(nid.eq.0)write(*,*)'PASSED mshchk ',s6
else
if(nid.eq.0)write(*,*)'FAILED mshchk ',s6
$ ,' rhs fail',ierr2,' neg jac ',ierr3
call rzero(dummy,nel*nxyz)
do e=1,nel
if(elist(e).gt.0) then
call cfill(dummy(1,1,1,e),lglel(e)*1.0,nxyz)
endif
enddo
endif
! find neighbors of bad elements
if(ifdump.AND.ierr.gt.0) then
do e=1,nel
if(elist(e).gt.0) elist(e)=1
enddo
nplt0=iglsum(elist,nel)
do inei=1,nnei
call rzero(tmp,nxyz*nel)
do e=1,nel
if(elist(e).gt.0) call rone(tmp(1,1,1,e),nxyz)
enddo
call dssum(tmp,lx1,ly1,lz1) ! use dssum to extend neighbor
do e=1,nel
if(vlmax(tmp(1,1,1,e),nxyz).gt.0) elist(e)=1
enddo
do e=1,nel
if(elist(e).gt.0) elist(e)=1
enddo
nplt1=iglsum(elist,nel)
if(nid.eq.0)write(*,*)'err exi',inei,nplt1
enddo
if(nid.eq.0)write(*,*)'err dmp',nnei,nplt0,nplt1
ifxyo=.true.
call my_outpost(elist,dummy,vy,vz,pr,t,'err')
endif
50 format(a,i3,i6,i6,i12,1p3e16.8)
return
end
c-----------------------------------------------------------------------
subroutine my_mshchk_vb(s6,ifdump)
implicit none
include 'SIZE'
include 'TOTAL'
character*6 s6
integer ierr,ierr1,ierr2,ierr3
integer e,nel,nxyz,iwk(lelt),elist(lelt),ic
integer inei,nnei,nplt0,nplt1,iglsum
real tmp(lx1,ly1,lz1,lelt), dummy(lx1,ly1,lz1,lelt), vlmax
logical ifdump
nel=nelv
nxyz=lx1*ly1*lz1
nnei=3 ! # neighbor for dumping files
ierr=0
call izero(elist,lelt)
! rhs
call izero(iwk,lelt)
call my_chkrhs(ierr2,iwk)
if(ierr2.gt.0) then
ierr=1
call my_iadd2(elist,iwk,nel)
do e=1,nel
if(iwk(e).gt.0) then
write(*,51)'Eerr rhs',istep,nid,lglel(e),int(bc(5,5,e,1))
$ ,xm1(1,1,1,e),ym1(1,1,1,e),zm1(1,1,1,e)
$ ,' ',(cbc(ic,e,1),ic=1,6)
endif
enddo
endif
! jac
call izero(iwk,lelt)
call my_chkjac(ierr3,iwk)
if(ierr3.gt.0) then
ierr=1
call my_iadd2(elist,iwk,nel)
do e=1,nel
if(iwk(e).gt.0) then
write(*,51)'Eerr jac',istep,nid,lglel(e),int(bc(5,5,e,1))
$ ,xm1(1,1,1,e),ym1(1,1,1,e),zm1(1,1,1,e)
$ ,' ',(cbc(ic,e,1),ic=1,6)
endif
enddo
endif
! summary
if(ierr.eq.0) then
if(nid.eq.0)write(*,*)'PASSED mshchk ',s6
else
if(nid.eq.0)write(*,*)'FAILED mshchk ',s6
$ ,' rhs fail',ierr2,' neg jac ',ierr3
call rzero(dummy,nel*nxyz)
do e=1,nel
if(elist(e).gt.0) then
call cfill(dummy(1,1,1,e),lglel(e)*1.0,nxyz)
endif
enddo
endif
! find neighbors of bad elements
if(ifdump.AND.ierr.gt.0) then
do e=1,nel
if(elist(e).gt.0) elist(e)=1
enddo
nplt0=iglsum(elist,nel)
do inei=1,nnei
call rzero(tmp,nxyz*nel)
do e=1,nel
if(elist(e).gt.0) call rone(tmp(1,1,1,e),nxyz)
enddo
call dssum(tmp,lx1,ly1,lz1)
do e=1,nel
if(vlmax(tmp(1,1,1,e),nxyz).gt.0) elist(e)=1
enddo
do e=1,nel
if(elist(e).gt.0) elist(e)=1
enddo
nplt1=iglsum(elist,nel)
if(nid.eq.0)write(*,*)'err exi',inei,nplt1
enddo
if(nid.eq.0)write(*,*)'err dmp',nnei,nplt0,nplt1
ifxyo=.true.
call my_outpost(elist,vx,vy,dummy,pr,t,'err')
endif
50 format(a,i3,i6,i6,i12,1p3e16.8)
51 format(a,i3,i6,i12,i3,1p3e16.8,a1,6a3)
return
end
c-----------------------------------------------------------------------
subroutine my_chkjac(ierr,elist)
c Check negative Jacobian
c Note, Nek only requires Jac in the same sign of V(1)
c At here, we want Jac>0 for all vertices
implicit none
include 'SIZE'
integer e,nel,nxyz,ierr,iglsum,elist(lelt)
real Jac(lx1*ly1*lz1,lelt),vlmin
nel=nelv
nxyz=lx1*ly1*lz1
ierr=0
call izero(elist,nel)
if(ldim.eq.2)return
call chk_comp_Jac0(Jac)
do e=1,nel
if(vlmin(Jac(1,e),nxyz).lt.0.0) then
elist(e)=1
ierr=ierr+1
endif
enddo
ierr=iglsum(ierr,1)
return
end
c-----------------------------------------------------------------------
subroutine my_chkrhs(ierr,elist)
c Check right-handness
implicit none
include 'SIZE' ! nelv
include 'GEOM' ! xm1
include 'INPUT'! xc
include 'SCRCT'! xyz
integer ierr,ie,nel,j,ivtx,iglsum,elist(lelt)
real VOLUM0,v1,v2,v3,v4,v5,v6,v7,v8
integer indx(8) ! local usage, todo: include TOPOL?
data indx /1,2,4,3,5,6,8,7/
nel=nelv
ierr=0
call izero(elist,nel)
if(ldim.eq.2)return
call xm1toxc
call rzero(xyz,24*nel)
do ie=1,nel
do j=1,8
ivtx = indx(j)
xyz(1,ivtx,ie) = xc(j,ie)
xyz(2,ivtx,ie) = yc(j,ie)
xyz(3,ivtx,ie) = zc(j,ie)
enddo
V1= VOLUM0(XYZ(1,2,IE),XYZ(1,3,IE),XYZ(1,5,IE),XYZ(1,1,IE))
V2= VOLUM0(XYZ(1,4,IE),XYZ(1,1,IE),XYZ(1,6,IE),XYZ(1,2,IE))
V3= VOLUM0(XYZ(1,1,IE),XYZ(1,4,IE),XYZ(1,7,IE),XYZ(1,3,IE))
V4= VOLUM0(XYZ(1,3,IE),XYZ(1,2,IE),XYZ(1,8,IE),XYZ(1,4,IE))
V5=-VOLUM0(XYZ(1,6,IE),XYZ(1,7,IE),XYZ(1,1,IE),XYZ(1,5,IE))
V6=-VOLUM0(XYZ(1,8,IE),XYZ(1,5,IE),XYZ(1,2,IE),XYZ(1,6,IE))
V7=-VOLUM0(XYZ(1,5,IE),XYZ(1,8,IE),XYZ(1,3,IE),XYZ(1,7,IE))
V8=-VOLUM0(XYZ(1,7,IE),XYZ(1,6,IE),XYZ(1,4,IE),XYZ(1,8,IE))
if (V1.LE.0.0.OR.V2.LE.0.0.OR.
$ V3.LE.0.0.OR.V4.LE.0.0.OR.
$ V5.LE.0.0.OR.V6.LE.0.0.OR.
$ V7.LE.0.0.OR.V8.LE.0.0 ) then
elist(ie)=1
ierr=ierr+1
endif
enddo
ierr=iglsum(ierr,1)
return
end
c-----------------------------------------------------------------------
subroutine chk_comp_Jac0(Jac)
implicit none
include 'SIZE'
integer nel,ntot,i
real Jac(lx1*ly1*lz1*lelt),J(lx1*ly1*lz1*lelt,ldim*ldim)
nel=nelv
ntot=lx1*ly1*lz1*nel
call xyzrst(J(1,1),J(1,2),J(1,3),J(1,4),J(1,5),J(1,6)
$ ,J(1,7),J(1,8),J(1,9),.false.) ! mesh 1
do i=1,ntot
Jac(i) = J(i,1)*J(i,5)*J(i,9)
$ + J(i,7)*J(i,2)*J(i,6)
$ + J(i,4)*J(i,8)*J(i,3)
$ - J(i,1)*J(i,8)*J(i,6)
$ - J(i,4)*J(i,2)*J(i,9)
$ - J(i,7)*J(i,5)*J(i,3)
enddo
return
end
c-----------------------------------------------------------------------
subroutine xm1toxc ! xc: vertices
implicit none
include 'SIZE' ! nelv
include 'GEOM' ! xm1
include 'INPUT'! xc
integer ie,nel
nel=nelv
do ie=1,nel
xc(1,ie) = XM1(1 ,1 ,1 ,ie)
xc(2,ie) = XM1(lx1,1 ,1 ,ie)
xc(3,ie) = XM1(lx1,ly1,1 ,ie)
xc(4,ie) = XM1(1 ,ly1,1 ,ie)
xc(5,ie) = XM1(1 ,1 ,lz1,ie)
xc(6,ie) = XM1(lx1,1 ,lz1,ie)
xc(7,ie) = XM1(lx1,ly1,lz1,ie)
xc(8,ie) = XM1(1 ,ly1,lz1,ie)
! y
yc(1,ie) = YM1(1 ,1 ,1 ,ie)
yc(2,ie) = YM1(lx1,1 ,1 ,ie)
yc(3,ie) = YM1(lx1,ly1,1 ,ie)
yc(4,ie) = YM1(1 ,ly1,1 ,ie)
yc(5,ie) = YM1(1 ,1 ,lz1,ie)
yc(6,ie) = YM1(lx1,1 ,lz1,ie)
yc(7,ie) = YM1(lx1,ly1,lz1,ie)
yc(8,ie) = YM1(1 ,ly1,lz1,ie)
! z
zc(1,ie) = ZM1(1 ,1 ,1 ,ie)
zc(2,ie) = ZM1(lx1,1 ,1 ,ie)
zc(3,ie) = ZM1(lx1,ly1,1 ,ie)
zc(4,ie) = ZM1(1 ,ly1,1 ,ie)
zc(5,ie) = ZM1(1 ,1 ,lz1,ie)
zc(6,ie) = ZM1(lx1,1 ,lz1,ie)
zc(7,ie) = ZM1(lx1,ly1,lz1,ie)
zc(8,ie) = ZM1(1 ,ly1,lz1,ie)
enddo
return
end
c-----------------------------------------------------------------------
subroutine xm1toxyz
implicit none
include 'SIZE' ! nelv
include 'GEOM' ! xm1
include 'INPUT'! xc
include 'SCRCT'! xyz
integer ierr,ie,nel,j,ivtx
integer indx(8) ! local usage, todo: include TOPOL?
data indx /1,2,4,3,5,6,8,7/
nel=nelv
call xm1toxc
call rzero(xyz,24*nel)
do ie=1,nel
do j=1,8
ivtx = indx(j)
xyz(1,ivtx,ie) = xc(j,ie)
xyz(2,ivtx,ie) = yc(j,ie)
xyz(3,ivtx,ie) = zc(j,ie)
enddo
enddo
return
end
c-----------------------------------------------------------------------
subroutine my_iadd2(i1,i2,n)
DIMENSION I1(1),I2(1)
C
DO 10 I=1,N
I1(I)=I1(I)+I2(I)
10 CONTINUE
return
END
c-----------------------------------------------------------------------
subroutine my_prt_lglel
include 'SIZE'
include 'INPUT'
include 'PARALLEL'
include 'SCRCT'
include 'SOLN'
include 'TSTEP'
include 'CTIMER'
logical ifverbm
C Output the processor-element map:
ifverbm=.false.
if (loglevel .gt. 2) ifverbm=.true.
loglevel_bak=loglevel
loglevel=3
ifverbm=.true.
if(ifverbm) then
idum = 1
if(nid.eq.0) then
N8 = min(8,nelt)
write(6 ,1310) node-1,(lglel(ie),ie=1,n8)
if (NELT.GT.8) write(6 ,1315) (lglel(ie),ie=9,NELT)
if (NELT.GT.8) write(6 ,*)' '
DO inid=1,NP-1
mtype = inid
call csend(mtype,idum,4,inid,0) ! handshake
call crecv(mtype,inelt,4) ! nelt of other cpus
N8 = min(8,inelt)
ENDDO
1310 FORMAT(' RANK',I6,' IEG',8I12)
1315 FORMAT(' ',6X,' ',8I12)
else
mtype = nid
call crecv(mtype,idum,4) ! hand-shake
call csend(mtype,nelt,4,0,0) ! nelt
if (loglevel .gt. 2) then
N8 = min(8,nelt)
write(6 ,1310) node-1,(lglel(ie),ie=1,n8)
if (NELT.GT.8) write(6 ,1315) (lglel(ie),ie=9,NELT)
if (NELT.GT.8) write(6 ,*)' '
endif
endif
endif
loglevel=loglevel_bak
return
end
c----------------------------------------------------------------------
subroutine my_dump_mesh(s3,imode)
include 'SIZE'
include 'SOLN'
include 'GEOM'
include 'INPUT'
character*3 s3
integer imid
call my_mshchk('bfrdmp',.true.) ! last chk before dump
call domain_size(xmn,xmx,ymn,ymx,zmn,zmx)
if (nio.eq.0) then
write(6,18) xmn,ymn,zmn
write(6,19) xmx,ymx,zmx
18 format('new xyz min ',5g13.5)
19 format('new xyz max ',5g13.5)
endif
call my_mesh_metric
ifxyo = .true.
call outpost(vx,vy,vz,pr,t,s3)
if (abs(imode).eq.1) then
imid=1
call gen_rea(imid)
if (nid.eq.0) write(6,*)'done genrea'
call gen_re2(imid) ! FIXME, bug for large E (in bc)
if (nid.eq.0) write(6,*)'done genre2'
elseif (abs(imode).eq.2) then
call my_gen_re2
if (nid.eq.0) write(6,*)'done genre2 (xyz only)'
call gen_re2(mid) ! FIXME, bug for large E (in bc) ! TODO for dbg
if (nid.eq.0) write(6,*)'done genre2'
endif
return
end
c-----------------------------------------------------------------------
subroutine my_gen_re2 ! <case>_xyz.re2
implicit none
include 'SIZE'
include 'TOTAL'
integer imid
character*80 hdr
real*4 test
data test / 6.54321 /
integer ierr
character*32 ffname
character*132 datfle
character*1 datfle1(132)
equivalence (datfle,datfle1)
integer lfname,ltrunc
logical ifdat
ierr=0
if(nid.eq.0) then
lfname = ltrunc(reafle,132) - 4
call blank (datfle,132)
call chcopy(datfle,reafle,lfname)
call chcopy(datfle1(lfname+1),'_xyz.re2',8)
write(*,*)'my_gen_re2: ',wdsize,datfle
call byte_open(datfle, ierr)