| 1 | subroutine meshsmoother
|
|---|
| 2 | include 'SIZE'
|
|---|
| 3 | include 'TOTAL'
|
|---|
| 4 |
|
|---|
| 5 | parameter(ndfsbc=1) !number of boundary conditions
|
|---|
| 6 | character*3 dfsbc(ndfsbc)
|
|---|
| 7 | save dfsbc
|
|---|
| 8 | data dfsbc /'W '/ !BCs listed here
|
|---|
| 9 |
|
|---|
| 10 | idftyp = 0 !distance function - 0 -> exponential, 1-> tanh
|
|---|
| 11 | alpha = 15. !Input for wall distance function
|
|---|
| 12 | beta = 0.1 !
|
|---|
| 13 |
|
|---|
| 14 | nouter = 50 !total loops around laplacian and optimizer smoothing
|
|---|
| 15 | nlap = 20 !number of laplacian iterations in each loop
|
|---|
| 16 | nopt = 20 !number of optimization iterations in each loop
|
|---|
| 17 |
|
|---|
| 18 | mtyp = 1 !metric type
|
|---|
| 19 |
|
|---|
| 20 | call smoothmesh(mtyp,nouter,nlap,nopt,ndfsbc,dfsbc,idftyp,alpha,
|
|---|
| 21 | $ beta)
|
|---|
| 22 |
|
|---|
| 23 | return
|
|---|
| 24 | end
|
|---|
| 25 | c-----------------------------------------------------------------------
|
|---|
| 26 | subroutine smoothmesh(mtyp,nouter,nlap,nopt,nbc,dcbc,
|
|---|
| 27 | $ idftyp,alpha,beta)
|
|---|
| 28 | include 'SIZE'
|
|---|
| 29 | include 'TOTAL'
|
|---|
| 30 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 31 | parameter (nxyzc=lxc*lyc*lzc,nxyz=lx1*ly1*lz1)
|
|---|
| 32 |
|
|---|
| 33 | real xmc(lxc,lyc,lzc,lelv),
|
|---|
| 34 | $ ymc(lxc,lyc,lzc,lelv),
|
|---|
| 35 | $ zmc(lxc,lyc,lzc,lelv),
|
|---|
| 36 | $ mltc(lxc,lyc,lzc,lelv)
|
|---|
| 37 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 38 |
|
|---|
| 39 | real x8(2**ldim,lelv*(2**ldim)),
|
|---|
| 40 | $ y8(2**ldim,lelv*(2**ldim)),
|
|---|
| 41 | $ z8(2**ldim,lelv*(2**ldim)),
|
|---|
| 42 | $ nodmask(2**ldim,lelv*(2**ldim)),
|
|---|
| 43 | $ dis((2**ldim)*lelv*(2**ldim)),
|
|---|
| 44 | $ mlt((2**ldim)*lelv*(2**ldim))
|
|---|
| 45 | common /coarsemesh8/ x8,y8,z8
|
|---|
| 46 |
|
|---|
| 47 | real dx(nxyzc,lelv),dy(nxyzc,lelv),dz(nxyzc,lelv),
|
|---|
| 48 | $ xbp(nxyzc,lelv),ybp(nxyzc,lelv),zbp(nxyzc,lelv)
|
|---|
| 49 | common / msmbackup / dx,dy,dz,xbp,ybp,zbp
|
|---|
| 50 |
|
|---|
| 51 | real dxm(nxyz,lelv),dym(nxyz,lelv),dzm(nxyz,lelv),
|
|---|
| 52 | $ dxn(nxyz,lelv),dyn(nxyz,lelv),dzn(nxyz,lelv)
|
|---|
| 53 |
|
|---|
| 54 | integer opt,gshl,gshlc,lapinv,optinv
|
|---|
| 55 | real alpha,beta,etstart, etend,f1sav,f1,f1pre
|
|---|
| 56 | character*3 dcbc(nbc)
|
|---|
| 57 | c -------------------
|
|---|
| 58 | c -------------------
|
|---|
| 59 | ccc INITIALIZE VARIABLES
|
|---|
| 60 | lapinv = 0
|
|---|
| 61 | optinv = 0
|
|---|
| 62 |
|
|---|
| 63 | if (nid.eq.0.and.loglevel.ge.5)
|
|---|
| 64 | $ write(6,*) 'SMOOTHER-check original mesh'
|
|---|
| 65 | call fix_geom
|
|---|
| 66 | c Copy original mesh from xm1 to dxm
|
|---|
| 67 | call copy(dxm,xm1,lx1*ly1*lz1*nelv)
|
|---|
| 68 | call copy(dym,ym1,lx1*ly1*lz1*nelv)
|
|---|
| 69 | if (ldim.eq.3) call copy(dzm,zm1,lx1*ly1*lz1*nelv)
|
|---|
| 70 |
|
|---|
| 71 | ccc Generate interpolators for going to lx1 = 3
|
|---|
| 72 | call gen_int_lx1_to_3
|
|---|
| 73 |
|
|---|
| 74 | ccc COPY THE ORIGINAL MESH TO dx,dy,dz vectors
|
|---|
| 75 | call copy(dx,xmc,lxc*lyc*lzc*nelv)
|
|---|
| 76 | call copy(dy,ymc,lxc*lyc*lzc*nelv)
|
|---|
| 77 | if (ldim.eq.3) call copy(dz,zmc,lxc*lyc*lzc*nelv)
|
|---|
| 78 | ccc COPY THE ORIGINAL MESH FOR BACKUP IF MESH BECOMES INVERTED
|
|---|
| 79 | call copy(xbp,xmc,lxc*lyc*lzc*nelt)
|
|---|
| 80 | call copy(ybp,ymc,lxc*lyc*lzc*nelt)
|
|---|
| 81 | if (ldim.eq.3) call copy(zbp,zmc,lxc*lyc*lzc*nelt)
|
|---|
| 82 |
|
|---|
| 83 | call xmtox8(xmc,x8)
|
|---|
| 84 | call xmtox8(ymc,y8)
|
|---|
| 85 | if (ldim.eq.3) call xmtox8(zmc,z8)
|
|---|
| 86 |
|
|---|
| 87 | ccc CREATE MASK
|
|---|
| 88 | call genmask(nodmask,mlt,gshl,mltc,gshlc)
|
|---|
| 89 |
|
|---|
| 90 | ccc CONSTRUCT WEIGHT FUNCTION
|
|---|
| 91 | call disfun(dis,idftyp,alpha,beta,dcbc,nbc)
|
|---|
| 92 |
|
|---|
| 93 | etstart=dnekclock()
|
|---|
| 94 | ccc GET INITIAL ENERGY
|
|---|
| 95 | call getglobsum(x8,y8,z8,mlt,gshl,2**ldim,1,f1sav)
|
|---|
| 96 | if (nid.eq.0)
|
|---|
| 97 | $ write(6,'(A,1p1e13.6)') 'SMOOTHER-initial energy',f1sav
|
|---|
| 98 | f1 = f1sav
|
|---|
| 99 |
|
|---|
| 100 | ccc START SMOOTHING HERE
|
|---|
| 101 | do j=1,nouter
|
|---|
| 102 | f1pre = f1
|
|---|
| 103 | if (nid.eq.0.and.loglevel.ge.2)
|
|---|
| 104 | $ write(6,'(A,I5)') 'SMOOTHER-iteration',j
|
|---|
| 105 | mtyp = 1 !if mtyp = 1, jacobian, 2 = l^2, 3 - scale jacobian
|
|---|
| 106 | if (nlap.gt.0) then
|
|---|
| 107 | call fastlap(nlap,nodmask,mlt,gshl,dis,lapinv,mltc,gshlc)
|
|---|
| 108 | if (lapinv.eq.1) nlap = 0
|
|---|
| 109 | if (lapinv.eq.1) call restbackmesh
|
|---|
| 110 | endif
|
|---|
| 111 | if (nopt.gt.0) then
|
|---|
| 112 | call optCG(nopt,nodmask,mlt,gshl,dis,mtyp,optinv,mltc,gshlc)
|
|---|
| 113 | if (optinv.eq.1) call restbackmesh !xbp->xm1 and xbp->x8
|
|---|
| 114 | if (optinv.eq.1) nopt = 0
|
|---|
| 115 | endif
|
|---|
| 116 |
|
|---|
| 117 | call getglobsum(x8,y8,z8,mlt,gshl,2**ldim,mtyp,f1)
|
|---|
| 118 |
|
|---|
| 119 | call xmtox8(xmc,x8)
|
|---|
| 120 | call xmtox8(ymc,y8)
|
|---|
| 121 | if (ldim.eq.3) call xmtox8(zmc,z8)
|
|---|
| 122 |
|
|---|
| 123 | call xmctoxm1
|
|---|
| 124 |
|
|---|
| 125 | if (nid.eq.0) write(6,'(A,I7,1p1e13.5)') 'SMOOTHER-energy',j,f1
|
|---|
| 126 | if (nid.eq.0.and.loglevel.ge.2) write(6,*) 'loop complete',j
|
|---|
| 127 | if (f1.ge.f1pre) goto 5001 !mesh didn't improve since last iteration
|
|---|
| 128 | if (nopt.eq.0.and.nlap.eq.0) goto 5001 !no iterations left to do
|
|---|
| 129 | enddo
|
|---|
| 130 | 5001 continue
|
|---|
| 131 |
|
|---|
| 132 | etend=dnekclock()
|
|---|
| 133 | ccc RESTORE BOUNDARY LAYER
|
|---|
| 134 | call restbndrlayer(dx,dy,dz,dis,mltc,gshlc) !dx,dy,dz is now actually smooth-coarse
|
|---|
| 135 | call fix_geom
|
|---|
| 136 |
|
|---|
| 137 | ccc Solve the Laplace's equation to morph the interior mesh to match
|
|---|
| 138 | ccc the actual surface mesh
|
|---|
| 139 | call opsub3(dxn,dyn,dzn,dxm,dym,dzm,xm1,ym1,zm1)
|
|---|
| 140 | call my_mv_mesh(dxn,dyn,dzn)
|
|---|
| 141 | call opadd2(xm1,ym1,zm1,dxn,dyn,dzn)
|
|---|
| 142 |
|
|---|
| 143 | call getglobsum(x8,y8,z8,mlt,gshl,2**ldim,mtyp,f1)
|
|---|
| 144 | if (nid.eq.0) then
|
|---|
| 145 | write(6,'(A,1p1e13.6)') 'SMOOTHER-initial energy ',f1sav
|
|---|
| 146 | write(6,'(A,1p1e13.6)') 'SMOOTHER-final energy ',f1
|
|---|
| 147 | write(6,'(A,f5.2)') 'SMOOTHER-improve % ',100.*(f1sav-f1)/f1sav
|
|---|
| 148 | write(6,'(A,1p1e13.6)') 'SMOOTHER-time taken ',etend-etstart
|
|---|
| 149 | endif
|
|---|
| 150 |
|
|---|
| 151 | call geom_reset(1) ! recompute Jacobians, etc.
|
|---|
| 152 |
|
|---|
| 153 | ccc OUTPUT THE MESH
|
|---|
| 154 | c call gen_rea_full(1) !output the rea for smooth mesh
|
|---|
| 155 |
|
|---|
| 156 | return
|
|---|
| 157 | end
|
|---|
| 158 | c-----------------------------------------------------------------------
|
|---|
| 159 | subroutine gen_int_lx1_to_3 !interpolate mesh from lx1 to lx1=3
|
|---|
| 160 | include 'SIZE'
|
|---|
| 161 | include 'TOTAL'
|
|---|
| 162 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 163 | real dxmc(3,3),dymc(3,3),dzmc(3,3)
|
|---|
| 164 | real dxtmc(3,3),dytmc(3,3),dztmc(3,3)
|
|---|
| 165 | common /coarseders/ dxmc,dymc,dzmc,dxtmc,dytmc,dztmc
|
|---|
| 166 | real ixtmc3(lx1,3),ixmc3(3,lx1)
|
|---|
| 167 | real ixtmf3(3,lx1),ixmf3(lx1,3)
|
|---|
| 168 | common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
|
|---|
| 169 | real zgmc(3),wxmc(3) !points and weights
|
|---|
| 170 | real zgmf(lx1),wxmf(lx1) !points and weights
|
|---|
| 171 | real w(lx1,lx1)
|
|---|
| 172 |
|
|---|
| 173 | c Generate GLL points and weight
|
|---|
| 174 | call zwgll(zgmc,wxmc,3)
|
|---|
| 175 | call zwgll(zgmf,wxmf,lx1)
|
|---|
| 176 | c Generate derivative matrices with transpose operators
|
|---|
| 177 | call dgll(dxmc,dxtmc,zgmc,3,3)
|
|---|
| 178 | call dgll(dymc,dytmc,zgmc,3,3)
|
|---|
| 179 | call rzero(dzmc,3*3)
|
|---|
| 180 | call rzero(dztmc,3*3)
|
|---|
| 181 | if (ldim.eq.3) call dgll(dzmc,dztmc,zgmc,3,3)
|
|---|
| 182 | c Generate interpolation matrices
|
|---|
| 183 | call igllm (ixmc3,ixtmc3,zgmf,zgmc,nx1,3,nx1,3)
|
|---|
| 184 | call igllm (ixmf3,ixtmf3,zgmc,zgmf,3,nx1,3,nx1)
|
|---|
| 185 |
|
|---|
| 186 | c Interpolate mesh to lx1=3
|
|---|
| 187 | call xm1toxmc
|
|---|
| 188 |
|
|---|
| 189 | return
|
|---|
| 190 | end
|
|---|
| 191 | c-----------------------------------------------------------------------
|
|---|
| 192 | subroutine int_fine_to_coarse_2d(u1,u2,n1,n2)
|
|---|
| 193 | include 'SIZE'
|
|---|
| 194 | c Interpolate u1 to u2 using interpolation matrices in12 and in12^T
|
|---|
| 195 | c u1 is size n1xn1 and u2 is size n2xn2
|
|---|
| 196 | real ixtmc3(lx1,3),ixmc3(3,lx1)
|
|---|
| 197 | real ixtmf3(3,lx1),ixmf3(lx1,3)
|
|---|
| 198 | common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
|
|---|
| 199 | real u1(n1,n1),u2(n2,n2)
|
|---|
| 200 | real w(20,20)
|
|---|
| 201 |
|
|---|
| 202 | if (lx1.gt.20) write(6,*) 'SMOOTHER-increase size of work array
|
|---|
| 203 | $ in int_fine_to_coarse and related routines '
|
|---|
| 204 | if (lx1.gt.20) call exitt
|
|---|
| 205 |
|
|---|
| 206 | call mxm(ixmc3,n2,u1,n1,w,n1)
|
|---|
| 207 | call mxm(w,n2,ixtmc3,n1,u2,n2)
|
|---|
| 208 |
|
|---|
| 209 | return
|
|---|
| 210 | end
|
|---|
| 211 | c-----------------------------------------------------------------------
|
|---|
| 212 | subroutine int_fine_to_coarse_3d(u1,u2,n1,n2)
|
|---|
| 213 | include 'SIZE'
|
|---|
| 214 | c Interpolate u1 to u2 using interpolation matrices in12 and in12^T
|
|---|
| 215 | c u1 is size n1xn1 and u2 is size n2xn2
|
|---|
| 216 | real ixtmc3(lx1,3),ixmc3(3,lx1)
|
|---|
| 217 | real ixtmf3(3,lx1),ixmf3(lx1,3)
|
|---|
| 218 | common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
|
|---|
| 219 | real u1(n1,n1,n1),u2(n2,n2,n2)
|
|---|
| 220 | real w(20*20*20)
|
|---|
| 221 | real v(20*20*20)
|
|---|
| 222 |
|
|---|
| 223 | if (lx1.gt.20) write(6,*) 'SMOOTHER-increase size of work array
|
|---|
| 224 | $ in int_fine_to_coarse and related routines '
|
|---|
| 225 | if (lx1.gt.20) call exitt
|
|---|
| 226 |
|
|---|
| 227 | mm = n1*n1
|
|---|
| 228 | mn = n1*n2
|
|---|
| 229 | nn = n2*n2
|
|---|
| 230 |
|
|---|
| 231 | call mxm(ixmc3,n2,u1,n1,v ,mm)
|
|---|
| 232 | iv=1
|
|---|
| 233 | iw=1
|
|---|
| 234 | do k=1,n1
|
|---|
| 235 | call mxm(v(iv),n2,ixtmc3,n1,w(iw),n2)
|
|---|
| 236 | iv = iv+mn
|
|---|
| 237 | iw = iw+nn
|
|---|
| 238 | enddo
|
|---|
| 239 | call mxm(w,nn,ixtmc3,n1,u2,n2)
|
|---|
| 240 |
|
|---|
| 241 | return
|
|---|
| 242 | end
|
|---|
| 243 | c-----------------------------------------------------------------------
|
|---|
| 244 | subroutine int_coarse_to_fine_2d(u1,u2,n1,n2)
|
|---|
| 245 | include 'SIZE'
|
|---|
| 246 | c Interpolate u1 to u2 using interpolation matrices in12 and in12^T
|
|---|
| 247 | c u1 is size n1xn1 and u2 is size n2xn2
|
|---|
| 248 | c u1 is fine, u2 is coarse
|
|---|
| 249 | real ixtmc3(lx1,3),ixmc3(3,lx1)
|
|---|
| 250 | real ixtmf3(3,lx1),ixmf3(lx1,3)
|
|---|
| 251 | common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
|
|---|
| 252 | real zgmc(3),wxmc(3) !points and weights
|
|---|
| 253 | real zgmf(lx1),wxmf(lx1) !points and weights
|
|---|
| 254 | common /coarswz/ zgmc,wxmc,zgmf,wxmf
|
|---|
| 255 | real u1(n1,n1),u2(n2,n2)
|
|---|
| 256 | real w(20,20)
|
|---|
| 257 |
|
|---|
| 258 | if (lx1.gt.20) write(6,*) 'SMOOTHER-increase size of work array
|
|---|
| 259 | $ in int_coarse_to_fine and related routines '
|
|---|
| 260 | if (lx1.gt.20) call exitt
|
|---|
| 261 |
|
|---|
| 262 | call mxm(ixmf3,n1,u2,n2,w,n2)
|
|---|
| 263 | call mxm(w,n1,ixtmf3,n2,u1,n1)
|
|---|
| 264 |
|
|---|
| 265 | return
|
|---|
| 266 | end
|
|---|
| 267 | c-----------------------------------------------------------------------
|
|---|
| 268 | subroutine int_coarse_to_fine_3d(u1,u2,n1,n2)
|
|---|
| 269 | include 'SIZE'
|
|---|
| 270 | c Interpolate u1 to u2 using interpolation matrices in12 and in12^T
|
|---|
| 271 | c u1 is size n1xn1 and u2 is size n2xn2
|
|---|
| 272 | real ixtmc3(lx1,3),ixmc3(3,lx1)
|
|---|
| 273 | real ixtmf3(3,lx1),ixmf3(lx1,3)
|
|---|
| 274 | common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
|
|---|
| 275 | real u1(n1,n1,n1),u2(n2,n2,n2)
|
|---|
| 276 | real w(20*20*20)
|
|---|
| 277 | real v(20*20*20)
|
|---|
| 278 |
|
|---|
| 279 | if (lx1.gt.20) write(6,*) 'SMOOTHER-increase size of work array
|
|---|
| 280 | $ in int_coarse_to_fine and related routines '
|
|---|
| 281 | if (lx1.gt.20) call exitt
|
|---|
| 282 |
|
|---|
| 283 | mm = n2*n2
|
|---|
| 284 | mn = n2*n1
|
|---|
| 285 | nn = n1*n1
|
|---|
| 286 |
|
|---|
| 287 | call mxm(ixmf3,n1,u2,n2,v ,mm)
|
|---|
| 288 | iv=1
|
|---|
| 289 | iw=1
|
|---|
| 290 | do k=1,n2
|
|---|
| 291 | call mxm(v(iv),n1,ixtmf3,n2,w(iw),n1)
|
|---|
| 292 | iv = iv+mn
|
|---|
| 293 | iw = iw+nn
|
|---|
| 294 | enddo
|
|---|
| 295 | call mxm(w,nn,ixtmf3,n2,u1,n1)
|
|---|
| 296 |
|
|---|
| 297 | return
|
|---|
| 298 | end
|
|---|
| 299 | c-----------------------------------------------------------------------
|
|---|
| 300 | subroutine genbackupmesh !xm1->xbp backup the mesh
|
|---|
| 301 | include 'SIZE'
|
|---|
| 302 | include 'TOTAL'
|
|---|
| 303 | ! backsup the mesh
|
|---|
| 304 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 305 | real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
|
|---|
| 306 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 307 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 308 |
|
|---|
| 309 | parameter(lt = lxc*lyc*lzc*lelv)
|
|---|
| 310 | real dx(lt),dy(lt),dz(lt),xbp(lt),ybp(lt),zbp(lt)
|
|---|
| 311 | common / msmbackup / dx,dy,dz,xbp,ybp,zbp
|
|---|
| 312 |
|
|---|
| 313 | call copy(xbp,xmc,lxc*lyc*lzc*nelt)
|
|---|
| 314 | call copy(ybp,ymc,lxc*lyc*lzc*nelt)
|
|---|
| 315 | if (ldim.eq.3) call copy(zbp,zmc,lxc*lyc*lzc*nelt)
|
|---|
| 316 |
|
|---|
| 317 | return
|
|---|
| 318 | end
|
|---|
| 319 | c-----------------------------------------------------------------------
|
|---|
| 320 | subroutine restbackmesh !xbp->xm1 and xbp->x8
|
|---|
| 321 | include 'SIZE'
|
|---|
| 322 | include 'TOTAL'
|
|---|
| 323 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 324 | real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
|
|---|
| 325 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 326 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 327 |
|
|---|
| 328 | real x8(2,2,ldim-1,lelv*(2**ldim)),y8(2,2,ldim-1,lelv*(2**ldim))
|
|---|
| 329 | real z8(2,2,ldim-1,lelv*(2**ldim))
|
|---|
| 330 | common /coarsemesh8/ x8,y8,z8
|
|---|
| 331 |
|
|---|
| 332 | parameter(lt = lxc*lyc*lzc*lelv)
|
|---|
| 333 | real dx(lt),dy(lt),dz(lt),xbp(lt),ybp(lt),zbp(lt)
|
|---|
| 334 | common / msmbackup / dx,dy,dz,xbp,ybp,zbp
|
|---|
| 335 |
|
|---|
| 336 | call copy(xmc,xbp,lxc*lyc*lzc*nelt)
|
|---|
| 337 | call copy(ymc,ybp,lxc*lyc*lzc*nelt)
|
|---|
| 338 | if (ldim.eq.3) call copy(zmc,zbp,lxc*lyc*lzc*nelt)
|
|---|
| 339 | call xmtox8(xmc,x8)
|
|---|
| 340 | call xmtox8(ymc,y8)
|
|---|
| 341 | if (ldim.eq.3) call xmtox8(zmc,z8)
|
|---|
| 342 |
|
|---|
| 343 | return
|
|---|
| 344 | end
|
|---|
| 345 | c-----------------------------------------------------------------------
|
|---|
| 346 | subroutine optCG(itmax,nodmask,mlt,gshl,dis,opt,optinv,mltc,gshlc)
|
|---|
| 347 | include 'SIZE'
|
|---|
| 348 | include 'TOTAL'
|
|---|
| 349 |
|
|---|
| 350 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 351 | parameter (nxyzc=lxc*lyc*lzc,nxyz=lx1*ly1*lz1)
|
|---|
| 352 |
|
|---|
| 353 | real xmc(lxc,lyc,lzc,lelv),
|
|---|
| 354 | $ ymc(lxc,lyc,lzc,lelv),
|
|---|
| 355 | $ zmc(lxc,lyc,lzc,lelv),
|
|---|
| 356 | $ mltc(lxc,lyc,lzc,lelv)
|
|---|
| 357 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 358 |
|
|---|
| 359 | real x8(2**ldim,lelv*(2**ldim)),
|
|---|
| 360 | $ y8(2**ldim,lelv*(2**ldim)),
|
|---|
| 361 | $ z8(2**ldim,lelv*(2**ldim)),
|
|---|
| 362 | $ nodmask(2**ldim,lelv*(2**ldim)),
|
|---|
| 363 | $ dis(2**ldim,lelv*(2**ldim)),
|
|---|
| 364 | $ mlt(2**ldim,lelv*(2**ldim)),
|
|---|
| 365 | $ dx(2**ldim,lelv*(2**ldim)),
|
|---|
| 366 | $ dy(2**ldim,lelv*(2**ldim)),
|
|---|
| 367 | $ dz(2**ldim,lelv*(2**ldim)),
|
|---|
| 368 | $ dfdx(2**ldim,lelv*(2**ldim),ldim),
|
|---|
| 369 | $ dfdxn(2**ldim,lelv*(2**ldim),ldim),
|
|---|
| 370 | $ sk(2**ldim,lelv*(2**ldim),ldim)
|
|---|
| 371 | common /coarsemesh8/ x8,y8,z8
|
|---|
| 372 |
|
|---|
| 373 | integer iter,gshl,siz,eg,opt,optinv,kerr,gshlc
|
|---|
| 374 | real f2,lambda,num,den,scale2
|
|---|
| 375 | real alpha,bk,num1,den1,wrk1
|
|---|
| 376 | real hval(2**ldim,lelv*(2**ldim))
|
|---|
| 377 |
|
|---|
| 378 | optinv = 0 !initialize to 0
|
|---|
| 379 |
|
|---|
| 380 | if (nid.eq.0.and.loglevel.ge.2)
|
|---|
| 381 | $ write(6,*) 'Optimization loop started'
|
|---|
| 382 |
|
|---|
| 383 | if (opt.eq.1) siz = 2**ldim !for jacobian
|
|---|
| 384 | if (opt.eq.2) siz = 4+8*(ldim-2) !for len
|
|---|
| 385 | n1 = 2**ldim
|
|---|
| 386 | n2 = nelv*(2**ldim)*n1
|
|---|
| 387 |
|
|---|
| 388 | call get_nodscale(hval,x8,y8,z8,gshl)
|
|---|
| 389 | c if (nid.eq.0.and.loglevel.ge.4)
|
|---|
| 390 | c $ write(6,'(A,1p1e13.5)') 'perturbation amount is ',hval
|
|---|
| 391 |
|
|---|
| 392 | call gradf(f2,dfdx,x8,y8,z8,mlt,gshl,siz,opt,hval)
|
|---|
| 393 | do i=1,ldim
|
|---|
| 394 | call copy(sk(1,1,i),dfdx(1,1,i),n2)
|
|---|
| 395 | call cmult(sk(1,1,i),-1.,n2)
|
|---|
| 396 | enddo
|
|---|
| 397 |
|
|---|
| 398 | do iter=1,itmax
|
|---|
| 399 | c do line search at this point
|
|---|
| 400 | call
|
|---|
| 401 | $ dolsalpha(x8,y8,z8,sk,alpha,siz,opt,mlt,gshl,nodmask,iter)
|
|---|
| 402 |
|
|---|
| 403 | if (alpha.eq.0) goto 5002
|
|---|
| 404 |
|
|---|
| 405 | call col4(dx,sk(1,1,1),nodmask,dis,n2)
|
|---|
| 406 | call col4(dy,sk(1,1,2),nodmask,dis,n2)
|
|---|
| 407 | if (ldim.eq.3) call col4(dz,sk(1,1,ldim),nodmask,dis,n2)
|
|---|
| 408 |
|
|---|
| 409 | call add2s2(x8,dx,alpha,n2)
|
|---|
| 410 | call add2s2(y8,dy,alpha,n2)
|
|---|
| 411 | if (ldim.eq.3) call add2s2(z8,dz,alpha,n2)
|
|---|
| 412 |
|
|---|
| 413 | c Gradf at new positions
|
|---|
| 414 | call gradf(f2,dfdxn,x8,y8,z8,mlt,gshl,siz,opt,hval)
|
|---|
| 415 |
|
|---|
| 416 | c get Bk = dfdxn'*dfdxn / (dfdx'*dfdx)
|
|---|
| 417 | num1 = 0.
|
|---|
| 418 | den1 = 0.
|
|---|
| 419 | do k=1,ldim
|
|---|
| 420 | do j=1,nelv*(2**ldim)
|
|---|
| 421 | do i=1,2**ldim
|
|---|
| 422 | num1 = dfdxn(i,j,k)*mlt(i,j)*dfdxn(i,j,k) + num1
|
|---|
| 423 | den1 = dfdx(i,j,k)*mlt(i,j)*dfdx(i,j,k) + den1
|
|---|
| 424 | enddo
|
|---|
| 425 | enddo
|
|---|
| 426 | enddo
|
|---|
| 427 | call gop(num1,wrk1,'+ ',1)
|
|---|
| 428 | call gop(den1,wrk1,'+ ',1)
|
|---|
| 429 | bk = num1/den1
|
|---|
| 430 |
|
|---|
| 431 | do k=1,ldim
|
|---|
| 432 | do j=1,nelv*(2**ldim)
|
|---|
| 433 | do i=1,2**ldim
|
|---|
| 434 | sk(i,j,k) = -dfdxn(i,j,k) + bk*sk(i,j,k)
|
|---|
| 435 | dfdx(i,j,k) = dfdxn(i,j,k)
|
|---|
| 436 | enddo
|
|---|
| 437 | enddo
|
|---|
| 438 | enddo
|
|---|
| 439 |
|
|---|
| 440 | if (mod(iter,5).eq.0.or.iter.eq.itmax) then
|
|---|
| 441 | call x8toxm(xmc,x8)
|
|---|
| 442 | call x8toxm(ymc,y8)
|
|---|
| 443 | if (ldim.eq.3) call x8toxm(zmc,z8)
|
|---|
| 444 | call fixcurs(mltc,gshlc)
|
|---|
| 445 | call xmtox8(xmc,x8)
|
|---|
| 446 | call xmtox8(ymc,y8)
|
|---|
| 447 | if (ldim.eq.3) call xmtox8(zmc,z8)
|
|---|
| 448 | endif
|
|---|
| 449 | enddo
|
|---|
| 450 |
|
|---|
| 451 | 5002 continue
|
|---|
| 452 | if (alpha.eq.0) then
|
|---|
| 453 | optinv = 2
|
|---|
| 454 | call x8toxm(xmc,x8)
|
|---|
| 455 | call x8toxm(ymc,y8)
|
|---|
| 456 | if (ldim.eq.3) call x8toxm(zmc,z8)
|
|---|
| 457 | call fixcurs(mltc,gshlc)
|
|---|
| 458 | call xmtox8(xmc,x8)
|
|---|
| 459 | call xmtox8(ymc,y8)
|
|---|
| 460 | if (ldim.eq.3) call xmtox8(zmc,z8)
|
|---|
| 461 | endif
|
|---|
| 462 |
|
|---|
| 463 | call glmapm1chkinv(kerr)
|
|---|
| 464 | if (kerr.gt.0) then
|
|---|
| 465 | optinv = 1 !Terminate smoother
|
|---|
| 466 | if (nid.eq.0) write(6,*) 'Terminating optimizer'
|
|---|
| 467 | else
|
|---|
| 468 | call genbackupmesh
|
|---|
| 469 | endif
|
|---|
| 470 |
|
|---|
| 471 | return
|
|---|
| 472 | end
|
|---|
| 473 | c-----------------------------------------------------------------------
|
|---|
| 474 | subroutine
|
|---|
| 475 | $ dolsalpha(xe,ye,ze,sk,alpha,siz,opt,mlt,gshl,nodmask,iter)
|
|---|
| 476 | include 'SIZE'
|
|---|
| 477 | include 'TOTAL'
|
|---|
| 478 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 479 | real xmc(lxc,lyc,lzc,lelv),
|
|---|
| 480 | $ ymc(lxc,lyc,lzc,lelv),
|
|---|
| 481 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 482 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 483 | real xe(2**ldim,lelv*(2**ldim)),
|
|---|
| 484 | $ ye(2**ldim,lelv*(2**ldim)),
|
|---|
| 485 | $ ze(2**ldim,lelv*(2**ldim)),
|
|---|
| 486 | $ dx(2**ldim,lelv*(2**ldim)),
|
|---|
| 487 | $ dy(2**ldim,lelv*(2**ldim)),
|
|---|
| 488 | $ dz(2**ldim,lelv*(2**ldim)),
|
|---|
| 489 | $ xs(2**ldim,lelv*(2**ldim)),
|
|---|
| 490 | $ ys(2**ldim,lelv*(2**ldim)),
|
|---|
| 491 | $ zs(2**ldim,lelv*(2**ldim)),
|
|---|
| 492 | $ nodmask(2**ldim,lelv*(2**ldim)),
|
|---|
| 493 | $ mlt(2**ldim,lelv*(2**ldim)),
|
|---|
| 494 | $ sk((2**ldim),lelv*(2**ldim),ldim),
|
|---|
| 495 | $ jac(2**ldim,lelv*(2**ldim))
|
|---|
| 496 | real alpha,wrk1,f1,f2
|
|---|
| 497 | integer siz,opt,z,gshl,invflag
|
|---|
| 498 | integer n1,n2,chk1,chk2,tstop,maxcount,iter
|
|---|
| 499 |
|
|---|
| 500 | common /lsinfo / lssteps
|
|---|
| 501 | real lssteps(5)
|
|---|
| 502 | integer icalld
|
|---|
| 503 | save icalld
|
|---|
| 504 | data icalld /0/
|
|---|
| 505 | integer idx
|
|---|
| 506 | save idx
|
|---|
| 507 | data idx /0/
|
|---|
| 508 |
|
|---|
| 509 | if (icalld.eq.0) then
|
|---|
| 510 | call rzero(lssteps,5)
|
|---|
| 511 | icalld = 1
|
|---|
| 512 | idx = 0
|
|---|
| 513 | endif
|
|---|
| 514 |
|
|---|
| 515 | if (idx.lt.6) then
|
|---|
| 516 | alpha = 1.0
|
|---|
| 517 | else
|
|---|
| 518 | dumsum = vlsum(lssteps,5)
|
|---|
| 519 | alpha = 2.*dumsum
|
|---|
| 520 | endif
|
|---|
| 521 | alphasav = alpha
|
|---|
| 522 | maxcount = 20
|
|---|
| 523 | chk1 = 0
|
|---|
| 524 | chk2 = 0
|
|---|
| 525 | tstop = 0
|
|---|
| 526 | n1 = 2**ldim
|
|---|
| 527 | n2 = nelv*(2**ldim)*n1
|
|---|
| 528 |
|
|---|
| 529 | call copy(xs(1,1),xe(1,1),n2)
|
|---|
| 530 | call copy(ys(1,1),ye(1,1),n2)
|
|---|
| 531 | if (ldim.eq.3) call copy(zs(1,1),ze(1,1),n2)
|
|---|
| 532 |
|
|---|
| 533 | call getglobsum(xs,ys,zs,mlt,gshl,siz,opt,f1)
|
|---|
| 534 | z = 0
|
|---|
| 535 | do while (tstop.eq.0.and.z.le.maxcount)
|
|---|
| 536 | z = z+1
|
|---|
| 537 | call col3c(dx,sk(1,1,1),nodmask,alpha,n2)
|
|---|
| 538 | call col3c(dy,sk(1,1,2),nodmask,alpha,n2)
|
|---|
| 539 | if (ldim.eq.3) call col3c(dz,sk(1,1,ldim),nodmask,alpha,n2)
|
|---|
| 540 |
|
|---|
| 541 | call add3(xs,xe,dx,n2)
|
|---|
| 542 | call add3(ys,ye,dy,n2)
|
|---|
| 543 | if (ldim.eq.3) call add3(zs,ze,dz,n2)
|
|---|
| 544 |
|
|---|
| 545 | call getglobsum(xs,ys,zs,mlt,gshl,siz,opt,f2)
|
|---|
| 546 |
|
|---|
| 547 | if (f2.lt.f1) then
|
|---|
| 548 | chk1 = 1
|
|---|
| 549 | endif
|
|---|
| 550 | call gop(chk1,wrk1,'m ',1)
|
|---|
| 551 |
|
|---|
| 552 | chk2 = 0
|
|---|
| 553 | if (chk1.eq.1) then
|
|---|
| 554 | call x8toxm(xmc,xs)
|
|---|
| 555 | call x8toxm(ymc,ys)
|
|---|
| 556 | if (ldim.eq.3) call x8toxm(zmc,zs)
|
|---|
| 557 | call glmapm1chkinv(invflag)
|
|---|
| 558 | if (invflag.eq.0) chk2 = 1
|
|---|
| 559 | endif
|
|---|
| 560 | call gop(chk2,wrk1,'m ',1)
|
|---|
| 561 |
|
|---|
| 562 | if (chk1.eq.1.and.chk2.eq.1.) then
|
|---|
| 563 | alphasav = alpha
|
|---|
| 564 | tstop = 1.
|
|---|
| 565 | else
|
|---|
| 566 | chk1 = 0.
|
|---|
| 567 | chk2 = 0.
|
|---|
| 568 | alpha = alpha/(10.)
|
|---|
| 569 | endif
|
|---|
| 570 | enddo
|
|---|
| 571 |
|
|---|
| 572 | if (idx.lt.6) then
|
|---|
| 573 | idx = idx+1
|
|---|
| 574 | lssteps(idx) = alphasav
|
|---|
| 575 | else
|
|---|
| 576 | lssteps(1) = lssteps(2)
|
|---|
| 577 | lssteps(2) = lssteps(3)
|
|---|
| 578 | lssteps(3) = lssteps(4)
|
|---|
| 579 | lssteps(4) = lssteps(5)
|
|---|
| 580 | lssteps(5) = alphasav
|
|---|
| 581 | endif
|
|---|
| 582 |
|
|---|
| 583 | if (tstop.eq.1) then
|
|---|
| 584 | alpha = alphasav
|
|---|
| 585 | if (nid.eq.0.and.loglevel.ge.3) write(6,101) iter,f2
|
|---|
| 586 | 101 format(i5,' glob_phi ',1p1e13.6)
|
|---|
| 587 | c write(6,*) z,'number of ls steps'
|
|---|
| 588 | else
|
|---|
| 589 | alpha = 0.
|
|---|
| 590 | if (nid.eq.0.and.loglevel.ge.4)
|
|---|
| 591 | $ write(6,*) 'SMOOTHER-line-search alpha set to 0'
|
|---|
| 592 | endif
|
|---|
| 593 |
|
|---|
| 594 | return
|
|---|
| 595 | end
|
|---|
| 596 | c-----------------------------------------------------------------------
|
|---|
| 597 | subroutine getglobsum(xe,ye,ze,mlt,gshl,siz,opt,fl)
|
|---|
| 598 | include 'SIZE'
|
|---|
| 599 | include 'TOTAL'
|
|---|
| 600 | real xe(2**ldim,lelv*(2**ldim))
|
|---|
| 601 | real ye(2**ldim,lelv*(2**ldim))
|
|---|
| 602 | real ze(2**ldim,lelv*(2**ldim))
|
|---|
| 603 | real mlt(2**ldim,lelv*(2**ldim))
|
|---|
| 604 | integer gshl,siz,e,opt
|
|---|
| 605 | real f1,fl,wrk1
|
|---|
| 606 |
|
|---|
| 607 | fl = 0.
|
|---|
| 608 | do e=1,nelv*(2**ldim)
|
|---|
| 609 | if (opt.eq.1)
|
|---|
| 610 | $ call get_jac(f1,xe(1,e),ye(1,e),ze(1,e),siz,e,1)
|
|---|
| 611 | if (opt.eq.2) call get_len(f1,xe(1,e),ye(1,e),ze(1,e),siz)
|
|---|
| 612 | fl = fl+f1
|
|---|
| 613 | enddo
|
|---|
| 614 | call gop(fl,wrk1,'+ ',1)
|
|---|
| 615 |
|
|---|
| 616 | return
|
|---|
| 617 | end
|
|---|
| 618 | c-----------------------------------------------------------------------
|
|---|
| 619 | subroutine disfun(dis,funtype,alpha,beta,dcbc,nbc)
|
|---|
| 620 | include 'SIZE'
|
|---|
| 621 | include 'TOTAL'
|
|---|
| 622 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 623 | real dd1(lx1*ly1*lz1*lelv),
|
|---|
| 624 | $ dd2(lx1*ly1*lz1*lelv),
|
|---|
| 625 | $ ddd(lx1*ly1*lz1*lelv),
|
|---|
| 626 | $ ddc(lxc*lyc*lzc,lelv),
|
|---|
| 627 | $ dis((2**ldim)*lelv*(2**ldim))
|
|---|
| 628 | integer i,nbc,j,funtype
|
|---|
| 629 | character*3 dcbc(nbc)
|
|---|
| 630 | real dscale,dmax,alpha,beta,dum2
|
|---|
| 631 |
|
|---|
| 632 | if (nid.eq.0.and.loglevel.ge.5)
|
|---|
| 633 | $ write(6,*) 'calculate distance function'
|
|---|
| 634 |
|
|---|
| 635 | call rone (dd1,lx1*ly1*lz1*nelv)
|
|---|
| 636 | call sm_cheap_dist(dd1,1,dcbc(1)) ! Distance function scaling
|
|---|
| 637 |
|
|---|
| 638 | if (nbc.ge.2) then
|
|---|
| 639 | do i=2,nbc
|
|---|
| 640 | call rone (dd2,lx1*ly1*lz1*nelv)
|
|---|
| 641 | call sm_cheap_dist(dd2,1,dcbc(i)) ! Distance function scaling
|
|---|
| 642 | do j=1,lx1*ly1*lz1*nelv
|
|---|
| 643 | dd1(j) = min(dd1(j),dd2(j))
|
|---|
| 644 | enddo
|
|---|
| 645 | enddo
|
|---|
| 646 | endif
|
|---|
| 647 |
|
|---|
| 648 | dmax = glamax(dd1,lx1*ly1*lz1*nelv)
|
|---|
| 649 | dscale = 1./dmax
|
|---|
| 650 | call cmult(dd1,dscale,lx1*ly1*lz1*nelv) !normalized 0 to 1
|
|---|
| 651 | c call outpost(dd1,vy,vz,pr,t,' ')
|
|---|
| 652 |
|
|---|
| 653 | nxyz = lx1*ly1*lz1
|
|---|
| 654 | if (ldim.eq.2) then
|
|---|
| 655 | do i=1,nelv
|
|---|
| 656 | call int_fine_to_coarse_2d(dd1((i-1)*nxyz+1),ddc(1,i),lx1,lxc)
|
|---|
| 657 | enddo
|
|---|
| 658 | else
|
|---|
| 659 | do i=1,nelv
|
|---|
| 660 | call int_fine_to_coarse_3d(dd1((i-1)*nxyz+1),ddc(1,i),lx1,lxc)
|
|---|
| 661 | enddo
|
|---|
| 662 | endif
|
|---|
| 663 | call xmtox8(ddc,dis)
|
|---|
| 664 | c
|
|---|
| 665 | if (funtype.eq.0) then
|
|---|
| 666 | do i=1,(2**ldim)*nelv*(2**ldim)
|
|---|
| 667 | dis(i) = (1-EXP(-dis(i)/beta))
|
|---|
| 668 | enddo
|
|---|
| 669 | elseif (funtype.eq.1) then
|
|---|
| 670 | do i=1,(2**ldim)*nelv*(2**ldim)
|
|---|
| 671 | dis(i) = 0.5*(tanh(alpha*(dis(i)-beta))+1)
|
|---|
| 672 | enddo
|
|---|
| 673 | else
|
|---|
| 674 | call exitti('Please set the funtype to 0 or 1$',funtype)
|
|---|
| 675 | endif
|
|---|
| 676 |
|
|---|
| 677 | if (nid.eq.0.and.loglevel.ge.5)
|
|---|
| 678 | $ write(6,'(1p1e13.4,A)') dmax,'max disfun'
|
|---|
| 679 |
|
|---|
| 680 | return
|
|---|
| 681 | end
|
|---|
| 682 | c-----------------------------------------------------------------------
|
|---|
| 683 | subroutine restbndrlayer(dx,dy,dz,dis,mltc,gshlc)
|
|---|
| 684 | include 'SIZE'
|
|---|
| 685 | include 'TOTAL'
|
|---|
| 686 | c dis*smoothmesh + (1-dis)*original mesh
|
|---|
| 687 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 688 | real xmc(lxc,lyc,lzc,lelv),
|
|---|
| 689 | $ ymc(lxc,lyc,lzc,lelv),
|
|---|
| 690 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 691 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 692 |
|
|---|
| 693 | real dx(lxc*lyc*lzc,lelv),
|
|---|
| 694 | $ dy(lxc*lyc*lzc,lelv),
|
|---|
| 695 | $ dz(lxc*lyc*lzc,lelv),
|
|---|
| 696 | $ dis((2**ldim)*lelv*(2**ldim)),
|
|---|
| 697 | $ dis2(lxc,lyc,lzc,lelv),
|
|---|
| 698 | $ mltc(lxc,lyc,lzc,lelv)
|
|---|
| 699 |
|
|---|
| 700 | integer gshlc
|
|---|
| 701 |
|
|---|
| 702 | call x8toxm(dis2,dis)
|
|---|
| 703 | dum2 = glamax(dis2,lxc*lyc*lzc*nelv)
|
|---|
| 704 | dum2 = 1./dum2
|
|---|
| 705 | call cmult(dis2,dum2,lxc*lyc*lzc*nelv)
|
|---|
| 706 |
|
|---|
| 707 | call col2(xmc,dis2,lxc*lyc*lzc*nelv) !w*xs
|
|---|
| 708 | call col2(ymc,dis2,lxc*lyc*lzc*nelv)
|
|---|
| 709 | if (ndim.eq.3) call col2(zmc,dis2,lxc*lyc*lzc*nelv)
|
|---|
| 710 |
|
|---|
| 711 | call add2(xmc,dx,lxc*lyc*lzc*nelv) !w*xs+xo
|
|---|
| 712 | call add2(ymc,dy,lxc*lyc*lzc*nelv)
|
|---|
| 713 | if (ndim.eq.3) call add2(zmc,dz,lxc*lyc*lzc*nelv)
|
|---|
| 714 |
|
|---|
| 715 | call col2(dx,dis2,lxc*lyc*lzc*nelv) !w*xo
|
|---|
| 716 | call col2(dy,dis2,lxc*lyc*lzc*nelv)
|
|---|
| 717 | if (ndim.eq.3) call col2(dz,dis2,lxc*lyc*lzc*nelv)
|
|---|
| 718 |
|
|---|
| 719 | call sub2(xmc,dx,lxc*lyc*lzc*nelv) !w*xs+xo-w*xo=w*xs+xo(1-w)
|
|---|
| 720 | call sub2(ymc,dy,lxc*lyc*lzc*nelv) !this is now store in xmc..zmc
|
|---|
| 721 | if (ndim.eq.3) call sub2(zmc,dz,lxc*lyc*lzc*nelv)
|
|---|
| 722 |
|
|---|
| 723 | call sub2(dx,xmc,lxc*lyc*lzc*nelv) !dx=xo-xmc
|
|---|
| 724 | call sub2(dy,ymc,lxc*lyc*lzc*nelv)
|
|---|
| 725 | if (ndim.eq.3) call sub2(dz,zmc,lxc*lyc*lzc*nelv)
|
|---|
| 726 |
|
|---|
| 727 | call fixcurs(mltc,gshlc)
|
|---|
| 728 | call xmctoxm1
|
|---|
| 729 |
|
|---|
| 730 | return
|
|---|
| 731 | end
|
|---|
| 732 | c-----------------------------------------------------------------------
|
|---|
| 733 | subroutine masklayers(nodmask,nlayers)
|
|---|
| 734 | include 'SIZE'
|
|---|
| 735 | include 'TOTAL'
|
|---|
| 736 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 737 | real nodmask(2**ldim,lelv*(2**ldim))
|
|---|
| 738 | real d(lx1*ly1*lz1,lelv)
|
|---|
| 739 | real vcmask(lxc,lyc,lzc,lelv)
|
|---|
| 740 | integer e,f,i,j,nlayers,valm,c
|
|---|
| 741 |
|
|---|
| 742 | character*3 tcbc(5)
|
|---|
| 743 | save tcbc
|
|---|
| 744 | data tcbc /'W ','v ','V ','mv ','MV ' /
|
|---|
| 745 |
|
|---|
| 746 | call x8toxm(v1mask,nodmask)
|
|---|
| 747 |
|
|---|
| 748 | call rone(d,nx1*ny1*nz1*nelv)
|
|---|
| 749 | do e=1,nelv
|
|---|
| 750 | do f=1,2*ldim
|
|---|
| 751 | do c=1,5
|
|---|
| 752 | if(cbc(f,e,1).eq.tcbc(c)) call facev(d,e,f,zero,nx1,ny1,nz1)
|
|---|
| 753 | enddo
|
|---|
| 754 | enddo
|
|---|
| 755 | enddo
|
|---|
| 756 | call dsop(d,'mul',nx1,ny1,nz1)
|
|---|
| 757 |
|
|---|
| 758 | do i=1,nlayers
|
|---|
| 759 | call dsop(d,'mul',nx1,ny1,nz1)
|
|---|
| 760 | do e=1,nelv
|
|---|
| 761 | val = glamin(d(1,e),lx1**ldim)
|
|---|
| 762 | call cmult(d(1,e),val,lx1**ldim)
|
|---|
| 763 | enddo
|
|---|
| 764 | enddo
|
|---|
| 765 |
|
|---|
| 766 | call dsop(d,'mul',nx1,ny1,nz1)
|
|---|
| 767 | call col2(v1mask,d,lx1*ly1*lz1*nelv)
|
|---|
| 768 | do i=1,nelv
|
|---|
| 769 | call int_fine_to_coarse_2d(v1mask(1,1,1,e),vcmask(1,1,1,e),lx1,3)
|
|---|
| 770 | enddo
|
|---|
| 771 | call xmtox8(vcmask,nodmask)
|
|---|
| 772 |
|
|---|
| 773 | return
|
|---|
| 774 | end
|
|---|
| 775 | c-----------------------------------------------------------------------
|
|---|
| 776 | subroutine fixcurs(mltc,gshlc)
|
|---|
| 777 | include 'SIZE'
|
|---|
| 778 | include 'TOTAL'
|
|---|
| 779 | c This routine fixes the curved edges post smoothing
|
|---|
| 780 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 781 | real xmc(lxc,lyc,lzc,lelv),
|
|---|
| 782 | $ ymc(lxc,lyc,lzc,lelv),
|
|---|
| 783 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 784 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 785 |
|
|---|
| 786 | common /ctmp0/ zg(3)
|
|---|
| 787 | real w1(lxc*lyc*lzc*lelv),w2(lxc*lyc*lzc*lelv)
|
|---|
| 788 | real vcmask(lxc,lyc,lzc,lelv)
|
|---|
| 789 | integer f,e,nedge,n,nfaces,gshlc
|
|---|
| 790 | integer ke(3,12)
|
|---|
| 791 | real xyz(3,3)
|
|---|
| 792 | real xd(lxc*lyc*lzc),
|
|---|
| 793 | $ yd(lxc*lyc*lzc),
|
|---|
| 794 | $ zd(lxc*lyc*lzc),
|
|---|
| 795 | $ mltc(lxc*lyc*lzc*lelv)
|
|---|
| 796 | real nxv,nyv,nzv
|
|---|
| 797 |
|
|---|
| 798 | save ke
|
|---|
| 799 | data ke / 1, 2, 3, 3, 6, 9, 9, 8, 7, 7, 4, 1
|
|---|
| 800 | $ , 19,20,21, 21,24,27, 27,26,25, 25,22,19
|
|---|
| 801 | $ , 1,10,19, 3,12,21, 9,18,27, 7,16,25 /
|
|---|
| 802 | integer kin(7)
|
|---|
| 803 |
|
|---|
| 804 | c START BY LOOPING OVER EACH ELEMENT AND THEN OVER EACH EDGE
|
|---|
| 805 | nedge = 4 + 8*(ldim-2)
|
|---|
| 806 | nfaces = 2*ldim
|
|---|
| 807 |
|
|---|
| 808 | n = lxc*lyc*lzc*nelv
|
|---|
| 809 | call rone (vcmask,n)
|
|---|
| 810 | do e=1,nelv
|
|---|
| 811 | do f=1,nfaces
|
|---|
| 812 | if (cbc(f,e,1).ne.'E '.and.cbc(f,e,1).ne.' ')
|
|---|
| 813 | $ call facev (vcmask,e,f,0.0,lxc,lyc,lzc)
|
|---|
| 814 | enddo
|
|---|
| 815 | enddo
|
|---|
| 816 | call fgslib_gs_op(gshlc,vcmask,1,2,0)
|
|---|
| 817 |
|
|---|
| 818 | do e=1,nelv
|
|---|
| 819 | do j=1,nedge
|
|---|
| 820 | call rzero(xyz,9)
|
|---|
| 821 | do i=1,3
|
|---|
| 822 | xyz(1,i) = xmc(ke(i,j),1,1,e)
|
|---|
| 823 | xyz(2,i) = ymc(ke(i,j),1,1,e)
|
|---|
| 824 | if (ldim.eq.3) xyz(3,i) = zmc(ke(i,j),1,1,e)
|
|---|
| 825 | enddo
|
|---|
| 826 |
|
|---|
| 827 | if (vcmask(ke(2,j),1,1,e).gt.0.5) then
|
|---|
| 828 | call fixedgs(xyz,nxv,nyv,nzv)
|
|---|
| 829 | xmc(ke(2,j),1,1,e) = nxv
|
|---|
| 830 | ymc(ke(2,j),1,1,e) = nyv
|
|---|
| 831 | if (ldim.eq.3) zmc(ke(2,j),1,1,e) = nzv
|
|---|
| 832 | endif
|
|---|
| 833 | enddo
|
|---|
| 834 |
|
|---|
| 835 | c fix face center and body center
|
|---|
| 836 | zg(1) = -1.
|
|---|
| 837 | zg(2) = 0.
|
|---|
| 838 | zg(3) = 1.
|
|---|
| 839 | call copy(xd,xmc(1,1,1,e),lxc*lyc*lzc)
|
|---|
| 840 | call copy(yd,ymc(1,1,1,e),lxc*lyc*lzc)
|
|---|
| 841 | if (ldim.eq.3) call copy(zd,zmc(1,1,1,e),lxc*lyc*lzc)
|
|---|
| 842 | call gh_face_extend(xd,zg,3,2,w1,w2)
|
|---|
| 843 | call gh_face_extend(yd,zg,3,2,w1,w2)
|
|---|
| 844 | if (ldim.eq.3) call gh_face_extend(zd,zg,3,2,w1,w2)
|
|---|
| 845 | kin(1) = 5
|
|---|
| 846 | kin(2) = 11
|
|---|
| 847 | kin(3) = 13
|
|---|
| 848 | kin(4) = 15
|
|---|
| 849 | kin(5) = 17
|
|---|
| 850 | kin(6) = 23
|
|---|
| 851 | kin(7) = 14
|
|---|
| 852 | do i=1,(ldim-2)*6+1
|
|---|
| 853 | xmc(kin(i),1,1,e) = xd(kin(i))
|
|---|
| 854 | ymc(kin(i),1,1,e) = yd(kin(i))
|
|---|
| 855 | if (ldim.eq.3) zmc(kin(i),1,1,e) = zd(kin(i))
|
|---|
| 856 | enddo
|
|---|
| 857 | enddo
|
|---|
| 858 |
|
|---|
| 859 | return
|
|---|
| 860 | end
|
|---|
| 861 | c-----------------------------------------------------------------------
|
|---|
| 862 | subroutine fixedgs(xyz,nx,ny,nz)
|
|---|
| 863 | real xyz(3,0:2) ! Coordinates
|
|---|
| 864 | real x0(3),x1(3),x2(3),v2(3),v1(3),l02,l01,n1i,u1(3),u2(3)
|
|---|
| 865 | real nx,ny,nz
|
|---|
| 866 | real h,tol,h2
|
|---|
| 867 |
|
|---|
| 868 | call copy(x0,xyz(1,0),3)
|
|---|
| 869 | call copy(x1,xyz(1,1),3)
|
|---|
| 870 | call copy(x2,xyz(1,2),3)
|
|---|
| 871 | ! ^ x2
|
|---|
| 872 | call sub3(v1,x1,x0,3) ! v1=x1-x0 / \ v2 p1 = x0+v2*dot
|
|---|
| 873 | call sub3(v2,x2,x0,3) ! v2=x2-x0 x1 .<-- x0
|
|---|
| 874 | ! v1
|
|---|
| 875 | l02 = vlsc2(v2,v2,3) ! || x2-x0 ||
|
|---|
| 876 | if (l02.gt.0) then
|
|---|
| 877 | l02 = sqrt(l02)
|
|---|
| 878 | scl = 1./l02
|
|---|
| 879 | call cmult2(u2,v2,scl,3) ! Unit tangent
|
|---|
| 880 | endif
|
|---|
| 881 |
|
|---|
| 882 | l01 = vlsc2(v1,v1,3) ! || x1-x0 ||
|
|---|
| 883 | if (l01.gt.0) then
|
|---|
| 884 | l01 = sqrt(l01)
|
|---|
| 885 | scl = 1./l01
|
|---|
| 886 | call cmult2(u1,v1,scl,3) ! Unit tangent
|
|---|
| 887 | endif
|
|---|
| 888 |
|
|---|
| 889 | dot = vlsc2(v1,u2,3)
|
|---|
| 890 |
|
|---|
| 891 | if (dot.le.0) then
|
|---|
| 892 | write(6,*) 'SMOOTHER-ERROR 1 IN SHIFT - YOU SHOULD ABORT'
|
|---|
| 893 | return
|
|---|
| 894 | elseif (dot.gt.l02) then
|
|---|
| 895 | write(6,*) 'SMOOTHER-ERROR 2 IN SHIFT - YOU SHOULD ABORT'
|
|---|
| 896 | return
|
|---|
| 897 | endif
|
|---|
| 898 | h = 0.
|
|---|
| 899 | h2 = 0.
|
|---|
| 900 | tol = 1.e-8
|
|---|
| 901 | do i=1,3
|
|---|
| 902 | p1i = x0(i) + u2(i)*dot ! Projection of x1 onto [x0,x2]
|
|---|
| 903 | n1i = x1(i) - p1i ! Normal vector
|
|---|
| 904 | h = h + n1i**2
|
|---|
| 905 | h2 = h2+ (x2(i)-x0(i))**2
|
|---|
| 906 | xmi = 0.5*(x0(i)+x2(i))
|
|---|
| 907 | x1(i) = xmi + n1i ! X1 point shifted to be centered at midpoint
|
|---|
| 908 | enddo
|
|---|
| 909 | if (h.le.h2*tol) then
|
|---|
| 910 | x1(1) = 0.5*(x0(1)+x2(1))
|
|---|
| 911 | x1(2) = 0.5*(x0(2)+x2(2))
|
|---|
| 912 | x1(3) = 0.5*(x0(3)+x2(3))
|
|---|
| 913 | endif
|
|---|
| 914 |
|
|---|
| 915 |
|
|---|
| 916 | nx = x1(1)
|
|---|
| 917 | ny = x1(2)
|
|---|
| 918 | nz = x1(3)
|
|---|
| 919 |
|
|---|
| 920 | return
|
|---|
| 921 | end
|
|---|
| 922 | c-----------------------------------------------------------------------
|
|---|
| 923 | subroutine gradf(f2,dfdx,x8,y8,z8,mlt,gshl,siz,opt,h)
|
|---|
| 924 | include 'SIZE'
|
|---|
| 925 | include 'TOTAL'
|
|---|
| 926 | real x8(2**ldim,lelv*(2**ldim)),
|
|---|
| 927 | $ y8(2**ldim,lelv*(2**ldim)),
|
|---|
| 928 | $ z8(2**ldim,lelv*(2**ldim)),
|
|---|
| 929 | $ mlt(2**ldim,lelv*(2**ldim)),
|
|---|
| 930 | $ dfdx(2**ldim,lelv*(2**ldim),ldim)
|
|---|
| 931 | real xt(2**ldim),
|
|---|
| 932 | $ yt(2**ldim),
|
|---|
| 933 | $ zt(2**ldim)
|
|---|
| 934 | integer siz,opt,vertex,gshl,e,eg,e0,f
|
|---|
| 935 | real par(siz)
|
|---|
| 936 | real f1,fl,gl,f2,h(2**ldim,lelv*(2**ldim))
|
|---|
| 937 |
|
|---|
| 938 | f1 = 0
|
|---|
| 939 | do e=1,nelv*(2**ldim)
|
|---|
| 940 | if (opt.eq.1)
|
|---|
| 941 | $ call get_jac(fl,x8(1,e),y8(1,e),z8(1,e),siz,e,0)
|
|---|
| 942 | if (opt.eq.2) call get_len(fl,x8(1,e),y8(1,e),z8(1,e),siz)
|
|---|
| 943 | f1 = f1+fl
|
|---|
| 944 | call copy(xt,x8(1,e),2**ldim)
|
|---|
| 945 | call copy(yt,y8(1,e),2**ldim)
|
|---|
| 946 | call copy(zt,z8(1,e),2**ldim)
|
|---|
| 947 | do j=1,2**ldim
|
|---|
| 948 | xt(j) = x8(j,e)+h(j,e)
|
|---|
| 949 | if (opt.eq.1) call get_jac(gl,xt,yt,zt,siz,e,1)
|
|---|
| 950 | if (opt.eq.2) call get_len(gl,xt,yt,zt,siz)
|
|---|
| 951 | xt(j) = x8(j,e)-h(j,e)
|
|---|
| 952 | if (opt.eq.1) call get_jac(fl,xt,yt,zt,siz,e,1)
|
|---|
| 953 | if (opt.eq.2) call get_len(fl,xt,yt,zt,siz)
|
|---|
| 954 | xt(j) = x8(j,e)
|
|---|
| 955 | dfdx(j,e,1) = (gl-fl)/(2.*h(j,e))
|
|---|
| 956 |
|
|---|
| 957 | yt(j) = y8(j,e)+h(j,e)
|
|---|
| 958 | if (opt.eq.1) call get_jac(gl,xt,yt,zt,siz,e,2)
|
|---|
| 959 | if (opt.eq.2) call get_len(gl,xt,yt,zt,siz)
|
|---|
| 960 | yt(j) = y8(j,e)-h(j,e)
|
|---|
| 961 | if (opt.eq.1) call get_jac(fl,xt,yt,zt,siz,e,2)
|
|---|
| 962 | if (opt.eq.2) call get_len(fl,xt,yt,zt,siz)
|
|---|
| 963 | yt(j) = y8(j,e)
|
|---|
| 964 | dfdx(j,e,2) = (gl-fl)/(2.*h(j,e))
|
|---|
| 965 |
|
|---|
| 966 | if (ldim.eq.3) then
|
|---|
| 967 | zt(j) = z8(j,e)+h(j,e)
|
|---|
| 968 | if (opt.eq.1) call get_jac(gl,xt,yt,zt,siz,e,3)
|
|---|
| 969 | if (opt.eq.2) call get_len(gl,xt,yt,zt,siz)
|
|---|
| 970 | zt(j) = z8(j,e)-h(j,e)
|
|---|
| 971 | if (opt.eq.1) call get_jac(fl,xt,yt,zt,siz,e,3)
|
|---|
| 972 | if (opt.eq.2) call get_len(fl,xt,yt,zt,siz)
|
|---|
| 973 | zt(j) = z8(j,e)
|
|---|
| 974 | dfdx(j,e,ldim) = (gl-fl)/(2.*h(j,e))
|
|---|
| 975 | endif
|
|---|
| 976 | enddo
|
|---|
| 977 | enddo
|
|---|
| 978 |
|
|---|
| 979 | call fgslib_gs_op(gshl,dfdx(1,1,1),1,1,0)
|
|---|
| 980 | call fgslib_gs_op(gshl,dfdx(1,1,2),1,1,0)
|
|---|
| 981 | if (ldim.eq.3) call fgslib_gs_op(gshl,dfdx(1,1,ldim),1,1,0)
|
|---|
| 982 |
|
|---|
| 983 | f2 = glsum(f1,1)
|
|---|
| 984 |
|
|---|
| 985 | return
|
|---|
| 986 | end
|
|---|
| 987 | c-----------------------------------------------------------------------
|
|---|
| 988 | subroutine get_jac(val,x,y,z,siz,el,dire)
|
|---|
| 989 | include 'SIZE'
|
|---|
| 990 | include 'TOTAL'
|
|---|
| 991 | c This routine compiles the jacobian matrix and then does the
|
|---|
| 992 | c frobenius norm of the matrix times that of the inverse
|
|---|
| 993 | integer i,j,k,siz,el,dire
|
|---|
| 994 |
|
|---|
| 995 | real par(siz),det(siz)
|
|---|
| 996 | real jm(ldim,ldim),jin(ldim,ldim),frn(ldim**2)
|
|---|
| 997 | real x(2**ldim),y(2**ldim),z(2**ldim),jac(2**ldim)
|
|---|
| 998 | real fr1,fr2,sum1,dumc,val
|
|---|
| 999 |
|
|---|
| 1000 | integer bzindx(24),czindx(24)
|
|---|
| 1001 | integer penalty
|
|---|
| 1002 | SAVE bzindx
|
|---|
| 1003 | DATA bzindx / 2,3,5, 1,4,6, 4,1,7, 3,2,8,
|
|---|
| 1004 | $ 6,7,1, 5,8,2, 8,5,3, 7,6,4 /
|
|---|
| 1005 | c bzindx tells which node is node connected to in r,s,t direction
|
|---|
| 1006 | c example: node 1 is connected to 2,3,5; 2 to 1,4,6 and so on
|
|---|
| 1007 | SAVE czindx
|
|---|
| 1008 | DATA czindx / 1,1,1, -1,1,1, 1,-1,1, -1,-1,1,
|
|---|
| 1009 | $ 1,1,-1, -1,1,-1, 1,-1,-1, -1,-1,-1 /
|
|---|
| 1010 | c this tells which node is to be subtracted.
|
|---|
| 1011 |
|
|---|
| 1012 | if (siz.ne.2**ldim) then
|
|---|
| 1013 | write(6,*) 'siz value wrong input into get_jac'
|
|---|
| 1014 | write(6,*) siz,'sizval',2**ldim,'should be this'
|
|---|
| 1015 | call exitt
|
|---|
| 1016 | endif
|
|---|
| 1017 |
|
|---|
| 1018 | do i=1,2**ldim
|
|---|
| 1019 | ind1 = (i-1)*3 !tells where to look in b/czindx
|
|---|
| 1020 | do j=1,ldim
|
|---|
| 1021 | ind2 = ind1+j
|
|---|
| 1022 | jm(1,j) = 0.5*czindx(ind2)*(x(bzindx(ind2))-x(i))
|
|---|
| 1023 | jm(2,j) = 0.5*czindx(ind2)*(y(bzindx(ind2))-y(i))
|
|---|
| 1024 | if (ldim.eq.3)
|
|---|
| 1025 | $ jm(ldim,j) = 0.5*czindx(ind2)*(z(bzindx(ind2))-z(i))
|
|---|
| 1026 | enddo
|
|---|
| 1027 | c jacobian matrix has been calculated at this point.
|
|---|
| 1028 | c now calculate determinant
|
|---|
| 1029 | if (ldim.eq.3) then
|
|---|
| 1030 | jac(i)=jm(1,1)*(jm(2,2)*jm(ldim,ldim)-jm(2,ldim)*jm(ldim,2))+
|
|---|
| 1031 | $ jm(1,ldim)*(jm(2,1)*jm(ldim,2)-jm(2,2)*jm(ldim,1)) +
|
|---|
| 1032 | $ jm(1,2)*(jm(2,ldim)*jm(ldim,1)-jm(2,1)*jm(ldim,ldim))
|
|---|
| 1033 | else
|
|---|
| 1034 | jac(i)=jm(1,1)*jm(2,2)-jm(2,1)*jm(1,2)
|
|---|
| 1035 | endif
|
|---|
| 1036 | c calculate inverse matrix
|
|---|
| 1037 | if (ldim.eq.3) then
|
|---|
| 1038 | jin(1,1) = jm(2,2)*jm(ldim,ldim)-jm(ldim,2)*jm(2,ldim)
|
|---|
| 1039 | jin(2,1) = jm(2,ldim)*jm(ldim,1)-jm(ldim,ldim)*jm(2,1)
|
|---|
| 1040 | jin(ldim,1) = jm(2,1)*jm(ldim,2)-jm(ldim,1)*jm(2,2)
|
|---|
| 1041 |
|
|---|
| 1042 | jin(1,2) = jm(1,ldim)*jm(ldim,2)-jm(ldim,ldim)*jm(1,2)
|
|---|
| 1043 | jin(2,2) = jm(1,1)*jm(ldim,ldim)-jm(ldim,1)*jm(1,ldim)
|
|---|
| 1044 | jin(ldim,2) = jm(1,2)*jm(ldim,1)-jm(ldim,2)*jm(1,1)
|
|---|
| 1045 |
|
|---|
| 1046 | jin(1,ldim) = jm(1,2)*jm(2,ldim)-jm(2,2)*jm(1,ldim)
|
|---|
| 1047 | jin(2,ldim) = jm(1,ldim)*jm(2,1)-jm(2,ldim)*jm(1,1)
|
|---|
| 1048 | jin(ldim,ldim) = jm(1,1)*jm(2,2)-jm(2,1)*jm(1,2)
|
|---|
| 1049 | else
|
|---|
| 1050 | jin(1,1) = jm(2,2)
|
|---|
| 1051 | jin(2,1) = -jm(2,1)
|
|---|
| 1052 |
|
|---|
| 1053 | jin(1,2) = -jm(1,2)
|
|---|
| 1054 | jin(2,2) = jm(1,1)
|
|---|
| 1055 | endif
|
|---|
| 1056 |
|
|---|
| 1057 | dumc = 1/jac(i)
|
|---|
| 1058 | call cmult(jin,dumc,ldim**2) !scale inverse by inv det
|
|---|
| 1059 |
|
|---|
| 1060 | call rzero(frn,ldim**2)
|
|---|
| 1061 | call col3(frn,jm,jm,ldim**2) !square the entries
|
|---|
| 1062 | sum1 = vlsum(frn,ldim**2)
|
|---|
| 1063 | fr1 = SQRT(sum1) !squareroot
|
|---|
| 1064 |
|
|---|
| 1065 | call rzero(frn,ldim**2)
|
|---|
| 1066 | call col3(frn,jin,jin,ldim**2)
|
|---|
| 1067 | sum1 = vlsum(frn,ldim**2)
|
|---|
| 1068 | fr2 = SQRT(sum1)
|
|---|
| 1069 |
|
|---|
| 1070 | par(i) = fr1*fr2
|
|---|
| 1071 | par(i) = (par(i)/ldim)**2
|
|---|
| 1072 |
|
|---|
| 1073 | enddo
|
|---|
| 1074 |
|
|---|
| 1075 | val = vlsum(par,siz)
|
|---|
| 1076 | val = val/siz !normalize to 1 for each element
|
|---|
| 1077 | val=abs(val)
|
|---|
| 1078 |
|
|---|
| 1079 | return
|
|---|
| 1080 | end
|
|---|
| 1081 | c-----------------------------------------------------------------------
|
|---|
| 1082 | subroutine get_len(val,x,y,z,siz)
|
|---|
| 1083 | include 'SIZE'
|
|---|
| 1084 | include 'TOTAL'
|
|---|
| 1085 | integer siz
|
|---|
| 1086 | real x(2**ldim),y(2**ldim),z(2**ldim)
|
|---|
| 1087 | real par(siz),xm,ym,zm,val
|
|---|
| 1088 |
|
|---|
| 1089 | integer azindx(24)
|
|---|
| 1090 | SAVE azindx
|
|---|
| 1091 | DATA azindx / 1,2 ,2,4, 3,4, 1,3, 2,6, 6,8,
|
|---|
| 1092 | $ 4,8, 5,6, 7,8, 5,7, 1,5, 3,7 /
|
|---|
| 1093 | c azindx tells what nodes make up each edge
|
|---|
| 1094 |
|
|---|
| 1095 | do i=1,siz
|
|---|
| 1096 | ind = (i-1)*2
|
|---|
| 1097 | xm = x(azindx(ind+1))-x(azindx(ind+2))
|
|---|
| 1098 | ym = y(azindx(ind+1))-y(azindx(ind+2))
|
|---|
| 1099 | zm = z(azindx(ind+1))-z(azindx(ind+2))
|
|---|
| 1100 | par(i) = sqrt(xm*xm+ym*ym+zm*zm)
|
|---|
| 1101 | par(i) = par(i)**2
|
|---|
| 1102 | enddo
|
|---|
| 1103 |
|
|---|
| 1104 | val = vlsum(par,siz)
|
|---|
| 1105 |
|
|---|
| 1106 | return
|
|---|
| 1107 | end
|
|---|
| 1108 | c-----------------------------------------------------------------------
|
|---|
| 1109 | subroutine get_nodscale(scalek,x8,y8,z8,gshl)
|
|---|
| 1110 | include 'SIZE'
|
|---|
| 1111 | include 'TOTAL'
|
|---|
| 1112 | real x8(2**ldim,lelv*(2**ldim)),
|
|---|
| 1113 | $ y8(2**ldim,lelv*(2**ldim)),
|
|---|
| 1114 | $ z8(2**ldim,lelv*(2**ldim))
|
|---|
| 1115 | real curval,xm,ym,zm,work1(1),dum,wrk1
|
|---|
| 1116 | c real scalek
|
|---|
| 1117 | real scalek(2**ldim,lelv*(2**ldim))
|
|---|
| 1118 | integer bzindx(24),e,gshl
|
|---|
| 1119 | DATA bzindx / 2,3,5, 1,4,6, 4,1,7, 3,2,8,
|
|---|
| 1120 | $ 6,7,1, 5,8,2, 8,5,3, 7,6,4 /
|
|---|
| 1121 | c bzindx tells what node is connected to what node
|
|---|
| 1122 |
|
|---|
| 1123 | n1 = nelv*(2**ldim)
|
|---|
| 1124 | n2 = n1*(2**ldim)
|
|---|
| 1125 | curval = 1e+10
|
|---|
| 1126 | call rone(scalek,n2)
|
|---|
| 1127 | call cmult(scalek,curval,n2)
|
|---|
| 1128 |
|
|---|
| 1129 | do e=1,nelv*(2**ldim)
|
|---|
| 1130 | do i=1,2**ldim
|
|---|
| 1131 | curval = 1e+10
|
|---|
| 1132 | ind = (i-1)*3
|
|---|
| 1133 | do j=1,ldim
|
|---|
| 1134 | xm = x8(i,e)-x8(bzindx(ind+j),e)
|
|---|
| 1135 | ym = y8(i,e)-y8(bzindx(ind+j),e)
|
|---|
| 1136 | zm = 0.
|
|---|
| 1137 | if (ldim.eq.3) zm = z8(i,e)-z8(bzindx(ind+j),e)
|
|---|
| 1138 | dum = sqrt(xm*xm+ym*ym+zm*zm)
|
|---|
| 1139 | curval = min(dum,curval)
|
|---|
| 1140 | enddo
|
|---|
| 1141 | scalek(i,e) = curval
|
|---|
| 1142 | enddo
|
|---|
| 1143 | enddo
|
|---|
| 1144 | c
|
|---|
| 1145 | fac = 1.e-2
|
|---|
| 1146 | call cmult(scalek,fac,n2)
|
|---|
| 1147 | call fgslib_gs_op(gshl,scalek,1,3,0)
|
|---|
| 1148 |
|
|---|
| 1149 | return
|
|---|
| 1150 | end
|
|---|
| 1151 | c-----------------------------------------------------------------------
|
|---|
| 1152 | subroutine genmask(nodmask,mlt,gshl,mltc,gshlc)
|
|---|
| 1153 | include 'SIZE'
|
|---|
| 1154 | include 'TOTAL'
|
|---|
| 1155 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 1156 | integer gs_handle
|
|---|
| 1157 | integer vertex(1)
|
|---|
| 1158 | common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
|
|---|
| 1159 | common /ivrtx/ vertex
|
|---|
| 1160 | integer*8 glo_num(lxc*lyc*lzc,lelv),ngv
|
|---|
| 1161 | integer*8 glo_numk2(2**ldim,lelv*(2**ldim))
|
|---|
| 1162 |
|
|---|
| 1163 | real mlt(2**ldim,lelv*(2**ldim)),
|
|---|
| 1164 | $ nodmask(2**ldim,lelv*(2**ldim))
|
|---|
| 1165 | real mltc(lxc*lyc*lzc*lelv),
|
|---|
| 1166 | $ wrkmask(lxc*lyc*lzc*lelv)
|
|---|
| 1167 | integer gshlc,gshl,e,eg,e0,f
|
|---|
| 1168 |
|
|---|
| 1169 | integer zindx(64)
|
|---|
| 1170 | SAVE zindx
|
|---|
| 1171 | DATA zindx / 1, 2, 4, 5, 10, 11, 13, 14,
|
|---|
| 1172 | $ 2, 3, 5, 6, 11, 12, 14, 15,
|
|---|
| 1173 | $ 4, 5, 7, 8, 13, 14, 16, 17,
|
|---|
| 1174 | $ 5, 6, 8, 9, 14, 15, 17, 18,
|
|---|
| 1175 | $ 10, 11, 13, 14, 19, 20, 22, 23,
|
|---|
| 1176 | $ 11, 12, 14, 15, 20, 21, 23, 24,
|
|---|
| 1177 | $ 13, 14, 16, 17, 22, 23, 25, 26,
|
|---|
| 1178 | $ 14, 15, 17, 18, 23, 24, 26, 27 /
|
|---|
| 1179 |
|
|---|
| 1180 |
|
|---|
| 1181 | call setupds_center(gs_handle,lxc,lyc,lzc,nelv,
|
|---|
| 1182 | $ nelgv,vertex,glo_num)
|
|---|
| 1183 |
|
|---|
| 1184 | c setup the mask first so that it can be distribute as well
|
|---|
| 1185 | zero = 0.
|
|---|
| 1186 | call rone(wrkmask,lxc*lyc*lzc*nelv)
|
|---|
| 1187 | do e=1,nelv
|
|---|
| 1188 | do f=1,2*ldim
|
|---|
| 1189 | if(cbc(f,e,1).ne.'E '.and.cbc(f,e,1).ne.' ')then
|
|---|
| 1190 | call facev(wrkmask,e,f,zero,lxc,lyc,lzc)
|
|---|
| 1191 | endif
|
|---|
| 1192 | enddo
|
|---|
| 1193 | enddo
|
|---|
| 1194 |
|
|---|
| 1195 | n = (lxc*lyc*lzc)*nelv
|
|---|
| 1196 | call rone(mltc,n)
|
|---|
| 1197 | call fgslib_gs_setup(gshlc,glo_num,n,nekcomm,np)
|
|---|
| 1198 | call fgslib_gs_op(gshlc,mltc,1,1,0)
|
|---|
| 1199 |
|
|---|
| 1200 | call fgslib_gs_op(gshlc,wrkmask,1,2,0)
|
|---|
| 1201 | call xmtox8(wrkmask,nodmask)
|
|---|
| 1202 |
|
|---|
| 1203 | n = 0
|
|---|
| 1204 | do e=1,nelv !each element
|
|---|
| 1205 | do j=1,2**ldim !broken down into 4 quads/8hexes
|
|---|
| 1206 | n = n+1
|
|---|
| 1207 | ind1 = (j-1)*8
|
|---|
| 1208 | do k=1,2**ldim
|
|---|
| 1209 | glo_numk2(k,n) = glo_num(zindx(ind1+k),e)
|
|---|
| 1210 | enddo
|
|---|
| 1211 | enddo
|
|---|
| 1212 | enddo
|
|---|
| 1213 |
|
|---|
| 1214 | n = (2**ldim)*(nelv*(2**ldim))
|
|---|
| 1215 | call fgslib_gs_setup(gshl,glo_numk2,n,nekcomm,np)
|
|---|
| 1216 | call rone(mlt,n)
|
|---|
| 1217 | call fgslib_gs_op(gshl,mlt,1,1,0) ! '+'
|
|---|
| 1218 | xmlt = glmax(mlt,n)
|
|---|
| 1219 | call invcol1(mlt,n)
|
|---|
| 1220 |
|
|---|
| 1221 | call fgslib_gs_op(gshl,nodmask,1,2,0)
|
|---|
| 1222 |
|
|---|
| 1223 | return
|
|---|
| 1224 | end
|
|---|
| 1225 | c-----------------------------------------------------------------------
|
|---|
| 1226 | subroutine setupds_center(gs_handle,nx,ny,nz,nel,melg,
|
|---|
| 1227 | $ vertex,glo_num)
|
|---|
| 1228 | include 'SIZE'
|
|---|
| 1229 | include 'INPUT'
|
|---|
| 1230 | include 'PARALLEL'
|
|---|
| 1231 | include 'NONCON'
|
|---|
| 1232 | integer gs_handle
|
|---|
| 1233 | integer vertex(1)
|
|---|
| 1234 | integer*8 glo_num(1),ngv
|
|---|
| 1235 | common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
|
|---|
| 1236 |
|
|---|
| 1237 | n = nx*ny*nz*nel
|
|---|
| 1238 | call set_vert(glo_num,ngv,nx,nel,vertex,.true.)
|
|---|
| 1239 | call fgslib_gs_setup(gs_handle,glo_num,n,nekcomm,mp)
|
|---|
| 1240 | return
|
|---|
| 1241 | end
|
|---|
| 1242 | c-----------------------------------------------------------------------
|
|---|
| 1243 | subroutine fastlap(iter,nodmask,mlt,gshl,dis,lapinv,mltc,gshlc)
|
|---|
| 1244 | include 'SIZE'
|
|---|
| 1245 | include 'TOTAL'
|
|---|
| 1246 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 1247 | real xmc(lxc,lyc,lzc,lelv),
|
|---|
| 1248 | $ ymc(lxc,lyc,lzc,lelv),
|
|---|
| 1249 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 1250 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 1251 |
|
|---|
| 1252 | real x8(2,2,ldim-1,lelv*(2**ldim)),
|
|---|
| 1253 | $ y8(2,2,ldim-1,lelv*(2**ldim)),
|
|---|
| 1254 | $ z8(2,2,ldim-1,lelv*(2**ldim))
|
|---|
| 1255 | common /coarsemesh8/ x8,y8,z8
|
|---|
| 1256 |
|
|---|
| 1257 | real dxx(2**ldim,lelv*(2**ldim)),
|
|---|
| 1258 | $ dyy(2**ldim,lelv*(2**ldim)),
|
|---|
| 1259 | $ dzz(2**ldim,lelv*(2**ldim)),
|
|---|
| 1260 | $ dis(2**ldim,lelv*(2**ldim)),
|
|---|
| 1261 | $ mlt(2**ldim,lelv*(2**ldim)),
|
|---|
| 1262 | $ nodmask(2**ldim,lelv*(2**ldim)),
|
|---|
| 1263 | $ mltc(lxc*lyc*lzc,lelv)
|
|---|
| 1264 | integer z,e,f,gshl,lapinv,kerr,gshlc
|
|---|
| 1265 | real xbar,ybar,zbar,sfac
|
|---|
| 1266 |
|
|---|
| 1267 | if (nid.eq.0.and.loglevel.ge.2)
|
|---|
| 1268 | $ write(6,*) 'laplacian smoothing started'
|
|---|
| 1269 |
|
|---|
| 1270 | n = 2**ldim
|
|---|
| 1271 | sfac = 0.01
|
|---|
| 1272 | do k=1,iter
|
|---|
| 1273 | do e=1,nelv*(2**ldim)
|
|---|
| 1274 | xbar = vlsum(x8(1,1,1,e),n)/(n)
|
|---|
| 1275 | ybar = vlsum(y8(1,1,1,e),n)/(n)
|
|---|
| 1276 | if (ldim.eq.3) zbar = vlsum(z8(1,1,1,e),n)/(n)
|
|---|
| 1277 | do i=1,n
|
|---|
| 1278 | dxx(i,e) = sfac*(xbar-x8(i,1,1,e))
|
|---|
| 1279 | dyy(i,e) = sfac*(ybar-y8(i,1,1,e))
|
|---|
| 1280 | if (ldim.eq.3) dzz(i,e) = sfac*(zbar-z8(i,1,1,e))
|
|---|
| 1281 | enddo
|
|---|
| 1282 | enddo
|
|---|
| 1283 |
|
|---|
| 1284 | n2 = nelv*(2**ldim)*(n)
|
|---|
| 1285 |
|
|---|
| 1286 | call col2(dxx,nodmask,n2)
|
|---|
| 1287 | call col2(dxx,dis,n2)
|
|---|
| 1288 | call dsavg_general(dxx,mlt,gshl)
|
|---|
| 1289 | call add2(x8,dxx,n2)
|
|---|
| 1290 |
|
|---|
| 1291 | call col2(dyy,nodmask,n2)
|
|---|
| 1292 | call col2(dyy,dis,n2)
|
|---|
| 1293 | call dsavg_general(dyy,mlt,gshl)
|
|---|
| 1294 | call add2(y8,dyy,n2)
|
|---|
| 1295 |
|
|---|
| 1296 | if (ldim.eq.3) then
|
|---|
| 1297 | call col2(dzz,nodmask,n2)
|
|---|
| 1298 | call col2(dzz,dis,n2)
|
|---|
| 1299 | call dsavg_general(dzz,mlt,gshl)
|
|---|
| 1300 | call add2(z8,dzz,n2)
|
|---|
| 1301 | endif
|
|---|
| 1302 |
|
|---|
| 1303 | dxm =glamax(dxx,n2)
|
|---|
| 1304 | dym =glamax(dyy,n2)
|
|---|
| 1305 | if (ldim.eq.3) then
|
|---|
| 1306 | dzm =glamax(dzz,n2)
|
|---|
| 1307 | else
|
|---|
| 1308 | dzm = 0.
|
|---|
| 1309 | endif
|
|---|
| 1310 | if (nio.eq.0.and.loglevel.ge.3) write(6,1) k,dxm,dym,dzm
|
|---|
| 1311 | 1 format(i5,1p3e12.3,' dxm')
|
|---|
| 1312 |
|
|---|
| 1313 | enddo
|
|---|
| 1314 |
|
|---|
| 1315 | call x8toxm(xmc,x8)
|
|---|
| 1316 | call x8toxm(ymc,y8)
|
|---|
| 1317 | if (ldim.eq.3) call x8toxm(zmc,z8)
|
|---|
| 1318 | call fixcurs(mltc,gshlc)
|
|---|
| 1319 | call xmtox8(xmc,x8)
|
|---|
| 1320 | call xmtox8(ymc,y8)
|
|---|
| 1321 | if (ldim.eq.3) call xmtox8(zmc,z8)
|
|---|
| 1322 |
|
|---|
| 1323 | call glmapm1chkinv(kerr)
|
|---|
| 1324 | if (kerr.gt.0.) then
|
|---|
| 1325 | lapinv = 1 !Terminate Laplacian smoothing
|
|---|
| 1326 | if (nid.eq.0.and.loglevel.ge.2) write(6,*) 'terminating Laplacian'
|
|---|
| 1327 | else
|
|---|
| 1328 | call genbackupmesh
|
|---|
| 1329 | endif
|
|---|
| 1330 |
|
|---|
| 1331 | return
|
|---|
| 1332 | end
|
|---|
| 1333 | c-----------------------------------------------------------------------
|
|---|
| 1334 | subroutine dsavg_general(u,mlt,gshl)
|
|---|
| 1335 | include 'SIZE'
|
|---|
| 1336 | include 'TOTAL'
|
|---|
| 1337 | integer gshl
|
|---|
| 1338 | real u(2**ldim,lelv*(2**ldim))
|
|---|
| 1339 | real mlt(2**ldim,lelv*(2**ldim))
|
|---|
| 1340 |
|
|---|
| 1341 | call fgslib_gs_op(gshl,u,1,1,0) ! '+'
|
|---|
| 1342 | call col2(u,mlt,(2**ldim)*nelv*(2**ldim))
|
|---|
| 1343 |
|
|---|
| 1344 | return
|
|---|
| 1345 | end
|
|---|
| 1346 | c-----------------------------------------------------------------------
|
|---|
| 1347 | subroutine xmtox8(xd,x8)
|
|---|
| 1348 | include 'SIZE'
|
|---|
| 1349 | include 'TOTAL'
|
|---|
| 1350 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 1351 | real x8(2**ldim,lelv*(2**ldim))
|
|---|
| 1352 | real xd(lxc*lyc*lzc,lelv)
|
|---|
| 1353 | integer e,k,n,j,ind1
|
|---|
| 1354 | integer zindx(64)
|
|---|
| 1355 | SAVE zindx
|
|---|
| 1356 | DATA zindx / 1, 2, 4, 5, 10, 11, 13, 14,
|
|---|
| 1357 | $ 2, 3, 5, 6, 11, 12, 14, 15,
|
|---|
| 1358 | $ 4, 5, 7, 8, 13, 14, 16, 17,
|
|---|
| 1359 | $ 5, 6, 8, 9, 14, 15, 17, 18,
|
|---|
| 1360 | $ 10, 11, 13, 14, 19, 20, 22, 23,
|
|---|
| 1361 | $ 11, 12, 14, 15, 20, 21, 23, 24,
|
|---|
| 1362 | $ 13, 14, 16, 17, 22, 23, 25, 26,
|
|---|
| 1363 | $ 14, 15, 17, 18, 23, 24, 26, 27 /
|
|---|
| 1364 |
|
|---|
| 1365 | n = 0
|
|---|
| 1366 | do e=1,nelv !each element
|
|---|
| 1367 | do j=1,2**ldim !broken down into 4 quads/8hexes
|
|---|
| 1368 | n = n+1
|
|---|
| 1369 | ind1 = (j-1)*8
|
|---|
| 1370 | do k=1,2**ldim
|
|---|
| 1371 | x8(k,n) = xd(zindx(ind1+k),e)
|
|---|
| 1372 | enddo
|
|---|
| 1373 | enddo
|
|---|
| 1374 | enddo
|
|---|
| 1375 |
|
|---|
| 1376 | return
|
|---|
| 1377 | end
|
|---|
| 1378 | c-----------------------------------------------------------------------
|
|---|
| 1379 | subroutine x8toxm(xd,x8)
|
|---|
| 1380 | include 'SIZE'
|
|---|
| 1381 | include 'TOTAL'
|
|---|
| 1382 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 1383 | real x8(2**ldim,lelv*(2**ldim))
|
|---|
| 1384 | real xd(lxc*lyc*lzc,lelv)
|
|---|
| 1385 | integer zindx((2**3)**2)
|
|---|
| 1386 | integer e,k,n,j,ind1
|
|---|
| 1387 | DATA zindx / 1, 2, 4, 5, 10, 11, 13, 14,
|
|---|
| 1388 | $ 2, 3, 5, 6, 11, 12, 14, 15,
|
|---|
| 1389 | $ 4, 5, 7, 8, 13, 14, 16, 17,
|
|---|
| 1390 | $ 5, 6, 8, 9, 14, 15, 17, 18,
|
|---|
| 1391 | $ 10, 11, 13, 14, 19, 20, 22, 23,
|
|---|
| 1392 | $ 11, 12, 14, 15, 20, 21, 23, 24,
|
|---|
| 1393 | $ 13, 14, 16, 17, 22, 23, 25, 26,
|
|---|
| 1394 | $ 14, 15, 17, 18, 23, 24, 26, 27 /
|
|---|
| 1395 |
|
|---|
| 1396 |
|
|---|
| 1397 | n = 0
|
|---|
| 1398 | do e=1,nelv !each element
|
|---|
| 1399 | do j=1,2**ldim !broken down into 4 quads/8hexes
|
|---|
| 1400 | n = n+1
|
|---|
| 1401 | ind1 = (j-1)*8
|
|---|
| 1402 | do k=1,2**ldim
|
|---|
| 1403 | xd(zindx(ind1+k),e) = x8(k,n)
|
|---|
| 1404 | enddo
|
|---|
| 1405 | enddo
|
|---|
| 1406 | enddo
|
|---|
| 1407 |
|
|---|
| 1408 | return
|
|---|
| 1409 | end
|
|---|
| 1410 | c-----------------------------------------------------------------------
|
|---|
| 1411 | subroutine col3c(a,b,c,d,n)
|
|---|
| 1412 | real a(1),b(1),c(1),d
|
|---|
| 1413 |
|
|---|
| 1414 | do i=1,n
|
|---|
| 1415 | a(i)=b(i)*c(i)*d
|
|---|
| 1416 | enddo
|
|---|
| 1417 |
|
|---|
| 1418 | return
|
|---|
| 1419 | end
|
|---|
| 1420 | c-----------------------------------------------------------------------
|
|---|
| 1421 | subroutine glmapm1chkinv(kerr)
|
|---|
| 1422 | INCLUDE 'SIZE'
|
|---|
| 1423 | INCLUDE 'GEOM'
|
|---|
| 1424 | INCLUDE 'INPUT'
|
|---|
| 1425 | INCLUDE 'SOLN'
|
|---|
| 1426 | C
|
|---|
| 1427 | C Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
|
|---|
| 1428 | C share the same array structure in Scratch Common /SCRNS/.
|
|---|
| 1429 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 1430 | parameter (nxc=3,nyc=3,nzc=1+(ldim-2)*(nxc-1))
|
|---|
| 1431 | real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
|
|---|
| 1432 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 1433 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 1434 | C
|
|---|
| 1435 | real XRMc(LXc,LYc,LZc,LELv)
|
|---|
| 1436 | $ , YRMc(LXc,LYc,LZc,LELv)
|
|---|
| 1437 | $ , XSMc(LXc,LYc,LZc,LELv)
|
|---|
| 1438 | $ , YSMc(LXc,LYc,LZc,LELv)
|
|---|
| 1439 | $ , XTMc(LXc,LYc,LZc,LELv)
|
|---|
| 1440 | $ , YTMc(LXc,LYc,LZc,LELv)
|
|---|
| 1441 | $ , ZRMc(LXc,LYc,LZc,LELv)
|
|---|
| 1442 | real ZSMc(LXc,LYc,LZc,LELv)
|
|---|
| 1443 | $ , ZTMc(LXc,LYc,LZc,LELv)
|
|---|
| 1444 | $ , jacmc(lxc,lyc,lzc,lelv)
|
|---|
| 1445 | real rxmc(lxc,lyc,lzc,lelv)
|
|---|
| 1446 | $ , rymc(lxc,lyc,lzc,lelv)
|
|---|
| 1447 | $ , rzmc(lxc,lyc,lzc,lelv)
|
|---|
| 1448 | $ , sxmc(lxc,lyc,lzc,lelv)
|
|---|
| 1449 | $ , symc(lxc,lyc,lzc,lelv)
|
|---|
| 1450 | $ , szmc(lxc,lyc,lzc,lelv)
|
|---|
| 1451 | $ , txmc(lxc,lyc,lzc,lelv)
|
|---|
| 1452 | $ , tymc(lxc,lyc,lzc,lelv)
|
|---|
| 1453 | $ , tzmc(lxc,lyc,lzc,lelv)
|
|---|
| 1454 | integer kerr,ierr
|
|---|
| 1455 | C
|
|---|
| 1456 | NXYc = NXc*NYc
|
|---|
| 1457 | NYZc = NYc*NZc
|
|---|
| 1458 | NXYZc = NXc*NYc*NZc
|
|---|
| 1459 | NTOTc = NXYZc*NELv
|
|---|
| 1460 | C
|
|---|
| 1461 | CALL XYZRSTc (XRMc,YRMc,ZRMc,XSMc,YSMc,ZSMc,XTMc,YTMc,ZTMc,
|
|---|
| 1462 | $ IFAXIS)
|
|---|
| 1463 |
|
|---|
| 1464 |
|
|---|
| 1465 | C
|
|---|
| 1466 | IF (NDIM.EQ.2) THEN
|
|---|
| 1467 | CALL RZERO (JACMc,NTOTc)
|
|---|
| 1468 | CALL ADDCOL3 (JACMc,XRMc,YSMc,NTOTc)
|
|---|
| 1469 | CALL SUBCOL3 (JACMc,XSMc,YRMc,NTOTc)
|
|---|
| 1470 | CALL COPY (RXMc,YSMc,NTOTc)
|
|---|
| 1471 | CALL COPY (RYMc,XSMc,NTOTc)
|
|---|
| 1472 | CALL CHSIGN (RYMc,NTOTc)
|
|---|
| 1473 | CALL COPY (SXMc,YRMc,NTOTc)
|
|---|
| 1474 | CALL CHSIGN (SXMc,NTOTc)
|
|---|
| 1475 | CALL COPY (SYMc,XRMc,NTOTc)
|
|---|
| 1476 | CALL RZERO (RZMc,NTOTc)
|
|---|
| 1477 | CALL RZERO (SZMc,NTOTc)
|
|---|
| 1478 | CALL RONE (TZMc,NTOTc)
|
|---|
| 1479 | ELSE
|
|---|
| 1480 | CALL RZERO (JACMc,NTOTc)
|
|---|
| 1481 | CALL ADDCOL4 (JACMc,XRMc,YSMc,ZTMc,NTOTc)
|
|---|
| 1482 | CALL ADDCOL4 (JACMc,XTMc,YRMc,ZSMc,NTOTc)
|
|---|
| 1483 | CALL ADDCOL4 (JACMc,XSMc,YTMc,ZRMc,NTOTc)
|
|---|
| 1484 | CALL SUBCOL4 (JACMc,XRMc,YTMc,ZSMc,NTOTc)
|
|---|
| 1485 | CALL SUBCOL4 (JACMc,XSMc,YRMc,ZTMc,NTOTc)
|
|---|
| 1486 | CALL SUBCOL4 (JACMc,XTMc,YSMc,ZRMc,NTOTc)
|
|---|
| 1487 | CALL ASCOL5 (RXMc,YSMc,ZTMc,YTMc,ZSMc,NTOTc)
|
|---|
| 1488 | CALL ASCOL5 (RYMc,XTMc,ZSMc,XSMc,ZTMc,NTOTc)
|
|---|
| 1489 | CALL ASCOL5 (RZMc,XSMc,YTMc,XTMc,YSMc,NTOTc)
|
|---|
| 1490 | CALL ASCOL5 (SXMc,YTMc,ZRMc,YRMc,ZTMc,NTOTc)
|
|---|
| 1491 | CALL ASCOL5 (SYMc,XRMc,ZTMc,XTMc,ZRMc,NTOTc)
|
|---|
| 1492 | CALL ASCOL5 (SZMc,XTMc,YRMc,XRMc,YTMc,NTOTc)
|
|---|
| 1493 | CALL ASCOL5 (TXMc,YRMc,ZSMc,YSMc,ZRMc,NTOTc)
|
|---|
| 1494 | CALL ASCOL5 (TYMc,XSMc,ZRMc,XRMc,ZSMc,NTOTc)
|
|---|
| 1495 | CALL ASCOL5 (TZMc,XRMc,YSMc,XSMc,YRMc,NTOTc)
|
|---|
| 1496 | ENDIF
|
|---|
| 1497 | C
|
|---|
| 1498 | kerr = 0
|
|---|
| 1499 | DO 500 ie=1,NELv
|
|---|
| 1500 | CALL CHKJACINV(JACMc(1,1,1,ie),NXYZc,ie,xmc(1,1,1,ie),
|
|---|
| 1501 | $ ymc(1,1,1,ie),zmc(1,1,1,ie),ndim,ierr)
|
|---|
| 1502 | if (ierr.ne.0) kerr = kerr+1
|
|---|
| 1503 | 500 CONTINUE
|
|---|
| 1504 | kerr = iglsum(kerr,1)
|
|---|
| 1505 |
|
|---|
| 1506 | RETURN
|
|---|
| 1507 | END
|
|---|
| 1508 | C-----------------------------------------------------------------------
|
|---|
| 1509 | subroutine chkjacinv(jac,n,iel,X,Y,Z,ND,IERR)
|
|---|
| 1510 | c
|
|---|
| 1511 | include 'SIZE'
|
|---|
| 1512 | include 'PARALLEL'
|
|---|
| 1513 | C
|
|---|
| 1514 | C Check the array JAC for a change in sign.
|
|---|
| 1515 | C
|
|---|
| 1516 | REAL JAC(N),x(1),y(1),z(1)
|
|---|
| 1517 | integer ierr
|
|---|
| 1518 | c
|
|---|
| 1519 | ierr = 0
|
|---|
| 1520 | SIGN = JAC(1)
|
|---|
| 1521 | DO 100 I=2,N
|
|---|
| 1522 | IF (SIGN*JAC(I).LE.0.0) THEN
|
|---|
| 1523 | ierr = 1
|
|---|
| 1524 | ENDIF
|
|---|
| 1525 | 100 CONTINUE
|
|---|
| 1526 |
|
|---|
| 1527 | RETURN
|
|---|
| 1528 | END
|
|---|
| 1529 | C-----------------------------------------------------------------------
|
|---|
| 1530 | subroutine xyzrstc(xrmc,yrmc,zrmc,xsmc,ysmc,zsmc,
|
|---|
| 1531 | $ XTMc,YTMc,ZTMc,IFAXIS)
|
|---|
| 1532 | C-----------------------------------------------------------------------
|
|---|
| 1533 | C
|
|---|
| 1534 | C Compute global-to-local derivatives on mesh 1.
|
|---|
| 1535 | C
|
|---|
| 1536 | C-----------------------------------------------------------------------
|
|---|
| 1537 | INCLUDE 'SIZE'
|
|---|
| 1538 | C
|
|---|
| 1539 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 1540 | DIMENSION XRMc(LXc,LYc,LZc,1),YRMc(LXc,LYc,LZc,1)
|
|---|
| 1541 | $ , ZRMc(LXc,LYc,LZc,1),XSMc(LXc,LYc,LZc,1)
|
|---|
| 1542 | $ , YSMc(LXc,LYc,LZc,1),ZSMc(LXc,LYc,LZc,1)
|
|---|
| 1543 | $ , XTMc(LXc,LYc,LZc,1),YTMc(LXc,LYc,LZc,1)
|
|---|
| 1544 | $ , ZTMc(LXc,LYc,LZc,1)
|
|---|
| 1545 | LOGICAL IFAXIS
|
|---|
| 1546 | real dxmc(3,3),dymc(3,3),dzmc(3,3)
|
|---|
| 1547 | real dxtmc(3,3),dytmc(3,3),dztmc(3,3)
|
|---|
| 1548 | common /coarseders/ dxmc,dymc,dzmc,dxtmc,dytmc,dztmc
|
|---|
| 1549 | real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
|
|---|
| 1550 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 1551 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 1552 | C
|
|---|
| 1553 | NXYc=lxc*lyc
|
|---|
| 1554 | NYZc=lyc*lzc
|
|---|
| 1555 | C
|
|---|
| 1556 | DO 100 IEL=1,NELv
|
|---|
| 1557 | C
|
|---|
| 1558 | CALL MXM (DXMc,lxc,XMc(1,1,1,IEL),lxc,XRMc(1,1,1,IEL),NYZc)
|
|---|
| 1559 | CALL MXM (DXMc,lxc,YMc(1,1,1,IEL),lxc,YRMc(1,1,1,IEL),NYZc)
|
|---|
| 1560 | CALL MXM (DXMc,lxc,ZMc(1,1,1,IEL),lxc,ZRMc(1,1,1,IEL),NYZc)
|
|---|
| 1561 | C
|
|---|
| 1562 | DO 10 IZ=1,lzc
|
|---|
| 1563 | CALL MXM (XMc(1,1,IZ,IEL),lxc,DYTMc,lyc,XSMc(1,1,IZ,IEL),lyc)
|
|---|
| 1564 | CALL MXM (YMc(1,1,IZ,IEL),lxc,DYTMc,lyc,YSMc(1,1,IZ,IEL),lyc)
|
|---|
| 1565 | CALL MXM (ZMc(1,1,IZ,IEL),lxc,DYTMc,lyc,ZSMc(1,1,IZ,IEL),lyc)
|
|---|
| 1566 | 10 CONTINUE
|
|---|
| 1567 | C
|
|---|
| 1568 | IF (ldim.EQ.3) THEN
|
|---|
| 1569 | CALL MXM (XMc(1,1,1,IEL),NXYc,DZTMc,lzc,XTMc(1,1,1,IEL),lzc)
|
|---|
| 1570 | CALL MXM (YMc(1,1,1,IEL),NXYc,DZTMc,lzc,YTMc(1,1,1,IEL),lzc)
|
|---|
| 1571 | CALL MXM (ZMc(1,1,1,IEL),NXYc,DZTMc,lzc,ZTMc(1,1,1,IEL),lzc)
|
|---|
| 1572 | ELSE
|
|---|
| 1573 | CALL RZERO (XTMc(1,1,1,IEL),NXYc)
|
|---|
| 1574 | CALL RZERO (YTMc(1,1,1,IEL),NXYc)
|
|---|
| 1575 | CALL RONE (ZTMc(1,1,1,IEL),NXYc)
|
|---|
| 1576 | ENDIF
|
|---|
| 1577 | C
|
|---|
| 1578 | 100 CONTINUE
|
|---|
| 1579 | C
|
|---|
| 1580 | RETURN
|
|---|
| 1581 | END
|
|---|
| 1582 | C-----------------------------------------------------------------------
|
|---|
| 1583 | subroutine xm1toxmc
|
|---|
| 1584 | include 'SIZE'
|
|---|
| 1585 | include 'TOTAL'
|
|---|
| 1586 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 1587 | real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
|
|---|
| 1588 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 1589 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 1590 | integer e
|
|---|
| 1591 |
|
|---|
| 1592 | if (ldim.eq.2) then
|
|---|
| 1593 | do e=1,nelv
|
|---|
| 1594 | call int_fine_to_coarse_2d(xm1(1,1,1,e),xmc(1,1,1,e),lx1,3)
|
|---|
| 1595 | call int_fine_to_coarse_2d(ym1(1,1,1,e),ymc(1,1,1,e),lx1,3)
|
|---|
| 1596 | enddo
|
|---|
| 1597 | else
|
|---|
| 1598 | do e=1,nelv
|
|---|
| 1599 | call int_fine_to_coarse_3d(xm1(1,1,1,e),xmc(1,1,1,e),lx1,3)
|
|---|
| 1600 | call int_fine_to_coarse_3d(ym1(1,1,1,e),ymc(1,1,1,e),lx1,3)
|
|---|
| 1601 | call int_fine_to_coarse_3d(zm1(1,1,1,e),zmc(1,1,1,e),lx1,3)
|
|---|
| 1602 | enddo
|
|---|
| 1603 | endif
|
|---|
| 1604 |
|
|---|
| 1605 | RETURN
|
|---|
| 1606 | END
|
|---|
| 1607 | C-----------------------------------------------------------------------
|
|---|
| 1608 | subroutine xmctoxm1
|
|---|
| 1609 | include 'SIZE'
|
|---|
| 1610 | include 'TOTAL'
|
|---|
| 1611 | parameter (lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
|
|---|
| 1612 | real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
|
|---|
| 1613 | $ zmc(lxc,lyc,lzc,lelv)
|
|---|
| 1614 | common /coarsemesh/ xmc,ymc,zmc
|
|---|
| 1615 | integer e
|
|---|
| 1616 |
|
|---|
| 1617 | if (ldim.eq.2) then
|
|---|
| 1618 | do e=1,nelv
|
|---|
| 1619 | call int_coarse_to_fine_2d(xm1(1,1,1,e),xmc(1,1,1,e),lx1,3)
|
|---|
| 1620 | call int_coarse_to_fine_2d(ym1(1,1,1,e),ymc(1,1,1,e),lx1,3)
|
|---|
| 1621 | enddo
|
|---|
| 1622 | else
|
|---|
| 1623 | do e=1,nelv
|
|---|
| 1624 | call int_coarse_to_fine_3d(xm1(1,1,1,e),xmc(1,1,1,e),lx1,3)
|
|---|
| 1625 | call int_coarse_to_fine_3d(ym1(1,1,1,e),ymc(1,1,1,e),lx1,3)
|
|---|
| 1626 | call int_coarse_to_fine_3d(zm1(1,1,1,e),zmc(1,1,1,e),lx1,3)
|
|---|
| 1627 | enddo
|
|---|
| 1628 | endif
|
|---|
| 1629 |
|
|---|
| 1630 | RETURN
|
|---|
| 1631 | END
|
|---|
| 1632 | C-----------------------------------------------------------------------
|
|---|
| 1633 | subroutine my_mv_mesh(dx,dy,dz)
|
|---|
| 1634 | include 'SIZE'
|
|---|
| 1635 | include 'TOTAL'
|
|---|
| 1636 |
|
|---|
| 1637 | parameter (lt = lx1*ly1*lz1*lelv)
|
|---|
| 1638 | common /mrthoi/ napprx(2),nappry(2),napprz(2)
|
|---|
| 1639 | common /mrthov/ apprx(lt,0:mxprev)
|
|---|
| 1640 | $ , appry(lt,0:mxprev)
|
|---|
| 1641 | $ , apprz(lt,0:mxprev)
|
|---|
| 1642 | common /mstuff/ d(lt),h1(lt),h2(lt),mask(lt)
|
|---|
| 1643 | real mask
|
|---|
| 1644 | real mask2(lx1,ly1,lz1,lelv)
|
|---|
| 1645 | real dx(1),dy(1),dz(1)
|
|---|
| 1646 |
|
|---|
| 1647 | integer e,f
|
|---|
| 1648 |
|
|---|
| 1649 | integer icalld
|
|---|
| 1650 | save icalld
|
|---|
| 1651 | data icalld /0/
|
|---|
| 1652 |
|
|---|
| 1653 | n = nx1*ny1*nz1*nelv
|
|---|
| 1654 | nface = 2*ndim
|
|---|
| 1655 |
|
|---|
| 1656 | napprx(1)=0
|
|---|
| 1657 | nappry(1)=0
|
|---|
| 1658 | napprz(1)=0
|
|---|
| 1659 |
|
|---|
| 1660 | nxz = nx1*nz1
|
|---|
| 1661 | nxyz = nx1*ny1*nz1
|
|---|
| 1662 |
|
|---|
| 1663 | volbl = 0. ! Volume of elements in boundary layer
|
|---|
| 1664 | do e=1,nelv
|
|---|
| 1665 | do f=1,nface
|
|---|
| 1666 | if (cbc(f,e,1).eq.'W ') then
|
|---|
| 1667 | srfbl = srfbl + vlsum(area(1,1,f,e),nxz )
|
|---|
| 1668 | volbl = volbl + vlsum(bm1 (1,1,1,e),nxyz)
|
|---|
| 1669 | endif
|
|---|
| 1670 | enddo
|
|---|
| 1671 | enddo
|
|---|
| 1672 | srfbl = glsum(srfbl,1) ! Sum over all processors
|
|---|
| 1673 | volbl = glsum(volbl,1)
|
|---|
| 1674 |
|
|---|
| 1675 | delta = volbl / srfbl ! Avg thickness of b.l. elements
|
|---|
| 1676 |
|
|---|
| 1677 | call rone (h1,n)
|
|---|
| 1678 | call rzero(h2,n)
|
|---|
| 1679 | call sm_cheap_dist(d,1,'W ')
|
|---|
| 1680 |
|
|---|
| 1681 | deltap = 5*delta ! Protected b.l. thickness
|
|---|
| 1682 | do i=1,n
|
|---|
| 1683 | arg = -(d(i)/deltap)**2
|
|---|
| 1684 | h1(i) = h1(i) + 9*exp(arg)
|
|---|
| 1685 | enddo
|
|---|
| 1686 |
|
|---|
| 1687 | call rone(mask,n)
|
|---|
| 1688 | do e=1,nelv
|
|---|
| 1689 | do f=1,nface
|
|---|
| 1690 | if (cbc(f,e,1).ne.'E '.and.cbc(f,e,1).ne.' ')
|
|---|
| 1691 | $ call facev(mask,e,f,0.,nx1,ny1,nz1)
|
|---|
| 1692 | enddo
|
|---|
| 1693 | enddo
|
|---|
| 1694 | call dsop(mask,'mul',nx1,ny1,nz1)
|
|---|
| 1695 | call copy(mask2,mask,n)
|
|---|
| 1696 | call chsign(mask2,n)
|
|---|
| 1697 | call cadd(mask2,1.,n)
|
|---|
| 1698 | call col2(dx,mask2,n)
|
|---|
| 1699 | call col2(dy,mask2,n)
|
|---|
| 1700 | call col2(dz,mask2,n)
|
|---|
| 1701 | tol = -1.e-6
|
|---|
| 1702 | call laplaceh('mshx',dx,h1,h2,mask,vmult,1,tol,500,apprx,napprx)
|
|---|
| 1703 | call laplaceh('mshy',dy,h1,h2,mask,vmult,1,tol,500,appry,nappry)
|
|---|
| 1704 | if (if3d)
|
|---|
| 1705 | $ call laplaceh('mshz',dz,h1,h2,mask,vmult,1,tol,500,apprz,napprz)
|
|---|
| 1706 |
|
|---|
| 1707 | return
|
|---|
| 1708 | end
|
|---|
| 1709 | c-----------------------------------------------------------------------
|
|---|
| 1710 | subroutine laplaceh
|
|---|
| 1711 | $ (name,u,h1,h2,mask,mult,ifld,tli,maxi,approx,napprox)
|
|---|
| 1712 | include 'SIZE'
|
|---|
| 1713 | include 'TOTAL'
|
|---|
| 1714 | include 'CTIMER'
|
|---|
| 1715 | c
|
|---|
| 1716 | character*4 name
|
|---|
| 1717 | real u(1),h1(1),h2(1),mask(1),mult(1),approx (1)
|
|---|
| 1718 | integer napprox(1)
|
|---|
| 1719 |
|
|---|
| 1720 | parameter (lt=lx1*ly1*lz1*lelv)
|
|---|
| 1721 | common /scruz/ r (lt),ub(lt)
|
|---|
| 1722 |
|
|---|
| 1723 | logical ifstdh
|
|---|
| 1724 | character*4 cname
|
|---|
| 1725 | character*6 name6
|
|---|
| 1726 |
|
|---|
| 1727 | logical ifwt,ifvec
|
|---|
| 1728 |
|
|---|
| 1729 | call chcopy(cname,name,4)
|
|---|
| 1730 | call capit (cname,4)
|
|---|
| 1731 |
|
|---|
| 1732 | call blank (name6,6)
|
|---|
| 1733 | call chcopy(name6,name,4)
|
|---|
| 1734 | ifwt = .true.
|
|---|
| 1735 | ifvec = .false.
|
|---|
| 1736 | isd = 1
|
|---|
| 1737 | imsh = 1
|
|---|
| 1738 | nel = nelfld(ifld)
|
|---|
| 1739 |
|
|---|
| 1740 | n = nx1*ny1*nz1*nel
|
|---|
| 1741 |
|
|---|
| 1742 | call copy (ub,u,n) ! ub = u on boundary
|
|---|
| 1743 | call dsavg(ub) ! Make certain ub is in H1
|
|---|
| 1744 | ! _
|
|---|
| 1745 | call axhelm (r,ub,h1,h2,1,1) ! r = A*ub
|
|---|
| 1746 |
|
|---|
| 1747 | do i=1,n ! _
|
|---|
| 1748 | r(i)=-r(i)*mask(i) ! r = -M*A*ub
|
|---|
| 1749 | enddo
|
|---|
| 1750 |
|
|---|
| 1751 | call dssum (r,nx1,ny1,nz1) ! dssum rhs
|
|---|
| 1752 |
|
|---|
| 1753 | call project1
|
|---|
| 1754 | $ (r,n,approx,napprox,h1,h2,mask,mult,ifwt,ifvec,name6)
|
|---|
| 1755 |
|
|---|
| 1756 | tol = -abs(tli)
|
|---|
| 1757 | if (nel.eq.nelv) then
|
|---|
| 1758 | call hmhzpf (name,u,r,h1,h2,mask,mult,imsh,tol,maxi,isd,binvm1)
|
|---|
| 1759 | else
|
|---|
| 1760 | call hmhzpf (name,u,r,h1,h2,mask,mult,imsh,tol,maxi,isd,bintm1)
|
|---|
| 1761 | endif
|
|---|
| 1762 |
|
|---|
| 1763 | call project2
|
|---|
| 1764 | $ (u,n,approx,napprox,h1,h2,mask,mult,ifwt,ifvec,name6)
|
|---|
| 1765 |
|
|---|
| 1766 | call add2(u,ub,n)
|
|---|
| 1767 |
|
|---|
| 1768 | return
|
|---|
| 1769 | end
|
|---|
| 1770 | c-----------------------------------------------------------------------
|
|---|
| 1771 | subroutine sm_cheap_dist(d,ifld,b)
|
|---|
| 1772 | c this is a copy of the routine in navier5.f, modified to be less
|
|---|
| 1773 | c verbose
|
|---|
| 1774 |
|
|---|
| 1775 |
|
|---|
| 1776 | include 'SIZE'
|
|---|
| 1777 | include 'GEOM' ! Coordinates
|
|---|
| 1778 | include 'INPUT' ! cbc()
|
|---|
| 1779 | include 'TSTEP' ! nelfld
|
|---|
| 1780 | include 'PARALLEL' ! gather-scatter handle for field "ifld"
|
|---|
| 1781 |
|
|---|
| 1782 | real d(lx1,ly1,lz1,lelt)
|
|---|
| 1783 |
|
|---|
| 1784 | character*3 b ! Boundary condition of interest
|
|---|
| 1785 |
|
|---|
| 1786 | integer e,eg,f
|
|---|
| 1787 |
|
|---|
| 1788 | nel = nelfld(ifld)
|
|---|
| 1789 | n = lx1*ly1*lz1*nel
|
|---|
| 1790 |
|
|---|
| 1791 | call domain_size(xmin,xmax,ymin,ymax,zmin,zmax)
|
|---|
| 1792 |
|
|---|
| 1793 | xmn = min(xmin,ymin)
|
|---|
| 1794 | xmx = max(xmax,ymax)
|
|---|
| 1795 | if (if3d) xmn = min(xmn ,zmin)
|
|---|
| 1796 | if (if3d) xmx = max(xmx ,zmax)
|
|---|
| 1797 |
|
|---|
| 1798 | big = 10*(xmx-xmn)
|
|---|
| 1799 | call cfill(d,big,n)
|
|---|
| 1800 |
|
|---|
| 1801 | nface = 2*ldim
|
|---|
| 1802 | do e=1,nel ! Set d=0 on walls
|
|---|
| 1803 | do f=1,nface
|
|---|
| 1804 | if (cbc(f,e,ifld).eq.b) call facev(d,e,f,0.,lx1,ly1,lz1)
|
|---|
| 1805 | enddo
|
|---|
| 1806 | enddo
|
|---|
| 1807 |
|
|---|
| 1808 | do ipass=1,10000
|
|---|
| 1809 | dmax = 0
|
|---|
| 1810 | nchange = 0
|
|---|
| 1811 | do e=1,nel
|
|---|
| 1812 | do k=1,lz1
|
|---|
| 1813 | do j=1,ly1
|
|---|
| 1814 | do i=1,lx1
|
|---|
| 1815 | i0=max( 1,i-1)
|
|---|
| 1816 | j0=max( 1,j-1)
|
|---|
| 1817 | k0=max( 1,k-1)
|
|---|
| 1818 | i1=min(lx1,i+1)
|
|---|
| 1819 | j1=min(ly1,j+1)
|
|---|
| 1820 | k1=min(lz1,k+1)
|
|---|
| 1821 | do kk=k0,k1
|
|---|
| 1822 | do jj=j0,j1
|
|---|
| 1823 | do ii=i0,i1
|
|---|
| 1824 |
|
|---|
| 1825 | if (if3d) then
|
|---|
| 1826 | dtmp = d(ii,jj,kk,e) + dist3d(
|
|---|
| 1827 | $ xm1(ii,jj,kk,e),ym1(ii,jj,kk,e),zm1(ii,jj,kk,e)
|
|---|
| 1828 | $ ,xm1(i ,j ,k ,e),ym1(i ,j ,k ,e),zm1(i ,j ,k ,e))
|
|---|
| 1829 | else
|
|---|
| 1830 | dtmp = d(ii,jj,kk,e) + dist2d(
|
|---|
| 1831 | $ xm1(ii,jj,kk,e),ym1(ii,jj,kk,e)
|
|---|
| 1832 | $ ,xm1(i ,j ,k ,e),ym1(i ,j ,k ,e))
|
|---|
| 1833 | endif
|
|---|
| 1834 |
|
|---|
| 1835 | if (dtmp.lt.d(i,j,k,e)) then
|
|---|
| 1836 | d(i,j,k,e) = dtmp
|
|---|
| 1837 | nchange = nchange+1
|
|---|
| 1838 | dmax = max(dmax,d(i,j,k,e))
|
|---|
| 1839 | endif
|
|---|
| 1840 | enddo
|
|---|
| 1841 | enddo
|
|---|
| 1842 | enddo
|
|---|
| 1843 |
|
|---|
| 1844 | enddo
|
|---|
| 1845 | enddo
|
|---|
| 1846 | enddo
|
|---|
| 1847 | enddo
|
|---|
| 1848 | call fgslib_gs_op(gsh_fld(ifld),d,1,3,0) ! min over all elements
|
|---|
| 1849 | nchange = iglsum(nchange,1)
|
|---|
| 1850 | dmax = glmax(dmax,1)
|
|---|
| 1851 | if (nio.eq.0.and.loglevel.ge.3) write(6,1) ipass,nchange,dmax,b
|
|---|
| 1852 | 1 format(i9,i12,1pe12.4,' max distance b: ',a3)
|
|---|
| 1853 | if (nchange.eq.0) goto 1000
|
|---|
| 1854 | enddo
|
|---|
| 1855 | 1000 return
|
|---|
| 1856 | end
|
|---|
| 1857 | c-----------------------------------------------------------------------
|
|---|
