| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | c
|
|---|
| 3 | c Some relevant parameters
|
|---|
| 4 | c
|
|---|
| 5 | c param(41):
|
|---|
| 6 | c 0 - use additive SEMG
|
|---|
| 7 | c 1 - use hybrid SEMG (not yet working... but coming soon!)
|
|---|
| 8 | c
|
|---|
| 9 | c param(42): navier0.f, fasts.f
|
|---|
| 10 | c 0 - use GMRES for iterative solver, use non-symmetric weighting
|
|---|
| 11 | c 1 - use PCG for iterative solver, do not use weighting
|
|---|
| 12 | c
|
|---|
| 13 | c param(43): uzawa_gmres.f, navier6.f
|
|---|
| 14 | c 0 - use additive multilevel scheme (requires param(42).eq.0)
|
|---|
| 15 | c 1 - use original 2 level scheme
|
|---|
| 16 | c
|
|---|
| 17 | c param(44): fast3d.f, navier6.f
|
|---|
| 18 | c 0 - base top-level additive Schwarz on restrictions of E
|
|---|
| 19 | c 1 - base top-level additive Schwarz on restrictions of A
|
|---|
| 20 | c
|
|---|
| 21 | c----------------------------------------------------------------------
|
|---|
| 22 | subroutine hsmg_setup()
|
|---|
| 23 | include 'SIZE'
|
|---|
| 24 | include 'INPUT'
|
|---|
| 25 | include 'PARALLEL'
|
|---|
| 26 | include 'HSMG'
|
|---|
| 27 | include 'SEMHAT'
|
|---|
| 28 | include 'TSTEP'
|
|---|
| 29 |
|
|---|
| 30 | integer nf,nc,nr
|
|---|
| 31 | integer nx,ny,nz
|
|---|
| 32 |
|
|---|
| 33 | mg_fld = 1
|
|---|
| 34 | if (ifield.gt.1) mg_fld = 2
|
|---|
| 35 | if (ifield.eq.1) call hsmg_index_0 ! initialize index sets
|
|---|
| 36 |
|
|---|
| 37 | call hsmg_setup_mg_nx ! set nx values for each level of multigrid
|
|---|
| 38 | call hsmg_setup_semhat ! set spectral element hat matrices
|
|---|
| 39 | call hsmg_setup_intp
|
|---|
| 40 | call hsmg_setup_dssum ! set direct stiffness summation handles
|
|---|
| 41 | call hsmg_setup_wtmask ! set restriction weight matrices and bc masks
|
|---|
| 42 | call hsmg_setup_fdm ! set up fast diagonalization method
|
|---|
| 43 | call hsmg_setup_schwarz_wt(.false.)
|
|---|
| 44 | call hsmg_setup_solve ! set up the solver
|
|---|
| 45 | c call hsmg_setup_dbg
|
|---|
| 46 |
|
|---|
| 47 | return
|
|---|
| 48 | end
|
|---|
| 49 | c----------------------------------------------------------------------
|
|---|
| 50 | subroutine hsmg_setup_semhat
|
|---|
| 51 | include 'SIZE'
|
|---|
| 52 | include 'INPUT'
|
|---|
| 53 | include 'HSMG'
|
|---|
| 54 | include 'SEMHAT'
|
|---|
| 55 | integer n,l
|
|---|
| 56 | c generate the SEM hat matrices for each level
|
|---|
| 57 | c top level
|
|---|
| 58 | n = mg_nx(mg_lmax)
|
|---|
| 59 | call semhat(ah,bh,ch,dh,zh,dph,jph,bgl,zglhat,dgl,jgl,n,wh)
|
|---|
| 60 | call copy(mg_zh(1,mg_lmax),zglhat,n-1) !top level based on gl points
|
|---|
| 61 | mg_nh(mg_lmax)=n-1
|
|---|
| 62 | mg_nhz(mg_lmax)=n-1
|
|---|
| 63 | if(.not.if3d) mg_nhz(mg_lmax)=1
|
|---|
| 64 | c lower levels
|
|---|
| 65 | do l=1,mg_lmax-1
|
|---|
| 66 | n = mg_nx(l)
|
|---|
| 67 | if(n.gt.1) then
|
|---|
| 68 | call semhat(ah,bh,ch,dh,zh,dph,jph,bgl,zglhat,dgl,jgl,n,wh)
|
|---|
| 69 | call copy(mg_ah(1,l),ah,(n+1)*(n+1))
|
|---|
| 70 | call copy(mg_bh(1,l),bh,n+1)
|
|---|
| 71 | call copy(mg_dh(1,l),dh,(n+1)*(n+1))
|
|---|
| 72 | call transpose(mg_dht(1,l),n+1,dh,n+1)
|
|---|
| 73 | call copy(mg_zh(1,l),zh,n+1)
|
|---|
| 74 | else
|
|---|
| 75 | mg_zh(1,l) = -1.
|
|---|
| 76 | mg_zh(2,l) = 1.
|
|---|
| 77 | endif
|
|---|
| 78 | mg_nh(l)=n+1
|
|---|
| 79 | mg_nhz(l)=mg_nz(l)+1
|
|---|
| 80 | enddo
|
|---|
| 81 | end
|
|---|
| 82 | c----------------------------------------------------------------------
|
|---|
| 83 | subroutine hsmg_setup_intp
|
|---|
| 84 | include 'SIZE'
|
|---|
| 85 | include 'HSMG'
|
|---|
| 86 | include 'SEMHAT'
|
|---|
| 87 | integer l,nf,nc
|
|---|
| 88 |
|
|---|
| 89 | do l=1,mg_lmax-1
|
|---|
| 90 |
|
|---|
| 91 | nf=mg_nh(l+1)
|
|---|
| 92 | nc=mg_nh(l)
|
|---|
| 93 |
|
|---|
| 94 | ! Standard multigrid coarse-to-fine interpolation
|
|---|
| 95 | call hsmg_setup_intpm(
|
|---|
| 96 | $ mg_jh(1,l),mg_zh(1,l+1),mg_zh(1,l),nf,nc)
|
|---|
| 97 | call transpose(mg_jht(1,l),nc,mg_jh(1,l),nf)
|
|---|
| 98 |
|
|---|
| 99 | ! Fine-to-coarse interpolation for variable-coefficient operators
|
|---|
| 100 | call hsmg_setup_intpm(
|
|---|
| 101 | $ mg_jhfc(1,l),mg_zh(1,l),mg_zh(1,l+1),nc,nf)
|
|---|
| 102 | call transpose(mg_jhfct(1,l),nf,mg_jhfc(1,l),nc)
|
|---|
| 103 | c call outmat(mg_jhfc(1,l),nc,nf,'MG_JHFC',l)
|
|---|
| 104 |
|
|---|
| 105 | enddo
|
|---|
| 106 | end
|
|---|
| 107 | c----------------------------------------------------------------------
|
|---|
| 108 | subroutine hsmg_setup_intpm(jh,zf,zc,nf,nc)
|
|---|
| 109 | integer nf,nc
|
|---|
| 110 | real jh(nf,nc),zf(1),zc(1)
|
|---|
| 111 | include 'SIZE'
|
|---|
| 112 | real w(2*lx1+2)
|
|---|
| 113 | do i=1,nf
|
|---|
| 114 | call fd_weights_full(zf(i),zc,nc-1,1,w)
|
|---|
| 115 | do j=1,nc
|
|---|
| 116 | jh(i,j)=w(j)
|
|---|
| 117 | enddo
|
|---|
| 118 | enddo
|
|---|
| 119 | return
|
|---|
| 120 | end
|
|---|
| 121 | c----------------------------------------------------------------------
|
|---|
| 122 | subroutine hsmg_setup_dssum
|
|---|
| 123 | include 'SIZE'
|
|---|
| 124 | include 'INPUT'
|
|---|
| 125 | include 'PARALLEL'
|
|---|
| 126 | include 'HSMG'
|
|---|
| 127 | parameter (lxyz=(lx1+2)*(ly1+2)*(lz1+2))
|
|---|
| 128 | common /c_is1/ glo_num(lxyz*lelv)
|
|---|
| 129 | common /ivrtx/ vertex ((2**ldim)*lelt)
|
|---|
| 130 |
|
|---|
| 131 | integer*8 glo_num
|
|---|
| 132 | integer vertex
|
|---|
| 133 | integer nx,ny,nz
|
|---|
| 134 | integer l
|
|---|
| 135 |
|
|---|
| 136 | c set up direct stiffness summation for each level
|
|---|
| 137 | ncrnr = 2**ldim
|
|---|
| 138 | call get_vert
|
|---|
| 139 |
|
|---|
| 140 | c++ write(6,*) mg_fld,' mgfld in hsmg_setup_dssum'
|
|---|
| 141 |
|
|---|
| 142 | do l=1,mg_lmax-1
|
|---|
| 143 | nx=mg_nh(l)
|
|---|
| 144 | ny=mg_nh(l)
|
|---|
| 145 | nz=mg_nhz(l)
|
|---|
| 146 | call setupds(mg_gsh_handle(l,mg_fld),nx,ny,nz
|
|---|
| 147 | $ ,nelv,nelgv,vertex,glo_num)
|
|---|
| 148 | nx=nx+2
|
|---|
| 149 | ny=ny+2
|
|---|
| 150 | nz=nz+2
|
|---|
| 151 | if(.not.if3d) nz=1
|
|---|
| 152 | call setupds(mg_gsh_schwarz_handle(l,mg_fld),nx,ny,nz
|
|---|
| 153 | $ ,nelv,nelgv,vertex,glo_num)
|
|---|
| 154 | enddo
|
|---|
| 155 | end
|
|---|
| 156 | c----------------------------------------------------------------------
|
|---|
| 157 | subroutine h1mg_setup_wtmask
|
|---|
| 158 | include 'SIZE'
|
|---|
| 159 | include 'HSMG'
|
|---|
| 160 | integer i,l
|
|---|
| 161 | i = mg_mask_index(mg_lmax,mg_fld-1)
|
|---|
| 162 | do l=1,mg_lmax
|
|---|
| 163 | mg_rstr_wt_index(l,mg_fld)=i
|
|---|
| 164 | mg_mask_index (l,mg_fld)=i
|
|---|
| 165 | i=i+mg_nh(l)*mg_nhz(l)*2*ldim*nelv
|
|---|
| 166 | if(i .gt. lmgs*lmg_rwt*2*ldim*lelv) then
|
|---|
| 167 | itmp = i/(2*ldim*lelv)
|
|---|
| 168 | write(6,*) 'parameter lmg_rwt too small',i,itmp,lmg_rwt
|
|---|
| 169 | call exitt
|
|---|
| 170 | endif
|
|---|
| 171 | call hsmg_setup_rstr_wt(
|
|---|
| 172 | $ mg_rstr_wt(mg_rstr_wt_index(l,mg_fld))
|
|---|
| 173 | $ ,mg_nh(l),mg_nh(l),mg_nhz(l),l,mg_work)
|
|---|
| 174 | call hsmg_setup_mask(
|
|---|
| 175 | $ mg_mask(mg_mask_index(l,mg_fld))
|
|---|
| 176 | $ ,mg_nh(l),mg_nh(l),mg_nhz(l),l,mg_work)
|
|---|
| 177 | enddo
|
|---|
| 178 | mg_mask_index(l,mg_fld)=i
|
|---|
| 179 | end
|
|---|
| 180 | c----------------------------------------------------------------------
|
|---|
| 181 | subroutine hsmg_setup_wtmask
|
|---|
| 182 | include 'SIZE'
|
|---|
| 183 | include 'HSMG'
|
|---|
| 184 | integer i,l
|
|---|
| 185 | i = mg_mask_index(mg_lmax,mg_fld-1)
|
|---|
| 186 | do l=1,mg_lmax-1
|
|---|
| 187 | mg_rstr_wt_index(l,mg_fld)=i
|
|---|
| 188 | mg_mask_index (l,mg_fld)=i
|
|---|
| 189 | i=i+mg_nh(l)*mg_nhz(l)*2*ldim*nelv
|
|---|
| 190 | if(i .gt. lmgs*lmg_rwt*2*ldim*lelv) then
|
|---|
| 191 | itmp = i/(2*ldim*lelv)
|
|---|
| 192 | write(6,*) 'parameter lmg_rwt too small',i,itmp,lmg_rwt
|
|---|
| 193 | call exitt
|
|---|
| 194 | endif
|
|---|
| 195 | call hsmg_setup_rstr_wt(
|
|---|
| 196 | $ mg_rstr_wt(mg_rstr_wt_index(l,mg_fld))
|
|---|
| 197 | $ ,mg_nh(l),mg_nh(l),mg_nhz(l),l,mg_work)
|
|---|
| 198 | call hsmg_setup_mask(
|
|---|
| 199 | $ mg_mask(mg_mask_index(l,mg_fld))
|
|---|
| 200 | $ ,mg_nh(l),mg_nh(l),mg_nhz(l),l,mg_work)
|
|---|
| 201 | enddo
|
|---|
| 202 | mg_mask_index(l,mg_fld)=i
|
|---|
| 203 | end
|
|---|
| 204 | c----------------------------------------------------------------------
|
|---|
| 205 | subroutine hsmg_intp(uf,uc,l) ! l is coarse level
|
|---|
| 206 | real uf(1),uc(1)
|
|---|
| 207 | integer l
|
|---|
| 208 | include 'SIZE'
|
|---|
| 209 | include 'HSMG'
|
|---|
| 210 | call hsmg_tnsr(uf,mg_nh(l+1),uc,mg_nh(l),mg_jh(1,l),mg_jht(1,l))
|
|---|
| 211 | return
|
|---|
| 212 | end
|
|---|
| 213 | c----------------------------------------------------------------------
|
|---|
| 214 | subroutine hsmg_rstr(uc,uf,l) ! l is coarse level
|
|---|
| 215 | real uf(1),uc(1)
|
|---|
| 216 | integer l
|
|---|
| 217 | include 'SIZE'
|
|---|
| 218 | include 'HSMG'
|
|---|
| 219 | if(l.ne.mg_lmax-1)
|
|---|
| 220 | $ call hsmg_do_wt(uf,mg_rstr_wt(mg_rstr_wt_index(l+1,mg_fld))
|
|---|
| 221 | $ ,mg_nh(l+1),mg_nh(l+1),mg_nhz(l+1))
|
|---|
| 222 | call hsmg_tnsr(uc,mg_nh(l),uf,mg_nh(l+1),mg_jht(1,l),mg_jh(1,l))
|
|---|
| 223 | call hsmg_dssum(uc,l)
|
|---|
| 224 | return
|
|---|
| 225 | end
|
|---|
| 226 | c----------------------------------------------------------------------
|
|---|
| 227 | subroutine hsmg_rstr_no_dssum(uc,uf,l) ! l is coarse level
|
|---|
| 228 | real uf(1),uc(1)
|
|---|
| 229 | integer l
|
|---|
| 230 | include 'SIZE'
|
|---|
| 231 | include 'HSMG'
|
|---|
| 232 | if(l.ne.mg_lmax-1)
|
|---|
| 233 | $ call hsmg_do_wt(uf,mg_rstr_wt(mg_rstr_wt_index(l+1,mg_fld))
|
|---|
| 234 | $ ,mg_nh(l+1),mg_nh(l+1),mg_nhz(l+1))
|
|---|
| 235 | call hsmg_tnsr(uc,mg_nh(l),uf,mg_nh(l+1),mg_jht(1,l),mg_jh(1,l))
|
|---|
| 236 | return
|
|---|
| 237 | end
|
|---|
| 238 | c----------------------------------------------------------------------
|
|---|
| 239 | c computes
|
|---|
| 240 | c v = [A (x) A] u or
|
|---|
| 241 | c v = [A (x) A (x) A] u
|
|---|
| 242 | subroutine hsmg_tnsr(v,nv,u,nu,A,At)
|
|---|
| 243 | integer nv,nu
|
|---|
| 244 | real v(1),u(1),A(1),At(1)
|
|---|
| 245 | include 'SIZE'
|
|---|
| 246 | include 'INPUT'
|
|---|
| 247 | if (.not. if3d) then
|
|---|
| 248 | call hsmg_tnsr2d(v,nv,u,nu,A,At)
|
|---|
| 249 | else
|
|---|
| 250 | call hsmg_tnsr3d(v,nv,u,nu,A,At,At)
|
|---|
| 251 | endif
|
|---|
| 252 | return
|
|---|
| 253 | end
|
|---|
| 254 | c----------------------------------------------------------------------
|
|---|
| 255 | c computes
|
|---|
| 256 | c T
|
|---|
| 257 | c v = A u B
|
|---|
| 258 | subroutine hsmg_tnsr2d(v,nv,u,nu,A,Bt)
|
|---|
| 259 | integer nv,nu
|
|---|
| 260 | real v(nv*nv,nelv),u(nu*nu,nelv),A(1),Bt(1)
|
|---|
| 261 | include 'SIZE'
|
|---|
| 262 | common /hsmgw/ work((lx1+2)*(lx1+2))
|
|---|
| 263 | integer ie
|
|---|
| 264 | do ie=1,nelv
|
|---|
| 265 | call mxm(A,nv,u(1,ie),nu,work,nu)
|
|---|
| 266 | call mxm(work,nv,Bt,nu,v(1,ie),nv)
|
|---|
| 267 | enddo
|
|---|
| 268 | return
|
|---|
| 269 | end
|
|---|
| 270 | c----------------------------------------------------------------------
|
|---|
| 271 | c computes
|
|---|
| 272 | c
|
|---|
| 273 | c v = [C (x) B (x) A] u
|
|---|
| 274 | subroutine hsmg_tnsr3d(v,nv,u,nu,A,Bt,Ct)
|
|---|
| 275 | integer nv,nu
|
|---|
| 276 | real v(nv*nv*nv,nelv),u(nu*nu*nu,nelv),A(1),Bt(1),Ct(1)
|
|---|
| 277 | include 'SIZE'
|
|---|
| 278 | parameter (lwk=(lx1+2)*(ly1+2)*(lz1+2))
|
|---|
| 279 | common /hsmgw/ work(0:lwk-1),work2(0:lwk-1)
|
|---|
| 280 | integer ie, i
|
|---|
| 281 | do ie=1,nelv
|
|---|
| 282 | call mxm(A,nv,u(1,ie),nu,work,nu*nu)
|
|---|
| 283 | do i=0,nu-1
|
|---|
| 284 | call mxm(work(nv*nu*i),nv,Bt,nu,work2(nv*nv*i),nv)
|
|---|
| 285 | enddo
|
|---|
| 286 | call mxm(work2,nv*nv,Ct,nu,v(1,ie),nv)
|
|---|
| 287 | enddo
|
|---|
| 288 | return
|
|---|
| 289 | end
|
|---|
| 290 | c----------------------------------------------------------------------
|
|---|
| 291 | c computes
|
|---|
| 292 | c T
|
|---|
| 293 | c v = A u B
|
|---|
| 294 | subroutine hsmg_tnsr2d_el(v,nv,u,nu,A,Bt)
|
|---|
| 295 | integer nv,nu
|
|---|
| 296 | real v(nv*nv),u(nu*nu),A(1),Bt(1)
|
|---|
| 297 | include 'SIZE'
|
|---|
| 298 | common /hsmgw/ work((lx1+2)*(lx1+2))
|
|---|
| 299 | c
|
|---|
| 300 | call mxm(A,nv,u,nu,work,nu)
|
|---|
| 301 | call mxm(work,nv,Bt,nu,v,nv)
|
|---|
| 302 | c
|
|---|
| 303 | return
|
|---|
| 304 | end
|
|---|
| 305 | c----------------------------------------------------------------------
|
|---|
| 306 | c computes
|
|---|
| 307 | c
|
|---|
| 308 | c v = [C (x) B (x) A] u
|
|---|
| 309 | subroutine hsmg_tnsr3d_el(v,nv,u,nu,A,Bt,Ct)
|
|---|
| 310 | integer nv,nu
|
|---|
| 311 | real v(nv*nv*nv),u(nu*nu*nu),A(1),Bt(1),Ct(1)
|
|---|
| 312 | include 'SIZE'
|
|---|
| 313 | parameter (lwk=(lx1+2)*(ly1+2)*(lz1+2))
|
|---|
| 314 | common /hsmgw/ work(0:lwk-1),work2(0:lwk-1)
|
|---|
| 315 | integer i
|
|---|
| 316 | c
|
|---|
| 317 | call mxm(A,nv,u,nu,work,nu*nu)
|
|---|
| 318 | do i=0,nu-1
|
|---|
| 319 | call mxm(work(nv*nu*i),nv,Bt,nu,work2(nv*nv*i),nv)
|
|---|
| 320 | enddo
|
|---|
| 321 | call mxm(work2,nv*nv,Ct,nu,v,nv)
|
|---|
| 322 | c
|
|---|
| 323 | return
|
|---|
| 324 | end
|
|---|
| 325 | c----------------------------------------------------------------------
|
|---|
| 326 | subroutine hsmg_dssum(u,l)
|
|---|
| 327 | include 'SIZE'
|
|---|
| 328 | include 'HSMG'
|
|---|
| 329 | include 'CTIMER'
|
|---|
| 330 | real u(1)
|
|---|
| 331 |
|
|---|
| 332 | if (ifsync) call nekgsync()
|
|---|
| 333 | etime1=dnekclock()
|
|---|
| 334 |
|
|---|
| 335 | call fgslib_gs_op(mg_gsh_handle(l,mg_fld),u,1,1,0)
|
|---|
| 336 | tdadd =tdadd + dnekclock()-etime1
|
|---|
| 337 |
|
|---|
| 338 |
|
|---|
| 339 | return
|
|---|
| 340 | end
|
|---|
| 341 | c----------------------------------------------------------------------
|
|---|
| 342 | subroutine hsmg_dsprod(u,l)
|
|---|
| 343 | include 'SIZE'
|
|---|
| 344 | include 'HSMG'
|
|---|
| 345 | include 'CTIMER'
|
|---|
| 346 | real u(1)
|
|---|
| 347 |
|
|---|
| 348 | if (ifsync) call nekgsync()
|
|---|
| 349 |
|
|---|
| 350 | call fgslib_gs_op(mg_gsh_handle(l,mg_fld),u,1,2,0)
|
|---|
| 351 | return
|
|---|
| 352 | end
|
|---|
| 353 | c----------------------------------------------------------------------
|
|---|
| 354 | subroutine hsmg_schwarz_dssum(u,l)
|
|---|
| 355 | include 'SIZE'
|
|---|
| 356 | include 'HSMG'
|
|---|
| 357 | include 'CTIMER'
|
|---|
| 358 |
|
|---|
| 359 | if (ifsync) call nekgsync()
|
|---|
| 360 | etime1=dnekclock()
|
|---|
| 361 |
|
|---|
| 362 | call fgslib_gs_op(mg_gsh_schwarz_handle(l,mg_fld),u,1,1,0)
|
|---|
| 363 | tdadd =tdadd + dnekclock()-etime1
|
|---|
| 364 |
|
|---|
| 365 | return
|
|---|
| 366 | end
|
|---|
| 367 | c----------------------------------------------------------------------
|
|---|
| 368 | subroutine hsmg_extrude(arr1,l1,f1,arr2,l2,f2,nx,ny,nz)
|
|---|
| 369 | include 'SIZE'
|
|---|
| 370 | include 'INPUT'
|
|---|
| 371 | integer l1,l2,nx,ny,nz
|
|---|
| 372 | real arr1(nx,ny,nz,nelv),arr2(nx,ny,nz,nelv)
|
|---|
| 373 | real f1,f2
|
|---|
| 374 |
|
|---|
| 375 | integer i,j,k,ie,i0,i1
|
|---|
| 376 | i0=2
|
|---|
| 377 | i1=nx-1
|
|---|
| 378 |
|
|---|
| 379 | if(.not.if3d) then
|
|---|
| 380 | do ie=1,nelv
|
|---|
| 381 | do j=i0,i1
|
|---|
| 382 | arr1(l1+1 ,j,1,ie) = f1*arr1(l1+1 ,j,1,ie)
|
|---|
| 383 | $ +f2*arr2(l2+1 ,j,1,ie)
|
|---|
| 384 | arr1(nx-l1,j,1,ie) = f1*arr1(nx-l1,j,1,ie)
|
|---|
| 385 | $ +f2*arr2(nx-l2,j,1,ie)
|
|---|
| 386 | enddo
|
|---|
| 387 | do i=i0,i1
|
|---|
| 388 | arr1(i,l1+1 ,1,ie) = f1*arr1(i,l1+1 ,1,ie)
|
|---|
| 389 | $ +f2*arr2(i,l2+1 ,1,ie)
|
|---|
| 390 | arr1(i,ny-l1,1,ie) = f1*arr1(i,ny-l1,1,ie)
|
|---|
| 391 | $ +f2*arr2(i,nx-l2,1,ie)
|
|---|
| 392 | enddo
|
|---|
| 393 | enddo
|
|---|
| 394 | else
|
|---|
| 395 | do ie=1,nelv
|
|---|
| 396 | do k=i0,i1
|
|---|
| 397 | do j=i0,i1
|
|---|
| 398 | arr1(l1+1 ,j,k,ie) = f1*arr1(l1+1 ,j,k,ie)
|
|---|
| 399 | $ +f2*arr2(l2+1 ,j,k,ie)
|
|---|
| 400 | arr1(nx-l1,j,k,ie) = f1*arr1(nx-l1,j,k,ie)
|
|---|
| 401 | $ +f2*arr2(nx-l2,j,k,ie)
|
|---|
| 402 | enddo
|
|---|
| 403 | enddo
|
|---|
| 404 | do k=i0,i1
|
|---|
| 405 | do i=i0,i1
|
|---|
| 406 | arr1(i,l1+1 ,k,ie) = f1*arr1(i,l1+1 ,k,ie)
|
|---|
| 407 | $ +f2*arr2(i,l2+1 ,k,ie)
|
|---|
| 408 | arr1(i,nx-l1,k,ie) = f1*arr1(i,nx-l1,k,ie)
|
|---|
| 409 | $ +f2*arr2(i,nx-l2,k,ie)
|
|---|
| 410 | enddo
|
|---|
| 411 | enddo
|
|---|
| 412 | do j=i0,i1
|
|---|
| 413 | do i=i0,i1
|
|---|
| 414 | arr1(i,j,l1+1 ,ie) = f1*arr1(i,j,l1+1 ,ie)
|
|---|
| 415 | $ +f2*arr2(i,j,l2+1 ,ie)
|
|---|
| 416 | arr1(i,j,nx-l1,ie) = f1*arr1(i,j,nx-l1,ie)
|
|---|
| 417 | $ +f2*arr2(i,j,nx-l2,ie)
|
|---|
| 418 | enddo
|
|---|
| 419 | enddo
|
|---|
| 420 | enddo
|
|---|
| 421 | endif
|
|---|
| 422 | return
|
|---|
| 423 | end
|
|---|
| 424 | c----------------------------------------------------------------------
|
|---|
| 425 | subroutine h1mg_schwarz(e,r,sigma,l)
|
|---|
| 426 | include 'SIZE'
|
|---|
| 427 | include 'HSMG'
|
|---|
| 428 |
|
|---|
| 429 | real e(1),r(1)
|
|---|
| 430 |
|
|---|
| 431 | n = mg_h1_n(l,mg_fld)
|
|---|
| 432 |
|
|---|
| 433 | call h1mg_schwarz_part1 (e,r,l)
|
|---|
| 434 | call hsmg_schwarz_wt (e,l) ! e := W e
|
|---|
| 435 | call cmult (e,sigma,n) ! l l
|
|---|
| 436 |
|
|---|
| 437 | return
|
|---|
| 438 | end
|
|---|
| 439 | c----------------------------------------------------------------------
|
|---|
| 440 | subroutine h1mg_schwarz_part1 (e,r,l)
|
|---|
| 441 | include 'SIZE'
|
|---|
| 442 | include 'INPUT' ! if3d
|
|---|
| 443 | include 'TSTEP' ! ifield
|
|---|
| 444 | include 'HSMG'
|
|---|
| 445 |
|
|---|
| 446 | real e(1),r(1)
|
|---|
| 447 |
|
|---|
| 448 | integer enx,eny,enz,pm
|
|---|
| 449 |
|
|---|
| 450 | zero = 0
|
|---|
| 451 | one = 1
|
|---|
| 452 | onem = -1
|
|---|
| 453 |
|
|---|
| 454 | n = mg_h1_n (l,mg_fld)
|
|---|
| 455 | pm = p_mg_msk(l,mg_fld)
|
|---|
| 456 |
|
|---|
| 457 | call h1mg_mask (r,mg_imask(pm),nelfld(ifield)) ! Zero Dirichlet nodes
|
|---|
| 458 |
|
|---|
| 459 | if (if3d) then ! extended array
|
|---|
| 460 | call hsmg_schwarz_toext3d(mg_work,r,mg_nh(l))
|
|---|
| 461 | else
|
|---|
| 462 | call hsmg_schwarz_toext2d(mg_work,r,mg_nh(l))
|
|---|
| 463 | endif
|
|---|
| 464 |
|
|---|
| 465 | enx=mg_nh(l)+2
|
|---|
| 466 | eny=mg_nh(l)+2
|
|---|
| 467 | enz=mg_nh(l)+2
|
|---|
| 468 | if(.not.if3d) enz=1
|
|---|
| 469 | i = enx*eny*enz*nelv+1
|
|---|
| 470 |
|
|---|
| 471 | c exchange interior nodes
|
|---|
| 472 | call hsmg_extrude(mg_work,0,zero,mg_work,2,one,enx,eny,enz)
|
|---|
| 473 | call hsmg_schwarz_dssum(mg_work,l)
|
|---|
| 474 | call hsmg_extrude(mg_work,0,one ,mg_work,2,onem,enx,eny,enz)
|
|---|
| 475 |
|
|---|
| 476 | call hsmg_fdm(mg_work(i),mg_work,l) ! Do the local solves
|
|---|
| 477 |
|
|---|
| 478 | c Sum overlap region (border excluded)
|
|---|
| 479 | call hsmg_extrude(mg_work,0,zero,mg_work(i),0,one ,enx,eny,enz)
|
|---|
| 480 | call hsmg_schwarz_dssum(mg_work(i),l)
|
|---|
| 481 | call hsmg_extrude(mg_work(i),0,one ,mg_work,0,onem,enx,eny,enz)
|
|---|
| 482 | call hsmg_extrude(mg_work(i),2,one,mg_work(i),0,one,enx,eny,enz)
|
|---|
| 483 |
|
|---|
| 484 | if(.not.if3d) then ! Go back to regular size array
|
|---|
| 485 | call hsmg_schwarz_toreg2d(e,mg_work(i),mg_nh(l))
|
|---|
| 486 | else
|
|---|
| 487 | call hsmg_schwarz_toreg3d(e,mg_work(i),mg_nh(l))
|
|---|
| 488 | endif
|
|---|
| 489 |
|
|---|
| 490 | call hsmg_dssum(e,l) ! sum border nodes
|
|---|
| 491 | call h1mg_mask (e,mg_imask(pm),nelfld(ifield)) ! apply mask
|
|---|
| 492 |
|
|---|
| 493 | return
|
|---|
| 494 | end
|
|---|
| 495 | c----------------------------------------------------------------------
|
|---|
| 496 | subroutine hsmg_schwarz(e,r,l)
|
|---|
| 497 | include 'SIZE'
|
|---|
| 498 | include 'INPUT'
|
|---|
| 499 | include 'HSMG'
|
|---|
| 500 | real e(1),r(1)
|
|---|
| 501 | integer l
|
|---|
| 502 | integer enx,eny,enz
|
|---|
| 503 | integer i
|
|---|
| 504 |
|
|---|
| 505 | real zero,one,onem
|
|---|
| 506 | zero = 0
|
|---|
| 507 | one = 1
|
|---|
| 508 | onem = -1
|
|---|
| 509 |
|
|---|
| 510 | c apply mask (zeros Dirichlet nodes)
|
|---|
| 511 | !!!!! uncommenting
|
|---|
| 512 | call hsmg_do_wt(r,mg_mask(mg_mask_index(l,mg_fld)),
|
|---|
| 513 | $ mg_nh(l),mg_nh(l),mg_nhz(l))
|
|---|
| 514 |
|
|---|
| 515 | c go to extended size array (room for overlap)
|
|---|
| 516 | if (if3d) then
|
|---|
| 517 | call hsmg_schwarz_toext3d(mg_work,r,mg_nh(l))
|
|---|
| 518 | else
|
|---|
| 519 | call hsmg_schwarz_toext2d(mg_work,r,mg_nh(l))
|
|---|
| 520 | endif
|
|---|
| 521 |
|
|---|
| 522 | enx=mg_nh(l)+2
|
|---|
| 523 | eny=mg_nh(l)+2
|
|---|
| 524 | enz=mg_nh(l)+2
|
|---|
| 525 | if(.not.if3d) enz=1
|
|---|
| 526 | i = enx*eny*enz*nelv+1
|
|---|
| 527 |
|
|---|
| 528 | c exchange interior nodes
|
|---|
| 529 | call hsmg_extrude(mg_work,0,zero,mg_work,2,one,enx,eny,enz)
|
|---|
| 530 | call hsmg_schwarz_dssum(mg_work,l)
|
|---|
| 531 | call hsmg_extrude(mg_work,0,one ,mg_work,2,onem,enx,eny,enz)
|
|---|
| 532 |
|
|---|
| 533 | c do the local solves
|
|---|
| 534 | call hsmg_fdm(mg_work(i),mg_work,l)
|
|---|
| 535 | c sum overlap region (border excluded)
|
|---|
| 536 | call hsmg_extrude(mg_work,0,zero,mg_work(i),0,one ,enx,eny,enz)
|
|---|
| 537 | call hsmg_schwarz_dssum(mg_work(i),l)
|
|---|
| 538 | call hsmg_extrude(mg_work(i),0,one ,mg_work,0,onem,enx,eny,enz)
|
|---|
| 539 | call hsmg_extrude(mg_work(i),2,one,mg_work(i),0,one,enx,eny,enz)
|
|---|
| 540 | c go back to regular size array
|
|---|
| 541 | if(.not.if3d) then
|
|---|
| 542 | call hsmg_schwarz_toreg2d(e,mg_work(i),mg_nh(l))
|
|---|
| 543 | else
|
|---|
| 544 | call hsmg_schwarz_toreg3d(e,mg_work(i),mg_nh(l))
|
|---|
| 545 | endif
|
|---|
| 546 | c sum border nodes
|
|---|
| 547 | call hsmg_dssum(e,l)
|
|---|
| 548 | c apply mask (zeros Dirichlet nodes)
|
|---|
| 549 | !!!!!! changing r to e
|
|---|
| 550 | call hsmg_do_wt(e,mg_mask(mg_mask_index(l,mg_fld)),
|
|---|
| 551 | $ mg_nh(l),mg_nh(l),mg_nhz(l))
|
|---|
| 552 | return
|
|---|
| 553 | end
|
|---|
| 554 | c----------------------------------------------------------------------
|
|---|
| 555 | subroutine hsmg_schwarz_toext2d(a,b,n)
|
|---|
| 556 | include 'SIZE'
|
|---|
| 557 | integer n
|
|---|
| 558 | real a(0:n+1,0:n+1,nelv),b(n,n,nelv)
|
|---|
| 559 |
|
|---|
| 560 | integer i,j,ie
|
|---|
| 561 | c call rzero(a,(n+2)*(n+2)*nelv)
|
|---|
| 562 | do ie=1,nelv
|
|---|
| 563 | do i=0,n+1
|
|---|
| 564 | a(i,0,ie)=0.
|
|---|
| 565 | enddo
|
|---|
| 566 | do j=1,n
|
|---|
| 567 | a(0 ,j,ie)=0.
|
|---|
| 568 | do i=1,n
|
|---|
| 569 | a(i,j,ie)=b(i,j,ie)
|
|---|
| 570 | enddo
|
|---|
| 571 | a(n+1,j,ie)=0.
|
|---|
| 572 | enddo
|
|---|
| 573 | do i=0,n+1
|
|---|
| 574 | a(i,n+1,ie)=0.
|
|---|
| 575 | enddo
|
|---|
| 576 | enddo
|
|---|
| 577 | return
|
|---|
| 578 | end
|
|---|
| 579 | c----------------------------------------------------------------------
|
|---|
| 580 | subroutine hsmg_schwarz_toext3d(a,b,n)
|
|---|
| 581 | include 'SIZE'
|
|---|
| 582 | integer n
|
|---|
| 583 | real a(0:n+1,0:n+1,0:n+1,nelv),b(n,n,n,nelv)
|
|---|
| 584 |
|
|---|
| 585 | integer i,j,k,ie
|
|---|
| 586 | call rzero(a,(n+2)*(n+2)*(n+2)*nelv)
|
|---|
| 587 | do ie=1,nelv
|
|---|
| 588 | do k=1,n
|
|---|
| 589 | do j=1,n
|
|---|
| 590 | do i=1,n
|
|---|
| 591 | a(i,j,k,ie)=b(i,j,k,ie)
|
|---|
| 592 | enddo
|
|---|
| 593 | enddo
|
|---|
| 594 | enddo
|
|---|
| 595 | enddo
|
|---|
| 596 | return
|
|---|
| 597 | end
|
|---|
| 598 | c----------------------------------------------------------------------
|
|---|
| 599 | subroutine hsmg_schwarz_toreg2d(b,a,n)
|
|---|
| 600 | include 'SIZE'
|
|---|
| 601 | integer n
|
|---|
| 602 | real a(0:n+1,0:n+1,nelv),b(n,n,nelv)
|
|---|
| 603 |
|
|---|
| 604 | integer i,j,ie
|
|---|
| 605 | do ie=1,nelv
|
|---|
| 606 | do j=1,n
|
|---|
| 607 | do i=1,n
|
|---|
| 608 | b(i,j,ie)=a(i,j,ie)
|
|---|
| 609 | enddo
|
|---|
| 610 | enddo
|
|---|
| 611 | enddo
|
|---|
| 612 | return
|
|---|
| 613 | end
|
|---|
| 614 | c----------------------------------------------------------------------
|
|---|
| 615 | subroutine hsmg_schwarz_toreg3d(b,a,n)
|
|---|
| 616 | include 'SIZE'
|
|---|
| 617 | integer n
|
|---|
| 618 | real a(0:n+1,0:n+1,0:n+1,nelv),b(n,n,n,nelv)
|
|---|
| 619 |
|
|---|
| 620 | integer i,j,k,ie
|
|---|
| 621 | do ie=1,nelv
|
|---|
| 622 | do k=1,n
|
|---|
| 623 | do j=1,n
|
|---|
| 624 | do i=1,n
|
|---|
| 625 | b(i,j,k,ie)=a(i,j,k,ie)
|
|---|
| 626 | enddo
|
|---|
| 627 | enddo
|
|---|
| 628 | enddo
|
|---|
| 629 | enddo
|
|---|
| 630 | return
|
|---|
| 631 | end
|
|---|
| 632 | c----------------------------------------------------------------------
|
|---|
| 633 | subroutine h1mg_setup_fdm()
|
|---|
| 634 | include 'SIZE'
|
|---|
| 635 | include 'INPUT'
|
|---|
| 636 | include 'HSMG'
|
|---|
| 637 |
|
|---|
| 638 | integer l,i,j,nl
|
|---|
| 639 | i = mg_fast_s_index(mg_lmax,mg_fld-1)
|
|---|
| 640 | j = mg_fast_d_index(mg_lmax,mg_fld-1)
|
|---|
| 641 | do l=2,mg_lmax
|
|---|
| 642 | mg_fast_s_index(l,mg_fld)=i
|
|---|
| 643 | nl = mg_nh(l)+2
|
|---|
| 644 | i=i+nl*nl*2*ldim*nelv
|
|---|
| 645 | if(i .gt. lmg_fasts*2*ldim*lelv) then
|
|---|
| 646 | itmp = i/(2*ldim*lelv)
|
|---|
| 647 | write(6,*) 'lmg_fasts too small',i,itmp,lmg_fasts,l
|
|---|
| 648 | call exitt
|
|---|
| 649 | endif
|
|---|
| 650 | mg_fast_d_index(l,mg_fld)=j
|
|---|
| 651 | j=j+(nl**ldim)*nelv
|
|---|
| 652 | if(j .gt. lmg_fastd*lelv) then
|
|---|
| 653 | itmp = i/(2*ldim*lelv)
|
|---|
| 654 | write(6,*) 'lmg_fastd too small',i,itmp,lmg_fastd,l
|
|---|
| 655 | call exitt
|
|---|
| 656 | endif
|
|---|
| 657 | call hsmg_setup_fast(
|
|---|
| 658 | $ mg_fast_s(mg_fast_s_index(l,mg_fld))
|
|---|
| 659 | $ ,mg_fast_d(mg_fast_d_index(l,mg_fld))
|
|---|
| 660 | $ ,mg_nh(l)+2,mg_ah(1,l),mg_bh(1,l),mg_nx(l))
|
|---|
| 661 | enddo
|
|---|
| 662 | mg_fast_s_index(l,mg_fld)=i
|
|---|
| 663 | mg_fast_d_index(l,mg_fld)=j
|
|---|
| 664 | return
|
|---|
| 665 | end
|
|---|
| 666 | c----------------------------------------------------------------------
|
|---|
| 667 | subroutine hsmg_setup_fdm()
|
|---|
| 668 | include 'SIZE'
|
|---|
| 669 | include 'INPUT'
|
|---|
| 670 | include 'HSMG'
|
|---|
| 671 |
|
|---|
| 672 | integer l,i,j,nl
|
|---|
| 673 | i = mg_fast_s_index(mg_lmax,mg_fld-1)
|
|---|
| 674 | j = mg_fast_d_index(mg_lmax,mg_fld-1)
|
|---|
| 675 | do l=2,mg_lmax-1
|
|---|
| 676 | mg_fast_s_index(l,mg_fld)=i
|
|---|
| 677 | nl = mg_nh(l)+2
|
|---|
| 678 | i=i+nl*nl*2*ldim*nelv
|
|---|
| 679 | if(i .gt. lmg_fasts*2*ldim*lelv) then
|
|---|
| 680 | itmp = i/(2*ldim*lelv)
|
|---|
| 681 | write(6,*) 'lmg_fasts too small',i,itmp,lmg_fasts,l
|
|---|
| 682 | call exitt
|
|---|
| 683 | endif
|
|---|
| 684 | mg_fast_d_index(l,mg_fld)=j
|
|---|
| 685 | j=j+(nl**ldim)*nelv
|
|---|
| 686 | if(j .gt. lmg_fastd*lelv) then
|
|---|
| 687 | itmp = i/(2*ldim*lelv)
|
|---|
| 688 | write(6,*) 'lmg_fastd too small',i,itmp,lmg_fastd,l
|
|---|
| 689 | call exitt
|
|---|
| 690 | endif
|
|---|
| 691 | call hsmg_setup_fast(
|
|---|
| 692 | $ mg_fast_s(mg_fast_s_index(l,mg_fld))
|
|---|
| 693 | $ ,mg_fast_d(mg_fast_d_index(l,mg_fld))
|
|---|
| 694 | $ ,mg_nh(l)+2,mg_ah(1,l),mg_bh(1,l),mg_nx(l))
|
|---|
| 695 | enddo
|
|---|
| 696 | mg_fast_s_index(l,mg_fld)=i
|
|---|
| 697 | mg_fast_d_index(l,mg_fld)=j
|
|---|
| 698 | return
|
|---|
| 699 | end
|
|---|
| 700 | c----------------------------------------------------------------------
|
|---|
| 701 | subroutine hsmg_setup_fast(s,d,nl,ah,bh,n)
|
|---|
| 702 | include 'SIZE'
|
|---|
| 703 | include 'INPUT'
|
|---|
| 704 | include 'HSMG'
|
|---|
| 705 | real s(nl*nl,2,ldim,nelv)
|
|---|
| 706 | real d(nl**ldim,nelv)
|
|---|
| 707 | real ah(1),bh(1)
|
|---|
| 708 | common /ctmpf/ lr(2*lx1+4),ls(2*lx1+4),lt(2*lx1+4)
|
|---|
| 709 | $ , llr(lelt),lls(lelt),llt(lelt)
|
|---|
| 710 | $ , lmr(lelt),lms(lelt),lmt(lelt)
|
|---|
| 711 | $ , lrr(lelt),lrs(lelt),lrt(lelt)
|
|---|
| 712 | real lr ,ls ,lt
|
|---|
| 713 | real llr,lls,llt
|
|---|
| 714 | real lmr,lms,lmt
|
|---|
| 715 | real lrr,lrs,lrt
|
|---|
| 716 |
|
|---|
| 717 | integer i,j,k
|
|---|
| 718 | integer ie,il,nr,ns,nt
|
|---|
| 719 | integer lbr,rbr,lbs,rbs,lbt,rbt,two
|
|---|
| 720 | real eps,diag
|
|---|
| 721 |
|
|---|
| 722 | two = 2
|
|---|
| 723 | ierr = 0
|
|---|
| 724 | do ie=1,nelv
|
|---|
| 725 | call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,ie,two,ierr)
|
|---|
| 726 | nr=nl
|
|---|
| 727 | ns=nl
|
|---|
| 728 | nt=nl
|
|---|
| 729 | call hsmg_setup_fast1d(s(1,1,1,ie),lr,nr,lbr,rbr
|
|---|
| 730 | $ ,llr(ie),lmr(ie),lrr(ie),ah,bh,n,ie)
|
|---|
| 731 | call hsmg_setup_fast1d(s(1,1,2,ie),ls,ns,lbs,rbs
|
|---|
| 732 | $ ,lls(ie),lms(ie),lrs(ie),ah,bh,n,ie)
|
|---|
| 733 | if(if3d) call hsmg_setup_fast1d(s(1,1,3,ie),lt,nt,lbt,rbt
|
|---|
| 734 | $ ,llt(ie),lmt(ie),lrt(ie),ah,bh,n,ie)
|
|---|
| 735 | il=1
|
|---|
| 736 | if(.not.if3d) then
|
|---|
| 737 | eps = 1.e-5*(vlmax(lr(2),nr-2) + vlmax(ls(2),ns-2))
|
|---|
| 738 | do j=1,ns
|
|---|
| 739 | do i=1,nr
|
|---|
| 740 | diag = lr(i)+ls(j)
|
|---|
| 741 | if (diag.gt.eps) then
|
|---|
| 742 | d(il,ie) = 1.0/diag
|
|---|
| 743 | else
|
|---|
| 744 | c write(6,2) ie,'Reset Eig in hsmg setup fast:',i,j,l
|
|---|
| 745 | c $ ,eps,diag,lr(i),ls(j)
|
|---|
| 746 | 2 format(i6,1x,a21,3i5,1p4e12.4)
|
|---|
| 747 | d(il,ie) = 0.0
|
|---|
| 748 | endif
|
|---|
| 749 | il=il+1
|
|---|
| 750 | enddo
|
|---|
| 751 | enddo
|
|---|
| 752 | else
|
|---|
| 753 | eps = 1.e-5 * (vlmax(lr(2),nr-2)
|
|---|
| 754 | $ + vlmax(ls(2),ns-2) + vlmax(lt(2),nt-2))
|
|---|
| 755 | do k=1,nt
|
|---|
| 756 | do j=1,ns
|
|---|
| 757 | do i=1,nr
|
|---|
| 758 | diag = lr(i)+ls(j)+lt(k)
|
|---|
| 759 | if (diag.gt.eps) then
|
|---|
| 760 | d(il,ie) = 1.0/diag
|
|---|
| 761 | else
|
|---|
| 762 | c write(6,3) ie,'Reset Eig in hsmg setup fast:',i,j,k,l
|
|---|
| 763 | c $ ,eps,diag,lr(i),ls(j),lt(k)
|
|---|
| 764 | 3 format(i6,1x,a21,4i5,1p5e12.4)
|
|---|
| 765 | d(il,ie) = 0.0
|
|---|
| 766 | endif
|
|---|
| 767 | il=il+1
|
|---|
| 768 | enddo
|
|---|
| 769 | enddo
|
|---|
| 770 | enddo
|
|---|
| 771 | endif
|
|---|
| 772 | enddo
|
|---|
| 773 |
|
|---|
| 774 | ierrmx = iglmax(ierr,1)
|
|---|
| 775 | if (ierrmx.gt.0) then
|
|---|
| 776 | if (ierr.gt.0) write(6,*) nid,ierr,' BC FAIL'
|
|---|
| 777 | call exitti('A INVALID BC FOUND in genfast$',ierrmx)
|
|---|
| 778 | endif
|
|---|
| 779 |
|
|---|
| 780 | return
|
|---|
| 781 | end
|
|---|
| 782 | c----------------------------------------------------------------------
|
|---|
| 783 | subroutine hsmg_setup_fast1d(s,lam,nl,lbc,rbc,ll,lm,lr,ah,bh,n,ie)
|
|---|
| 784 | integer nl,lbc,rbc,n
|
|---|
| 785 | real s(nl,nl,2),lam(nl),ll,lm,lr
|
|---|
| 786 | real ah(0:n,0:n),bh(0:n)
|
|---|
| 787 |
|
|---|
| 788 | include 'SIZE'
|
|---|
| 789 | parameter(lxm=lx1+2)
|
|---|
| 790 | common /ctmp0/ b(2*lxm*lxm),w(2*lxm*lxm)
|
|---|
| 791 |
|
|---|
| 792 | call hsmg_setup_fast1d_a(s,lbc,rbc,ll,lm,lr,ah,n)
|
|---|
| 793 | call hsmg_setup_fast1d_b(b,lbc,rbc,ll,lm,lr,bh,n)
|
|---|
| 794 |
|
|---|
| 795 | c if (nid.eq.0) write(6,*) 'THIS is generalev call',nl,lbc
|
|---|
| 796 | call generalev(s,b,lam,nl,w)
|
|---|
| 797 | if(lbc.gt.0) call row_zero(s,nl,nl,1)
|
|---|
| 798 | if(lbc.eq.1) call row_zero(s,nl,nl,2)
|
|---|
| 799 | if(rbc.gt.0) call row_zero(s,nl,nl,nl)
|
|---|
| 800 | if(rbc.eq.1) call row_zero(s,nl,nl,nl-1)
|
|---|
| 801 |
|
|---|
| 802 | call transpose(s(1,1,2),nl,s,nl)
|
|---|
| 803 | return
|
|---|
| 804 | end
|
|---|
| 805 | c----------------------------------------------------------------------
|
|---|
| 806 | subroutine hsmg_setup_fast1d_a(a,lbc,rbc,ll,lm,lr,ah,n)
|
|---|
| 807 | integer lbc,rbc,n
|
|---|
| 808 | real a(0:n+2,0:n+2),ll,lm,lr
|
|---|
| 809 | real ah(0:n,0:n)
|
|---|
| 810 |
|
|---|
| 811 | real fac
|
|---|
| 812 | integer i,j,i0,i1
|
|---|
| 813 | i0=0
|
|---|
| 814 | if(lbc.eq.1) i0=1
|
|---|
| 815 | i1=n
|
|---|
| 816 | if(rbc.eq.1) i1=n-1
|
|---|
| 817 |
|
|---|
| 818 | call rzero(a,(n+3)*(n+3))
|
|---|
| 819 | fac = 2.0/lm
|
|---|
| 820 | a(1,1)=1.0
|
|---|
| 821 | a(n+1,n+1)=1.0
|
|---|
| 822 | do j=i0,i1
|
|---|
| 823 | do i=i0,i1
|
|---|
| 824 | a(i+1,j+1)=fac*ah(i,j)
|
|---|
| 825 | enddo
|
|---|
| 826 | enddo
|
|---|
| 827 | if(lbc.eq.0) then
|
|---|
| 828 | fac = 2.0/ll
|
|---|
| 829 | a(0,0)=fac*ah(n-1,n-1)
|
|---|
| 830 | a(1,0)=fac*ah(n ,n-1)
|
|---|
| 831 | a(0,1)=fac*ah(n-1,n )
|
|---|
| 832 | a(1,1)=a(1,1)+fac*ah(n ,n )
|
|---|
| 833 | else
|
|---|
| 834 | a(0,0)=1.0
|
|---|
| 835 | endif
|
|---|
| 836 | if(rbc.eq.0) then
|
|---|
| 837 | fac = 2.0/lr
|
|---|
| 838 | a(n+1,n+1)=a(n+1,n+1)+fac*ah(0,0)
|
|---|
| 839 | a(n+2,n+1)=fac*ah(1,0)
|
|---|
| 840 | a(n+1,n+2)=fac*ah(0,1)
|
|---|
| 841 | a(n+2,n+2)=fac*ah(1,1)
|
|---|
| 842 | else
|
|---|
| 843 | a(n+2,n+2)=1.0
|
|---|
| 844 | endif
|
|---|
| 845 | return
|
|---|
| 846 | end
|
|---|
| 847 | c----------------------------------------------------------------------
|
|---|
| 848 | subroutine hsmg_setup_fast1d_b(b,lbc,rbc,ll,lm,lr,bh,n)
|
|---|
| 849 | integer lbc,rbc,n
|
|---|
| 850 | real b(0:n+2,0:n+2),ll,lm,lr
|
|---|
| 851 | real bh(0:n)
|
|---|
| 852 |
|
|---|
| 853 | real fac
|
|---|
| 854 | integer i,j,i0,i1
|
|---|
| 855 | i0=0
|
|---|
| 856 | if(lbc.eq.1) i0=1
|
|---|
| 857 | i1=n
|
|---|
| 858 | if(rbc.eq.1) i1=n-1
|
|---|
| 859 |
|
|---|
| 860 | call rzero(b,(n+3)*(n+3))
|
|---|
| 861 | fac = 0.5*lm
|
|---|
| 862 | b(1,1)=1.0
|
|---|
| 863 | b(n+1,n+1)=1.0
|
|---|
| 864 | do i=i0,i1
|
|---|
| 865 | b(i+1,i+1)=fac*bh(i)
|
|---|
| 866 | enddo
|
|---|
| 867 | if(lbc.eq.0) then
|
|---|
| 868 | fac = 0.5*ll
|
|---|
| 869 | b(0,0)=fac*bh(n-1)
|
|---|
| 870 | b(1,1)=b(1,1)+fac*bh(n )
|
|---|
| 871 | else
|
|---|
| 872 | b(0,0)=1.0
|
|---|
| 873 | endif
|
|---|
| 874 | if(rbc.eq.0) then
|
|---|
| 875 | fac = 0.5*lr
|
|---|
| 876 | b(n+1,n+1)=b(n+1,n+1)+fac*bh(0)
|
|---|
| 877 | b(n+2,n+2)=fac*bh(1)
|
|---|
| 878 | else
|
|---|
| 879 | b(n+2,n+2)=1.0
|
|---|
| 880 | endif
|
|---|
| 881 | return
|
|---|
| 882 | end
|
|---|
| 883 | c----------------------------------------------------------------------
|
|---|
| 884 | c clobbers r
|
|---|
| 885 | subroutine hsmg_fdm(e,r,l)
|
|---|
| 886 | include 'SIZE'
|
|---|
| 887 | include 'INPUT'
|
|---|
| 888 | include 'HSMG'
|
|---|
| 889 | call hsmg_do_fast(e,r,
|
|---|
| 890 | $ mg_fast_s(mg_fast_s_index(l,mg_fld)),
|
|---|
| 891 | $ mg_fast_d(mg_fast_d_index(l,mg_fld)),
|
|---|
| 892 | $ mg_nh(l)+2)
|
|---|
| 893 | return
|
|---|
| 894 | end
|
|---|
| 895 | c----------------------------------------------------------------------
|
|---|
| 896 | c clobbers r
|
|---|
| 897 | subroutine hsmg_do_fast(e,r,s,d,nl)
|
|---|
| 898 | include 'SIZE'
|
|---|
| 899 | include 'INPUT'
|
|---|
| 900 | real e(nl**ldim,nelv)
|
|---|
| 901 | real r(nl**ldim,nelv)
|
|---|
| 902 | real s(nl*nl,2,ldim,nelv)
|
|---|
| 903 | real d(nl**ldim,nelv)
|
|---|
| 904 |
|
|---|
| 905 | integer ie,nn,i
|
|---|
| 906 | nn=nl**ldim
|
|---|
| 907 | if(.not.if3d) then
|
|---|
| 908 | do ie=1,nelv
|
|---|
| 909 | call hsmg_tnsr2d_el(e(1,ie),nl,r(1,ie),nl
|
|---|
| 910 | $ ,s(1,2,1,ie),s(1,1,2,ie))
|
|---|
| 911 | do i=1,nn
|
|---|
| 912 | r(i,ie)=d(i,ie)*e(i,ie)
|
|---|
| 913 | enddo
|
|---|
| 914 | call hsmg_tnsr2d_el(e(1,ie),nl,r(1,ie),nl
|
|---|
| 915 | $ ,s(1,1,1,ie),s(1,2,2,ie))
|
|---|
| 916 | enddo
|
|---|
| 917 | else
|
|---|
| 918 | do ie=1,nelv
|
|---|
| 919 | call hsmg_tnsr3d_el(e(1,ie),nl,r(1,ie),nl
|
|---|
| 920 | $ ,s(1,2,1,ie),s(1,1,2,ie),s(1,1,3,ie))
|
|---|
| 921 | do i=1,nn
|
|---|
| 922 | r(i,ie)=d(i,ie)*e(i,ie)
|
|---|
| 923 | enddo
|
|---|
| 924 | call hsmg_tnsr3d_el(e(1,ie),nl,r(1,ie),nl
|
|---|
| 925 | $ ,s(1,1,1,ie),s(1,2,2,ie),s(1,2,3,ie))
|
|---|
| 926 | enddo
|
|---|
| 927 | endif
|
|---|
| 928 | return
|
|---|
| 929 | end
|
|---|
| 930 | c----------------------------------------------------------------------
|
|---|
| 931 | c u = wt .* u
|
|---|
| 932 | subroutine hsmg_do_wt(u,wt,nx,ny,nz)
|
|---|
| 933 | include 'SIZE'
|
|---|
| 934 | include 'INPUT'
|
|---|
| 935 | integer nx,ny,nz
|
|---|
| 936 | real u(nx,ny,nz,nelv)
|
|---|
| 937 | real wt(nx,nz,2,ldim,nelv)
|
|---|
| 938 |
|
|---|
| 939 | integer e
|
|---|
| 940 |
|
|---|
| 941 | c if (nx.eq.2) then
|
|---|
| 942 | c do e=1,nelv
|
|---|
| 943 | c call outmat(wt(1,1,1,1,e),nx,nz,'wt 1-1',e)
|
|---|
| 944 | c call outmat(wt(1,1,2,1,e),nx,nz,'wt 2-1',e)
|
|---|
| 945 | c call outmat(wt(1,1,1,2,e),nx,nz,'wt 1-2',e)
|
|---|
| 946 | c call outmat(wt(1,1,2,2,e),nx,nz,'wt 2-2',e)
|
|---|
| 947 | c enddo
|
|---|
| 948 | c call exitti('hsmg_do_wt quit$',nelv)
|
|---|
| 949 | c endif
|
|---|
| 950 |
|
|---|
| 951 | if (.not. if3d) then
|
|---|
| 952 | do ie=1,nelv
|
|---|
| 953 | do j=1,ny
|
|---|
| 954 | u( 1,j,1,ie)=u( 1,j,1,ie)*wt(j,1,1,1,ie)
|
|---|
| 955 | u(nx,j,1,ie)=u(nx,j,1,ie)*wt(j,1,2,1,ie)
|
|---|
| 956 | enddo
|
|---|
| 957 | do i=2,nx-1
|
|---|
| 958 | u(i, 1,1,ie)=u(i, 1,1,ie)*wt(i,1,1,2,ie)
|
|---|
| 959 | u(i,ny,1,ie)=u(i,ny,1,ie)*wt(i,1,2,2,ie)
|
|---|
| 960 | enddo
|
|---|
| 961 | enddo
|
|---|
| 962 | else
|
|---|
| 963 | do ie=1,nelv
|
|---|
| 964 | do k=1,nz
|
|---|
| 965 | do j=1,ny
|
|---|
| 966 | u( 1,j,k,ie)=u( 1,j,k,ie)*wt(j,k,1,1,ie)
|
|---|
| 967 | u(nx,j,k,ie)=u(nx,j,k,ie)*wt(j,k,2,1,ie)
|
|---|
| 968 | enddo
|
|---|
| 969 | enddo
|
|---|
| 970 | do k=1,nz
|
|---|
| 971 | do i=2,nx-1
|
|---|
| 972 | u(i, 1,k,ie)=u(i, 1,k,ie)*wt(i,k,1,2,ie)
|
|---|
| 973 | u(i,ny,k,ie)=u(i,ny,k,ie)*wt(i,k,2,2,ie)
|
|---|
| 974 | enddo
|
|---|
| 975 | enddo
|
|---|
| 976 | do j=2,ny-1
|
|---|
| 977 | do i=2,nx-1
|
|---|
| 978 | u(i,j, 1,ie)=u(i,j, 1,ie)*wt(i,j,1,3,ie)
|
|---|
| 979 | u(i,j,nz,ie)=u(i,j,nz,ie)*wt(i,j,2,3,ie)
|
|---|
| 980 | enddo
|
|---|
| 981 | enddo
|
|---|
| 982 | enddo
|
|---|
| 983 | endif
|
|---|
| 984 | return
|
|---|
| 985 | end
|
|---|
| 986 | c----------------------------------------------------------------------
|
|---|
| 987 | subroutine hsmg_setup_rstr_wt(wt,nx,ny,nz,l,w)
|
|---|
| 988 | include 'SIZE'
|
|---|
| 989 | include 'INPUT'
|
|---|
| 990 | integer nx,ny,nz,l
|
|---|
| 991 | real w(nx,ny,nz,nelv)
|
|---|
| 992 | real wt(nx,nz,2,ldim,nelv)
|
|---|
| 993 |
|
|---|
| 994 | integer ie
|
|---|
| 995 | !init border nodes to 1
|
|---|
| 996 | call rzero(w,nx*ny*nz*nelv)
|
|---|
| 997 | c print *, 'Setup rstr wt: ',nx,ny,nz,nelv
|
|---|
| 998 | if (.not.if3d) then
|
|---|
| 999 | do ie=1,nelv
|
|---|
| 1000 | do i=1,nx
|
|---|
| 1001 | w(i,1,1,ie)=1.0
|
|---|
| 1002 | w(i,ny,1,ie)=1.0
|
|---|
| 1003 | enddo
|
|---|
| 1004 | do j=1,ny
|
|---|
| 1005 | w(1,j,1,ie)=1.0
|
|---|
| 1006 | w(nx,j,1,ie)=1.0
|
|---|
| 1007 | enddo
|
|---|
| 1008 | enddo
|
|---|
| 1009 | else
|
|---|
| 1010 | do ie=1,nelv
|
|---|
| 1011 | do j=1,ny
|
|---|
| 1012 | do i=1,nx
|
|---|
| 1013 | w(i,j,1,ie)=1.0
|
|---|
| 1014 | w(i,j,nz,ie)=1.0
|
|---|
| 1015 | enddo
|
|---|
| 1016 | enddo
|
|---|
| 1017 | do k=1,nz
|
|---|
| 1018 | do i=1,nx
|
|---|
| 1019 | w(i,1,k,ie)=1.0
|
|---|
| 1020 | w(i,ny,k,ie)=1.0
|
|---|
| 1021 | enddo
|
|---|
| 1022 | enddo
|
|---|
| 1023 | do k=1,nz
|
|---|
| 1024 | do j=1,ny
|
|---|
| 1025 | w(1,j,k,ie)=1.0
|
|---|
| 1026 | w(nx,j,k,ie)=1.0
|
|---|
| 1027 | enddo
|
|---|
| 1028 | enddo
|
|---|
| 1029 | enddo
|
|---|
| 1030 | endif
|
|---|
| 1031 | call hsmg_dssum(w,l)
|
|---|
| 1032 | !invert the count w to get the weight wt
|
|---|
| 1033 | if (.not. if3d) then
|
|---|
| 1034 | do ie=1,nelv
|
|---|
| 1035 | do j=1,ny
|
|---|
| 1036 | wt(j,1,1,1,ie)=1.0/w(1,j,1,ie)
|
|---|
| 1037 | wt(j,1,2,1,ie)=1.0/w(nx,j,1,ie)
|
|---|
| 1038 | enddo
|
|---|
| 1039 | do i=1,nx
|
|---|
| 1040 | wt(i,1,1,2,ie)=1.0/w(i,1,1,ie)
|
|---|
| 1041 | wt(i,1,2,2,ie)=1.0/w(i,ny,1,ie)
|
|---|
| 1042 | enddo
|
|---|
| 1043 | enddo
|
|---|
| 1044 | else
|
|---|
| 1045 | do ie=1,nelv
|
|---|
| 1046 | do k=1,nz
|
|---|
| 1047 | do j=1,ny
|
|---|
| 1048 | wt(j,k,1,1,ie)=1.0/w(1,j,k,ie)
|
|---|
| 1049 | wt(j,k,2,1,ie)=1.0/w(nx,j,k,ie)
|
|---|
| 1050 | enddo
|
|---|
| 1051 | enddo
|
|---|
| 1052 | do k=1,nz
|
|---|
| 1053 | do i=1,nx
|
|---|
| 1054 | wt(i,k,1,2,ie)=1.0/w(i,1,k,ie)
|
|---|
| 1055 | wt(i,k,2,2,ie)=1.0/w(i,ny,k,ie)
|
|---|
| 1056 | enddo
|
|---|
| 1057 | enddo
|
|---|
| 1058 | do j=1,ny
|
|---|
| 1059 | do i=1,nx
|
|---|
| 1060 | wt(i,j,1,3,ie)=1.0/w(i,j,1,ie)
|
|---|
| 1061 | wt(i,j,2,3,ie)=1.0/w(i,j,nz,ie)
|
|---|
| 1062 | enddo
|
|---|
| 1063 | enddo
|
|---|
| 1064 | enddo
|
|---|
| 1065 | endif
|
|---|
| 1066 | return
|
|---|
| 1067 | end
|
|---|
| 1068 | c----------------------------------------------------------------------
|
|---|
| 1069 | subroutine hsmg_setup_mask(wt,nx,ny,nz,l,w)
|
|---|
| 1070 | include 'SIZE'
|
|---|
| 1071 | include 'INPUT'
|
|---|
| 1072 | integer nx,ny,nz,l
|
|---|
| 1073 | real w(nx,ny,nz,nelv)
|
|---|
| 1074 | real wt(nx,nz,2,ldim,nelv)
|
|---|
| 1075 |
|
|---|
| 1076 | integer ie
|
|---|
| 1077 | integer lbr,rbr,lbs,rbs,lbt,rbt,two
|
|---|
| 1078 | c init everything to 1
|
|---|
| 1079 |
|
|---|
| 1080 | n = nx*ny*nz*nelv
|
|---|
| 1081 | call rone(w,n)
|
|---|
| 1082 |
|
|---|
| 1083 | c set dirichlet nodes to zero
|
|---|
| 1084 | ierr = 0
|
|---|
| 1085 | two = 2
|
|---|
| 1086 | do ie=1,nelv
|
|---|
| 1087 | call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,ie,two,ierr)
|
|---|
| 1088 | if (ierr.ne.0) then
|
|---|
| 1089 | ierr = -1
|
|---|
| 1090 | call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,ie,two,ierr)
|
|---|
| 1091 | endif
|
|---|
| 1092 |
|
|---|
| 1093 | if(lbr.eq.1) then
|
|---|
| 1094 | do k=1,nz
|
|---|
| 1095 | do j=1,ny
|
|---|
| 1096 | w(1,j,k,ie)=0.0
|
|---|
| 1097 | enddo
|
|---|
| 1098 | enddo
|
|---|
| 1099 | endif
|
|---|
| 1100 | if(rbr.eq.1) then
|
|---|
| 1101 | do k=1,nz
|
|---|
| 1102 | do j=1,ny
|
|---|
| 1103 | w(nx,j,k,ie)=0.0
|
|---|
| 1104 | enddo
|
|---|
| 1105 | enddo
|
|---|
| 1106 | endif
|
|---|
| 1107 | if(lbs.eq.1) then
|
|---|
| 1108 | do k=1,nz
|
|---|
| 1109 | do i=1,nx
|
|---|
| 1110 | w(i,1,k,ie)=0.0
|
|---|
| 1111 | enddo
|
|---|
| 1112 | enddo
|
|---|
| 1113 | endif
|
|---|
| 1114 | if(rbs.eq.1) then
|
|---|
| 1115 | do k=1,nz
|
|---|
| 1116 | do i=1,nx
|
|---|
| 1117 | w(i,ny,k,ie)=0.0
|
|---|
| 1118 | enddo
|
|---|
| 1119 | enddo
|
|---|
| 1120 | endif
|
|---|
| 1121 | if(if3d) then
|
|---|
| 1122 | if(lbt.eq.1) then
|
|---|
| 1123 | do j=1,ny
|
|---|
| 1124 | do i=1,nx
|
|---|
| 1125 | w(i,j,1,ie)=0.0
|
|---|
| 1126 | enddo
|
|---|
| 1127 | enddo
|
|---|
| 1128 | endif
|
|---|
| 1129 | if(rbt.eq.1) then
|
|---|
| 1130 | do j=1,ny
|
|---|
| 1131 | do i=1,nx
|
|---|
| 1132 | w(i,j,nz,ie)=0.0
|
|---|
| 1133 | enddo
|
|---|
| 1134 | enddo
|
|---|
| 1135 | endif
|
|---|
| 1136 | endif
|
|---|
| 1137 | enddo
|
|---|
| 1138 | c do direct stiffness multiply
|
|---|
| 1139 |
|
|---|
| 1140 | call hsmg_dsprod(w,l)
|
|---|
| 1141 |
|
|---|
| 1142 |
|
|---|
| 1143 | c store weight
|
|---|
| 1144 | if (.not. if3d) then
|
|---|
| 1145 | do ie=1,nelv
|
|---|
| 1146 | do j=1,ny
|
|---|
| 1147 | wt(j,1,1,1,ie)=w(1,j,1,ie)
|
|---|
| 1148 | wt(j,1,2,1,ie)=w(nx,j,1,ie)
|
|---|
| 1149 | enddo
|
|---|
| 1150 | do i=1,nx
|
|---|
| 1151 | wt(i,1,1,2,ie)=w(i,1,1,ie)
|
|---|
| 1152 | wt(i,1,2,2,ie)=w(i,ny,1,ie)
|
|---|
| 1153 | enddo
|
|---|
| 1154 | enddo
|
|---|
| 1155 | else
|
|---|
| 1156 | do ie=1,nelv
|
|---|
| 1157 | do k=1,nz
|
|---|
| 1158 | do j=1,ny
|
|---|
| 1159 | wt(j,k,1,1,ie)=w(1,j,k,ie)
|
|---|
| 1160 | wt(j,k,2,1,ie)=w(nx,j,k,ie)
|
|---|
| 1161 | enddo
|
|---|
| 1162 | enddo
|
|---|
| 1163 | do k=1,nz
|
|---|
| 1164 | do i=1,nx
|
|---|
| 1165 | wt(i,k,1,2,ie)=w(i,1,k,ie)
|
|---|
| 1166 | wt(i,k,2,2,ie)=w(i,ny,k,ie)
|
|---|
| 1167 | enddo
|
|---|
| 1168 | enddo
|
|---|
| 1169 | do k=1,nz
|
|---|
| 1170 | do j=1,ny
|
|---|
| 1171 | wt(j,k,1,3,ie)=w(i,j,1,ie)
|
|---|
| 1172 | wt(j,k,2,3,ie)=w(i,j,nz,ie)
|
|---|
| 1173 | enddo
|
|---|
| 1174 | enddo
|
|---|
| 1175 | enddo
|
|---|
| 1176 | endif
|
|---|
| 1177 |
|
|---|
| 1178 | ierrmx = iglmax(ierr,1)
|
|---|
| 1179 | if (ierrmx.gt.0) then
|
|---|
| 1180 | if (ierr.gt.0) write(6,*) nid,ierr,' BC FAIL'
|
|---|
| 1181 | call exitti('B INVALID BC FOUND in genfast$',ierrmx)
|
|---|
| 1182 | endif
|
|---|
| 1183 |
|
|---|
| 1184 | return
|
|---|
| 1185 | end
|
|---|
| 1186 | c----------------------------------------------------------------------
|
|---|
| 1187 | subroutine hsmg_setup_schwarz_wt(ifsqrt)
|
|---|
| 1188 | logical ifsqrt
|
|---|
| 1189 | include 'SIZE'
|
|---|
| 1190 | include 'INPUT'
|
|---|
| 1191 | include 'HSMG'
|
|---|
| 1192 |
|
|---|
| 1193 | integer l,i,nl,nlz
|
|---|
| 1194 |
|
|---|
| 1195 | i = mg_schwarz_wt_index(mg_lmax,mg_fld-1)
|
|---|
| 1196 | do l=2,mg_lmax-1
|
|---|
| 1197 | mg_schwarz_wt_index(l,mg_fld)=i
|
|---|
| 1198 | nl = mg_nh(l)
|
|---|
| 1199 | nlz = mg_nh(l)
|
|---|
| 1200 | if(.not.if3d) nlz=1
|
|---|
| 1201 | i=i+nl*nlz*4*ldim*nelv
|
|---|
| 1202 | if(i .gt. lmg_swt*4*ldim*lelv) then
|
|---|
| 1203 | itmp = i/(4*ldim*lelv)
|
|---|
| 1204 | write(6,*) 'lmg_swt too small',i,itmp,lmg_swt,l
|
|---|
| 1205 | call exitt
|
|---|
| 1206 | endif
|
|---|
| 1207 |
|
|---|
| 1208 | call h1mg_setup_schwarz_wt_1(
|
|---|
| 1209 | $ mg_schwarz_wt(mg_schwarz_wt_index(l,mg_fld)),l,ifsqrt)
|
|---|
| 1210 |
|
|---|
| 1211 | enddo
|
|---|
| 1212 | mg_schwarz_wt_index(l,mg_fld)=i
|
|---|
| 1213 |
|
|---|
| 1214 | return
|
|---|
| 1215 | end
|
|---|
| 1216 | c----------------------------------------------------------------------
|
|---|
| 1217 | subroutine h1mg_setup_schwarz_wt(ifsqrt)
|
|---|
| 1218 | logical ifsqrt
|
|---|
| 1219 | include 'SIZE'
|
|---|
| 1220 | include 'INPUT'
|
|---|
| 1221 | include 'HSMG'
|
|---|
| 1222 |
|
|---|
| 1223 | integer l,i,nl,nlz
|
|---|
| 1224 |
|
|---|
| 1225 | i = mg_schwarz_wt_index(mg_lmax,mg_fld-1)
|
|---|
| 1226 | do l=2,mg_lmax
|
|---|
| 1227 |
|
|---|
| 1228 | mg_schwarz_wt_index(l,mg_fld)=i
|
|---|
| 1229 | nl = mg_nh(l)
|
|---|
| 1230 | nlz = mg_nhz(l)
|
|---|
| 1231 | i = i+nl*nlz*4*ldim*nelv
|
|---|
| 1232 |
|
|---|
| 1233 | if (i .gt. lmg_swt*4*ldim*lelv) then
|
|---|
| 1234 | itmp = i/(4*ldim*lelv)
|
|---|
| 1235 | write(6,*) 'lmg_swt too small',i,itmp,lmg_swt,l
|
|---|
| 1236 | call exitt
|
|---|
| 1237 | endif
|
|---|
| 1238 |
|
|---|
| 1239 | call h1mg_setup_schwarz_wt_1(
|
|---|
| 1240 | $ mg_schwarz_wt(mg_schwarz_wt_index(l,mg_fld)),l,ifsqrt)
|
|---|
| 1241 |
|
|---|
| 1242 | enddo
|
|---|
| 1243 |
|
|---|
| 1244 | mg_schwarz_wt_index(l,mg_fld)=i
|
|---|
| 1245 |
|
|---|
| 1246 | return
|
|---|
| 1247 | end
|
|---|
| 1248 | c----------------------------------------------------------------------
|
|---|
| 1249 | subroutine hsmg_schwarz_wt(e,l)
|
|---|
| 1250 | include 'SIZE'
|
|---|
| 1251 | include 'INPUT'
|
|---|
| 1252 | include 'HSMG'
|
|---|
| 1253 |
|
|---|
| 1254 | if(.not.if3d) call hsmg_schwarz_wt2d(
|
|---|
| 1255 | $ e,mg_schwarz_wt(mg_schwarz_wt_index(l,mg_fld)),mg_nh(l))
|
|---|
| 1256 | if(if3d) call hsmg_schwarz_wt3d(
|
|---|
| 1257 | $ e,mg_schwarz_wt(mg_schwarz_wt_index(l,mg_fld)),mg_nh(l))
|
|---|
| 1258 | return
|
|---|
| 1259 | end
|
|---|
| 1260 | c----------------------------------------------------------------------
|
|---|
| 1261 | subroutine hsmg_schwarz_wt2d(e,wt,n)
|
|---|
| 1262 | include 'SIZE'
|
|---|
| 1263 | integer n
|
|---|
| 1264 | real e(n,n,nelv)
|
|---|
| 1265 | real wt(n,4,2,nelv)
|
|---|
| 1266 |
|
|---|
| 1267 | integer ie,i,j
|
|---|
| 1268 | do ie=1,nelv
|
|---|
| 1269 | do j=1,n
|
|---|
| 1270 | e(1 ,j,ie)=e(1 ,j,ie)*wt(j,1,1,ie)
|
|---|
| 1271 | e(2 ,j,ie)=e(2 ,j,ie)*wt(j,2,1,ie)
|
|---|
| 1272 | e(n-1,j,ie)=e(n-1,j,ie)*wt(j,3,1,ie)
|
|---|
| 1273 | e(n ,j,ie)=e(n ,j,ie)*wt(j,4,1,ie)
|
|---|
| 1274 | enddo
|
|---|
| 1275 | do i=3,n-2
|
|---|
| 1276 | e(i,1 ,ie)=e(i,1 ,ie)*wt(i,1,2,ie)
|
|---|
| 1277 | e(i,2 ,ie)=e(i,2 ,ie)*wt(i,2,2,ie)
|
|---|
| 1278 | e(i,n-1,ie)=e(i,n-1,ie)*wt(i,3,2,ie)
|
|---|
| 1279 | e(i,n ,ie)=e(i,n ,ie)*wt(i,4,2,ie)
|
|---|
| 1280 | enddo
|
|---|
| 1281 | enddo
|
|---|
| 1282 | return
|
|---|
| 1283 | end
|
|---|
| 1284 | c----------------------------------------------------------------------
|
|---|
| 1285 | subroutine hsmg_schwarz_wt3d(e,wt,n)
|
|---|
| 1286 | include 'SIZE'
|
|---|
| 1287 | integer n
|
|---|
| 1288 | real e(n,n,n,nelv)
|
|---|
| 1289 | real wt(n,n,4,3,nelv)
|
|---|
| 1290 |
|
|---|
| 1291 | integer ie,i,j,k
|
|---|
| 1292 | do ie=1,nelv
|
|---|
| 1293 | do k=1,n
|
|---|
| 1294 | do j=1,n
|
|---|
| 1295 | e(1 ,j,k,ie)=e(1 ,j,k,ie)*wt(j,k,1,1,ie)
|
|---|
| 1296 | e(2 ,j,k,ie)=e(2 ,j,k,ie)*wt(j,k,2,1,ie)
|
|---|
| 1297 | e(n-1,j,k,ie)=e(n-1,j,k,ie)*wt(j,k,3,1,ie)
|
|---|
| 1298 | e(n ,j,k,ie)=e(n ,j,k,ie)*wt(j,k,4,1,ie)
|
|---|
| 1299 | enddo
|
|---|
| 1300 | enddo
|
|---|
| 1301 | do k=1,n
|
|---|
| 1302 | do i=3,n-2
|
|---|
| 1303 | e(i,1 ,k,ie)=e(i,1 ,k,ie)*wt(i,k,1,2,ie)
|
|---|
| 1304 | e(i,2 ,k,ie)=e(i,2 ,k,ie)*wt(i,k,2,2,ie)
|
|---|
| 1305 | e(i,n-1,k,ie)=e(i,n-1,k,ie)*wt(i,k,3,2,ie)
|
|---|
| 1306 | e(i,n ,k,ie)=e(i,n ,k,ie)*wt(i,k,4,2,ie)
|
|---|
| 1307 | enddo
|
|---|
| 1308 | enddo
|
|---|
| 1309 | do j=3,n-2
|
|---|
| 1310 | do i=3,n-2
|
|---|
| 1311 | e(i,j,1 ,ie)=e(i,j,1 ,ie)*wt(i,j,1,3,ie)
|
|---|
| 1312 | e(i,j,2 ,ie)=e(i,j,2 ,ie)*wt(i,j,2,3,ie)
|
|---|
| 1313 | e(i,j,n-1,ie)=e(i,j,n-1,ie)*wt(i,j,3,3,ie)
|
|---|
| 1314 | e(i,j,n ,ie)=e(i,j,n ,ie)*wt(i,j,4,3,ie)
|
|---|
| 1315 | enddo
|
|---|
| 1316 | enddo
|
|---|
| 1317 | enddo
|
|---|
| 1318 | return
|
|---|
| 1319 | end
|
|---|
| 1320 | c----------------------------------------------------------------------
|
|---|
| 1321 | subroutine hsmg_coarse_solve(e,r)
|
|---|
| 1322 | include 'SIZE'
|
|---|
| 1323 | include 'DOMAIN'
|
|---|
| 1324 | include 'ESOLV'
|
|---|
| 1325 | include 'GEOM'
|
|---|
| 1326 | include 'SOLN'
|
|---|
| 1327 | include 'PARALLEL'
|
|---|
| 1328 | include 'HSMG'
|
|---|
| 1329 | include 'CTIMER'
|
|---|
| 1330 | include 'INPUT'
|
|---|
| 1331 | include 'TSTEP'
|
|---|
| 1332 | real e(1),r(1)
|
|---|
| 1333 | c
|
|---|
| 1334 | integer n_crs_tot
|
|---|
| 1335 | save n_crs_tot
|
|---|
| 1336 | data n_crs_tot /0/
|
|---|
| 1337 | c
|
|---|
| 1338 | if (icalld.eq.0) then ! timer info
|
|---|
| 1339 | ncrsl=0
|
|---|
| 1340 | tcrsl=0.0
|
|---|
| 1341 | endif
|
|---|
| 1342 | icalld = 1
|
|---|
| 1343 |
|
|---|
| 1344 | if (ifsync) call nekgsync()
|
|---|
| 1345 |
|
|---|
| 1346 | ncrsl = ncrsl + 1
|
|---|
| 1347 | etime1=dnekclock()
|
|---|
| 1348 |
|
|---|
| 1349 | call fgslib_crs_solve(xxth(ifield),e,r)
|
|---|
| 1350 |
|
|---|
| 1351 | tcrsl=tcrsl+dnekclock()-etime1
|
|---|
| 1352 |
|
|---|
| 1353 | return
|
|---|
| 1354 | end
|
|---|
| 1355 | c----------------------------------------------------------------------
|
|---|
| 1356 | subroutine hsmg_setup_solve
|
|---|
| 1357 | include 'SIZE'
|
|---|
| 1358 | include 'HSMG'
|
|---|
| 1359 |
|
|---|
| 1360 | integer l,i,nl,nlz
|
|---|
| 1361 | i = mg_solve_index(mg_lmax+1,mg_fld-1)
|
|---|
| 1362 | do l=1,mg_lmax
|
|---|
| 1363 | mg_solve_index(l,mg_fld)=i
|
|---|
| 1364 | i=i+mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
|
|---|
| 1365 | if(i .gt. lmg_solve*lelv) then
|
|---|
| 1366 | itmp = i/lelv
|
|---|
| 1367 | write(6,*) 'lmg_solve too small',i,itmp,lmg_solve,l
|
|---|
| 1368 | call exitt
|
|---|
| 1369 | endif
|
|---|
| 1370 | enddo
|
|---|
| 1371 | mg_solve_index(l,mg_fld)=i
|
|---|
| 1372 |
|
|---|
| 1373 | return
|
|---|
| 1374 | end
|
|---|
| 1375 | c----------------------------------------------------------------------
|
|---|
| 1376 | subroutine hsmg_solve(e,r)
|
|---|
| 1377 | real e(1),r(1)
|
|---|
| 1378 | include 'SIZE'
|
|---|
| 1379 | include 'HSMG'
|
|---|
| 1380 | include 'GEOM'
|
|---|
| 1381 | include 'INPUT'
|
|---|
| 1382 | include 'MASS'
|
|---|
| 1383 | include 'SOLN'
|
|---|
| 1384 | include 'TSTEP'
|
|---|
| 1385 | include 'CTIMER'
|
|---|
| 1386 | include 'PARALLEL'
|
|---|
| 1387 |
|
|---|
| 1388 | common /quick/ ecrs (2) ! quick work array
|
|---|
| 1389 | $ , ecrs2 (2) ! quick work array
|
|---|
| 1390 | c common /quick/ ecrs (lx2*ly2*lz2*lelv) ! quick work array
|
|---|
| 1391 | c $ , ecrs2 (lx2*ly2*lz2*lelv) ! quick work array
|
|---|
| 1392 |
|
|---|
| 1393 | common /scrhi/ h2inv (lx1,ly1,lz1,lelv)
|
|---|
| 1394 | common /scrvh/ h1 (lx1,ly1,lz1,lelv),
|
|---|
| 1395 | $ h2 (lx1,ly1,lz1,lelv)
|
|---|
| 1396 |
|
|---|
| 1397 | integer ilstep,iter
|
|---|
| 1398 | save ilstep,iter
|
|---|
| 1399 | data ilstep,iter /0,0/
|
|---|
| 1400 |
|
|---|
| 1401 | real rhoavg,copt(2),copw(2)
|
|---|
| 1402 | save rhoavg,copt1,copt2
|
|---|
| 1403 | data rhoavg,copt1,copt2 /3*1./ ! Default copt = 1 for additive
|
|---|
| 1404 |
|
|---|
| 1405 | integer l,nt
|
|---|
| 1406 | integer*8 ntotg,nxyz2
|
|---|
| 1407 |
|
|---|
| 1408 | logical if_hybrid
|
|---|
| 1409 |
|
|---|
| 1410 | mg_fld = 1
|
|---|
| 1411 | if (ifield.gt.1) mg_fld = 2
|
|---|
| 1412 |
|
|---|
| 1413 | if (istep.ne.ilstep) then
|
|---|
| 1414 | ilstep = istep
|
|---|
| 1415 | ntot1 = lx1*ly1*lz1*nelv
|
|---|
| 1416 | rhoavg = glsc2(vtrans,bm1,ntot1)/volvm1
|
|---|
| 1417 | endif
|
|---|
| 1418 |
|
|---|
| 1419 | n = lx2*ly2*lz2*nelv
|
|---|
| 1420 | c call copy(e,r,n)
|
|---|
| 1421 | c return
|
|---|
| 1422 |
|
|---|
| 1423 | if (icalld.eq.0) then
|
|---|
| 1424 |
|
|---|
| 1425 | tddsl=0.0
|
|---|
| 1426 | nddsl=0
|
|---|
| 1427 |
|
|---|
| 1428 | icalld = 1
|
|---|
| 1429 | taaaa = 0
|
|---|
| 1430 | tbbbb = 0
|
|---|
| 1431 | tcccc = 0
|
|---|
| 1432 | tdddd = 0
|
|---|
| 1433 | teeee = 0
|
|---|
| 1434 | endif
|
|---|
| 1435 |
|
|---|
| 1436 | nddsl = nddsl + 1
|
|---|
| 1437 | etime1 = dnekclock()
|
|---|
| 1438 |
|
|---|
| 1439 |
|
|---|
| 1440 | c n = lx2*ly2*lz2*nelv
|
|---|
| 1441 | c rmax = glmax(r,n)
|
|---|
| 1442 | c if (nid.eq.0) write(6,*) istep,n,rmax,' rmax1'
|
|---|
| 1443 |
|
|---|
| 1444 | iter = iter + 1
|
|---|
| 1445 |
|
|---|
| 1446 | l = mg_lmax
|
|---|
| 1447 | nt = mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
|
|---|
| 1448 | ! e := W M r
|
|---|
| 1449 | ! Schwarz
|
|---|
| 1450 | time_0 = dnekclock()
|
|---|
| 1451 | call local_solves_fdm(e,r)
|
|---|
| 1452 |
|
|---|
| 1453 | time_1 = dnekclock()
|
|---|
| 1454 |
|
|---|
| 1455 | c if (param(41).eq.1) if_hybrid = .true.
|
|---|
| 1456 | if_hybrid = .false.
|
|---|
| 1457 |
|
|---|
| 1458 | if (if_hybrid) then
|
|---|
| 1459 | ! w := E e
|
|---|
| 1460 | rbd1dt = rhoavg*bd(1)/dt ! Assumes constant density!!!
|
|---|
| 1461 | call cdabdtp(mg_work2,e,h1,h2,h2inv,1)
|
|---|
| 1462 | call cmult (mg_work2,rbd1dt,nt)
|
|---|
| 1463 | time_2 = dnekclock()
|
|---|
| 1464 | if (istep.eq.1) then
|
|---|
| 1465 | copt(1) = vlsc2(r ,mg_work2,nt)
|
|---|
| 1466 | copt(2) = vlsc2(mg_work2,mg_work2,nt)
|
|---|
| 1467 | call gop(copt,copw,'+ ', 2)
|
|---|
| 1468 | copt(1) = copt(1)/copt(2)
|
|---|
| 1469 | avg2 = 1./iter
|
|---|
| 1470 | avg1 = 1.-avg2
|
|---|
| 1471 | copt1 = avg1*copt1 + avg2*copt(1)
|
|---|
| 1472 | if(nio.eq.0)write(6,1)istep,iter,rbd1dt,copt(1),copt1,'cpt1'
|
|---|
| 1473 | 1 format(2i6,1p3e14.5,2x,a4)
|
|---|
| 1474 | endif
|
|---|
| 1475 | ! w := r - w
|
|---|
| 1476 | do i = 1,nt
|
|---|
| 1477 | mg_work2(i) = r(i) - copt1*mg_work2(i)
|
|---|
| 1478 | e (i) = copt1*e(i)
|
|---|
| 1479 | ecrs2 (i) = mg_work2(i)
|
|---|
| 1480 | enddo
|
|---|
| 1481 |
|
|---|
| 1482 | else ! Additive
|
|---|
| 1483 | ! w := r - w
|
|---|
| 1484 | do i = 1,nt
|
|---|
| 1485 | mg_work2(i) = r(i)
|
|---|
| 1486 | enddo
|
|---|
| 1487 | time_2 = dnekclock()
|
|---|
| 1488 | endif
|
|---|
| 1489 |
|
|---|
| 1490 | do l = mg_lmax-1,2,-1
|
|---|
| 1491 |
|
|---|
| 1492 | c rmax = glmax(mg_work2,nt)
|
|---|
| 1493 | c if (nid.eq.0) write(6,*) l,nt,rmax,' rmax2'
|
|---|
| 1494 |
|
|---|
| 1495 | nt = mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
|
|---|
| 1496 | ! T
|
|---|
| 1497 | ! r := J w
|
|---|
| 1498 | ! l
|
|---|
| 1499 | call hsmg_rstr(mg_solve_r(mg_solve_index(l,mg_fld)),mg_work2,l)
|
|---|
| 1500 |
|
|---|
| 1501 | ! w := r
|
|---|
| 1502 | ! l
|
|---|
| 1503 | call copy(mg_work2,mg_solve_r(mg_solve_index(l,mg_fld)),nt)
|
|---|
| 1504 | ! e := M w
|
|---|
| 1505 | ! l Schwarz
|
|---|
| 1506 | call hsmg_schwarz(
|
|---|
| 1507 | $ mg_solve_e(mg_solve_index(l,mg_fld)),mg_work2,l)
|
|---|
| 1508 |
|
|---|
| 1509 | ! e := W e
|
|---|
| 1510 | ! l l
|
|---|
| 1511 | call hsmg_schwarz_wt(mg_solve_e(mg_solve_index(l,mg_fld)),l)
|
|---|
| 1512 |
|
|---|
| 1513 | c call exitti('quit in mg$',l)
|
|---|
| 1514 |
|
|---|
| 1515 | ! w := r - w
|
|---|
| 1516 | ! l
|
|---|
| 1517 | do i = 0,nt-1
|
|---|
| 1518 | mg_work2(i+1) = mg_solve_r(mg_solve_index(l,mg_fld)+i)
|
|---|
| 1519 | $ !-alpha*mg_work2(i+1)
|
|---|
| 1520 | enddo
|
|---|
| 1521 | enddo
|
|---|
| 1522 |
|
|---|
| 1523 | call hsmg_rstr_no_dssum(
|
|---|
| 1524 | $ mg_solve_r(mg_solve_index(1,mg_fld)),mg_work2,1)
|
|---|
| 1525 |
|
|---|
| 1526 | nzw = ldim-1
|
|---|
| 1527 |
|
|---|
| 1528 | call hsmg_do_wt(mg_solve_r(mg_solve_index(1,mg_fld)),
|
|---|
| 1529 | $ mg_mask(mg_mask_index(1,mg_fld)),2,2,nzw)
|
|---|
| 1530 |
|
|---|
| 1531 | ! -1
|
|---|
| 1532 | ! e := A r
|
|---|
| 1533 | ! 1 1
|
|---|
| 1534 | call hsmg_coarse_solve(mg_solve_e(mg_solve_index(1,mg_fld)),
|
|---|
| 1535 | $ mg_solve_r(mg_solve_index(1,mg_fld)))
|
|---|
| 1536 |
|
|---|
| 1537 | call hsmg_do_wt(mg_solve_e(mg_solve_index(1,mg_fld)),
|
|---|
| 1538 | $ mg_mask(mg_mask_index(1,mg_fld)),2,2,nzw)
|
|---|
| 1539 | time_3 = dnekclock()
|
|---|
| 1540 | do l = 2,mg_lmax-1
|
|---|
| 1541 | nt = mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
|
|---|
| 1542 | ! w := J e
|
|---|
| 1543 | ! l-1
|
|---|
| 1544 | call hsmg_intp
|
|---|
| 1545 | $ (mg_work2,mg_solve_e(mg_solve_index(l-1,mg_fld)),l-1)
|
|---|
| 1546 |
|
|---|
| 1547 | ! e := e + w
|
|---|
| 1548 | ! l l
|
|---|
| 1549 | do i = 0,nt-1
|
|---|
| 1550 | mg_solve_e(mg_solve_index(l,mg_fld)+i) =
|
|---|
| 1551 | $ + mg_solve_e(mg_solve_index(l,mg_fld)+i) + mg_work2(i+1)
|
|---|
| 1552 | enddo
|
|---|
| 1553 | enddo
|
|---|
| 1554 | l = mg_lmax
|
|---|
| 1555 | nt = mg_nh(l)*mg_nh(l)*mg_nhz(l)*nelv
|
|---|
| 1556 | ! w := J e
|
|---|
| 1557 | ! m-1
|
|---|
| 1558 |
|
|---|
| 1559 | call hsmg_intp(mg_work2,
|
|---|
| 1560 | $ mg_solve_e(mg_solve_index(l-1,mg_fld)),l-1)
|
|---|
| 1561 |
|
|---|
| 1562 | if (if_hybrid.and.istep.eq.1) then
|
|---|
| 1563 | ! ecrs := E e_c
|
|---|
| 1564 | call cdabdtp(ecrs,mg_work2,h1,h2,h2inv,1)
|
|---|
| 1565 | call cmult (ecrs,rbd1dt,nt)
|
|---|
| 1566 | copt(1) = vlsc2(ecrs2,ecrs,nt)
|
|---|
| 1567 | copt(2) = vlsc2(ecrs ,ecrs,nt)
|
|---|
| 1568 | call gop(copt,copw,'+ ', 2)
|
|---|
| 1569 | copt(1) = copt(1)/copt(2)
|
|---|
| 1570 | avg2 = 1./iter
|
|---|
| 1571 | avg1 = 1.-avg2
|
|---|
| 1572 | copt2 = avg1*copt2 + avg2*copt(1)
|
|---|
| 1573 | if(nio.eq.0)write(6,1)istep,iter,rbd1dt,copt(1),copt2,'cpt2'
|
|---|
| 1574 | endif
|
|---|
| 1575 | ! e := e + w
|
|---|
| 1576 |
|
|---|
| 1577 | do i = 1,nt
|
|---|
| 1578 | e(i) = e(i) + copt2*mg_work2(i)
|
|---|
| 1579 | enddo
|
|---|
| 1580 | time_4 = dnekclock()
|
|---|
| 1581 | c print *, 'Did an MG iteration'
|
|---|
| 1582 | c
|
|---|
| 1583 | taaaa = taaaa + (time_1 - time_0)
|
|---|
| 1584 | tbbbb = tbbbb + (time_2 - time_1)
|
|---|
| 1585 | tcccc = tcccc + (time_3 - time_2)
|
|---|
| 1586 | tdddd = tdddd + (time_4 - time_3)
|
|---|
| 1587 | teeee = teeee + (time_4 - time_0)
|
|---|
| 1588 | c
|
|---|
| 1589 | c A typical time breakdown:
|
|---|
| 1590 | c
|
|---|
| 1591 | c 1.3540E+01 5.4390E+01 1.1440E+01 1.2199E+00 8.0590E+01 HSMG time
|
|---|
| 1592 | c
|
|---|
| 1593 | c ==> 54/80 = 67 % of preconditioner time is in residual evaluation!
|
|---|
| 1594 | c
|
|---|
| 1595 | call ortho (e)
|
|---|
| 1596 |
|
|---|
| 1597 | tddsl = tddsl + ( dnekclock()-etime1 )
|
|---|
| 1598 |
|
|---|
| 1599 |
|
|---|
| 1600 |
|
|---|
| 1601 | return
|
|---|
| 1602 | end
|
|---|
| 1603 | c----------------------------------------------------------------------
|
|---|
| 1604 | subroutine hsmg_setup_mg_nx()
|
|---|
| 1605 | include 'SIZE'
|
|---|
| 1606 | include 'INPUT'
|
|---|
| 1607 | include 'PARALLEL'
|
|---|
| 1608 | include 'HSMG'
|
|---|
| 1609 | include 'SEMHAT'
|
|---|
| 1610 | real w(lx1+2)
|
|---|
| 1611 | integer nf,nc,nr
|
|---|
| 1612 | integer nx,ny,nz
|
|---|
| 1613 |
|
|---|
| 1614 | integer mgn2(10)
|
|---|
| 1615 | save mgn2
|
|---|
| 1616 | data mgn2 / 1, 2, 2, 2, 2, 3, 3, 5, 5, 5/
|
|---|
| 1617 | c data mgn2 / 1, 2, 3, 4, 5, 6, 7, 8, 9, 0
|
|---|
| 1618 |
|
|---|
| 1619 | c if (param(82).eq.0) param(82)=2 ! nek default
|
|---|
| 1620 | c if (np.eq.1) param(82)=2 ! single proc. too slow
|
|---|
| 1621 | p82 = 2 ! potentially variable nxc
|
|---|
| 1622 | mg_lmax = 3
|
|---|
| 1623 | c mg_lmax = 4
|
|---|
| 1624 | if (lx1.eq.4) mg_lmax = 2
|
|---|
| 1625 | c if (param(79).ne.0) mg_lmax = param(79)
|
|---|
| 1626 | mglx1 = p82-1 !1
|
|---|
| 1627 | mg_nx(1) = mglx1
|
|---|
| 1628 | mg_ny(1) = mglx1
|
|---|
| 1629 | mg_nz(1) = mglx1
|
|---|
| 1630 | if (.not.if3d) mg_nz(1) = 0
|
|---|
| 1631 |
|
|---|
| 1632 | mglx2 = 2*(lx2/4) + 1
|
|---|
| 1633 | if (lx1.eq.5) mglx2 = 3
|
|---|
| 1634 | c if (lx1.eq.6) mglx2 = 3
|
|---|
| 1635 | if (lx1.le.10) mglx2 = mgn2(lx1)
|
|---|
| 1636 | if (lx1.eq.8) mglx2 = 4
|
|---|
| 1637 | if (lx1.eq.8) mglx2 = 3
|
|---|
| 1638 |
|
|---|
| 1639 | c mglx2 = min(3,mglx2)
|
|---|
| 1640 |
|
|---|
| 1641 |
|
|---|
| 1642 | mg_nx(2) = mglx2
|
|---|
| 1643 | mg_ny(2) = mglx2
|
|---|
| 1644 | mg_nz(2) = mglx2
|
|---|
| 1645 | if (.not.if3d) mg_nz(2) = 0
|
|---|
| 1646 |
|
|---|
| 1647 | mg_nx(3) = mglx2+1
|
|---|
| 1648 | mg_ny(3) = mglx2+1
|
|---|
| 1649 | mg_nz(3) = mglx2+1
|
|---|
| 1650 | if (.not.if3d) mg_nz(3) = 0
|
|---|
| 1651 |
|
|---|
| 1652 | mg_nx(mg_lmax) = lx1-1
|
|---|
| 1653 | mg_ny(mg_lmax) = ly1-1
|
|---|
| 1654 | mg_nz(mg_lmax) = lz1-1
|
|---|
| 1655 |
|
|---|
| 1656 | if (nio.eq.0) write(*,*) 'mg_nx:',(mg_nx(i),i=1,mg_lmax)
|
|---|
| 1657 | if (nio.eq.0) write(*,*) 'mg_ny:',(mg_ny(i),i=1,mg_lmax)
|
|---|
| 1658 | if (nio.eq.0) write(*,*) 'mg_nz:',(mg_nz(i),i=1,mg_lmax)
|
|---|
| 1659 |
|
|---|
| 1660 | return
|
|---|
| 1661 | end
|
|---|
| 1662 | c----------------------------------------------------------------------
|
|---|
| 1663 | subroutine hsmg_index_0 ! initialize index sets
|
|---|
| 1664 | include 'SIZE'
|
|---|
| 1665 | include 'HSMG'
|
|---|
| 1666 |
|
|---|
| 1667 | n = lmgn*(lmgs+1)
|
|---|
| 1668 |
|
|---|
| 1669 | call izero( mg_rstr_wt_index , n )
|
|---|
| 1670 | call izero( mg_mask_index , n )
|
|---|
| 1671 | call izero( mg_solve_index , n )
|
|---|
| 1672 | call izero( mg_fast_s_index , n )
|
|---|
| 1673 | call izero( mg_fast_d_index , n )
|
|---|
| 1674 | call izero( mg_schwarz_wt_index , n )
|
|---|
| 1675 |
|
|---|
| 1676 | return
|
|---|
| 1677 | end
|
|---|
| 1678 | c----------------------------------------------------------------------
|
|---|
| 1679 | subroutine outfldn (x,n,txt10,ichk) ! writes into unit=40+ifiled
|
|---|
| 1680 | INCLUDE 'SIZE'
|
|---|
| 1681 | INCLUDE 'TSTEP'
|
|---|
| 1682 | real x(n,n,1,lelt)
|
|---|
| 1683 | character*10 txt10
|
|---|
| 1684 | c
|
|---|
| 1685 | integer idum,e
|
|---|
| 1686 | save idum
|
|---|
| 1687 | data idum /3/
|
|---|
| 1688 | if (idum.lt.0) return
|
|---|
| 1689 | m = 40 + ifield ! unit #
|
|---|
| 1690 | c
|
|---|
| 1691 | C
|
|---|
| 1692 | mtot = n*n*nelv
|
|---|
| 1693 | if (n.gt.7.or.nelv.gt.16) return
|
|---|
| 1694 | xmin = glmin(x,mtot)
|
|---|
| 1695 | xmax = glmax(x,mtot)
|
|---|
| 1696 | c
|
|---|
| 1697 | rnel = nelv
|
|---|
| 1698 | snel = sqrt(rnel)+.1
|
|---|
| 1699 | ne = snel
|
|---|
| 1700 | ne1 = nelv-ne+1
|
|---|
| 1701 | do ie=ne1,1,-ne
|
|---|
| 1702 | l=ie-1
|
|---|
| 1703 | do k=1,1
|
|---|
| 1704 | if (ie.eq.ne1) write(m,116) txt10,k,ie,xmin,xmax,ichk,time
|
|---|
| 1705 | write(m,117)
|
|---|
| 1706 | do j=n,1,-1
|
|---|
| 1707 | if (n.eq.2) write(m,102) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1708 | if (n.eq.3) write(m,103) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1709 | if (n.eq.4) write(m,104) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1710 | if (n.eq.5) write(m,105) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1711 | if (n.eq.6) write(m,106) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1712 | if (n.eq.7) write(m,107) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1713 | if (n.eq.8) write(m,108) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1714 | enddo
|
|---|
| 1715 | enddo
|
|---|
| 1716 | enddo
|
|---|
| 1717 |
|
|---|
| 1718 | C
|
|---|
| 1719 | 102 FORMAT(4(2f9.5,2x))
|
|---|
| 1720 | 103 FORMAT(4(3f9.5,2x))
|
|---|
| 1721 | 104 FORMAT(4(4f7.3,2x))
|
|---|
| 1722 | 105 FORMAT(5f9.5,10x,5f9.5)
|
|---|
| 1723 | 106 FORMAT(6f9.5,5x,6f9.5)
|
|---|
| 1724 | 107 FORMAT(7f8.4,5x,7f8.4)
|
|---|
| 1725 | 108 FORMAT(8f8.4,4x,8f8.4)
|
|---|
| 1726 | c
|
|---|
| 1727 | 116 FORMAT( /,5X,' ^ ',/,
|
|---|
| 1728 | $ 5X,' Y | ',/,
|
|---|
| 1729 | $ 5X,' | ',A10,/,
|
|---|
| 1730 | $ 5X,' +----> ','Plane = ',I2,'/',I2,2x,2e12.4,/,
|
|---|
| 1731 | $ 5X,' X ','Step =',I9,f15.5)
|
|---|
| 1732 | 117 FORMAT(' ')
|
|---|
| 1733 | c
|
|---|
| 1734 | c if (ichk.eq.1.and.idum.gt.0) call checkit(idum)
|
|---|
| 1735 | return
|
|---|
| 1736 | end
|
|---|
| 1737 | c-----------------------------------------------------------------------
|
|---|
| 1738 | subroutine outfldn0 (x,n,txt10,ichk)
|
|---|
| 1739 | INCLUDE 'SIZE'
|
|---|
| 1740 | INCLUDE 'TSTEP'
|
|---|
| 1741 | real x(n,n,1,lelt)
|
|---|
| 1742 | character*10 txt10
|
|---|
| 1743 | c
|
|---|
| 1744 | integer idum,e
|
|---|
| 1745 | save idum
|
|---|
| 1746 | data idum /3/
|
|---|
| 1747 | if (idum.lt.0) return
|
|---|
| 1748 | c
|
|---|
| 1749 | C
|
|---|
| 1750 | mtot = n*n*nelv
|
|---|
| 1751 | if (n.gt.7.or.nelv.gt.16) return
|
|---|
| 1752 | xmin = glmin(x,mtot)
|
|---|
| 1753 | xmax = glmax(x,mtot)
|
|---|
| 1754 | c
|
|---|
| 1755 | rnel = nelv
|
|---|
| 1756 | snel = sqrt(rnel)+.1
|
|---|
| 1757 | ne = snel
|
|---|
| 1758 | ne1 = nelv-ne+1
|
|---|
| 1759 | do ie=ne1,1,-ne
|
|---|
| 1760 | l=ie-1
|
|---|
| 1761 | do k=1,1
|
|---|
| 1762 | if (ie.eq.ne1) write(6,116) txt10,k,ie,xmin,xmax,ichk,time
|
|---|
| 1763 | write(6,117)
|
|---|
| 1764 | do j=n,1,-1
|
|---|
| 1765 | if (n.eq.2) write(6,102) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1766 | if (n.eq.3) write(6,103) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1767 | if (n.eq.4) write(6,104) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1768 | if (n.eq.5) write(6,105) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1769 | if (n.eq.6) write(6,106) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1770 | if (n.eq.7) write(6,107) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1771 | if (n.eq.8) write(6,108) ((x(i,j,k,e+l),i=1,n),e=1,ne)
|
|---|
| 1772 | enddo
|
|---|
| 1773 | enddo
|
|---|
| 1774 | enddo
|
|---|
| 1775 |
|
|---|
| 1776 | C
|
|---|
| 1777 | 102 FORMAT(4(2f9.5,2x))
|
|---|
| 1778 | 103 FORMAT(4(3f9.5,2x))
|
|---|
| 1779 | 104 FORMAT(4(4f7.3,2x))
|
|---|
| 1780 | 105 FORMAT(5f9.5,10x,5f9.5)
|
|---|
| 1781 | 106 FORMAT(6f9.5,5x,6f9.5)
|
|---|
| 1782 | 107 FORMAT(7f8.4,5x,7f8.4)
|
|---|
| 1783 | 108 FORMAT(8f8.4,4x,8f8.4)
|
|---|
| 1784 | c
|
|---|
| 1785 | 116 FORMAT( /,5X,' ^ ',/,
|
|---|
| 1786 | $ 5X,' Y | ',/,
|
|---|
| 1787 | $ 5X,' | ',A10,/,
|
|---|
| 1788 | $ 5X,' +----> ','Plane = ',I2,'/',I2,2x,2e12.4,/,
|
|---|
| 1789 | $ 5X,' X ','Step =',I9,f15.5)
|
|---|
| 1790 | 117 FORMAT(' ')
|
|---|
| 1791 | c
|
|---|
| 1792 | c if (ichk.eq.1.and.idum.gt.0) call checkit(idum)
|
|---|
| 1793 | return
|
|---|
| 1794 | end
|
|---|
| 1795 | c-----------------------------------------------------------------------
|
|---|
| 1796 | subroutine outflda (x,n,txt10,ichk) ! writes into unit=p130+ifiled
|
|---|
| 1797 | INCLUDE 'SIZE' ! or into std. output for p130<9
|
|---|
| 1798 | INCLUDE 'TSTEP' ! truncated below eps=p131
|
|---|
| 1799 | INCLUDE 'INPUT' ! param(130)
|
|---|
| 1800 | real x(1)
|
|---|
| 1801 | character*10 txt10 ! note: n is not used
|
|---|
| 1802 | c parameter (eps=1.e-7)
|
|---|
| 1803 | C
|
|---|
| 1804 | p130 = param(130)
|
|---|
| 1805 | eps = param(131)
|
|---|
| 1806 | if (p130.le.0) return
|
|---|
| 1807 | m = 6
|
|---|
| 1808 | if (p130.gt.9) m = p130 + ifield
|
|---|
| 1809 |
|
|---|
| 1810 | ntot = lx1*ly1*lz1*nelfld(ifield)
|
|---|
| 1811 |
|
|---|
| 1812 | xmin = glmin(x,ntot)
|
|---|
| 1813 | xmax = glmax(x,ntot)
|
|---|
| 1814 | xavg = glsum(x,ntot)/ntot
|
|---|
| 1815 |
|
|---|
| 1816 | if (abs(xavg).lt.eps) xavg = 0. ! truncation
|
|---|
| 1817 |
|
|---|
| 1818 | if (nid.eq.0) write(m,10) txt10,ichk,ntot,xavg,xmin,xmax
|
|---|
| 1819 |
|
|---|
| 1820 | 10 format(3X,a10,2i8,' pts, avg,min,max = ',1p3g14.6)
|
|---|
| 1821 | c
|
|---|
| 1822 | return
|
|---|
| 1823 | end
|
|---|
| 1824 | c-----------------------------------------------------------------------
|
|---|
| 1825 | subroutine outfldan(x,n,txt10,ichk) ! writes x(1:n) into unit=p130+ifiled
|
|---|
| 1826 | INCLUDE 'SIZE' ! or into std. output for 0<p130<9
|
|---|
| 1827 | INCLUDE 'TSTEP' ! truncated below eps=p131
|
|---|
| 1828 | INCLUDE 'INPUT'
|
|---|
| 1829 | real x(1)
|
|---|
| 1830 | character*10 txt10
|
|---|
| 1831 | c parameter (eps=1.e-7)
|
|---|
| 1832 | C
|
|---|
| 1833 | p130 = param(130)
|
|---|
| 1834 | eps = param(131)
|
|---|
| 1835 | if (p130.le.0) return
|
|---|
| 1836 | m = 6
|
|---|
| 1837 | if (p130.gt.9) m = p130 + ifield
|
|---|
| 1838 |
|
|---|
| 1839 | ntot = n
|
|---|
| 1840 |
|
|---|
| 1841 | xmin = glmin(x,ntot)
|
|---|
| 1842 | xmax = glmax(x,ntot)
|
|---|
| 1843 | xavg = glsum(x,ntot)/ntot
|
|---|
| 1844 |
|
|---|
| 1845 | if (abs(xavg).lt.eps) xavg = 0. ! truncation
|
|---|
| 1846 |
|
|---|
| 1847 | if (nid.eq.0) write(m,10) txt10,ichk,ntot,xavg,xmin,xmax
|
|---|
| 1848 |
|
|---|
| 1849 | 10 format(3X,a10,2i8,' pts, avg,min,max = ',1p3g11.3)
|
|---|
| 1850 | c 10 format(3X,a10,2i8,' pts, avg,min,max = ',1p3g14.6)
|
|---|
| 1851 | c
|
|---|
| 1852 | return
|
|---|
| 1853 | end
|
|---|
| 1854 | c-----------------------------------------------------------------------
|
|---|
| 1855 | subroutine h1mg_solve(z,rhs,if_hybrid) ! Solve preconditioner: Mz=rhs
|
|---|
| 1856 | real z(1),rhs(1)
|
|---|
| 1857 |
|
|---|
| 1858 | c Assumes that preprocessing has been completed via h1mg_setup()
|
|---|
| 1859 |
|
|---|
| 1860 |
|
|---|
| 1861 | include 'SIZE'
|
|---|
| 1862 | include 'HSMG' ! Same array space as HSMG
|
|---|
| 1863 | include 'GEOM'
|
|---|
| 1864 | include 'INPUT'
|
|---|
| 1865 | include 'MASS'
|
|---|
| 1866 | include 'SOLN'
|
|---|
| 1867 | include 'TSTEP'
|
|---|
| 1868 | include 'CTIMER'
|
|---|
| 1869 | include 'PARALLEL'
|
|---|
| 1870 |
|
|---|
| 1871 | common /scrhi/ h2inv (lx1,ly1,lz1,lelv)
|
|---|
| 1872 | common /scrvh/ h1 (lx1,ly1,lz1,lelv),
|
|---|
| 1873 | $ h2 (lx1,ly1,lz1,lelv)
|
|---|
| 1874 | parameter (lt=lx1*ly1*lz1*lelt)
|
|---|
| 1875 | common /scrmg/ e(2*lt),w(lt),r(lt)
|
|---|
| 1876 | integer p_msk,p_b
|
|---|
| 1877 | logical if_hybrid
|
|---|
| 1878 |
|
|---|
| 1879 | c if_hybrid = .true. ! Control this from gmres, according
|
|---|
| 1880 | c if_hybrid = .false. ! to convergence efficiency
|
|---|
| 1881 |
|
|---|
| 1882 | nel = nelfld(ifield)
|
|---|
| 1883 |
|
|---|
| 1884 | op = 1. ! Coefficients for h1mg_ax
|
|---|
| 1885 | om = -1.
|
|---|
| 1886 | sigma = 1.
|
|---|
| 1887 | if (if_hybrid) sigma = 2./3.
|
|---|
| 1888 |
|
|---|
| 1889 | l = mg_h1_lmax
|
|---|
| 1890 | n = mg_h1_n(l,mg_fld)
|
|---|
| 1891 | is = 1 ! solve index
|
|---|
| 1892 |
|
|---|
| 1893 | call h1mg_schwarz(z,rhs,sigma,l) ! z := sigma W M rhs
|
|---|
| 1894 | ! Schwarz
|
|---|
| 1895 | call copy(r,rhs,n) ! r := rhs
|
|---|
| 1896 | if (if_hybrid) call h1mg_axm(r,z,op,om,l,w) ! r := rhs - A z
|
|---|
| 1897 | ! l
|
|---|
| 1898 |
|
|---|
| 1899 | do l = mg_h1_lmax-1,2,-1 ! DOWNWARD Leg of V-cycle
|
|---|
| 1900 | is = is + n
|
|---|
| 1901 | n = mg_h1_n(l,mg_fld)
|
|---|
| 1902 | ! T
|
|---|
| 1903 | call h1mg_rstr(r,l,.true.) ! r := J r
|
|---|
| 1904 | ! l l+1
|
|---|
| 1905 | ! OVERLAPPING Schwarz exchange and solve:
|
|---|
| 1906 | call h1mg_schwarz(e(is),r,sigma,l) ! e := sigma W M r
|
|---|
| 1907 | ! l Schwarz l
|
|---|
| 1908 |
|
|---|
| 1909 | if(if_hybrid)call h1mg_axm(r,e(is),op,om,l,w)! r := r - A e
|
|---|
| 1910 | ! l l
|
|---|
| 1911 | enddo
|
|---|
| 1912 | is = is+n
|
|---|
| 1913 | ! T
|
|---|
| 1914 | call h1mg_rstr(r,1,.false.) ! r := J r
|
|---|
| 1915 | ! l l+1
|
|---|
| 1916 | p_msk = p_mg_msk(l,mg_fld)
|
|---|
| 1917 | call h1mg_mask(r,mg_imask(p_msk),nel) ! -1
|
|---|
| 1918 | call hsmg_coarse_solve ( e(is) , r ) ! e := A r
|
|---|
| 1919 | call h1mg_mask(e(is),mg_imask(p_msk),nel) ! 1 1 1
|
|---|
| 1920 |
|
|---|
| 1921 | c nx = mg_nh(1)
|
|---|
| 1922 | c call outnxfld (e(is),nx,nelv,'ecrsb4',is)
|
|---|
| 1923 | c call h1mg_mask(e(is),mg_imask(p_msk),nel) ! 1 1 1
|
|---|
| 1924 | c call outnxfld (e(is),nx,nelv,'ecrsaf',is)
|
|---|
| 1925 | c call exitt
|
|---|
| 1926 |
|
|---|
| 1927 | do l = 2,mg_h1_lmax-1 ! UNWIND. No smoothing.
|
|---|
| 1928 | im = is
|
|---|
| 1929 | is = is - n
|
|---|
| 1930 | n = mg_h1_n(l,mg_fld)
|
|---|
| 1931 | call hsmg_intp (w,e(im),l-1) ! w := J e
|
|---|
| 1932 | i1=is-1 ! l-1
|
|---|
| 1933 | do i=1,n
|
|---|
| 1934 | e(i1+i) = e(i1+i) + w(i) ! e := e + w
|
|---|
| 1935 | enddo ! l l
|
|---|
| 1936 | enddo
|
|---|
| 1937 |
|
|---|
| 1938 | l = mg_h1_lmax
|
|---|
| 1939 | n = mg_h1_n(l,mg_fld)
|
|---|
| 1940 | im = is ! solve index
|
|---|
| 1941 | call hsmg_intp(w,e(im),l-1) ! w := J e
|
|---|
| 1942 | do i = 1,n ! l-1
|
|---|
| 1943 | z(i) = z(i) + w(i) ! z := z + w
|
|---|
| 1944 | enddo
|
|---|
| 1945 |
|
|---|
| 1946 | call dsavg(z) ! Emergency hack --- to ensure continuous z!
|
|---|
| 1947 |
|
|---|
| 1948 | return
|
|---|
| 1949 | end
|
|---|
| 1950 | c-----------------------------------------------------------------------
|
|---|
| 1951 | subroutine h1mg_axm(w,p,aw,ap,l,wk)
|
|---|
| 1952 | c
|
|---|
| 1953 | c w := aw*w + ap*H*p, level l, with mask and dssum
|
|---|
| 1954 | c
|
|---|
| 1955 | c Hu := div. h1 grad u + h2 u
|
|---|
| 1956 | c
|
|---|
| 1957 | c ~= h1 A u + h2 B u
|
|---|
| 1958 | c
|
|---|
| 1959 | c Here, we assume that pointers into g() and h1() and h2() have
|
|---|
| 1960 | c been established
|
|---|
| 1961 | c
|
|---|
| 1962 | include 'SIZE'
|
|---|
| 1963 | include 'HSMG'
|
|---|
| 1964 | include 'TSTEP' ! nelfld
|
|---|
| 1965 |
|
|---|
| 1966 | real w(1),p(1),wk(1)
|
|---|
| 1967 |
|
|---|
| 1968 | integer p_h1,p_h2,p_g,p_b,p_msk
|
|---|
| 1969 | logical ifh2
|
|---|
| 1970 |
|
|---|
| 1971 | p_h1 = p_mg_h1 (l,mg_fld)
|
|---|
| 1972 | p_h2 = p_mg_h2 (l,mg_fld)
|
|---|
| 1973 | p_g = p_mg_g (l,mg_fld)
|
|---|
| 1974 | p_b = p_mg_b (l,mg_fld)
|
|---|
| 1975 | p_msk = p_mg_msk (l,mg_fld)
|
|---|
| 1976 |
|
|---|
| 1977 | if (p_h1 .eq.0) call mg_set_h1 (p_h1 ,l)
|
|---|
| 1978 | if (p_h2 .eq.0) call mg_set_h2 (p_h2 ,l)
|
|---|
| 1979 | if (p_g .eq.0) call mg_set_gb (p_g,p_b,l)
|
|---|
| 1980 | if (p_msk.eq.0) call mg_set_msk (p_msk,l)
|
|---|
| 1981 |
|
|---|
| 1982 | ifh2 = .true.
|
|---|
| 1983 | if (p_h2.eq.0) ifh2 = .false. ! Need a much better mech.
|
|---|
| 1984 | ifh2 = .false.
|
|---|
| 1985 |
|
|---|
| 1986 | nx = mg_nh(l)
|
|---|
| 1987 | ny = mg_nh(l)
|
|---|
| 1988 | nz = mg_nhz(l)
|
|---|
| 1989 | ng = 3*ldim-3
|
|---|
| 1990 |
|
|---|
| 1991 |
|
|---|
| 1992 | call h1mg_axml (wk,p
|
|---|
| 1993 | $ ,mg_h1(p_h1),mg_h2(p_h2),nx,ny,nz,nelfld(ifield)
|
|---|
| 1994 | $ ,mg_g (p_g) , ng ,mg_b(p_b), mg_imask(p_msk),ifh2)
|
|---|
| 1995 |
|
|---|
| 1996 |
|
|---|
| 1997 | call hsmg_dssum (wk,l)
|
|---|
| 1998 |
|
|---|
| 1999 | n = nx*ny*nz*nelfld(ifield)
|
|---|
| 2000 | call add2sxy (w,aw,wk,ap,n)
|
|---|
| 2001 |
|
|---|
| 2002 | return
|
|---|
| 2003 | end
|
|---|
| 2004 | c-----------------------------------------------------------------------
|
|---|
| 2005 | subroutine h1mg_axml
|
|---|
| 2006 | $ (w,p,h1,h2,nx,ny,nz,nel,g,ng,b,mask,ifh2)
|
|---|
| 2007 | c
|
|---|
| 2008 | c w := aw*w + ap*H*p, level l, with mask and dssum
|
|---|
| 2009 | c
|
|---|
| 2010 | c Hu := div. h1 grad u + h2 u
|
|---|
| 2011 | c
|
|---|
| 2012 | c ~= h1 A u + h2 B u
|
|---|
| 2013 | c
|
|---|
| 2014 |
|
|---|
| 2015 | include 'SIZE'
|
|---|
| 2016 | include 'TOTAL'
|
|---|
| 2017 | include 'HSMG'
|
|---|
| 2018 |
|
|---|
| 2019 | real w (nx*ny*nz,nel), p (nx*ny*nz,nel)
|
|---|
| 2020 | $ , h1(nx*ny*nz,nel), h2(nx*ny*nz,nel)
|
|---|
| 2021 | $ , b (nx*ny*nz,nel), g (ng*nx*ny*nz,nel)
|
|---|
| 2022 | integer mask(1)
|
|---|
| 2023 |
|
|---|
| 2024 | logical ifh2
|
|---|
| 2025 |
|
|---|
| 2026 | parameter (lxyz=lx1*ly1*lz1)
|
|---|
| 2027 | common /ctmp0/ ur(lxyz),us(lxyz),ut(lxyz)
|
|---|
| 2028 |
|
|---|
| 2029 | integer e
|
|---|
| 2030 |
|
|---|
| 2031 | do e=1,nel
|
|---|
| 2032 |
|
|---|
| 2033 | call axe(w(1,e),p(1,e),h1(1,e),h2(1,e),g(1,e),ng,b(1,e)
|
|---|
| 2034 | $ ,nx,ny,nz,ur,us,ut,ifh2,ifrzer(e),e)
|
|---|
| 2035 |
|
|---|
| 2036 | im = mask(e)
|
|---|
| 2037 | call mg_mask_e(w,mask(im)) ! Zero out Dirichlet conditions
|
|---|
| 2038 |
|
|---|
| 2039 | enddo
|
|---|
| 2040 |
|
|---|
| 2041 | return
|
|---|
| 2042 | end
|
|---|
| 2043 | c-----------------------------------------------------------------------
|
|---|
| 2044 | subroutine h1mg_mask(w,mask,nel)
|
|---|
| 2045 | include 'SIZE'
|
|---|
| 2046 |
|
|---|
| 2047 | real w (1)
|
|---|
| 2048 | integer mask(1) ! Pointer to Dirichlet BCs
|
|---|
| 2049 | integer e
|
|---|
| 2050 |
|
|---|
| 2051 | do e=1,nel
|
|---|
| 2052 | im = mask(e)
|
|---|
| 2053 | call mg_mask_e(w,mask(im)) ! Zero out Dirichlet conditions
|
|---|
| 2054 | enddo
|
|---|
| 2055 |
|
|---|
| 2056 | return
|
|---|
| 2057 | end
|
|---|
| 2058 | c----------------------------------------------------------------------
|
|---|
| 2059 | subroutine mg_mask_e(w,mask) ! Zero out Dirichlet conditions
|
|---|
| 2060 | include 'SIZE'
|
|---|
| 2061 | real w(1)
|
|---|
| 2062 | integer mask(0:1)
|
|---|
| 2063 |
|
|---|
| 2064 | n=mask(0)
|
|---|
| 2065 | do i=1,n
|
|---|
| 2066 | c write(6,*) i,mask(i),n,' MG_MASK'
|
|---|
| 2067 | w(mask(i)) = 0.
|
|---|
| 2068 | enddo
|
|---|
| 2069 |
|
|---|
| 2070 | return
|
|---|
| 2071 | end
|
|---|
| 2072 | c-----------------------------------------------------------------------
|
|---|
| 2073 | subroutine axe
|
|---|
| 2074 | $ (w,p,h1,h2,g,ng,b,nx,ny,nz,ur,us,ut,ifh2,ifrz,e)
|
|---|
| 2075 |
|
|---|
| 2076 | include 'SIZE'
|
|---|
| 2077 | include 'INPUT' ! if3d
|
|---|
| 2078 | logical ifh2,ifrz
|
|---|
| 2079 |
|
|---|
| 2080 | real w (nx*ny*nz), p (nx*ny*nz)
|
|---|
| 2081 | $ , h1(nx*ny*nz), h2(nx*ny*nz)
|
|---|
| 2082 | $ , b (nx*ny*nz), g (ng,nx*ny*nz)
|
|---|
| 2083 | $ , ur(nx*ny*nz), us(nx*ny*nz), ut(nx*ny*nz)
|
|---|
| 2084 | integer e
|
|---|
| 2085 |
|
|---|
| 2086 | nxyz = nx*ny*nz
|
|---|
| 2087 |
|
|---|
| 2088 | call gradl_rst(ur,us,ut,p,nx,if3d)
|
|---|
| 2089 | c if (e.eq.1) then
|
|---|
| 2090 | c call outmat(p ,nx,ny,'ur A p',e)
|
|---|
| 2091 | c call outmat(ur,nx,ny,'ur A r',e)
|
|---|
| 2092 | c call outmat(us,nx,ny,'ur A s',e)
|
|---|
| 2093 | c endif
|
|---|
| 2094 |
|
|---|
| 2095 | if (if3d) then
|
|---|
| 2096 | do i=1,nxyz
|
|---|
| 2097 | wr = g(1,i)*ur(i) + g(4,i)*us(i) + g(5,i)*ut(i)
|
|---|
| 2098 | ws = g(4,i)*ur(i) + g(2,i)*us(i) + g(6,i)*ut(i)
|
|---|
| 2099 | wt = g(5,i)*ur(i) + g(6,i)*us(i) + g(3,i)*ut(i)
|
|---|
| 2100 | ur(i) = wr*h1(i)
|
|---|
| 2101 | us(i) = ws*h1(i)
|
|---|
| 2102 | ut(i) = wt*h1(i)
|
|---|
| 2103 | enddo
|
|---|
| 2104 | elseif (ifaxis) then
|
|---|
| 2105 | call exitti('Blame Paul for no gradl_rst support yet$',nx)
|
|---|
| 2106 | else
|
|---|
| 2107 | do i=1,nxyz
|
|---|
| 2108 | wr = g(1,i)*ur(i) + g(3,i)*us(i)
|
|---|
| 2109 | ws = g(3,i)*ur(i) + g(2,i)*us(i)
|
|---|
| 2110 | c write(6,3) i,ur(i),wr,g(1,i)/b(i),b(i)
|
|---|
| 2111 | c 3 format(i5,1p4e12.4,' wrws')
|
|---|
| 2112 | ur(i) = wr*h1(i)
|
|---|
| 2113 | us(i) = ws*h1(i)
|
|---|
| 2114 | enddo
|
|---|
| 2115 | endif
|
|---|
| 2116 |
|
|---|
| 2117 | call gradl_rst_t(w,ur,us,ut,nx,if3d)
|
|---|
| 2118 |
|
|---|
| 2119 | c if (e.eq.1) then
|
|---|
| 2120 | c call outmat(w ,nx,ny,'ur B w',e)
|
|---|
| 2121 | c call outmat(ur,nx,ny,'ur B r',e)
|
|---|
| 2122 | c call outmat(us,nx,ny,'ur B s',e)
|
|---|
| 2123 | c endif
|
|---|
| 2124 |
|
|---|
| 2125 | if (ifh2) then
|
|---|
| 2126 | do i=1,nxyz
|
|---|
| 2127 | w(i)=w(i)+h2(i)*b(i)*p(i)
|
|---|
| 2128 | enddo
|
|---|
| 2129 | endif
|
|---|
| 2130 |
|
|---|
| 2131 | return
|
|---|
| 2132 | end
|
|---|
| 2133 | c----------------------------------------------------------------------
|
|---|
| 2134 | subroutine hsmg_tnsr1(v,nv,nu,A,At)
|
|---|
| 2135 | c
|
|---|
| 2136 | c v = [A (x) A] u or
|
|---|
| 2137 | c v = [A (x) A (x) A] u
|
|---|
| 2138 | c
|
|---|
| 2139 | integer nv,nu
|
|---|
| 2140 | real v(1),A(1),At(1)
|
|---|
| 2141 | include 'SIZE'
|
|---|
| 2142 | include 'INPUT'
|
|---|
| 2143 | if (.not. if3d) then
|
|---|
| 2144 | call hsmg_tnsr1_2d(v,nv,nu,A,At)
|
|---|
| 2145 | else
|
|---|
| 2146 | call hsmg_tnsr1_3d(v,nv,nu,A,At,At)
|
|---|
| 2147 | endif
|
|---|
| 2148 | return
|
|---|
| 2149 | end
|
|---|
| 2150 | c-------------------------------------------------------T--------------
|
|---|
| 2151 | subroutine hsmg_tnsr1_2d(v,nv,nu,A,Bt) ! u = A u B
|
|---|
| 2152 | integer nv,nu
|
|---|
| 2153 | real v(1),A(1),Bt(1)
|
|---|
| 2154 | include 'SIZE'
|
|---|
| 2155 | common /hsmgw/ work(lx1*lx1)
|
|---|
| 2156 | integer e
|
|---|
| 2157 |
|
|---|
| 2158 | nv2 = nv*nv
|
|---|
| 2159 | nu2 = nu*nu
|
|---|
| 2160 |
|
|---|
| 2161 | if (nv.le.nu) then
|
|---|
| 2162 | iv=1
|
|---|
| 2163 | iu=1
|
|---|
| 2164 | do e=1,nelv
|
|---|
| 2165 | call mxm(A,nv,v(iu),nu,work,nu)
|
|---|
| 2166 | call mxm(work,nv,Bt,nu,v(iv),nv)
|
|---|
| 2167 | iv = iv + nv2
|
|---|
| 2168 | iu = iu + nu2
|
|---|
| 2169 | enddo
|
|---|
| 2170 | else
|
|---|
| 2171 | do e=nelv,1,-1
|
|---|
| 2172 | iu=1+nu2*(e-1)
|
|---|
| 2173 | iv=1+nv2*(e-1)
|
|---|
| 2174 | call mxm(A,nv,v(iu),nu,work,nu)
|
|---|
| 2175 | call mxm(work,nv,Bt,nu,v(iv),nv)
|
|---|
| 2176 | enddo
|
|---|
| 2177 | endif
|
|---|
| 2178 |
|
|---|
| 2179 | return
|
|---|
| 2180 | end
|
|---|
| 2181 | c----------------------------------------------------------------------
|
|---|
| 2182 | subroutine hsmg_tnsr1_3d(v,nv,nu,A,Bt,Ct) ! v = [C (x) B (x) A] u
|
|---|
| 2183 | integer nv,nu
|
|---|
| 2184 | real v(1),A(1),Bt(1),Ct(1)
|
|---|
| 2185 | include 'SIZE'
|
|---|
| 2186 | parameter (lwk=(lx1+2)*(ly1+2)*(lz1+2))
|
|---|
| 2187 | common /hsmgw/ work(0:lwk-1),work2(0:lwk-1)
|
|---|
| 2188 | integer e,e0,ee,es
|
|---|
| 2189 |
|
|---|
| 2190 | e0=1
|
|---|
| 2191 | es=1
|
|---|
| 2192 | ee=nelv
|
|---|
| 2193 |
|
|---|
| 2194 | if (nv.gt.nu) then
|
|---|
| 2195 | e0=nelv
|
|---|
| 2196 | es=-1
|
|---|
| 2197 | ee=1
|
|---|
| 2198 | endif
|
|---|
| 2199 |
|
|---|
| 2200 | nu3 = nu**3
|
|---|
| 2201 | nv3 = nv**3
|
|---|
| 2202 |
|
|---|
| 2203 | do e=e0,ee,es
|
|---|
| 2204 | iu = 1 + (e-1)*nu3
|
|---|
| 2205 | iv = 1 + (e-1)*nv3
|
|---|
| 2206 | call mxm(A,nv,v(iu),nu,work,nu*nu)
|
|---|
| 2207 | do i=0,nu-1
|
|---|
| 2208 | call mxm(work(nv*nu*i),nv,Bt,nu,work2(nv*nv*i),nv)
|
|---|
| 2209 | enddo
|
|---|
| 2210 | call mxm(work2,nv*nv,Ct,nu,v(iv),nv)
|
|---|
| 2211 | enddo
|
|---|
| 2212 |
|
|---|
| 2213 | return
|
|---|
| 2214 | end
|
|---|
| 2215 | c------------------------------------------ T -----------------------
|
|---|
| 2216 | subroutine h1mg_rstr(r,l,ifdssum) ! r =J r, l is coarse level
|
|---|
| 2217 | include 'SIZE'
|
|---|
| 2218 | include 'HSMG'
|
|---|
| 2219 | logical ifdssum
|
|---|
| 2220 |
|
|---|
| 2221 | real r(1)
|
|---|
| 2222 | integer l
|
|---|
| 2223 |
|
|---|
| 2224 | call hsmg_do_wt(r,mg_rstr_wt(mg_rstr_wt_index(l+1,mg_fld))
|
|---|
| 2225 | $ ,mg_nh(l+1),mg_nh(l+1),mg_nhz(l+1))
|
|---|
| 2226 |
|
|---|
| 2227 | call hsmg_tnsr1(r,mg_nh(l),mg_nh(l+1),mg_jht(1,l),mg_jh(1,l))
|
|---|
| 2228 |
|
|---|
| 2229 | if (ifdssum) call hsmg_dssum(r,l)
|
|---|
| 2230 |
|
|---|
| 2231 | return
|
|---|
| 2232 | end
|
|---|
| 2233 | c----------------------------------------------------------------------
|
|---|
| 2234 | subroutine h1mg_setup()
|
|---|
| 2235 | include 'SIZE'
|
|---|
| 2236 | include 'TOTAL'
|
|---|
| 2237 | include 'HSMG'
|
|---|
| 2238 |
|
|---|
| 2239 | common /scrhi/ h2inv (lx1,ly1,lz1,lelt)
|
|---|
| 2240 | common /scrvh/ h1 (lx1,ly1,lz1,lelt),
|
|---|
| 2241 | $ h2 (lx1,ly1,lz1,lelt)
|
|---|
| 2242 |
|
|---|
| 2243 | integer p_h1,p_h2,p_g,p_b,p_msk
|
|---|
| 2244 |
|
|---|
| 2245 |
|
|---|
| 2246 | param(59) = 1
|
|---|
| 2247 | call geom_reset(1) ! Recompute g1m1 etc. with deformed only
|
|---|
| 2248 |
|
|---|
| 2249 | n = lx1*ly1*lz1*nelt
|
|---|
| 2250 | call rone (h1 ,n)
|
|---|
| 2251 | call rzero(h2 ,n)
|
|---|
| 2252 | call rzero(h2inv,n)
|
|---|
| 2253 |
|
|---|
| 2254 | call h1mg_setup_mg_nx
|
|---|
| 2255 | call h1mg_setup_semhat ! SEM hat matrices for each level
|
|---|
| 2256 | call hsmg_setup_intp ! Interpolation operators
|
|---|
| 2257 | call h1mg_setup_dssum ! set direct stiffness summation handles
|
|---|
| 2258 | call h1mg_setup_wtmask ! set restriction weight matrices and bc masks
|
|---|
| 2259 | call h1mg_setup_fdm ! set up fast diagonalization method
|
|---|
| 2260 | call h1mg_setup_schwarz_wt(.false.)
|
|---|
| 2261 | call hsmg_setup_solve ! set up the solver
|
|---|
| 2262 |
|
|---|
| 2263 | l=mg_h1_lmax
|
|---|
| 2264 | call mg_set_h1 (p_h1 ,l)
|
|---|
| 2265 | call mg_set_h2 (p_h2 ,l)
|
|---|
| 2266 | call mg_set_gb (p_g,p_b,l)
|
|---|
| 2267 | call mg_set_msk (p_msk,l)
|
|---|
| 2268 |
|
|---|
| 2269 | return
|
|---|
| 2270 | end
|
|---|
| 2271 | c-----------------------------------------------------------------------
|
|---|
| 2272 | subroutine h1mg_setup_mg_nx()
|
|---|
| 2273 | include 'SIZE'
|
|---|
| 2274 | include 'INPUT'
|
|---|
| 2275 | include 'PARALLEL'
|
|---|
| 2276 | include 'HSMG'
|
|---|
| 2277 | include 'SEMHAT'
|
|---|
| 2278 | include 'TSTEP' ! nelfld
|
|---|
| 2279 | real w(lx1+2)
|
|---|
| 2280 | integer nf,nc,nr
|
|---|
| 2281 | integer nx,ny,nz
|
|---|
| 2282 |
|
|---|
| 2283 | integer mgn2(10)
|
|---|
| 2284 | save mgn2
|
|---|
| 2285 | data mgn2 / 1, 2, 2, 2, 2, 3, 3, 5, 5, 5/
|
|---|
| 2286 | c data mgn2 / 1, 2, 3, 4, 5, 6, 7, 8, 9, 0
|
|---|
| 2287 |
|
|---|
| 2288 | c if (param(82).eq.0) param(82)=2 ! nek default
|
|---|
| 2289 | c if (np.eq.1) param(82)=2 ! single proc. too slow
|
|---|
| 2290 | p82 = 2 ! potentially variable nxc
|
|---|
| 2291 | mg_h1_lmax = 3
|
|---|
| 2292 | c mg_h1_lmax = 4
|
|---|
| 2293 | if (lx1.eq.4) mg_h1_lmax = 2
|
|---|
| 2294 | c if (param(79).ne.0) mg_h1_lmax = param(79)
|
|---|
| 2295 | mg_lmax = mg_h1_lmax
|
|---|
| 2296 | mglx1 = p82-1 !1
|
|---|
| 2297 | mg_nx(1) = mglx1
|
|---|
| 2298 | mg_ny(1) = mglx1
|
|---|
| 2299 | mg_nz(1) = mglx1
|
|---|
| 2300 | if (.not.if3d) mg_nz(1) = 0
|
|---|
| 2301 |
|
|---|
| 2302 | mglx2 = 2*(lx2/4) + 1
|
|---|
| 2303 | if (lx1.eq.5) mglx2 = 3
|
|---|
| 2304 | c if (lx1.eq.6) mglx2 = 3
|
|---|
| 2305 | if (lx1.le.10) mglx2 = mgn2(lx1)
|
|---|
| 2306 | if (lx1.eq.8) mglx2 = 4
|
|---|
| 2307 | if (lx1.eq.8) mglx2 = 3
|
|---|
| 2308 |
|
|---|
| 2309 | mglx2 = min(3,mglx2) ! This choice seems best (9/24/12)
|
|---|
| 2310 |
|
|---|
| 2311 | mg_nx(2) = mglx2
|
|---|
| 2312 | mg_ny(2) = mglx2
|
|---|
| 2313 | mg_nz(2) = mglx2
|
|---|
| 2314 | if (.not.if3d) mg_nz(2) = 0
|
|---|
| 2315 |
|
|---|
| 2316 | mg_nx(3) = mglx2+1
|
|---|
| 2317 | mg_ny(3) = mglx2+1
|
|---|
| 2318 | mg_nz(3) = mglx2+1
|
|---|
| 2319 | if (.not.if3d) mg_nz(3) = 0
|
|---|
| 2320 |
|
|---|
| 2321 | mg_nx(mg_h1_lmax) = lx1-1
|
|---|
| 2322 | mg_ny(mg_h1_lmax) = ly1-1
|
|---|
| 2323 | mg_nz(mg_h1_lmax) = lz1-1
|
|---|
| 2324 |
|
|---|
| 2325 | if (nio.eq.0) write(*,*) 'h1_mg_nx:',(mg_nx(i),i=1,mg_h1_lmax)
|
|---|
| 2326 | if (nio.eq.0) write(*,*) 'h1_mg_ny:',(mg_ny(i),i=1,mg_h1_lmax)
|
|---|
| 2327 | if (nio.eq.0) write(*,*) 'h1_mg_nz:',(mg_nz(i),i=1,mg_h1_lmax)
|
|---|
| 2328 |
|
|---|
| 2329 | do ifld=1,ldimt1
|
|---|
| 2330 | do l=1,mg_lmax
|
|---|
| 2331 | mg_h1_n(l,ifld)=(mg_nx(l)+1)
|
|---|
| 2332 | $ *(mg_ny(l)+1)
|
|---|
| 2333 | $ *(mg_nz(l)+1)*nelfld(ifld)
|
|---|
| 2334 | enddo
|
|---|
| 2335 | enddo
|
|---|
| 2336 |
|
|---|
| 2337 | return
|
|---|
| 2338 | end
|
|---|
| 2339 | c----------------------------------------------------------------------
|
|---|
| 2340 | subroutine h1mg_setup_semhat ! SEM hat matrices for each level
|
|---|
| 2341 | include 'SIZE'
|
|---|
| 2342 | include 'INPUT'
|
|---|
| 2343 | include 'HSMG'
|
|---|
| 2344 | include 'SEMHAT'
|
|---|
| 2345 |
|
|---|
| 2346 | do l=1,mg_h1_lmax
|
|---|
| 2347 | n = mg_nx(l) ! polynomial order
|
|---|
| 2348 | call semhat(ah,bh,ch,dh,zh,dph,jph,bgl,zglhat,dgl,jgl,n,wh)
|
|---|
| 2349 | call copy(mg_ah(1,l),ah,(n+1)*(n+1))
|
|---|
| 2350 | call copy(mg_bh(1,l),bh,n+1)
|
|---|
| 2351 | call copy(mg_dh(1,l),dh,(n+1)*(n+1))
|
|---|
| 2352 | call transpose(mg_dht(1,l),n+1,dh,n+1)
|
|---|
| 2353 | call copy(mg_zh(1,l),zh,n+1)
|
|---|
| 2354 |
|
|---|
| 2355 | mg_nh(l)=n+1
|
|---|
| 2356 | mg_nhz(l)=mg_nz(l)+1
|
|---|
| 2357 |
|
|---|
| 2358 | enddo
|
|---|
| 2359 | end
|
|---|
| 2360 | c----------------------------------------------------------------------
|
|---|
| 2361 | subroutine h1mg_setup_dssum
|
|---|
| 2362 | include 'SIZE'
|
|---|
| 2363 | include 'INPUT'
|
|---|
| 2364 | include 'PARALLEL'
|
|---|
| 2365 | include 'HSMG'
|
|---|
| 2366 | parameter (lxyz=(lx1+2)*(ly1+2)*(lz1+2))
|
|---|
| 2367 | common /c_is1/ glo_num(lxyz*lelt)
|
|---|
| 2368 | common /ivrtx/ vertex ((2**ldim)*lelt)
|
|---|
| 2369 |
|
|---|
| 2370 | integer*8 glo_num
|
|---|
| 2371 | integer vertex
|
|---|
| 2372 | integer nx,ny,nz
|
|---|
| 2373 | integer l
|
|---|
| 2374 |
|
|---|
| 2375 | ncrnr = 2**ldim
|
|---|
| 2376 | call get_vert
|
|---|
| 2377 |
|
|---|
| 2378 |
|
|---|
| 2379 | do l=1,mg_lmax ! set up direct stiffness summation for each level
|
|---|
| 2380 | nx=mg_nh(l)
|
|---|
| 2381 | ny=mg_nh(l)
|
|---|
| 2382 | nz=mg_nhz(l)
|
|---|
| 2383 | call setupds(mg_gsh_handle(l,mg_fld),nx,ny,nz
|
|---|
| 2384 | $ ,nelv,nelgv,vertex,glo_num)
|
|---|
| 2385 | nx=nx+2
|
|---|
| 2386 | ny=ny+2
|
|---|
| 2387 | nz=nz+2
|
|---|
| 2388 | if(.not.if3d) nz=1
|
|---|
| 2389 | call setupds(mg_gsh_schwarz_handle(l,mg_fld),nx,ny,nz
|
|---|
| 2390 | $ ,nelv,nelgv,vertex,glo_num)
|
|---|
| 2391 | enddo
|
|---|
| 2392 |
|
|---|
| 2393 | return
|
|---|
| 2394 | end
|
|---|
| 2395 | c----------------------------------------------------------------------
|
|---|
| 2396 | subroutine mg_set_msk(p_msk ,l0)
|
|---|
| 2397 | include 'SIZE'
|
|---|
| 2398 | include 'HSMG'
|
|---|
| 2399 | include 'TSTEP'
|
|---|
| 2400 | integer p_msk
|
|---|
| 2401 |
|
|---|
| 2402 | l = mg_h1_lmax
|
|---|
| 2403 | p_mg_msk(l,mg_fld) = 0
|
|---|
| 2404 | n = mg_h1_n(l,mg_fld)
|
|---|
| 2405 |
|
|---|
| 2406 |
|
|---|
| 2407 | do l=mg_h1_lmax,1,-1
|
|---|
| 2408 | nx = mg_nh (l)
|
|---|
| 2409 | ny = mg_nh (l)
|
|---|
| 2410 | nz = mg_nhz (l)
|
|---|
| 2411 |
|
|---|
| 2412 | p_msk = p_mg_msk(l,mg_fld)
|
|---|
| 2413 |
|
|---|
| 2414 | call h1mg_setup_mask
|
|---|
| 2415 | $ (mg_imask(p_msk),nm,nx,ny,nz,nelfld(ifield),l,mg_work)
|
|---|
| 2416 |
|
|---|
| 2417 | if (l.gt.1) p_mg_msk(l-1,mg_fld)=p_mg_msk(l,mg_fld)+nm
|
|---|
| 2418 |
|
|---|
| 2419 | enddo
|
|---|
| 2420 |
|
|---|
| 2421 | p_msk = p_mg_msk(l0,mg_fld)
|
|---|
| 2422 |
|
|---|
| 2423 | return
|
|---|
| 2424 | end
|
|---|
| 2425 | c----------------------------------------------------------------------
|
|---|
| 2426 | subroutine h1mg_setup_mask(mask,nm,nx,ny,nz,nel,l,w)
|
|---|
| 2427 | include 'SIZE'
|
|---|
| 2428 | include 'INPUT' ! if3d
|
|---|
| 2429 |
|
|---|
| 2430 | integer mask(1) ! Pointer to Dirichlet BCs
|
|---|
| 2431 | integer nx,ny,nz,l
|
|---|
| 2432 | real w(nx,ny,nz,nel)
|
|---|
| 2433 |
|
|---|
| 2434 | integer e,count,ptr
|
|---|
| 2435 | integer lbr,rbr,lbs,rbs,lbt,rbt,two
|
|---|
| 2436 |
|
|---|
| 2437 | zero = 0
|
|---|
| 2438 | nxyz = nx*ny*nz
|
|---|
| 2439 | n = nx*ny*nz*nel
|
|---|
| 2440 |
|
|---|
| 2441 | call rone(w,n) ! Init everything to 1
|
|---|
| 2442 |
|
|---|
| 2443 | ierrmx = 0 ! BC verification
|
|---|
| 2444 | two = 2
|
|---|
| 2445 | do e=1,nel ! Set dirichlet nodes to zero
|
|---|
| 2446 |
|
|---|
| 2447 | call get_fast_bc(lbr,rbr,lbs,rbs,lbt,rbt,e,two,ierr)
|
|---|
| 2448 | c write(6,6) e,lbr,rbr,lbs,rbs,ierr,nx
|
|---|
| 2449 | c 6 format(i5,2x,4i3,2x,i2,3x,i5,' lbr,rbr,lbs')
|
|---|
| 2450 |
|
|---|
| 2451 | if (lbr.eq.1) call facev(w,e,4,zero,nx,ny,nz)
|
|---|
| 2452 | if (rbr.eq.1) call facev(w,e,2,zero,nx,ny,nz)
|
|---|
| 2453 | if (lbs.eq.1) call facev(w,e,1,zero,nx,ny,nz)
|
|---|
| 2454 | if (rbs.eq.1) call facev(w,e,3,zero,nx,ny,nz)
|
|---|
| 2455 | if (if3d) then
|
|---|
| 2456 | if (lbt.eq.1) call facev(w,e,5,zero,nx,ny,nz)
|
|---|
| 2457 | if (rbt.eq.1) call facev(w,e,6,zero,nx,ny,nz)
|
|---|
| 2458 | endif
|
|---|
| 2459 | ierrmx = max(ierrmx,ierr)
|
|---|
| 2460 | enddo
|
|---|
| 2461 |
|
|---|
| 2462 | call hsmg_dsprod(w,l) ! direct stiffness multiply
|
|---|
| 2463 |
|
|---|
| 2464 | c
|
|---|
| 2465 | c Prototypical mask layout, nel=5:
|
|---|
| 2466 | c
|
|---|
| 2467 | c e=1 ... 10
|
|---|
| 2468 | c 1 2 3 4 5 ... 10 | 11 12 13 14 | 15 | 16 |
|
|---|
| 2469 | c 11 15 16 ... | 3 p1 p2 p3 | 0 | 0 | ...
|
|---|
| 2470 | c ^
|
|---|
| 2471 | c |
|
|---|
| 2472 | c |_count for e=1
|
|---|
| 2473 | c
|
|---|
| 2474 |
|
|---|
| 2475 | nm = 1 ! store mask
|
|---|
| 2476 | do e=1,nel
|
|---|
| 2477 |
|
|---|
| 2478 | mask(e) = nel+nm
|
|---|
| 2479 | count = 0 ! # Dirchlet points on element e
|
|---|
| 2480 | ptr = mask(e)
|
|---|
| 2481 |
|
|---|
| 2482 | do i=1,nxyz
|
|---|
| 2483 | if (w(i,1,1,e).eq.0) then
|
|---|
| 2484 | nm = nm +1
|
|---|
| 2485 | count = count+1
|
|---|
| 2486 | ptr = ptr +1
|
|---|
| 2487 | mask(ptr) = i + nxyz*(e-1) ! where I mask on element e
|
|---|
| 2488 | endif
|
|---|
| 2489 | enddo
|
|---|
| 2490 |
|
|---|
| 2491 |
|
|---|
| 2492 | ptr = mask(e)
|
|---|
| 2493 | mask(ptr) = count
|
|---|
| 2494 |
|
|---|
| 2495 | nm = nm+1 ! bump pointer to hold next count
|
|---|
| 2496 |
|
|---|
| 2497 | enddo
|
|---|
| 2498 |
|
|---|
| 2499 | nm = nel + nm-1 ! Return total number of mask pointers/counters
|
|---|
| 2500 |
|
|---|
| 2501 | ierrmx = iglmax(ierrmx,1)
|
|---|
| 2502 | if (ierrmx.gt.0) then
|
|---|
| 2503 | if (ierr.gt.0) write(6,*) nid,ierr,' BC FAIL h1'
|
|---|
| 2504 | call exitti('D INVALID BC FOUND in h1mg_setup_mask$',ierrmx)
|
|---|
| 2505 | endif
|
|---|
| 2506 |
|
|---|
| 2507 | return
|
|---|
| 2508 | end
|
|---|
| 2509 | c----------------------------------------------------------------------
|
|---|
| 2510 | subroutine mg_set_h1 (p_h1 ,l0)
|
|---|
| 2511 | include 'SIZE'
|
|---|
| 2512 | include 'HSMG'
|
|---|
| 2513 | integer pf,pc
|
|---|
| 2514 |
|
|---|
| 2515 | c As a first pass, rely on the cheesy common-block interface to get h1
|
|---|
| 2516 |
|
|---|
| 2517 | common /scrvh/ h1 (lx1,ly1,lz1,lelv)
|
|---|
| 2518 | $ , h2 (lx1,ly1,lz1,lelv)
|
|---|
| 2519 | $ , h2inv (lx1,ly1,lz1,lelv)
|
|---|
| 2520 |
|
|---|
| 2521 | integer p_h1
|
|---|
| 2522 |
|
|---|
| 2523 | l = mg_h1_lmax
|
|---|
| 2524 | p_mg_h1(l,mg_fld) = 0
|
|---|
| 2525 | n = mg_h1_n(l,mg_fld)
|
|---|
| 2526 |
|
|---|
| 2527 | call copy (mg_h1,h1,n) ! Fine grid is just original h1
|
|---|
| 2528 |
|
|---|
| 2529 | nx = mg_nh(l)
|
|---|
| 2530 | ny = mg_nh(l)
|
|---|
| 2531 | nz = mg_nhz(l)
|
|---|
| 2532 |
|
|---|
| 2533 | do l=mg_h1_lmax-1,1,-1
|
|---|
| 2534 |
|
|---|
| 2535 | p_mg_h1(l,mg_fld) = p_mg_h1(l+1,mg_fld) + n
|
|---|
| 2536 | n = mg_h1_n(l ,mg_fld)
|
|---|
| 2537 |
|
|---|
| 2538 | pf = p_mg_h1(l+1,mg_fld)
|
|---|
| 2539 | pc = p_mg_h1(l ,mg_fld)
|
|---|
| 2540 |
|
|---|
| 2541 | call hsmg_intp_fc (mg_h1(pc),mg_h1(pf),l)
|
|---|
| 2542 |
|
|---|
| 2543 | enddo
|
|---|
| 2544 |
|
|---|
| 2545 | p_h1 = p_mg_h1(l0,mg_fld)
|
|---|
| 2546 |
|
|---|
| 2547 | return
|
|---|
| 2548 | end
|
|---|
| 2549 | c-----------------------------------------------------------------------
|
|---|
| 2550 | subroutine mg_set_h2 (p_h2 ,l0)
|
|---|
| 2551 | include 'SIZE'
|
|---|
| 2552 | include 'HSMG'
|
|---|
| 2553 |
|
|---|
| 2554 | c As a first pass, rely on the cheesy common-block interface to get h2
|
|---|
| 2555 |
|
|---|
| 2556 | common /scrvh/ h1 (lx1,ly1,lz1,lelv)
|
|---|
| 2557 | $ , h2 (lx1,ly1,lz1,lelv)
|
|---|
| 2558 | $ , h2inv (lx1,ly1,lz1,lelv)
|
|---|
| 2559 |
|
|---|
| 2560 | integer p_h2,pf,pc
|
|---|
| 2561 |
|
|---|
| 2562 | l = mg_h1_lmax
|
|---|
| 2563 | p_mg_h2(l,mg_fld) = 0
|
|---|
| 2564 | n = mg_h1_n(l,mg_fld)
|
|---|
| 2565 |
|
|---|
| 2566 | call copy (mg_h2,h2,n) ! Fine grid is just original h2
|
|---|
| 2567 |
|
|---|
| 2568 | nx = mg_nh(l)
|
|---|
| 2569 | ny = mg_nh(l)
|
|---|
| 2570 | nz = mg_nhz(l)
|
|---|
| 2571 |
|
|---|
| 2572 | do l=mg_h1_lmax-1,1,-1
|
|---|
| 2573 |
|
|---|
| 2574 | p_mg_h2(l,mg_fld) = p_mg_h2(l+1,mg_fld) + n
|
|---|
| 2575 | n = mg_h1_n(l ,mg_fld)
|
|---|
| 2576 |
|
|---|
| 2577 | pf = p_mg_h2(l+1,mg_fld)
|
|---|
| 2578 | pc = p_mg_h2(l ,mg_fld)
|
|---|
| 2579 |
|
|---|
| 2580 | call hsmg_intp_fc (mg_h2(pc),mg_h2(pf),l)
|
|---|
| 2581 |
|
|---|
| 2582 | enddo
|
|---|
| 2583 |
|
|---|
| 2584 | p_h2 = p_mg_h2(l0,mg_fld)
|
|---|
| 2585 |
|
|---|
| 2586 | return
|
|---|
| 2587 | end
|
|---|
| 2588 | c-----------------------------------------------------------------------
|
|---|
| 2589 | subroutine hsmg_intp_fc(uc,uf,l) ! l is coarse level
|
|---|
| 2590 |
|
|---|
| 2591 | include 'SIZE'
|
|---|
| 2592 | include 'HSMG'
|
|---|
| 2593 |
|
|---|
| 2594 | real uc(1),uf(1)
|
|---|
| 2595 |
|
|---|
| 2596 |
|
|---|
| 2597 | nc = mg_nh(l)
|
|---|
| 2598 | nf = mg_nh(l+1)
|
|---|
| 2599 | call hsmg_tnsr(uc,nc,uf,nf,mg_jhfc(1,l),mg_jhfct(1,l))
|
|---|
| 2600 |
|
|---|
| 2601 | return
|
|---|
| 2602 | end
|
|---|
| 2603 | c-----------------------------------------------------------------------
|
|---|
| 2604 | subroutine mg_intp_fc_e(uc,uf,nxc,nyc,nzc,nxf,nyf,nzf,e,l,w)
|
|---|
| 2605 | include 'SIZE'
|
|---|
| 2606 | include 'INPUT' ! if3d
|
|---|
| 2607 | include 'HSMG'
|
|---|
| 2608 |
|
|---|
| 2609 | real uf(nxf,nyf,nzf),uc(nxc,nyc,nzc),w(1)
|
|---|
| 2610 | integer e
|
|---|
| 2611 |
|
|---|
| 2612 | if (if3d) then
|
|---|
| 2613 |
|
|---|
| 2614 | n1=nxf*nyf
|
|---|
| 2615 | n2=nzf
|
|---|
| 2616 | n3=nzc
|
|---|
| 2617 | call mxm(uf,n1,mg_jhfct(1,l),n2,w,n3)
|
|---|
| 2618 |
|
|---|
| 2619 | lf=1 ! pointers into work array w()
|
|---|
| 2620 | lc=1 + n1*n3
|
|---|
| 2621 | lc0=lc
|
|---|
| 2622 |
|
|---|
| 2623 | n1=nxf
|
|---|
| 2624 | n2=nyf
|
|---|
| 2625 | n3=nyc
|
|---|
| 2626 |
|
|---|
| 2627 | do k=1,nzc
|
|---|
| 2628 | call mxm(w(lf),n1,mg_jhfct(1,l),n2,w(lc),n3)
|
|---|
| 2629 | lf = lf + n1*n2
|
|---|
| 2630 | lc = lc + n1*n3
|
|---|
| 2631 | enddo
|
|---|
| 2632 |
|
|---|
| 2633 | lf=lc0 ! Rewind fine pointer to start of coarse data
|
|---|
| 2634 | n1=nxc
|
|---|
| 2635 | n2=nxf
|
|---|
| 2636 | n3=nyc*nzc
|
|---|
| 2637 | call mxm(mg_jhfc(1,l),n1,w(lf),n2,uc,n3)
|
|---|
| 2638 |
|
|---|
| 2639 | else ! 2D
|
|---|
| 2640 |
|
|---|
| 2641 | n1=nxf
|
|---|
| 2642 | n2=nyf
|
|---|
| 2643 | n3=nyc
|
|---|
| 2644 | call mxm(uf,n1,mg_jhfct(1,l),n2,w,n3)
|
|---|
| 2645 |
|
|---|
| 2646 | n1=nxc
|
|---|
| 2647 | n2=nxf
|
|---|
| 2648 | n3=nyc
|
|---|
| 2649 | call mxm(mg_jhfc(1,l),n1,w,n2,uc,n3)
|
|---|
| 2650 |
|
|---|
| 2651 | endif
|
|---|
| 2652 |
|
|---|
| 2653 | return
|
|---|
| 2654 | end
|
|---|
| 2655 | c-----------------------------------------------------------------------
|
|---|
| 2656 | subroutine mg_intp_gfc_e(gc,gf,ng,nxc,nyc,nzc,nxf,nyf,nzf,e,l,w)
|
|---|
| 2657 | include 'SIZE'
|
|---|
| 2658 | include 'INPUT' ! if3d
|
|---|
| 2659 | include 'HSMG'
|
|---|
| 2660 |
|
|---|
| 2661 | real gf(ng,nxf,nyf,nzf),gc(ng,nxc,nyc,nzc),w(1)
|
|---|
| 2662 | integer e
|
|---|
| 2663 |
|
|---|
| 2664 |
|
|---|
| 2665 | if (if3d) then
|
|---|
| 2666 |
|
|---|
| 2667 | n1=ng*nxf*nyf
|
|---|
| 2668 | n2=nzf
|
|---|
| 2669 | n3=nzc
|
|---|
| 2670 | call mxm(gf,n1,mg_jhfct(1,l),n2,w,n3)
|
|---|
| 2671 |
|
|---|
| 2672 | lf=1 ! pointers into work array w()
|
|---|
| 2673 | lc=1 + n1*n3
|
|---|
| 2674 | lc0=lc
|
|---|
| 2675 |
|
|---|
| 2676 | n1=ng*nxf
|
|---|
| 2677 | n2=nyf
|
|---|
| 2678 | n3=nyc
|
|---|
| 2679 |
|
|---|
| 2680 | do k=1,nzc
|
|---|
| 2681 | call mxm(w(lf),n1,mg_jhfct(1,l),n2,w(lc),n3)
|
|---|
| 2682 | lf = lf + n1*n2
|
|---|
| 2683 | lc = lc + n1*n3
|
|---|
| 2684 | enddo
|
|---|
| 2685 |
|
|---|
| 2686 | lf=lc0 ! Rewind fine pointer to start of coarse data
|
|---|
| 2687 | n1=ng
|
|---|
| 2688 | n2=nxf
|
|---|
| 2689 | n3=nxc
|
|---|
| 2690 |
|
|---|
| 2691 | do k=1,nyc*nzc
|
|---|
| 2692 | call mxm(w(lf),n1,mg_jhfct(1,l),n2,gc(1,1,k,1),n3)
|
|---|
| 2693 | lf = lf + n1*n2
|
|---|
| 2694 | enddo
|
|---|
| 2695 |
|
|---|
| 2696 | else ! 2D
|
|---|
| 2697 |
|
|---|
| 2698 | n1=ng*nxf
|
|---|
| 2699 | n2=nyf
|
|---|
| 2700 | n3=nyc
|
|---|
| 2701 | call mxm(gf,n1,mg_jhfct(1,l),n2,w,n3)
|
|---|
| 2702 |
|
|---|
| 2703 | lf=1 ! pointers into work array w()
|
|---|
| 2704 |
|
|---|
| 2705 | n1=ng
|
|---|
| 2706 | n2=nxf
|
|---|
| 2707 | n3=nxc
|
|---|
| 2708 |
|
|---|
| 2709 | do k=1,nyc
|
|---|
| 2710 | call mxm(w(lf),n1,mg_jhfct(1,l),n2,gc(1,1,k,1),n3)
|
|---|
| 2711 | lf = lf + n1*n2
|
|---|
| 2712 | enddo
|
|---|
| 2713 |
|
|---|
| 2714 | endif
|
|---|
| 2715 |
|
|---|
| 2716 | return
|
|---|
| 2717 | end
|
|---|
| 2718 | c-----------------------------------------------------------------------
|
|---|
| 2719 | subroutine mg_scale_mass (b,g,wt,ng,nx,ny,nz,wk,ifinv)
|
|---|
| 2720 | include 'SIZE'
|
|---|
| 2721 | include 'INPUT' ! if3d
|
|---|
| 2722 | include 'HSMG'
|
|---|
| 2723 |
|
|---|
| 2724 | real b(1),g(ng,1),wt(1),wk(1)
|
|---|
| 2725 | logical ifinv
|
|---|
| 2726 |
|
|---|
| 2727 | common /ctmp0/ wi(2*lx1+4)
|
|---|
| 2728 |
|
|---|
| 2729 | n = nx*ny*nz
|
|---|
| 2730 |
|
|---|
| 2731 | if (nx.le.2*lx1) then
|
|---|
| 2732 | if (ifinv) then
|
|---|
| 2733 | call invers2(wi,wt,nx)
|
|---|
| 2734 | else
|
|---|
| 2735 | call copy(wi,wt,nx)
|
|---|
| 2736 | endif
|
|---|
| 2737 | else
|
|---|
| 2738 | call exitti('mg_scale_mass: wi too small$',nx)
|
|---|
| 2739 | endif
|
|---|
| 2740 |
|
|---|
| 2741 | if (if3d) then
|
|---|
| 2742 | l=0
|
|---|
| 2743 | do k=1,nz
|
|---|
| 2744 | do j=1,ny
|
|---|
| 2745 | wjk=wi(j)*wi(k)
|
|---|
| 2746 | do i=1,nx
|
|---|
| 2747 | l=l+1
|
|---|
| 2748 | wk(l) = wjk*wi(i)
|
|---|
| 2749 | enddo
|
|---|
| 2750 | enddo
|
|---|
| 2751 | enddo
|
|---|
| 2752 |
|
|---|
| 2753 | do k=1,n
|
|---|
| 2754 | b(k) = wk(k)*b(k)
|
|---|
| 2755 | g(1,k) = wk(k)*g(1,k)
|
|---|
| 2756 | g(2,k) = wk(k)*g(2,k)
|
|---|
| 2757 | g(3,k) = wk(k)*g(3,k)
|
|---|
| 2758 | g(4,k) = wk(k)*g(4,k)
|
|---|
| 2759 | g(5,k) = wk(k)*g(5,k)
|
|---|
| 2760 | g(6,k) = wk(k)*g(6,k)
|
|---|
| 2761 | enddo
|
|---|
| 2762 |
|
|---|
| 2763 | else ! 2D
|
|---|
| 2764 | l=0
|
|---|
| 2765 | do j=1,ny
|
|---|
| 2766 | do i=1,nx
|
|---|
| 2767 | l=l+1
|
|---|
| 2768 | wk(l) = wi(i)*wi(j)
|
|---|
| 2769 | enddo
|
|---|
| 2770 | enddo
|
|---|
| 2771 |
|
|---|
| 2772 | do k=1,n
|
|---|
| 2773 | b(k) = wk(k)*b(k)
|
|---|
| 2774 | g(1,k) = wk(k)*g(1,k)
|
|---|
| 2775 | g(2,k) = wk(k)*g(2,k)
|
|---|
| 2776 | g(3,k) = wk(k)*g(3,k)
|
|---|
| 2777 | enddo
|
|---|
| 2778 |
|
|---|
| 2779 | endif
|
|---|
| 2780 |
|
|---|
| 2781 | return
|
|---|
| 2782 | end
|
|---|
| 2783 | c-----------------------------------------------------------------------
|
|---|
| 2784 | subroutine mg_set_gb (p_g,p_b,l0)
|
|---|
| 2785 | include 'SIZE'
|
|---|
| 2786 | include 'HSMG'
|
|---|
| 2787 | include 'MASS' ! bm1
|
|---|
| 2788 | include 'TSTEP' ! nelfld
|
|---|
| 2789 |
|
|---|
| 2790 | integer p_g,p_b,e
|
|---|
| 2791 | common /ctmp1/ w(lx1*ly1*lz1*lelt*2)
|
|---|
| 2792 |
|
|---|
| 2793 | l = mg_h1_lmax
|
|---|
| 2794 | p_mg_b (l,mg_fld) = 0
|
|---|
| 2795 | p_mg_g (l,mg_fld) = 0
|
|---|
| 2796 | n = mg_h1_n(l,mg_fld)
|
|---|
| 2797 |
|
|---|
| 2798 |
|
|---|
| 2799 | ng = 3*(ldim-1) ! 3 or 6 elements to symm dxd tensor
|
|---|
| 2800 |
|
|---|
| 2801 | do l=mg_h1_lmax-1,1,-1
|
|---|
| 2802 |
|
|---|
| 2803 | p_mg_b (l,mg_fld) = p_mg_b (l+1,mg_fld) + n
|
|---|
| 2804 | p_mg_g (l,mg_fld) = p_mg_g (l+1,mg_fld) + n*ng
|
|---|
| 2805 | n = mg_h1_n(l ,mg_fld)
|
|---|
| 2806 |
|
|---|
| 2807 | enddo
|
|---|
| 2808 |
|
|---|
| 2809 | do e=1,nelfld(ifield)
|
|---|
| 2810 | do l=mg_h1_lmax,1,-1
|
|---|
| 2811 |
|
|---|
| 2812 | nx = mg_nh(l)
|
|---|
| 2813 | ny = mg_nh(l)
|
|---|
| 2814 | nz = mg_nhz(l)
|
|---|
| 2815 | nxyz = nx*ny*nz
|
|---|
| 2816 |
|
|---|
| 2817 | p_g = p_mg_g (l,mg_fld) + ng*nx*ny*nz*(e-1)
|
|---|
| 2818 | p_b = p_mg_b (l,mg_fld) + nx*ny*nz*(e-1)
|
|---|
| 2819 |
|
|---|
| 2820 | if (l.eq.mg_h1_lmax) then
|
|---|
| 2821 | call gxfer_e (mg_g(p_g) ,ng,e ) ! Fine grid=original G
|
|---|
| 2822 | call copy (mg_b(p_b) ,bm1(1,1,1,e),nxyz) ! Fine grid=original B
|
|---|
| 2823 | call mg_scale_mass ! Divide out Wghts
|
|---|
| 2824 | $ (mg_b(p_b),mg_g(p_g),mg_bh(1,l),ng,nx,ny,nz,w,.true.)
|
|---|
| 2825 | else
|
|---|
| 2826 |
|
|---|
| 2827 | c Generate G and B by interpolating their continous counterparts onto
|
|---|
| 2828 | c the coarse grid and collocating with coarse-grid quadrature weights
|
|---|
| 2829 |
|
|---|
| 2830 | call mg_intp_gfc_e
|
|---|
| 2831 | $ (mg_g(p_g),mg_g(l_g),ng,nx,ny,nz,nxl,nyl,nzl,e,l,w)
|
|---|
| 2832 |
|
|---|
| 2833 | call mg_intp_fc_e
|
|---|
| 2834 | $ (mg_b(p_b),mg_b(l_b) ,nx,ny,nz,nxl,nyl,nzl,e,l,w)
|
|---|
| 2835 |
|
|---|
| 2836 | call mg_scale_mass ! Reinstate weights
|
|---|
| 2837 | $ (mg_b(l_b),mg_g(l_g),mg_bh(1,l+1),ng,nxl,nyl,nzl,w,.false.)
|
|---|
| 2838 |
|
|---|
| 2839 | endif
|
|---|
| 2840 |
|
|---|
| 2841 | l_b = p_b
|
|---|
| 2842 | l_g = p_g
|
|---|
| 2843 |
|
|---|
| 2844 | nxl = nx
|
|---|
| 2845 | nyl = ny
|
|---|
| 2846 | nzl = nz
|
|---|
| 2847 |
|
|---|
| 2848 | enddo
|
|---|
| 2849 |
|
|---|
| 2850 | call mg_scale_mass ! Reinstate weights
|
|---|
| 2851 | $ (mg_b(l_b),mg_g(l_g),mg_bh(1,1),ng,nxl,nyl,nzl,w,.false.)
|
|---|
| 2852 |
|
|---|
| 2853 |
|
|---|
| 2854 | enddo
|
|---|
| 2855 |
|
|---|
| 2856 | p_b = p_mg_b (l0,mg_fld)
|
|---|
| 2857 | p_g = p_mg_g (l0,mg_fld)
|
|---|
| 2858 |
|
|---|
| 2859 | return
|
|---|
| 2860 | end
|
|---|
| 2861 | c-----------------------------------------------------------------------
|
|---|
| 2862 | subroutine gxfer_e (g,ng,e)
|
|---|
| 2863 | include 'SIZE'
|
|---|
| 2864 | include 'TOTAL'
|
|---|
| 2865 |
|
|---|
| 2866 | real g(ng,1)
|
|---|
| 2867 | integer e
|
|---|
| 2868 |
|
|---|
| 2869 | nxyz = lx1*ly1*lz1
|
|---|
| 2870 |
|
|---|
| 2871 | c ifdfrm(e) = .true. ! TOO LATE
|
|---|
| 2872 |
|
|---|
| 2873 | if (if3d) then
|
|---|
| 2874 | do i=1,nxyz
|
|---|
| 2875 | g(1,i) = g1m1(i,1,1,e)
|
|---|
| 2876 | g(2,i) = g2m1(i,1,1,e)
|
|---|
| 2877 | g(3,i) = g3m1(i,1,1,e)
|
|---|
| 2878 | g(4,i) = g4m1(i,1,1,e)
|
|---|
| 2879 | g(5,i) = g5m1(i,1,1,e)
|
|---|
| 2880 | g(6,i) = g6m1(i,1,1,e)
|
|---|
| 2881 | enddo
|
|---|
| 2882 | else
|
|---|
| 2883 | do i=1,nxyz
|
|---|
| 2884 | g(1,i) = g1m1(i,1,1,e)
|
|---|
| 2885 | g(2,i) = g2m1(i,1,1,e)
|
|---|
| 2886 | g(3,i) = g4m1(i,1,1,e)
|
|---|
| 2887 | enddo
|
|---|
| 2888 | endif
|
|---|
| 2889 |
|
|---|
| 2890 | return
|
|---|
| 2891 | end
|
|---|
| 2892 | c-----------------------------------------------------------------------
|
|---|
| 2893 | subroutine chkr(name3,ii)
|
|---|
| 2894 | include 'SIZE'
|
|---|
| 2895 | include 'TOTAL'
|
|---|
| 2896 | include 'HSMG'
|
|---|
| 2897 | character*3 name3
|
|---|
| 2898 |
|
|---|
| 2899 | write(6,*) mg_h1_lmax,ii,' ',name3,' CHKR'
|
|---|
| 2900 |
|
|---|
| 2901 | return
|
|---|
| 2902 | end
|
|---|
| 2903 | c-----------------------------------------------------------------------
|
|---|
| 2904 | subroutine outgmat(a,ng,nx,ny,name6,k,e)
|
|---|
| 2905 |
|
|---|
| 2906 | integer e
|
|---|
| 2907 | real a(ng,nx,ny)
|
|---|
| 2908 | common /ctmp0/ w(100000)
|
|---|
| 2909 | character*6 name6
|
|---|
| 2910 |
|
|---|
| 2911 | c do i=1,ng
|
|---|
| 2912 | do i=1,1
|
|---|
| 2913 | sum = 0.
|
|---|
| 2914 | do ii=1,nx*ny
|
|---|
| 2915 | w(ii)=a(i,ii,1)
|
|---|
| 2916 | sum = sum + a(i,ii,1)
|
|---|
| 2917 | enddo
|
|---|
| 2918 |
|
|---|
| 2919 | write(6,1) name6,i,k,e,nx,ny,ng,sum
|
|---|
| 2920 | 1 format(a6,6i5,f12.5,' outgmat')
|
|---|
| 2921 |
|
|---|
| 2922 | call outmatz(w,nx,ny,name6,i)
|
|---|
| 2923 | enddo
|
|---|
| 2924 |
|
|---|
| 2925 | return
|
|---|
| 2926 | end
|
|---|
| 2927 | c-----------------------------------------------------------------------
|
|---|
| 2928 | subroutine outmatz(a,m,n,name6,ie)
|
|---|
| 2929 | real a(m,n)
|
|---|
| 2930 | character*6 name6
|
|---|
| 2931 |
|
|---|
| 2932 | sum=0.
|
|---|
| 2933 | sua=0.
|
|---|
| 2934 | do i=1,m*n
|
|---|
| 2935 | sum=sum+ a(i,1)
|
|---|
| 2936 | sua=sua+abs(a(i,1))
|
|---|
| 2937 | enddo
|
|---|
| 2938 | sum=sum/(m*n)
|
|---|
| 2939 | sua=sua/(m*n)
|
|---|
| 2940 |
|
|---|
| 2941 | write(6,*)
|
|---|
| 2942 | write(6,1) ie,name6,m,n,sum,sua
|
|---|
| 2943 | 1 format(i8,' matrix: ',a6,2i5,1p2e12.4)
|
|---|
| 2944 |
|
|---|
| 2945 | n12 = min(m,12)
|
|---|
| 2946 | do j=m,1,-1
|
|---|
| 2947 | write(6,6) ie,name6,(a(i,j),i=1,n12)
|
|---|
| 2948 | enddo
|
|---|
| 2949 | 6 format(i3,1x,a6,12f9.5)
|
|---|
| 2950 | c write(6,*)
|
|---|
| 2951 | return
|
|---|
| 2952 | end
|
|---|
| 2953 | c-----------------------------------------------------------------------
|
|---|
| 2954 | subroutine h1mg_setup_schwarz_wt_2(wt,ie,n,work,ifsqrt)
|
|---|
| 2955 | include 'SIZE'
|
|---|
| 2956 | real wt(1),work(1)
|
|---|
| 2957 | logical ifsqrt
|
|---|
| 2958 |
|
|---|
| 2959 | if (ldim.eq.2) call h1mg_setup_schwarz_wt2d_2(wt,ie,n,work,ifsqrt)
|
|---|
| 2960 | if (ldim.eq.3) call h1mg_setup_schwarz_wt3d_2(wt,ie,n,work,ifsqrt)
|
|---|
| 2961 |
|
|---|
| 2962 | return
|
|---|
| 2963 | end
|
|---|
| 2964 | c----------------------------------------------------------------------
|
|---|
| 2965 | subroutine h1mg_setup_schwarz_wt2d_2(wt,ie,n,work,ifsqrt)
|
|---|
| 2966 | include 'SIZE'
|
|---|
| 2967 | logical ifsqrt
|
|---|
| 2968 | integer n
|
|---|
| 2969 | real wt(n,4,2,nelv)
|
|---|
| 2970 | real work(n,n)
|
|---|
| 2971 |
|
|---|
| 2972 | integer ie,i,j
|
|---|
| 2973 | do j=1,n
|
|---|
| 2974 | wt(j,1,1,ie)=1.0/work(1,j)
|
|---|
| 2975 | wt(j,2,1,ie)=1.0/work(2,j)
|
|---|
| 2976 | wt(j,3,1,ie)=1.0/work(n-1,j)
|
|---|
| 2977 | wt(j,4,1,ie)=1.0/work(n,j)
|
|---|
| 2978 | enddo
|
|---|
| 2979 | do i=1,n
|
|---|
| 2980 | wt(i,1,2,ie)=1.0/work(i,1)
|
|---|
| 2981 | wt(i,2,2,ie)=1.0/work(i,2)
|
|---|
| 2982 | wt(i,3,2,ie)=1.0/work(i,n-1)
|
|---|
| 2983 | wt(i,4,2,ie)=1.0/work(i,n)
|
|---|
| 2984 | enddo
|
|---|
| 2985 | if(ifsqrt) then
|
|---|
| 2986 | do ii=1,2
|
|---|
| 2987 | do j=1,4
|
|---|
| 2988 | do i=1,n
|
|---|
| 2989 | wt(i,j,ii,ie)=sqrt(wt(i,j,ii,ie))
|
|---|
| 2990 | enddo
|
|---|
| 2991 | enddo
|
|---|
| 2992 | enddo
|
|---|
| 2993 | endif
|
|---|
| 2994 |
|
|---|
| 2995 | return
|
|---|
| 2996 | end
|
|---|
| 2997 | c----------------------------------------------------------------------
|
|---|
| 2998 | subroutine h1mg_setup_schwarz_wt3d_2(wt,ie,n,work,ifsqrt)
|
|---|
| 2999 | include 'SIZE'
|
|---|
| 3000 | logical ifsqrt
|
|---|
| 3001 | integer n
|
|---|
| 3002 | real wt(n,n,4,3,nelv)
|
|---|
| 3003 | real work(n,n,n)
|
|---|
| 3004 |
|
|---|
| 3005 | integer ie,i,j,k
|
|---|
| 3006 | integer lbr,rbr,lbs,rbs,lbt,rbt
|
|---|
| 3007 |
|
|---|
| 3008 | ierr = 0
|
|---|
| 3009 | do k=1,n
|
|---|
| 3010 | do j=1,n
|
|---|
| 3011 | wt(j,k,1,1,ie)=1.0/work(1,j,k)
|
|---|
| 3012 | wt(j,k,2,1,ie)=1.0/work(2,j,k)
|
|---|
| 3013 | wt(j,k,3,1,ie)=1.0/work(n-1,j,k)
|
|---|
| 3014 | wt(j,k,4,1,ie)=1.0/work(n,j,k)
|
|---|
| 3015 | enddo
|
|---|
| 3016 | enddo
|
|---|
| 3017 | do k=1,n
|
|---|
| 3018 | do i=1,n
|
|---|
| 3019 | wt(i,k,1,2,ie)=1.0/work(i,1,k)
|
|---|
| 3020 | wt(i,k,2,2,ie)=1.0/work(i,2,k)
|
|---|
| 3021 | wt(i,k,3,2,ie)=1.0/work(i,n-1,k)
|
|---|
| 3022 | wt(i,k,4,2,ie)=1.0/work(i,n,k)
|
|---|
| 3023 | enddo
|
|---|
| 3024 | enddo
|
|---|
| 3025 | do j=1,n
|
|---|
| 3026 | do i=1,n
|
|---|
| 3027 | wt(i,j,1,3,ie)=1.0/work(i,j,1)
|
|---|
| 3028 | wt(i,j,2,3,ie)=1.0/work(i,j,2)
|
|---|
| 3029 | wt(i,j,3,3,ie)=1.0/work(i,j,n-1)
|
|---|
| 3030 | wt(i,j,4,3,ie)=1.0/work(i,j,n)
|
|---|
| 3031 | enddo
|
|---|
| 3032 | enddo
|
|---|
| 3033 | if(ifsqrt) then
|
|---|
| 3034 | do ii=1,3
|
|---|
| 3035 | do k=1,4
|
|---|
| 3036 | do j=1,4
|
|---|
| 3037 | do i=1,n
|
|---|
| 3038 | wt(i,j,k,ii,ie)=sqrt(wt(i,j,k,ii,ie))
|
|---|
| 3039 | enddo
|
|---|
| 3040 | enddo
|
|---|
| 3041 | enddo
|
|---|
| 3042 | enddo
|
|---|
| 3043 | endif
|
|---|
| 3044 |
|
|---|
| 3045 | return
|
|---|
| 3046 | end
|
|---|
| 3047 | c----------------------------------------------------------------------
|
|---|
| 3048 | subroutine h1mg_setup_schwarz_wt_1(wt,l,ifsqrt)
|
|---|
| 3049 | include 'SIZE'
|
|---|
| 3050 | include 'INPUT' ! if3d
|
|---|
| 3051 | include 'TSTEP' ! ifield
|
|---|
| 3052 | include 'HSMG'
|
|---|
| 3053 |
|
|---|
| 3054 | real wt(1),work(1)
|
|---|
| 3055 | logical ifsqrt
|
|---|
| 3056 |
|
|---|
| 3057 | integer enx,eny,enz,pm
|
|---|
| 3058 |
|
|---|
| 3059 | zero = 0
|
|---|
| 3060 | one = 1
|
|---|
| 3061 | onem = -1
|
|---|
| 3062 |
|
|---|
| 3063 | n = mg_h1_n (l,mg_fld)
|
|---|
| 3064 | pm = p_mg_msk(l,mg_fld)
|
|---|
| 3065 |
|
|---|
| 3066 | enx=mg_nh(l)+2
|
|---|
| 3067 | eny=mg_nh(l)+2
|
|---|
| 3068 | enz=mg_nh(l)+2
|
|---|
| 3069 | if(.not.if3d) enz=1
|
|---|
| 3070 | ns = enx*eny*enz*nelfld(ifield)
|
|---|
| 3071 | i = ns+1
|
|---|
| 3072 |
|
|---|
| 3073 | call rone(mg_work(i),ns)
|
|---|
| 3074 |
|
|---|
| 3075 | c Sum overlap region (border excluded)
|
|---|
| 3076 | call hsmg_extrude(mg_work,0,zero,mg_work(i),0,one ,enx,eny,enz)
|
|---|
| 3077 | call hsmg_schwarz_dssum(mg_work(i),l)
|
|---|
| 3078 | call hsmg_extrude(mg_work(i),0,one ,mg_work,0,onem,enx,eny,enz)
|
|---|
| 3079 | call hsmg_extrude(mg_work(i),2,one,mg_work(i),0,one,enx,eny,enz)
|
|---|
| 3080 |
|
|---|
| 3081 | if(.not.if3d) then ! Go back to regular size array
|
|---|
| 3082 | call hsmg_schwarz_toreg2d(mg_work,mg_work(i),mg_nh(l))
|
|---|
| 3083 | else
|
|---|
| 3084 | call hsmg_schwarz_toreg3d(mg_work,mg_work(i),mg_nh(l))
|
|---|
| 3085 | endif
|
|---|
| 3086 |
|
|---|
| 3087 | call hsmg_dssum(mg_work,l) ! sum border nodes
|
|---|
| 3088 |
|
|---|
| 3089 |
|
|---|
| 3090 | nx = mg_nh(l)
|
|---|
| 3091 | ny = mg_nh(l)
|
|---|
| 3092 | nz = mg_nh(l)
|
|---|
| 3093 | if (.not.if3d) nz=1
|
|---|
| 3094 | nxyz = nx*ny*nz
|
|---|
| 3095 | k = 1
|
|---|
| 3096 | do ie=1,nelfld(ifield)
|
|---|
| 3097 | c call outmat(mg_work(k),nx,ny,'NEW WT',ie)
|
|---|
| 3098 | call h1mg_setup_schwarz_wt_2(wt,ie,nx,mg_work(k),ifsqrt)
|
|---|
| 3099 | k = k+nxyz
|
|---|
| 3100 | enddo
|
|---|
| 3101 | c stop
|
|---|
| 3102 |
|
|---|
| 3103 | return
|
|---|
| 3104 | end
|
|---|
| 3105 | c----------------------------------------------------------------------
|
|---|
