source: CIVL/examples/fortran/nek5000/core/subs1.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100755
File size: 72.0 KB
Line 
1c-----------------------------------------------------------------------
2 subroutine cggosf (u1,u2,u3,r1,r2,r3,h1,h2,rmult,binv,
3 $ vol,tin,maxit,matmod)
4
5C Conjugate gradient iteration for solution of coupled
6C Helmholtz equations
7
8 include 'SIZE'
9 include 'TOTAL'
10 include 'DOMAIN'
11 include 'FDMH1'
12
13 common /screv/ dpc(lx1*ly1*lz1*lelt)
14 $ , p1 (lx1*ly1*lz1*lelt)
15 common /scrch/ p2 (lx1*ly1*lz1*lelt)
16 $ , p3 (lx1*ly1*lz1*lelt)
17 common /scrsl/ qq1(lx1*ly1*lz1*lelt)
18 $ , qq2(lx1*ly1*lz1*lelt)
19 $ , qq3(lx1*ly1*lz1*lelt)
20 common /scrmg/ pp1(lx1*ly1*lz1*lelt)
21 $ , pp2(lx1*ly1*lz1*lelt)
22 $ , pp3(lx1*ly1*lz1*lelt)
23 $ , wa (lx1*ly1*lz1*lelt)
24 real ap1(1),ap2(1),ap3(1)
25 equivalence (ap1,pp1),(ap2,pp2),(ap3,pp3)
26
27 common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
28 common /cprint/ ifprint
29 logical ifdfrm, iffast, ifh2, ifsolv, ifprint
30
31 real u1(1),u2(1),u3(1),
32 $ r1(1),r2(1),r3(1),h1(1),h2(1),rmult(1),binv(1)
33
34
35 logical iffdm,ifcrsl
36
37 iffdm = .true.
38 iffdm = .false.
39c ifcrsl = .true.
40 ifcrsl = .false.
41
42 nel = nelfld(ifield)
43 nxyz = lx1*ly1*lz1
44 n = nxyz*nel
45
46 if (istep.le.1.and.iffdm) call set_fdm_prec_h1A
47
48 tol = tin
49
50c overrule input tolerance
51 if (restol(ifield).ne.0) tol=restol(ifield)
52
53 if (ifcrsl) call set_up_h1_crs_strs(h1,h2,ifield,matmod)
54
55c if (nio.eq.0.and.istep.eq.1) write(6,6) ifield,tol,tin
56c 6 format(i3,1p2e12.4,' ifield, tol, tol_in')
57
58 if ( .not.ifsolv ) then ! Set logical flags
59 call setfast (h1,h2,imesh)
60 ifsolv = .true.
61 endif
62
63 call opdot (wa,r1,r2,r3,r1,r2,r3,n)
64 rbnorm = glsc3(wa,binv,rmult,n)
65 rbnorm = sqrt ( rbnorm / vol )
66 if (rbnorm .lt. tol**2) then
67 iter = 0
68 r0 = rbnorm
69c if ( .not.ifprint ) goto 9999
70 if (matmod.ge.0.and.nio.eq.0) write (6,3000)
71 $ istep,iter,rbnorm,r0,tol
72 if (matmod.lt.0.and.nio.eq.0) write (6,3010)
73 $ istep,iter,rbnorm,r0,tol
74 goto 9999
75 endif
76
77C Evaluate diagional pre-conidtioner for fluid solve
78 call setprec (dpc,h1,h2,imesh,1)
79 call setprec (wa ,h1,h2,imesh,2)
80 call add2 (dpc,wa,n)
81 if (ldim.eq.3) then
82 call setprec (wa,h1,h2,imesh,3)
83 call add2 (dpc,wa,n)
84 endif
85c call rone (dpc,n)
86c call copy (dpc,binv,n)
87
88 if (iffdm) then
89 call set_fdm_prec_h1b(dpc,h1,h2,nel)
90 call fdm_h1a (pp1,r1,dpc,nel,ktype(1,1,1),wa)
91 call fdm_h1a (pp2,r2,dpc,nel,ktype(1,1,2),wa)
92 call fdm_h1a (pp3,r3,dpc,nel,ktype(1,1,3),wa)
93 call rmask (pp1,pp2,pp3,nel)
94 call opdssum (pp1,pp2,pp3)
95 else
96 call col3 (pp1,dpc,r1,n)
97 call col3 (pp2,dpc,r2,n)
98 if (if3d) call col3 (pp3,dpc,r3,n)
99 endif
100 if (ifcrsl) then
101 call crs_strs(p1,p2,p3,r1,r2,r3)
102 call rmask (p1,p2,p3,nel)
103 else
104 call opzero(p1,p2,p3)
105 endif
106 call opadd2 (p1,p2,p3,pp1,pp2,pp3)
107 rpp1 = op_glsc2_wt(p1,p2,p3,r1,r2,r3,rmult)
108
109 maxit=200
110 do 1000 iter=1,maxit
111 call axhmsf (ap1,ap2,ap3,p1,p2,p3,h1,h2,matmod)
112 call rmask (ap1,ap2,ap3,nel)
113 call opdssum (ap1,ap2,ap3)
114 pap = op_glsc2_wt(p1,p2,p3,ap1,ap2,ap3,rmult)
115 alpha = rpp1 / pap
116
117 call opadds (u1,u2,u3,p1 ,p2 ,p3 , alpha,n,2)
118 call opadds (r1,r2,r3,ap1,ap2,ap3,-alpha,n,2)
119
120 call opdot (wa,r1,r2,r3,r1,r2,r3,n)
121 rbnorm = glsc3(wa,binv,rmult,n)
122 rbnorm = sqrt (rbnorm/vol)
123
124 if (iter.eq.1) r0 = rbnorm
125
126 if (rbnorm.lt.tol) then
127 ifin = iter
128 if (nio.eq.0) then
129 if (matmod.ge.0) write(6,3000) istep,ifin,rbnorm,r0,tol
130 if (matmod.lt.0) write(6,3010) istep,ifin,rbnorm,r0,tol
131 endif
132 goto 9999
133 endif
134
135 if (iffdm) then
136 call fdm_h1a (pp1,r1,dpc,nel,ktype(1,1,1),wa)
137 call fdm_h1a (pp2,r2,dpc,nel,ktype(1,1,2),wa)
138 call fdm_h1a (pp3,r3,dpc,nel,ktype(1,1,3),wa)
139 call rmask (pp1,pp2,pp3,nel)
140 call opdssum (pp1,pp2,pp3)
141 else
142 call col3 (pp1,dpc,r1,n)
143 call col3 (pp2,dpc,r2,n)
144 if (if3d) call col3 (pp3,dpc,r3,n)
145 endif
146
147 if (ifcrsl) then
148 call crs_strs(qq1,qq2,qq3,r1,r2,r3)
149 call rmask (qq1,qq2,qq3,nel)
150 call opadd2 (pp1,pp2,pp3,qq1,qq2,qq3)
151 endif
152
153 call opdot (wa,r1,r2,r3,pp1,pp2,pp3,n)
154
155 rpp2 = rpp1
156 rpp1 = glsc2(wa,rmult,n)
157 beta = rpp1/rpp2
158 call opadds (p1,p2,p3,pp1,pp2,pp3,beta,n,1)
159
160 1000 continue
161 if (matmod.ge.0.and.nio.eq.0) write (6,3001)
162 $ istep,iter,rbnorm,r0,tol
163 if (matmod.lt.0.and.nio.eq.0) write (6,3011)
164 $ istep,iter,rbnorm,r0,tol
165
166 9999 continue
167 ifsolv = .false.
168
169
170 3000 format(i11,' Helmh3 fluid ',I6,1p3E13.4)
171 3010 format(i11,' Helmh3 mesh ',I6,1p3E13.4)
172 3001 format(i11,' Helmh3 fluid unconverged! ',I6,1p3E13.4)
173 3011 format(i11,' Helmh3 mesh unconverged! ',I6,1p3E13.4)
174
175 return
176 end
177c-----------------------------------------------------------------------
178 subroutine setdt
179c
180c Set the new time step. All cases covered.
181c
182 include 'SIZE'
183 include 'SOLN'
184 include 'MVGEOM'
185 include 'INPUT'
186 include 'TSTEP'
187 include 'PARALLEL'
188
189 common /scruz/ cx(lx1*ly1*lz1*lelt)
190 $ , cy(lx1,ly1,lz1,lelt)
191 $ , cz(lx1,ly1,lz1,lelt)
192
193 common /cprint/ ifprint
194 logical ifprint
195 common /udxmax/ umax
196 REAL DTOLD
197 SAVE DTOLD
198 DATA DTOLD /0.0/
199 REAL DTOpf
200 SAVE DTOpf
201 DATA DTOpf /0.0/
202 logical iffxdt
203 save iffxdt
204 data iffxdt /.false./
205C
206
207 if (param(12).lt.0.or.iffxdt) then
208 iffxdt = .true.
209 param(12) = abs(param(12))
210 dt = param(12)
211 dtopf = dt
212 if (ifmvbd) then
213 call opsub3 (cx,cy,cz,vx,vy,vz,wx,wy,wz)
214 call compute_cfl(umax,cx,cy,cz,1.0)
215 else
216 call compute_cfl(umax,vx,vy,vz,1.0)
217 endif
218
219 goto 200
220 else IF (PARAM(84).NE.0.0) THEN
221 if (dtold.eq.0.0) then
222 dt =param(84)
223 dtold=param(84)
224 dtopf=param(84)
225 return
226 else
227 dtold=dt
228 dtopf=dt
229 dt=dtopf*param(85)
230 dt=min(dt,param(12))
231 endif
232 endif
233
234C
235C Find DT=DTCFL based on CFL-condition (if applicable)
236C
237 CALL SETDTC
238 DTCFL = DT
239C
240C Find DTFS based on surface tension (if applicable)
241C
242 CALL SETDTFS (DTFS)
243C
244C Select appropriate DT
245C
246 IF ((DT.EQ.0.).AND.(DTFS.GT.0.)) THEN
247 DT = DTFS
248 ELSEIF ((DT.GT.0.).AND.(DTFS.GT.0.)) THEN
249 DT = MIN(DT,DTFS)
250 ELSEIF ((DT.EQ.0.).AND.(DTFS.EQ.0.)) THEN
251 DT = 0.
252 IF (IFFLOW.AND.NID.EQ.0.AND.IFPRINT) THEN
253 WRITE (6,*) 'WARNING: CFL-condition & surface tension'
254 WRITE (6,*) ' are not applicable'
255 endif
256 ELSEIF ((DT.GT.0.).AND.(DTFS.EQ.0.)) THEN
257 DT = DT
258 ELSE
259 DT = 0.
260 IF (NIO.EQ.0) WRITE (6,*) 'WARNING: DT<0 or DTFS<0'
261 IF (NIO.EQ.0) WRITE (6,*) ' Reset DT '
262 endif
263C
264C Check DT against user-specified input, DTINIT=PARAM(12).
265C
266 IF ((DT.GT.0.).AND.(DTINIT.GT.0.)) THEN
267 DT = MIN(DT,DTINIT)
268 ELSEIF ((DT.EQ.0.).AND.(DTINIT.GT.0.)) THEN
269 DT = DTINIT
270 ELSEIF ((DT.GT.0.).AND.(DTINIT.EQ.0.)) THEN
271 DT = DT
272 ELSEIF (.not.iffxdt) THEN
273 DT = 0.001
274 IF(NIO.EQ.0)WRITE (6,*) 'WARNING: Set DT=0.001 (arbitrarily)'
275 endif
276C
277C Check if final time (user specified) has been reached.
278C
279 200 IF (TIME+DT .GE. FINTIM .AND. FINTIM.NE.0.0) THEN
280C Last step
281 LASTEP = 1
282 DT = FINTIM-TIME
283 IF (NIO.EQ.0) WRITE (6,*) 'Final time step = ',DT
284 endif
285C
286 COURNO = DT*UMAX
287 IF (NIO.EQ.0.AND.IFPRINT.AND.DT.NE.DTOLD)
288 $ WRITE (6,100) DT,DTCFL,DTFS,DTINIT
289 100 FORMAT(5X,'DT/DTCFL/DTFS/DTINIT',4E12.3)
290C
291C Put limits on how much DT can change.
292C
293 IF (DTOLD.NE.0.0 .AND. LASTEP.NE.1) THEN
294 DTMIN=0.8*DTOLD
295 DTMAX=1.2*DTOLD
296 DT = MIN(DTMAX,DT)
297 DT = MAX(DTMIN,DT)
298 endif
299 DTOLD=DT
300
301C IF (PARAM(84).NE.0.0) THEN
302C dt=dtopf*param(85)
303C dt=min(dt,param(12))
304C endif
305
306 if (iffxdt) dt=dtopf
307 COURNO = DT*UMAX
308
309! synchronize time step for multiple sessions
310 if (ifneknek) dt = glmin_ms(dt,1)
311c
312 if (iffxdt.and.abs(courno).gt.10.*abs(ctarg)) then
313 if (nid.eq.0) write(6,*) 'CFL, Ctarg!',courno,ctarg
314 call emerxit
315 endif
316
317
318 return
319 end
320C
321C--------------------------------------------------------
322C
323 subroutine cvgnlps (ifconv)
324C----------------------------------------------------------------------
325C
326C Check convergence for non-linear passisve scalar solver.
327C Relevant for solving heat transport problems with radiation b.c.
328C
329C----------------------------------------------------------------------
330 include 'SIZE'
331 include 'INPUT'
332 include 'TSTEP'
333 LOGICAL IFCONV
334C
335 IF (IFNONL(IFIELD)) THEN
336 IFCONV = .FALSE.
337 ELSE
338 IFCONV = .TRUE.
339 return
340 endif
341C
342 TNORM1 = TNRMH1(IFIELD-1)
343 CALL UNORM
344 TNORM2 = TNRMH1(IFIELD-1)
345 EPS = ABS((TNORM2-TNORM1)/TNORM2)
346 IF (EPS .LT. TOLNL) IFCONV = .TRUE.
347C
348 return
349 end
350C
351 subroutine unorm
352C---------------------------------------------------------------------
353C
354C Norm calculation.
355C
356C---------------------------------------------------------------------
357 include 'SIZE'
358 include 'SOLN'
359 include 'TSTEP'
360C
361 IF (IFIELD.EQ.1) THEN
362C
363C Compute norms of the velocity.
364C Compute time mean (L2) of the inverse of the time step.
365C Compute L2 in time, H1 in space of the velocity.
366C
367 CALL NORMVC (VNRMH1,VNRMSM,VNRML2,VNRML8,VX,VY,VZ)
368 IF (ISTEP.EQ.0) return
369 IF (ISTEP.EQ.1) THEN
370 DTINVM = 1./DT
371 VMEAN = VNRML8
372 ELSE
373 tden = time
374 if (time.le.0) tden = abs(time)+1.e-9
375 arg = ((TIME-DT)*DTINVM**2+1./DT)/tden
376 if (arg.gt.0) DTINVM = SQRT(arg)
377 arg = ((TIME-DT)*VMEAN**2+DT*VNRMH1**2)/tden
378 if (arg.gt.0) VMEAN = SQRT(arg)
379 endif
380 ELSE
381C
382C Compute norms of a passive scalar
383C
384 CALL NORMSC (TNRMH1(IFIELD-1),TNRMSM(IFIELD-1),
385 $ TNRML2(IFIELD-1),TNRML8(IFIELD-1),
386 $ T(1,1,1,1,IFIELD-1),IMESH)
387 TMEAN(IFIELD-1) = 0.
388 endif
389C
390 return
391 end
392C
393 subroutine chktmg (tol,res,w1,w2,mult,mask,imesh)
394C-------------------------------------------------------------------
395C
396C Check that the tolerances are not too small for the MG-solver.
397C Important when calling the MG-solver (Gauss-Lobatto mesh).
398C Note: direct stiffness summation
399C
400C-------------------------------------------------------------------
401 include 'SIZE'
402 include 'INPUT'
403 include 'MASS'
404 include 'EIGEN'
405C
406 REAL RES (LX1,LY1,LZ1,1)
407 REAL W1 (LX1,LY1,LZ1,1)
408 REAL W2 (LX1,LY1,LZ1,1)
409 REAL MULT (LX1,LY1,LZ1,1)
410 REAL MASK (LX1,LY1,LZ1,1)
411C
412C Single or double precision???
413C
414 DELTA = 1.E-9
415 X = 1.+DELTA
416 Y = 1.
417 DIFF = ABS(X-Y)
418 IF (DIFF.EQ.0.) EPS = 1.E-6*EIGGA/EIGAA
419 IF (DIFF.GT.0.) EPS = 1.E-13*EIGGA/EIGAA
420C
421 IF (IMESH.EQ.1) NL = NELV
422 IF (IMESH.EQ.2) NL = NELT
423 NTOT1 = lx1*ly1*lz1*NL
424 CALL COPY (W1,RES,NTOT1)
425C
426 CALL DSSUM (W1,lx1,ly1,lz1)
427C
428 IF (IMESH.EQ.1) THEN
429 CALL COL3 (W2,BINVM1,W1,NTOT1)
430 RINIT = SQRT(GLSC3 (W2,W1,MULT,NTOT1)/VOLVM1)
431 ELSE
432 CALL COL3 (W2,BINTM1,W1,NTOT1)
433 RINIT = SQRT(GLSC3 (W2,W1,MULT,NTOT1)/VOLTM1)
434 endif
435 RMIN = EPS*RINIT
436 IF (TOL.LT.RMIN) THEN
437 TOLOLD=TOL
438 TOL = RMIN
439 IF (NIO.EQ.0)
440 $ WRITE (6,*) 'New MG-tolerance (RINIT*epsm*cond) = ',TOL,TOLOLD
441 endif
442C
443 CALL RONE (W1,NTOT1)
444 BCNEU1 = GLSC3(W1,MASK,MULT,NTOT1)
445 BCNEU2 = GLSC3(W1,W1 ,MULT,NTOT1)
446 BCTEST = ABS(BCNEU1-BCNEU2)
447 IF (BCTEST .LT. .1) THEN
448 OTR = GLSUM (RES,NTOT1)
449 TOLMIN = ABS(OTR)*EIGGA/EIGAA
450 IF (TOL .LT. TOLMIN) THEN
451 TOLOLD = TOL
452 TOL = TOLMIN
453 IF (NIO.EQ.0)
454 $ WRITE (6,*) 'New MG-tolerance (OTR) = ',TOL,TOLOLD
455 endif
456 endif
457C
458 return
459 end
460C
461C
462 subroutine setdtc
463C--------------------------------------------------------------
464C
465C Compute new timestep based on CFL-condition
466C
467C--------------------------------------------------------------
468 include 'SIZE'
469 include 'GEOM'
470 include 'MVGEOM'
471 include 'MASS'
472 include 'INPUT'
473 include 'SOLN'
474 include 'TSTEP'
475C
476 common /ctmp1/ u(lx1,ly1,lz1,lelv)
477 $ , v(lx1,ly1,lz1,lelv)
478 $ , w(lx1,ly1,lz1,lelv)
479 common /ctmp0/ x(lx1,ly1,lz1,lelv)
480 $ , r(lx1,ly1,lz1,lelv)
481 common /udxmax/ umax
482
483 common /scruz/ cx(lx1*ly1*lz1*lelv)
484 $ , cy(lx1,ly1,lz1,lelv)
485 $ , cz(lx1,ly1,lz1,lelv)
486C
487C
488 REAL VCOUR
489 SAVE VCOUR
490C
491 INTEGER IFIRST
492 SAVE IFIRST
493 DATA IFIRST/0/
494C
495C
496C Steady state => all done
497C
498C
499 IF (.NOT.IFTRAN) THEN
500 IFIRST=1
501 LASTEP=1
502 return
503 endif
504C
505 irst = param(46)
506 if (irst.gt.0) ifirst=1
507C
508C First time around
509C
510 IF (IFIRST.EQ.0) THEN
511 DT = DTINIT
512 IF (IFFLOW) THEN
513 IFIELD = 1
514 CALL BCDIRVC (VX,VY,VZ,v1mask,v2mask,v3mask)
515 endif
516 endif
517 IFIRST=IFIRST+1
518C
519 DTOLD = DT
520C
521C Convection ?
522C
523C Don't enforce Courant condition if there is no convection.
524C
525C
526 ICONV=0
527 IF (IFFLOW .AND. IFNAV) ICONV=1
528 IF (IFWCNO) ICONV=1
529 IF (IFHEAT) THEN
530 DO 10 IPSCAL=0,NPSCAL
531 IF (IFADVC(IPSCAL+2)) ICONV=1
532 10 CONTINUE
533 endif
534
535 IF (ICONV.EQ.0) THEN
536 DT=0.
537 return
538 endif
539C
540C
541C Find Courant and Umax
542C
543C
544 NTOT = lx1*ly1*lz1*NELV
545 NTOTL = LX1*LY1*LZ1*LELV
546 NTOTD = NTOTL*ldim
547 COLD = COURNO
548 CMAX = 1.2*CTARG
549 CMIN = 0.8*CTARG
550
551 if (ifmvbd) then
552 call opsub3 (cx,cy,cz,vx,vy,vz,wx,wy,wz)
553 call compute_cfl(umax,cx,cy,cz,1.0)
554 else
555 call compute_cfl(umax,vx,vy,vz,1.0)
556 endif
557
558c if (nio.eq.0) write(6,1) istep,time,umax,cmax
559c 1 format(i9,1p3e12.4,' cmax')
560
561C Zero DT
562
563 IF (DT .EQ. 0.0) THEN
564
565 IF (UMAX .NE. 0.0) THEN
566 DT = CTARG/UMAX
567 VCOUR = UMAX
568 ELSEIF (IFFLOW) THEN
569C
570C We'll use the body force to predict max velocity
571C
572 CALL SETPROP
573 IFIELD = 1
574C
575 CALL MAKEUF
576 CALL OPDSSUM (BFX,BFY,BFZ)
577 CALL OPCOLV (BFX,BFY,BFZ,BINVM1)
578 FMAX=0.0
579 CALL RZERO (U,NTOTD)
580 DO 600 I=1,NTOT
581 U(I,1,1,1) = ABS(BFX(I,1,1,1))
582 V(I,1,1,1) = ABS(BFY(I,1,1,1))
583 W(I,1,1,1) = ABS(BFZ(I,1,1,1))
584 600 CONTINUE
585 FMAX = GLMAX (U,NTOTD)
586 DENSITY = AVTRAN(1)
587 AMAX = FMAX/DENSITY
588 DXCHAR = SQRT( (XM1(1,1,1,1)-XM1(2,1,1,1))**2 +
589 $ (YM1(1,1,1,1)-YM1(2,1,1,1))**2 +
590 $ (ZM1(1,1,1,1)-ZM1(2,1,1,1))**2 )
591 DXCHAR = GLMIN (dxchar,1)
592 IF (AMAX.NE.0.) THEN
593 DT = SQRT(CTARG*DXCHAR/AMAX)
594 ELSE
595 IF (NID.EQ.0)
596 $ WRITE (6,*) 'CFL: Zero velocity and body force'
597 DT = 0.0
598 return
599 endif
600 ELSEIF (IFWCNO) THEN
601 IF (NID.EQ.0)
602 $ WRITE (6,*) ' Stefan problem with no fluid flow'
603 DT = 0.0
604 return
605 endif
606C
607 ELSEIF ((DT.GT.0.0).AND.(UMAX.NE.0.0)) THEN
608C
609C
610C Nonzero DT & nonzero velocity
611C
612C
613 COURNO = DT*UMAX
614 VOLD = VCOUR
615 VCOUR = UMAX
616 IF (IFIRST.EQ.1) THEN
617 COLD = COURNO
618 VOLD = VCOUR
619 endif
620 CPRED = 2.*COURNO-COLD
621C
622C Change DT if it is too big or if it is too small
623C
624c if (nid.eq.0)
625c $write(6,917) dt,umax,vold,vcour,cpred,cmax,courno,cmin
626c 917 format(' dt',4f9.5,4f10.6)
627 IF(COURNO.GT.CMAX .OR. CPRED.GT.CMAX .OR. COURNO.LT.CMIN) THEN
628C
629 A=(VCOUR-VOLD)/DT
630 B=VCOUR
631C -C IS Target Courant number
632 C=-CTARG
633 DISCR=B**2-4*A*C
634 DTOLD=DT
635 IF(DISCR.LE.0.0)THEN
636c if (nio.eq.0)
637c $ PRINT*,'Problem calculating new DT Discriminant=',discr
638 DT=DT*(CTARG/COURNO)
639C IF(DT.GT.DTOLD) DT=DTOLD
640 ELSE IF(ABS((VCOUR-VOLD)/VCOUR).LT.0.001)THEN
641C Easy: same v as before (LINEARIZED)
642 DT=DT*(CTARG/COURNO)
643c if (nid.eq.0)
644c $write(6,918) dt,dthi,dtlow,discr,a,b,c
645c 918 format(' d2',4f9.5,4f10.6)
646 ELSE
647 DTLOW=(-B+SQRT(DISCR) )/(2.0*A)
648 DTHI =(-B-SQRT(DISCR) )/(2.0*A)
649 IF(DTHI .GT. 0.0 .AND. DTLOW .GT. 0.0)THEN
650 DT = MIN (DTHI,DTLOW)
651c if (nid.eq.0)
652c $write(6,919) dt,dthi,dtlow,discr,a,b,c
653c 919 format(' d3',4f9.5,4f10.6)
654 ELSE IF(DTHI .LE. 0.0 .AND. DTLOW .LE. 0.0)THEN
655c PRINT*,'DTLOW,DTHI',DTLOW,DTHI
656c PRINT*,'WARNING: Abnormal DT from CFL-condition'
657c PRINT*,' Keep going'
658 DT=DT*(CTARG/COURNO)
659 ELSE
660C Normal case; 1 positive root, one negative root
661 DT = MAX (DTHI,DTLOW)
662c if (nid.eq.0)
663c $write(6,929) dt,dthi,dtlow,discr,a,b,c
664c 929 format(' d4',4f9.5,4f10.6)
665 endif
666 endif
667C We'll increase gradually-- make it the geometric mean between
668c if (nid.eq.0)
669c $write(6,939) dt,dtold
670c 939 format(' d5',4f9.5,4f10.6)
671 IF (DTOLD/DT .LT. 0.2) DT = DTOLD*5
672 endif
673C
674 endif
675C
676 return
677 end
678C
679 subroutine setdtfs (dtfs)
680C
681 include 'SIZE'
682 include 'INPUT'
683 include 'SOLN'
684 include 'GEOM'
685 include 'TSTEP'
686 common /ctmp0/ stc(lx1,ly1,lz1),sigst(lx1,ly1),dtst(lx1,ly1)
687 character cb*3,cb2*2
688C
689C Applicable?
690C
691 IF (.NOT.IFSURT) THEN
692 DTFS = 0.0
693 return
694 endif
695C
696 NFACE = 2*ldim
697 NXZ1 = lx1*lz1
698 NTOT1 = lx1*ly1*lz1*NELV
699 DTFS = 1.1E+10
700C
701C Kludge : for debugging purpose Fudge Factor comes in
702C as PARAM(45)
703C
704 FACTOR = PARAM(45)
705 IF (FACTOR .LE. 0.0) FACTOR=1.0
706 FACTOR = FACTOR / SQRT( PI**3 )
707C
708 IF (ISTEP.EQ.1) CALL SETPROP
709C
710 IFIELD = 1
711 RHOMIN = GLMIN( VTRANS(1,1,1,1,IFIELD),NTOT1 )
712C
713 DO 100 IEL=1,NELV
714 DO 100 IFC=1,NFACE
715 CB =CBC(IFC,IEL,IFIELD)
716 CB2=CBC(IFC,IEL,IFIELD)
717 IF (CB2.NE.'MS' .AND. CB2.NE.'ms') GOTO 100
718 IF (CB2(1:1).EQ.'M') THEN
719 BC4 = BC(4,IFC,IEL,IFIELD)
720 CALL CFILL (SIGST,BC4,NXZ1)
721 ELSE
722 CALL FACEIS (CB,STC,IEL,IFC,lx1,ly1,lz1)
723 CALL FACEXS (SIGST,STC,IFC,0)
724 endif
725 SIGMAX = VLMAX (SIGST,NXZ1)
726 IF (SIGMAX.LE.0.0) GOTO 100
727 RHOSIG = SQRT( RHOMIN/SIGMAX )
728 IF (ldim.EQ.2) THEN
729 CALL CDXMIN2 (DTST,RHOSIG,IEL,IFC,IFAXIS)
730 ELSE
731 CALL CDXMIN3 (DTST,RHOSIG,IEL,IFC)
732 endif
733 DTMIN = VLMIN( DTST,NXZ1 )
734 DTFS = MIN( DTFS,DTMIN )
735 100 CONTINUE
736 DTFS = GLMIN( dtfs,1 )
737C
738 IF (DTFS .GT. 1.E+10) DTFS = 0.0
739 DTFS = DTFS * FACTOR
740C
741 IF (DTFS.EQ.0.0) THEN
742 IF (ISTEP.EQ.1.AND.NIO.EQ.0) THEN
743 WRITE (6,*) ' Warning - zero surface-tension may results in'
744 WRITE (6,*) ' instability of free-surface update'
745 endif
746 endif
747C
748 return
749 end
750 subroutine cdxmin2 (dtst,rhosig,iel,ifc,ifaxis)
751C
752 include 'SIZE'
753 include 'GEOM'
754 include 'DXYZ'
755 common /delrst/ drst(lx1),drsti(lx1)
756 common /ctmp0/ xfm1(lx1),yfm1(lx1),t1xf(lx1),t1yf(lx1)
757 DIMENSION DTST(LX1,1)
758 LOGICAL IFAXIS
759C
760 DELTA = 1.E-9
761 X = 1.+DELTA
762 Y = 1.
763 DIFF = ABS(X-Y)
764 IF (DIFF.EQ.0.) EPS = 1.E-6
765 IF (DIFF.GT.0.) EPS = 1.E-13
766C
767 CALL FACEC2 (XFM1,YFM1,XM1(1,1,1,IEL),YM1(1,1,1,IEL),IFC)
768C
769 IF (IFC.EQ.1 .OR. IFC.EQ.3) THEN
770 CALL MXM (DXM1,lx1,XFM1,lx1,T1XF,1)
771 CALL MXM (DXM1,lx1,YFM1,lx1,T1YF,1)
772 ELSE
773 IF (IFAXIS) CALL SETAXDY ( IFRZER(IEL) )
774 CALL MXM (DYM1,ly1,XFM1,ly1,T1XF,1)
775 CALL MXM (DYM1,ly1,YFM1,ly1,T1YF,1)
776 endif
777C
778 IF (IFAXIS) THEN
779 DO 100 IX=1,lx1
780 IF (YFM1(IX) .LT. EPS) THEN
781 DTST(IX,1) = 1.e+10
782 ELSE
783 XJ = SQRT( T1XF(IX)**2 + T1YF(IX)**2 )*DRST(IX)
784 DTST(IX,1) = RHOSIG * SQRT( XJ**3 ) * YFM1(IX)
785 endif
786 100 CONTINUE
787 ELSE
788 DO 200 IX=1,lx1
789 XJ = SQRT( T1XF(IX)**2 + T1YF(IX)**2 )*DRST(IX)
790 DTST(IX,1) = RHOSIG * SQRT( XJ**3 )
791 200 CONTINUE
792 endif
793C
794 return
795 end
796 subroutine cdxmin3 (dtst,rhosig,iel,ifc)
797C
798 include 'SIZE'
799 include 'GEOM'
800 include 'DXYZ'
801 common /delrst/ drst(lx1),drsti(lx1)
802 common /ctmp0/ xfm1(lx1,ly1),yfm1(lx1,ly1),zfm1(lx1,ly1)
803 common /ctmp1/ drm1(lx1,lx1),drtm1(lx1,ly1)
804 $ , dsm1(lx1,lx1),dstm1(lx1,ly1)
805 common /scrmg/ xrm1(lx1,ly1),yrm1(lx1,ly1),zrm1(lx1,ly1)
806 $ , xsm1(lx1,ly1),ysm1(lx1,ly1),zsm1(lx1,ly1)
807 dimension dtst(lx1,ly1)
808C
809 call facexv (xfm1,yfm1,zfm1,xm1(1,1,1,iel),ym1(1,1,1,iel),
810 $ zm1(1,1,1,iel),ifc,0)
811 call setdrs (drm1,drtm1,dsm1,dstm1,ifc)
812C
813 CALL MXM (DRM1,lx1, XFM1,lx1,XRM1,ly1)
814 CALL MXM (DRM1,lx1, YFM1,lx1,YRM1,ly1)
815 CALL MXM (DRM1,lx1, ZFM1,lx1,ZRM1,ly1)
816 CALL MXM (XFM1,lx1,DSTM1,ly1,XSM1,ly1)
817 CALL MXM (YFM1,lx1,DSTM1,ly1,YSM1,ly1)
818 CALL MXM (ZFM1,lx1,DSTM1,ly1,ZSM1,ly1)
819C
820 DO 100 IX=1,lx1
821 DO 100 IY=1,ly1
822 DELR = XRM1(IX,IY)**2 + YRM1(IX,IY)**2 + ZRM1(IX,IY)**2
823 DELS = XSM1(IX,IY)**2 + YSM1(IX,IY)**2 + ZSM1(IX,IY)**2
824 DELR = SQRT( DELR )*DRST(IX)
825 DELS = SQRT( DELS )*DRST(IY)
826 XJ = MIN( DELR,DELS )
827 DTST(IX,IY) = RHOSIG * SQRT( XJ**3 )
828 100 CONTINUE
829C
830 return
831 end
832C
833 FUNCTION FACDOT(A,B,IFACE1)
834C
835C Take the dot product of A and B on the surface IFACE1 of element IE.
836C
837C IFACE1 is in the preprocessor notation
838C IFACE is the dssum notation.
839C 5 Jan 1989 15:12:22 PFF
840C
841 include 'SIZE'
842 include 'TOPOL'
843 DIMENSION A(LX1,LY1,LZ1),B(LX1,LY1)
844C
845C Set up counters
846C
847 CALL DSSET(lx1,ly1,lz1)
848 IFACE = EFACE1(IFACE1)
849 JS1 = SKPDAT(1,IFACE)
850 JF1 = SKPDAT(2,IFACE)
851 JSKIP1 = SKPDAT(3,IFACE)
852 JS2 = SKPDAT(4,IFACE)
853 JF2 = SKPDAT(5,IFACE)
854 JSKIP2 = SKPDAT(6,IFACE)
855C
856 SUM=0.0
857 I = 0
858 DO 100 J2=JS2,JF2,JSKIP2
859 DO 100 J1=JS1,JF1,JSKIP1
860 I = I+1
861 SUM = SUM + A(J1,J2,1)*B(I,1)
862 100 CONTINUE
863C
864 FACDOT = SUM
865C
866 return
867 end
868C
869 subroutine fcaver(xaver,a,iel,iface1)
870C------------------------------------------------------------------------
871C
872C Compute the average of A over the face IFACE1 in element IEL.
873C
874C A is a (NX,NY,NZ) data structure
875C IFACE1 is in the preprocessor notation
876C IFACE is the dssum notation.
877C------------------------------------------------------------------------
878 include 'SIZE'
879 include 'GEOM'
880 include 'TOPOL'
881 REAL A(LX1,LY1,LZ1,1)
882C
883 FCAREA = 0.
884 XAVER = 0.
885C
886C Set up counters
887C
888 CALL DSSET(lx1,ly1,lz1)
889 IFACE = EFACE1(IFACE1)
890 JS1 = SKPDAT(1,IFACE)
891 JF1 = SKPDAT(2,IFACE)
892 JSKIP1 = SKPDAT(3,IFACE)
893 JS2 = SKPDAT(4,IFACE)
894 JF2 = SKPDAT(5,IFACE)
895 JSKIP2 = SKPDAT(6,IFACE)
896C
897 I = 0
898 DO 100 J2=JS2,JF2,JSKIP2
899 DO 100 J1=JS1,JF1,JSKIP1
900 I = I+1
901 FCAREA = FCAREA+AREA(I,1,IFACE1,IEL)
902 XAVER = XAVER +AREA(I,1,IFACE1,IEL)*A(J1,J2,1,IEL)
903 100 CONTINUE
904C
905 XAVER = XAVER/FCAREA
906 return
907 end
908 subroutine faccl2(a,b,iface1)
909C
910C Collocate B with A on the surface IFACE1 of element IE.
911C
912C A is a (NX,NY,NZ) data structure
913C B is a (NX,NY,IFACE) data structure
914C IFACE1 is in the preprocessor notation
915C IFACE is the dssum notation.
916C 5 Jan 1989 15:12:22 PFF
917C
918 include 'SIZE'
919 include 'TOPOL'
920 DIMENSION A(LX1,LY1,LZ1),B(LX1,LY1)
921C
922C Set up counters
923C
924 CALL DSSET(lx1,ly1,lz1)
925 IFACE = EFACE1(IFACE1)
926 JS1 = SKPDAT(1,IFACE)
927 JF1 = SKPDAT(2,IFACE)
928 JSKIP1 = SKPDAT(3,IFACE)
929 JS2 = SKPDAT(4,IFACE)
930 JF2 = SKPDAT(5,IFACE)
931 JSKIP2 = SKPDAT(6,IFACE)
932C
933 I = 0
934 DO 100 J2=JS2,JF2,JSKIP2
935 DO 100 J1=JS1,JF1,JSKIP1
936 I = I+1
937 A(J1,J2,1) = A(J1,J2,1)*B(I,1)
938 100 CONTINUE
939C
940 return
941 end
942C
943 subroutine faccl3(a,b,c,iface1)
944C
945C Collocate B with A on the surface IFACE1 of element IE.
946C
947C A is a (NX,NY,NZ) data structure
948C B is a (NX,NY,IFACE) data structure
949C IFACE1 is in the preprocessor notation
950C IFACE is the dssum notation.
951C 5 Jan 1989 15:12:22 PFF
952C
953 include 'SIZE'
954 include 'TOPOL'
955 DIMENSION A(LX1,LY1,LZ1),B(LX1,LY1,LZ1),C(LX1,LY1)
956C
957C Set up counters
958C
959 CALL DSSET(lx1,ly1,lz1)
960 IFACE = EFACE1(IFACE1)
961 JS1 = SKPDAT(1,IFACE)
962 JF1 = SKPDAT(2,IFACE)
963 JSKIP1 = SKPDAT(3,IFACE)
964 JS2 = SKPDAT(4,IFACE)
965 JF2 = SKPDAT(5,IFACE)
966 JSKIP2 = SKPDAT(6,IFACE)
967C
968 I = 0
969 DO 100 J2=JS2,JF2,JSKIP2
970 DO 100 J1=JS1,JF1,JSKIP1
971 I = I+1
972 A(J1,J2,1) = B(J1,J2,1)*C(I,1)
973 100 CONTINUE
974C
975 return
976 end
977 subroutine faddcl3(a,b,c,iface1)
978C
979C Collocate B with C and add to A on the surface IFACE1 of element IE.
980C
981C A is a (NX,NY,NZ) data structure
982C B is a (NX,NY,NZ) data structure
983C C is a (NX,NY,IFACE) data structure
984C IFACE1 is in the preprocessor notation
985C IFACE is the dssum notation.
986C 29 Jan 1990 18:00 PST PFF
987C
988 include 'SIZE'
989 include 'TOPOL'
990 DIMENSION A(LX1,LY1,LZ1),B(LX1,LY1,LZ1),C(LX1,LY1)
991C
992C Set up counters
993C
994 CALL DSSET(lx1,ly1,lz1)
995 IFACE = EFACE1(IFACE1)
996 JS1 = SKPDAT(1,IFACE)
997 JF1 = SKPDAT(2,IFACE)
998 JSKIP1 = SKPDAT(3,IFACE)
999 JS2 = SKPDAT(4,IFACE)
1000 JF2 = SKPDAT(5,IFACE)
1001 JSKIP2 = SKPDAT(6,IFACE)
1002C
1003 I = 0
1004 DO 100 J2=JS2,JF2,JSKIP2
1005 DO 100 J1=JS1,JF1,JSKIP1
1006 I = I+1
1007 A(J1,J2,1) = A(J1,J2,1) + B(J1,J2,1)*C(I,1)
1008 100 CONTINUE
1009C
1010 return
1011 end
1012c-----------------------------------------------------------------------
1013 subroutine sethlm (h1,h2,intloc)
1014
1015c Set the variable property arrays H1 and H2
1016c in the Helmholtz equation.
1017c (associated with variable IFIELD)
1018c INTLOC = integration type
1019
1020 include 'SIZE'
1021 include 'INPUT'
1022 include 'SOLN'
1023 include 'TSTEP'
1024
1025 real h1(1),h2(1)
1026
1027 nel = nelfld(ifield)
1028 ntot1 = lx1*ly1*lz1*nel
1029
1030 if (iftran) then
1031 dtbd = bd(1)/dt
1032 call copy (h1,vdiff (1,1,1,1,ifield),ntot1)
1033 if (intloc.eq.0) then
1034 call rzero (h2,ntot1)
1035 else
1036 call cmult2 (h2,vtrans(1,1,1,1,ifield),dtbd,ntot1)
1037 endif
1038
1039c if (ifield.eq.1 .and. ifanls) then ! this should be replaced
1040c const = 2. ! with a correct stress
1041c call cmult (h1,const,ntot1) ! formulation
1042c endif
1043
1044 ELSE
1045 CALL COPY (H1,VDIFF (1,1,1,1,IFIELD),NTOT1)
1046 CALL RZERO (H2,NTOT1)
1047 endif
1048
1049 return
1050 end
1051c-----------------------------------------------------------------------
1052 subroutine nekuvp (iel)
1053C------------------------------------------------------------------
1054C
1055C Generate user-specified material properties
1056C
1057C------------------------------------------------------------------
1058 include 'SIZE'
1059 include 'INPUT'
1060 include 'SOLN'
1061 include 'TSTEP'
1062 include 'PARALLEL'
1063 include 'NEKUSE'
1064 ielg = lglel(iel)
1065c IF (IFSTRS .AND. IFIELD.EQ.1) CALL STNRINV ! don't call! pff, 2007
1066 DO 10 K=1,lz1
1067 DO 10 J=1,ly1
1068 DO 10 I=1,lx1
1069 if (optlevel.le.2) CALL NEKASGN (I,J,K,IEL)
1070 CALL USERVP (I,J,K,IELG)
1071 VDIFF (I,J,K,IEL,IFIELD) = UDIFF
1072 VTRANS(I,J,K,IEL,IFIELD) = UTRANS
1073 10 CONTINUE
1074 return
1075 end
1076C
1077 subroutine diagnos
1078 return
1079 end
1080
1081c-----------------------------------------------------------------------
1082 subroutine setsolv
1083 include 'SIZE'
1084
1085 common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
1086 logical ifdfrm, iffast, ifh2, ifsolv
1087
1088 ifsolv = .false.
1089
1090 return
1091 end
1092c-----------------------------------------------------------------------
1093 subroutine hmhzsf (name,u1,u2,u3,r1,r2,r3,h1,h2,
1094 $ rmask1,rmask2,rmask3,rmult,
1095 $ tol,maxit,matmod)
1096
1097c Solve coupled Helmholtz equations (stress formulation)
1098
1099
1100 include 'SIZE'
1101 include 'INPUT'
1102 include 'MASS'
1103 include 'SOLN' ! For outpost diagnostic call
1104 include 'TSTEP'
1105 include 'ORTHOSTRS'
1106 include 'CTIMER'
1107
1108 real u1(1),u2(1),u3(1),r1(1),r2(1),r3(1),h1(1),h2(1)
1109 real rmask1(1),rmask2(1),rmask3(1),rmult(1)
1110 character name*4
1111
1112#ifdef TIMER
1113 nhmhz = nhmhz + 1
1114 etime1 = dnekclock()
1115#endif
1116
1117 nel = nelfld(ifield)
1118 vol = volfld(ifield)
1119 n = lx1*ly1*lz1*nel
1120
1121 napproxstrs(1) = 0
1122 iproj = 0
1123 if (ifprojfld(ifield)) iproj = param(94)
1124 if (iproj.gt.0.and.istep.ge.iproj) napproxstrs(1)=param(93)
1125 napproxstrs(1)=min(napproxstrs(1),mxprev)
1126
1127 call rmask (r1,r2,r3,nel)
1128 call opdssum (r1,r2,r3)
1129 call rzero3 (u1,u2,u3,n)
1130
1131 if (imesh.eq.1) then
1132 call chktcgs (r1,r2,r3,rmask1,rmask2,rmask3,rmult,binvm1
1133 $ ,vol,tol,nel)
1134
1135 call strs_project_a(r1,r2,r3,h1,h2,rmult,ifield,ierr,matmod)
1136
1137 call cggosf (u1,u2,u3,r1,r2,r3,h1,h2,rmult,binvm1
1138 $ ,vol,tol,maxit,matmod)
1139
1140 call strs_project_b(u1,u2,u3,h1,h2,rmult,ifield,ierr)
1141
1142 else
1143
1144 call chktcgs (r1,r2,r3,rmask1,rmask2,rmask3,rmult,bintm1
1145 $ ,vol,tol,nel)
1146 call cggosf (u1,u2,u3,r1,r2,r3,h1,h2,rmult,bintm1
1147 $ ,vol,tol,maxit,matmod)
1148
1149 endif
1150
1151#ifdef TIMER
1152 thmhz=thmhz+(dnekclock()-etime1)
1153#endif
1154
1155 return
1156 end
1157
1158 subroutine chktcgs (r1,r2,r3,rmask1,rmask2,rmask3,rmult,binv,
1159 $ vol,tol,nel)
1160C-------------------------------------------------------------------
1161C
1162C Check that the tolerances are not too small for the CG-solver.
1163C Important when calling the CG-solver (Gauss-Lobatto mesh) with
1164C zero Neumann b.c.
1165C
1166C-------------------------------------------------------------------
1167 include 'SIZE'
1168 include 'INPUT'
1169 include 'MASS'
1170 include 'EIGEN'
1171 common /cprint/ ifprint
1172 logical ifprint
1173 common /ctmp0/ wa (lx1,ly1,lz1,lelt)
1174C
1175 dimension r1 (lx1,ly1,lz1,1)
1176 $ , r2 (lx1,ly1,lz1,1)
1177 $ , r3 (lx1,ly1,lz1,1)
1178 $ , rmask1(lx1,ly1,lz1,1)
1179 $ , rmask2(lx1,ly1,lz1,1)
1180 $ , rmask3(lx1,ly1,lz1,1)
1181 $ , rmult (lx1,ly1,lz1,1)
1182 $ , binv (lx1,ly1,lz1,1)
1183C
1184 NTOT1 = lx1*ly1*lz1*NEL
1185C
1186 IF (EIGAA .NE. 0.0) THEN
1187 ACONDNO = EIGGA/EIGAA
1188 ELSE
1189 ACONDNO = 10.0
1190 endif
1191C
1192C Check Single or double precision
1193C
1194 DELTA = 1.0E-9
1195 X = 1.0 + DELTA
1196 Y = 1.0
1197 DIFF = ABS(X - Y)
1198 IF (DIFF .EQ. 0.0) EPS = 1.0E-6
1199 IF (DIFF .GT. 0.0) EPS = 1.0E-13
1200C
1201 CALL OPDOT (WA,R1,R2,R3,R1,R2,R3,NTOT1)
1202 RINIT = GLSC3(WA,BINV,RMULT,NTOT1)
1203 RINIT = SQRT (RINIT/VOL)
1204 RMIN = EPS*RINIT
1205C
1206 IF (TOL.LT.RMIN) THEN
1207 TOLOLD = TOL
1208 TOL = RMIN
1209c IF (NIO.EQ.0 .AND. IFPRINT)
1210c $ WRITE(6,*)'New CG1(stress)-tolerance (RINIT*epsm) = ',TOL,TOLOLD
1211 endif
1212C
1213 IF (ldim.EQ.2) THEN
1214 CALL ADD3 (WA,RMASK1,RMASK2,NTOT1)
1215 ELSE
1216 CALL ADD4 (WA,RMASK1,RMASK2,RMASK3,NTOT1)
1217 endif
1218 BCNEU1 = GLSC2 (WA,RMULT,NTOT1)
1219 BCNEU2 = (ldim) * GLSUM(RMULT,NTOT1)
1220 BCTEST = ABS(BCNEU1 - BCNEU2)
1221 IF (BCTEST .LT. 0.1) THEN
1222 IF (ldim.EQ.2) THEN
1223 CALL ADD3 (WA,R1,R2,NTOT1)
1224 ELSE
1225 CALL ADD4 (WA,R1,R2,R3,NTOT1)
1226 endif
1227 OTR = GLSC2(WA,RMULT,NTOT1) / ( ldim )
1228 TOLMIN = ABS(OTR) * ACONDNO
1229 IF (TOL .LT. TOLMIN) THEN
1230 TOLOLD = TOL
1231 TOL = TOLMIN
1232c IF (NIO.EQ.0)
1233c $ WRITE (6,*) 'New CG1(stress)-tolerance (OTR) = ',TOL,TOLOLD
1234 endif
1235 endif
1236C
1237 return
1238 end
1239c-----------------------------------------------------------------------
1240 subroutine axhmsf (au1,au2,au3,u1,u2,u3,h1,h2,matmod)
1241
1242C Compute the coupled Helmholtz matrix-vector products
1243
1244C Fluid (MATMOD .GE. 0) : Hij Uj = Aij*Uj + H2*B*Ui
1245C Solid (MATMOD .LT. 0) : Hij Uj = Kij*Uj
1246C
1247C-----------------------------------------------------------------------
1248 include 'SIZE'
1249 include 'INPUT'
1250 include 'GEOM'
1251 include 'MASS'
1252 include 'TSTEP'
1253 include 'CTIMER'
1254
1255 common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
1256 logical ifdfrm, iffast, ifh2, ifsolv
1257C
1258 dimension au1(lx1,ly1,lz1,1)
1259 $ , au2(lx1,ly1,lz1,1)
1260 $ , au3(lx1,ly1,lz1,1)
1261 $ , u1 (lx1*ly1*lz1*1)
1262 $ , u2 (lx1*ly1*lz1*1)
1263 $ , u3 (lx1*ly1*lz1*1)
1264 $ , h1 (lx1,ly1,lz1,1)
1265 $ , h2 (lx1,ly1,lz1,1)
1266
1267 naxhm = naxhm + 1
1268 etime1 = dnekclock()
1269
1270 nel = nelfld(ifield)
1271 ntot1 = lx1*ly1*lz1*nel
1272
1273c if (ifaxis.and.ifsplit) call exitti(
1274c $'Axisymmetric stress w/PnPn not yet supported.$',istep)
1275
1276c icase = 1 --- axsf_fast (no axisymmetry)
1277c icase = 2 --- stress formulation and supports axisymmetry
1278c icase = 3 --- 3 separate axhelm calls
1279
1280 icase = 1 ! Fast mode for stress
1281 if (ifaxis) icase=2 ! Slow for stress, but supports axisymmetry
1282 if (matmod.lt.0) icase=2 ! Elasticity case
1283c if (matmod.lt.0) icase=3 ! Block-diagonal Axhelm
1284 if (matmod.lt.0) icase=1 ! Elasticity case (faster, 7/28/17,pff)
1285 if (.not.ifstrs) icase=3 ! Block-diagonal Axhelm
1286
1287 if (icase.eq.1) then
1288
1289 call axsf_fast(au1,au2,au3,u1,u2,u3,h1,h2,ifield)
1290
1291 elseif (icase.eq.3) then
1292
1293 call axhelm(au1,u1,h1,h2,1,1)
1294 call axhelm(au2,u2,h1,h2,1,2)
1295 if (if3d) call axhelm(au3,u3,h1,h2,1,3)
1296
1297 else ! calculate coupled Aij Uj products
1298
1299 if ( .not.ifsolv ) call setfast (h1,h2,imesh)
1300
1301 call stnrate (u1,u2,u3,nel,matmod)
1302 call stress (h1,h2,nel,matmod,ifaxis)
1303 call aijuj (au1,au2,au3,nel,ifaxis)
1304
1305 if (ifh2 .and. matmod.ge.0) then ! add Helmholtz contributions
1306 call addcol4 (au1,bm1,h2,u1,ntot1)
1307 call addcol4 (au2,bm1,h2,u2,ntot1)
1308 if (ldim.eq.3) call addcol4 (au3,bm1,h2,u3,ntot1)
1309 endif
1310
1311 endif
1312
1313 taxhm=taxhm+(dnekclock()-etime1)
1314
1315 return
1316 end
1317c-----------------------------------------------------------------------
1318 subroutine stnrate (u1,u2,u3,nel,matmod)
1319C
1320C Compute strainrates
1321C
1322C CAUTION : Stresses and strainrates share the same scratch commons
1323C
1324 include 'SIZE'
1325 include 'INPUT'
1326 include 'GEOM'
1327 include 'TSTEP'
1328
1329 common /ctmp0/ exz(lx1*ly1*lz1*lelt)
1330 $ , eyz(lx1*ly1*lz1*lelt)
1331 common /ctmp1/ exx(lx1*ly1*lz1*lelt)
1332 $ , exy(lx1*ly1*lz1*lelt)
1333 $ , eyy(lx1*ly1*lz1*lelt)
1334 $ , ezz(lx1*ly1*lz1*lelt)
1335c
1336 dimension u1(lx1,ly1,lz1,1)
1337 $ , u2(lx1,ly1,lz1,1)
1338 $ , u3(lx1,ly1,lz1,1)
1339C
1340 NTOT1 = lx1*ly1*lz1*NEL
1341
1342 CALL RZERO3 (EXX,EYY,EZZ,NTOT1)
1343 CALL RZERO3 (EXY,EXZ,EYZ,NTOT1)
1344
1345 CALL UXYZ (U1,EXX,EXY,EXZ,NEL)
1346 CALL UXYZ (U2,EXY,EYY,EYZ,NEL)
1347 IF (ldim.EQ.3) CALL UXYZ (U3,EXZ,EYZ,EZZ,NEL)
1348
1349 CALL INVCOL2 (EXX,JACM1,NTOT1)
1350 CALL INVCOL2 (EXY,JACM1,NTOT1)
1351 CALL INVCOL2 (EYY,JACM1,NTOT1)
1352
1353 IF (IFAXIS) CALL AXIEZZ (U2,EYY,EZZ,NEL)
1354C
1355 IF (ldim.EQ.3) THEN
1356 CALL INVCOL2 (EXZ,JACM1,NTOT1)
1357 CALL INVCOL2 (EYZ,JACM1,NTOT1)
1358 CALL INVCOL2 (EZZ,JACM1,NTOT1)
1359 endif
1360C
1361 return
1362 end
1363c-----------------------------------------------------------------------
1364 subroutine stress (h1,h2,nel,matmod,ifaxis)
1365C
1366C MATMOD.GE.0 Fluid material models
1367C MATMOD.LT.0 Solid material models
1368C
1369C CAUTION : Stresses and strainrates share the same scratch commons
1370C
1371 include 'SIZE'
1372 common /ctmp1/ txx(lx1,ly1,lz1,lelt)
1373 $ , txy(lx1,ly1,lz1,lelt)
1374 $ , tyy(lx1,ly1,lz1,lelt)
1375 $ , tzz(lx1,ly1,lz1,lelt)
1376 common /ctmp0/ txz(lx1,ly1,lz1,lelt)
1377 $ , tyz(lx1,ly1,lz1,lelt)
1378 common /scrsf/ t11(lx1,ly1,lz1,lelt)
1379 $ , t22(lx1,ly1,lz1,lelt)
1380 $ , t33(lx1,ly1,lz1,lelt)
1381 $ , hii(lx1,ly1,lz1,lelt)
1382C
1383 DIMENSION H1(LX1,LY1,LZ1,1),H2(LX1,LY1,LZ1,1)
1384 LOGICAL IFAXIS
1385
1386 NTOT1 = lx1*ly1*lz1*NEL
1387
1388 IF (MATMOD.EQ.0) THEN
1389
1390C Newtonian fluids
1391
1392 CONST = 2.0
1393 CALL CMULT2 (HII,H1,CONST,NTOT1)
1394 CALL COL2 (TXX,HII,NTOT1)
1395 CALL COL2 (TXY,H1 ,NTOT1)
1396 CALL COL2 (TYY,HII,NTOT1)
1397 IF (IFAXIS .OR. ldim.EQ.3) CALL COL2 (TZZ,HII,NTOT1)
1398 IF (ldim.EQ.3) THEN
1399 CALL COL2 (TXZ,H1 ,NTOT1)
1400 CALL COL2 (TYZ,H1 ,NTOT1)
1401 endif
1402C
1403 ELSEIF (MATMOD.EQ.-1) THEN
1404C
1405C Elastic solids
1406C
1407 CONST = 2.0
1408 CALL ADD3S (HII,H1,H2,CONST,NTOT1)
1409 CALL COPY (T11,TXX,NTOT1)
1410 CALL COPY (T22,TYY,NTOT1)
1411 CALL COL3 (TXX,HII,T11,NTOT1)
1412 CALL ADDCOL3 (TXX,H1 ,T22,NTOT1)
1413 CALL COL3 (TYY,H1 ,T11,NTOT1)
1414 CALL ADDCOL3 (TYY,HII,T22,NTOT1)
1415 CALL COL2 (TXY,H2 ,NTOT1)
1416 IF (IFAXIS .OR. ldim.EQ.3) THEN
1417 CALL COPY (T33,TZZ,NTOT1)
1418 CALL COL3 (TZZ,H1 ,T11,NTOT1)
1419 CALL ADDCOL3 (TZZ,H1 ,T22,NTOT1)
1420 CALL ADDCOL3 (TZZ,HII,T33,NTOT1)
1421 CALL ADDCOL3 (TXX,H1 ,T33,NTOT1)
1422 CALL ADDCOL3 (TYY,H1 ,T33,NTOT1)
1423 endif
1424 IF (ldim.EQ.3) THEN
1425 CALL COL2 (TXZ,H2 ,NTOT1)
1426 CALL COL2 (TYZ,H2 ,NTOT1)
1427 endif
1428C
1429 endif
1430C
1431 return
1432 end
1433c-----------------------------------------------------------------------
1434 subroutine aijuj (au1,au2,au3,nel,ifaxis)
1435C
1436 include 'SIZE'
1437 common /ctmp1/ txx(lx1,ly1,lz1,lelt)
1438 $ , txy(lx1,ly1,lz1,lelt)
1439 $ , tyy(lx1,ly1,lz1,lelt)
1440 $ , tzz(lx1,ly1,lz1,lelt)
1441 common /ctmp0/ txz(lx1,ly1,lz1,lelt)
1442 $ , tyz(lx1,ly1,lz1,lelt)
1443C
1444 DIMENSION AU1(LX1,LY1,LZ1,1)
1445 $ , AU2(LX1,LY1,LZ1,1)
1446 $ , AU3(LX1,LY1,LZ1,1)
1447 LOGICAL IFAXIS
1448C
1449 CALL TTXYZ (AU1,TXX,TXY,TXZ,NEL)
1450 CALL TTXYZ (AU2,TXY,TYY,TYZ,NEL)
1451 IF (IFAXIS) CALL AXITZZ (AU2,TZZ,NEL)
1452 IF (ldim.EQ.3) CALL TTXYZ (AU3,TXZ,TYZ,TZZ,NEL)
1453C
1454 return
1455 end
1456c-----------------------------------------------------------------------
1457 subroutine uxyz (u,ex,ey,ez,nel)
1458C
1459 include 'SIZE'
1460 include 'GEOM'
1461 common /scrsf/ ur(lx1,ly1,lz1,lelt)
1462 $ , us(lx1,ly1,lz1,lelt)
1463 $ , ut(lx1,ly1,lz1,lelt)
1464c
1465 dimension u (lx1,ly1,lz1,1)
1466 $ , ex(lx1,ly1,lz1,1)
1467 $ , ey(lx1,ly1,lz1,1)
1468 $ , ez(lx1,ly1,lz1,1)
1469C
1470 NTOT1 = lx1*ly1*lz1*NEL
1471C
1472 CALL URST (U,UR,US,UT,NEL)
1473C
1474 CALL ADDCOL3 (EX,RXM1,UR,NTOT1)
1475 CALL ADDCOL3 (EX,SXM1,US,NTOT1)
1476 CALL ADDCOL3 (EY,RYM1,UR,NTOT1)
1477 CALL ADDCOL3 (EY,SYM1,US,NTOT1)
1478C
1479 IF (ldim.EQ.3) THEN
1480 CALL ADDCOL3 (EZ,RZM1,UR,NTOT1)
1481 CALL ADDCOL3 (EZ,SZM1,US,NTOT1)
1482 CALL ADDCOL3 (EZ,TZM1,UT,NTOT1)
1483 CALL ADDCOL3 (EX,TXM1,UT,NTOT1)
1484 CALL ADDCOL3 (EY,TYM1,UT,NTOT1)
1485 endif
1486C
1487 return
1488 end
1489 subroutine urst (u,ur,us,ut,nel)
1490C
1491 include 'SIZE'
1492 include 'GEOM'
1493 include 'INPUT'
1494C
1495 DIMENSION U (LX1,LY1,LZ1,1)
1496 $ , UR(LX1,LY1,LZ1,1)
1497 $ , US(LX1,LY1,LZ1,1)
1498 $ , UT(LX1,LY1,LZ1,1)
1499C
1500 DO 100 IEL=1,NEL
1501 IF (IFAXIS) CALL SETAXDY ( IFRZER(IEL) )
1502 CALL DDRST (U (1,1,1,IEL),UR(1,1,1,IEL),
1503 $ US(1,1,1,IEL),UT(1,1,1,IEL))
1504 100 CONTINUE
1505C
1506 return
1507 end
1508 subroutine ddrst (u,ur,us,ut)
1509C
1510 include 'SIZE'
1511 include 'DXYZ'
1512C
1513 DIMENSION U (LX1,LY1,LZ1)
1514 $ , UR(LX1,LY1,LZ1)
1515 $ , US(LX1,LY1,LZ1)
1516 $ , UT(LX1,LY1,LZ1)
1517C
1518 NXY1 = lx1*ly1
1519 NYZ1 = ly1*lz1
1520C
1521 CALL MXM (DXM1,lx1,U,lx1,UR,NYZ1)
1522 IF (ldim.EQ.2) THEN
1523 CALL MXM (U,lx1,DYTM1,ly1,US,ly1)
1524 ELSE
1525 DO 10 IZ=1,lz1
1526 CALL MXM (U(1,1,IZ),lx1,DYTM1,ly1,US(1,1,IZ),ly1)
1527 10 CONTINUE
1528 CALL MXM (U,NXY1,DZTM1,lz1,UT,lz1)
1529 endif
1530C
1531 return
1532 end
1533 subroutine axiezz (u2,eyy,ezz,nel)
1534C
1535 include 'SIZE'
1536 include 'GEOM'
1537C
1538 DIMENSION U2 (LX1,LY1,LZ1,1)
1539 $ , EYY(LX1,LY1,LZ1,1)
1540 $ , EZZ(LX1,LY1,LZ1,1)
1541C
1542 NXYZ1 = lx1*ly1*lz1
1543C
1544 DO 100 IEL=1,NEL
1545 IF ( IFRZER(IEL) ) THEN
1546 DO 200 IX=1,lx1
1547 EZZ(IX, 1,1,IEL) = EYY(IX,1,1,IEL)
1548 DO 200 IY=2,ly1
1549 EZZ(IX,IY,1,IEL) = U2(IX,IY,1,IEL) / YM1(IX,IY,1,IEL)
1550 200 CONTINUE
1551 ELSE
1552 CALL INVCOL3 (EZZ(1,1,1,IEL),U2(1,1,1,IEL),YM1(1,1,1,IEL),
1553 $ NXYZ1)
1554 endif
1555 100 CONTINUE
1556C
1557 return
1558 end
1559c-----------------------------------------------------------------------
1560 subroutine flush_io
1561 return
1562 end
1563c-----------------------------------------------------------------------
1564 subroutine fcsum2(xsum,asum,x,e,f)
1565c
1566c Compute the weighted sum of X over face f of element e
1567c
1568c x is an (NX,NY,NZ) data structure
1569c f is in the preprocessor notation
1570c
1571c xsum is sum (X*area)
1572c asum is sum (area)
1573
1574
1575 include 'SIZE'
1576 include 'GEOM'
1577 include 'TOPOL'
1578 real x(lx1,ly1,lz1,1)
1579 integer e,f,fd
1580
1581 asum = 0.
1582 xsum = 0.
1583
1584c Set up counters ; fd is the dssum notation.
1585 call dsset(lx1,ly1,lz1)
1586 fd = eface1(f)
1587 js1 = skpdat(1,fd)
1588 jf1 = skpdat(2,fd)
1589 jskip1 = skpdat(3,fd)
1590 js2 = skpdat(4,fd)
1591 jf2 = skpdat(5,fd)
1592 jskip2 = skpdat(6,fd)
1593
1594 i = 0
1595 do j2=js2,jf2,jskip2
1596 do j1=js1,jf1,jskip1
1597 i = i+1
1598 xsum = xsum+area(i,1,f,e)*x(j1,j2,1,e)
1599 asum = asum+area(i,1,f,e)
1600 enddo
1601 enddo
1602
1603 return
1604 end
1605c-----------------------------------------------------------------------
1606 function surf_mean(u,ifld,bc_in,ierr)
1607
1608 include 'SIZE'
1609 include 'TOTAL'
1610
1611 real u(1)
1612
1613 integer e,f
1614 character*3 bc_in
1615
1616 usum = 0
1617 asum = 0
1618
1619 nface = 2*ldim
1620 do e=1,nelfld(ifld)
1621 do f=1,nface
1622 if (cbc(f,e,ifld).eq.bc_in) then
1623 call fcsum2(usum_f,asum_f,u,e,f)
1624 usum = usum + usum_f
1625 asum = asum + asum_f
1626 endif
1627 enddo
1628 enddo
1629
1630 usum = glsum(usum,1) ! sum across processors
1631 asum = glsum(asum,1)
1632
1633 surf_mean = usum
1634 ierr = 1
1635
1636 if (asum.gt.0) surf_mean = usum/asum
1637 if (asum.gt.0) ierr = 0
1638
1639 return
1640 end
1641c-----------------------------------------------------------------------
1642 subroutine fdm_h1a(z,r,d,nel,kt,rr)
1643 include 'SIZE'
1644 include 'TOTAL'
1645c
1646 common /ctmp0/ w(lx1,ly1,lz1)
1647c
1648 include 'FDMH1'
1649c
1650c Overlapping Schwarz, FDM based
1651c
1652 real z(lx1,ly1,lz1,1)
1653 real r(lx1,ly1,lz1,1)
1654 real d(lx1,ly1,lz1,1)
1655 real rr(lx1,ly1,lz1,1)
1656
1657 integer kt(lelt,3)
1658
1659 integer icalld
1660 save icalld
1661 data icalld /0/
1662
1663 n1 = lx1
1664 n2 = lx1*lx1
1665 n3 = lx1*lx1*lx1
1666 n = lx1*ly1*lz1*nel
1667
1668 if (ifbhalf) then
1669 call col3(rr,r,bhalf,n)
1670 else
1671 call copy(rr,r,n)
1672 endif
1673 icalld = icalld+1
1674c
1675c
1676 do ie=1,nel
1677 if (if3d) then
1678c Transfer to wave space:
1679 call mxm(fdst(1,kt(ie,1)),n1,rr(1,1,1,ie),n1,w,n2)
1680 do iz=1,n1
1681 call mxm(w(1,1,iz),n1,fds (1,kt(ie,2)),n1,z(1,1,iz,ie),n1)
1682 enddo
1683 call mxm(z(1,1,1,ie),n2,fds (1,kt(ie,3)),n1,w,n1)
1684c
1685c fdsolve:
1686c
1687 call col2(w,d(1,1,1,ie),n3)
1688c
1689c Transfer to physical space:
1690c
1691 call mxm(w,n2,fdst(1,kt(ie,3)),n1,z(1,1,1,ie),n1)
1692 do iz=1,n1
1693 call mxm(z(1,1,iz,ie),n1,fdst(1,kt(ie,2)),n1,w(1,1,iz),n1)
1694 enddo
1695 call mxm(fds (1,kt(ie,1)),n1,w,n1,z(1,1,1,ie),n2)
1696c
1697 else
1698c Transfer to wave space:
1699 call mxm(fdst(1,kt(ie,1)),n1,rr(1,1,1,ie),n1,w,n1)
1700 call mxm(w,n1,fds (1,kt(ie,2)),n1,z(1,1,1,ie),n1)
1701c
1702c fdsolve:
1703c
1704 call col2(z(1,1,1,ie),d(1,1,1,ie),n2)
1705c
1706c Transfer to physical space:
1707c
1708 call mxm(z(1,1,1,ie),n1,fdst(1,kt(ie,2)),n1,w,n1)
1709 call mxm(fds (1,kt(ie,1)),n1,w,n1,z(1,1,1,ie),n1)
1710c
1711 endif
1712 enddo
1713
1714 if (ifbhalf) call col2(z,bhalf,n)
1715
1716 call dssum(z,lx1,ly1,lz1)
1717
1718 return
1719 end
1720c-----------------------------------------------------------------------
1721 subroutine set_vert_strs(glo_num,ngv,nx,nel,vertex,ifcenter)
1722
1723c Given global array, vertex, pointing to hex vertices, set up
1724c a new array of global pointers for an nx^ldim set of elements.
1725
1726 include 'SIZE'
1727 include 'INPUT'
1728
1729 integer*8 glo_num(1),ngv
1730 integer vertex(1),nx
1731 logical ifcenter
1732
1733 common /ctmp0/ gnum(lx1*ly1*lz1*lelt)
1734 integer*8 gnum
1735
1736 call set_vert(gnum,ngv,nx,nel,vertex,ifcenter)
1737
1738 n=nel*nx*nx
1739 if (if3d) n=n*nx
1740
1741c Pack vertex ids component-wise (u,v,w_1; u,v,w_2; etc.)
1742
1743 k=0
1744 do i=1,n
1745
1746 do j=1,ldim
1747 glo_num(k+j) = ldim*(gnum(i)-1) + j
1748 enddo
1749 k=k+ldim
1750
1751 enddo
1752
1753 return
1754 end
1755c-----------------------------------------------------------------------
1756 subroutine get_strs_mask(mask,nxc,nzc,nel)
1757 include 'SIZE'
1758 include 'INPUT'
1759 include 'SOLN'
1760
1761 real mask(ldim,nxc,nxc,nzc,nel)
1762 integer e
1763
1764
1765 if (if3d) then
1766 do e=1,nel
1767 do kc=1,nxc
1768 do jc=1,nxc
1769 do ic=1,nxc
1770 k = 1 + ((lx1-1)*(kc-1))/(nxc-1)
1771 j = 1 + ((lx1-1)*(jc-1))/(nxc-1)
1772 i = 1 + ((lx1-1)*(ic-1))/(nxc-1)
1773 mask(1,ic,jc,kc,e) = v1mask(i,j,k,e)
1774 mask(2,ic,jc,kc,e) = v2mask(i,j,k,e)
1775 mask(3,ic,jc,kc,e) = v3mask(i,j,k,e)
1776 enddo
1777 enddo
1778 enddo
1779 enddo
1780 else
1781 do e=1,nel
1782 do jc=1,nxc
1783 do ic=1,nxc
1784 j = 1 + ((lx1-1)*(jc-1))/(nxc-1)
1785 i = 1 + ((lx1-1)*(ic-1))/(nxc-1)
1786 mask(1,ic,jc,1,e) = v1mask(i,j,1,e)
1787 mask(2,ic,jc,1,e) = v2mask(i,j,1,e)
1788 enddo
1789 enddo
1790 enddo
1791 endif
1792
1793 return
1794 end
1795c-----------------------------------------------------------------------
1796 subroutine axstrs(a1,a2,a3,p1,p2,p3,h1,h2,matmod,nel)
1797 real a1(1),a2(1),a3(1),p1(1),p2(1),p3(1)
1798
1799 call axhmsf (a1,a2,a3,p1,p2,p3,h1,h2,matmod)
1800 call rmask (a1,a2,a3,nel)
1801 call opdssum (a1,a2,a3)
1802
1803 return
1804 end
1805c-----------------------------------------------------------------------
1806 subroutine axstrs_nds(a1,a2,a3,p1,p2,p3,h1,h2,matmod,nel)
1807 real a1(1),a2(1),a3(1),p1(1),p2(1),p3(1)
1808
1809 call rmask (p1,p2,p3,nel)
1810 call axhmsf (a1,a2,a3,p1,p2,p3,h1,h2,matmod)
1811 call rmask (a1,a2,a3,nel)
1812c call opdssum (a1,a2,a3)
1813
1814 return
1815 end
1816c-----------------------------------------------------------------------
1817 subroutine get_local_crs_galerkin_strs(a,ncl,nxc,h1,h2,matmod)
1818
1819c This routine generates Nelv submatrices of order ncl using
1820c Galerkin projection
1821
1822 include 'SIZE'
1823 include 'INPUT'
1824 include 'TSTEP'
1825
1826
1827 real a(ldim,ncl,ldim,ncl,1),h1(1),h2(1)
1828c real a(ncl,ldim,ncl,ldim,1),h1(1),h2(1)
1829
1830 common /scrcr2/ a1(lx1*ly1*lz1,lelt),w1(lx1*ly1*lz1,lelt)
1831 $ , a2(lx1*ly1*lz1,lelt),w2(lx1*ly1*lz1,lelt)
1832 common /scrcr3/ a3(lx1*ly1*lz1,lelt),w3(lx1*ly1*lz1,lelt)
1833 $ , b (lx1*ly1*lz1,8)
1834
1835 integer e
1836
1837 nel = nelfld(ifield)
1838
1839 do j=1,ncl ! bi- or tri-linear interpolant, ONLY
1840 call gen_crs_basis(b(1,j),j)
1841 enddo
1842
1843 nxyz = lx1*ly1*lz1
1844 do k = 1,ldim
1845 call opzero(w1,w2,w3)
1846 call opzero(a1,a2,a3)
1847
1848 do j = 1,ncl
1849
1850 do e = 1,nel
1851 if (k.eq.1) call copy(w1(1,e),b(1,j),nxyz)
1852 if (k.eq.2) call copy(w2(1,e),b(1,j),nxyz)
1853 if (k.eq.3) call copy(w3(1,e),b(1,j),nxyz)
1854 enddo
1855
1856 call axstrs_nds(a1,a2,a3,w1,w2,w3,h1,h2,matmod,nel)
1857
1858 do e = 1,nel
1859 do i = 1,ncl
1860
1861 a(1,i,k,j,e)=vlsc2(b(1,i),a1(1,e),nxyz) ! bi^T * A^e * bj
1862 a(2,i,k,j,e)=vlsc2(b(1,i),a2(1,e),nxyz) ! bi^T * A^e * bj
1863 if (if3d)
1864 $ a(3,i,k,j,e)=vlsc2(b(1,i),a3(1,e),nxyz) ! bi^T * A^e * bj
1865
1866 enddo
1867 enddo
1868
1869 enddo
1870 enddo
1871
1872 return
1873 end
1874c-----------------------------------------------------------------------
1875 subroutine crs_strs(u1,u2,u3,v1,v2,v3)
1876c Given an input vector v, this generates the H1 coarse-grid solution
1877 include 'SIZE'
1878 include 'DOMAIN'
1879 include 'INPUT'
1880 include 'GEOM'
1881 include 'SOLN'
1882 include 'PARALLEL'
1883 include 'CTIMER'
1884 include 'TSTEP'
1885
1886 real u1(1),u2(1),u3(1),v1(1),v2(1),v3(1)
1887
1888 common /scrpr3/ uc1(lcr*lelt),uc2(lcr*lelt),uc3(lcr*lelt)
1889 common /scrpr2/ vc1(lcr*lelt),vc2(lcr*lelt),vc3(lcr*lelt)
1890
1891 integer icalld1
1892 save icalld1
1893 data icalld1 /0/
1894
1895
1896 if (icalld1.eq.0) then ! timer info
1897 ncrsl=0
1898 tcrsl=0.0
1899 icalld1=1
1900 endif
1901 ncrsl = ncrsl + 1
1902
1903 n = nelv*lx1*ly1*lz1
1904 m = nelv*lcr
1905
1906 call map_f_to_c_h1_bilin(uc1,v1) ! additive Schwarz
1907 call map_f_to_c_h1_bilin(uc2,v2) ! additive Schwarz
1908 if (if3d) call map_f_to_c_h1_bilin(uc3,v3) ! additive Schwarz
1909
1910 k=0
1911 if (if3d) then
1912 do i=1,m
1913 vc1(k+1)=uc1(i)
1914 vc1(k+2)=uc2(i)
1915 vc1(k+3)=uc3(i)
1916 k=k+3
1917 enddo
1918 else
1919 do i=1,m
1920 vc1(k+1)=uc1(i)
1921 vc1(k+2)=uc2(i)
1922 k=k+2
1923 enddo
1924 endif
1925
1926 etime1=dnekclock()
1927 call fgslib_crs_solve(xxth_strs,uc1,vc1)
1928 tcrsl=tcrsl+dnekclock()-etime1
1929
1930 k=0
1931 if (if3d) then
1932 do i=1,m
1933 vc1(i)=uc1(k+1)
1934 vc2(i)=uc1(k+2)
1935 vc3(i)=uc1(k+3)
1936 k=k+3
1937 enddo
1938 else
1939 do i=1,m
1940 vc1(i)=uc1(k+1)
1941 vc2(i)=uc1(k+2)
1942 k=k+2
1943 enddo
1944 endif
1945 call map_c_to_f_h1_bilin(u1,vc1)
1946 call map_c_to_f_h1_bilin(u2,vc2)
1947 if (if3d) call map_c_to_f_h1_bilin(u3,vc3)
1948
1949 return
1950 end
1951c-----------------------------------------------------------------------
1952 subroutine set_up_h1_crs_strs(h1,h2,ifld,matmod)
1953
1954 include 'SIZE'
1955 include 'GEOM'
1956 include 'DOMAIN'
1957 include 'INPUT'
1958 include 'PARALLEL'
1959 include 'TSTEP'
1960 common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
1961
1962 common /ivrtx/ vertex ((2**ldim)*lelt)
1963 integer vertex
1964
1965 integer null_space,e
1966
1967 character*3 cb
1968 common /scrxxti/ ia(ldim*ldim*lcr*lcr*lelv)
1969 $ , ja(ldim*ldim*lcr*lcr*lelv)
1970
1971 parameter (lcc=2**ldim)
1972 common /scrcr1/ a(ldim*ldim*lcc*lcc*lelt)
1973 real mask(ldim,lcr,lelv)
1974 equivalence (mask,a)
1975
1976 integer*8 ngv
1977
1978 t0 = dnekclock()
1979
1980 nel = nelfld(ifld)
1981
1982 nxc = 2
1983 nzc = 1
1984 if (if3d) nzc=nxc
1985
1986 ncr = nxc**ldim
1987 mcr = ldim*ncr ! Dimension of unassembled sub-block
1988 n = mcr*nel
1989
1990c Set SEM_to_GLOB index
1991 call get_vertex
1992 call set_vert_strs(se_to_gcrs,ngv,nxc,nel,vertex,.true.)
1993
1994
1995c Set mask for full array
1996 call get_strs_mask (mask,nxc,nzc,nel) ! Set mask
1997 call set_jl_crs_mask (n,mask,se_to_gcrs)
1998
1999
2000c Setup local SEM-based Neumann operators (for now, just full...)
2001 call get_local_crs_galerkin_strs(a,ncr,nxc,h1,h2,matmod)
2002 call set_mat_ij(ia,ja,mcr,nel)
2003 nnz=mcr*mcr*nel
2004
2005 null_space=0
2006
2007 t1 = dnekclock()-t0
2008 if (nio.eq.0)
2009 $ write(6,*) 'start:: setup h1 coarse grid ',t1, ' sec'
2010 $ ,nnz,mcr,ncr,nel
2011
2012 do i=1,nnz
2013 write(44,44) ia(i),ja(i),a(i)
2014 enddo
2015 44 format(2i9,1pe22.13)
2016c stop
2017
2018 imode = param(40)
2019 call fgslib_crs_setup(xxth_strs,imode,nekcomm,mp,n,se_to_gcrs,
2020 $ nnz,ia,ja,a,null_space)
2021
2022 t0 = dnekclock()-t0
2023 if (nio.eq.0) then
2024 write(6,*) 'done :: setup h1 coarse grid ',t0, ' sec',xxth_strs
2025 write(6,*) ' '
2026 endif
2027
2028 return
2029 end
2030c-----------------------------------------------------------------------
2031 subroutine axsf_e_3d(au,av,aw,u,v,w,h1,h2,ur,e)
2032c
2033c du_i
2034c Compute the gradient tensor G_ij := ---- , for element e
2035c du_j
2036c
2037 include 'SIZE'
2038 include 'TOTAL'
2039
2040 real au(1),av(1),aw(1),u(1),v(1),w(1),h1(1),h2(1)
2041
2042 parameter (l=lx1*ly1*lz1)
2043 real ur(l,ldim,ldim)
2044
2045 integer e,p
2046
2047 p = lx1-1 ! Polynomial degree
2048 n = lx1*ly1*lz1
2049
2050
2051 call local_grad3(ur(1,1,1),ur(1,2,1),ur(1,3,1),u,p,1,dxm1,dxtm1)
2052 call local_grad3(ur(1,1,2),ur(1,2,2),ur(1,3,2),v,p,1,dxm1,dxtm1)
2053 call local_grad3(ur(1,1,3),ur(1,2,3),ur(1,3,3),w,p,1,dxm1,dxtm1)
2054
2055 do i=1,n
2056
2057 u1 = ur(i,1,1)*rxm1(i,1,1,e) + ur(i,2,1)*sxm1(i,1,1,e)
2058 $ + ur(i,3,1)*txm1(i,1,1,e)
2059 u2 = ur(i,1,1)*rym1(i,1,1,e) + ur(i,2,1)*sym1(i,1,1,e)
2060 $ + ur(i,3,1)*tym1(i,1,1,e)
2061 u3 = ur(i,1,1)*rzm1(i,1,1,e) + ur(i,2,1)*szm1(i,1,1,e)
2062 $ + ur(i,3,1)*tzm1(i,1,1,e)
2063
2064 v1 = ur(i,1,2)*rxm1(i,1,1,e) + ur(i,2,2)*sxm1(i,1,1,e)
2065 $ + ur(i,3,2)*txm1(i,1,1,e)
2066 v2 = ur(i,1,2)*rym1(i,1,1,e) + ur(i,2,2)*sym1(i,1,1,e)
2067 $ + ur(i,3,2)*tym1(i,1,1,e)
2068 v3 = ur(i,1,2)*rzm1(i,1,1,e) + ur(i,2,2)*szm1(i,1,1,e)
2069 $ + ur(i,3,2)*tzm1(i,1,1,e)
2070
2071 w1 = ur(i,1,3)*rxm1(i,1,1,e) + ur(i,2,3)*sxm1(i,1,1,e)
2072 $ + ur(i,3,3)*txm1(i,1,1,e)
2073 w2 = ur(i,1,3)*rym1(i,1,1,e) + ur(i,2,3)*sym1(i,1,1,e)
2074 $ + ur(i,3,3)*tym1(i,1,1,e)
2075 w3 = ur(i,1,3)*rzm1(i,1,1,e) + ur(i,2,3)*szm1(i,1,1,e)
2076 $ + ur(i,3,3)*tzm1(i,1,1,e)
2077
2078 dj = h1(i)*w3m1(i,1,1)*jacmi(i,e)
2079 s11 = dj*(u1 + u1)! S_ij
2080 s12 = dj*(u2 + v1)
2081 s13 = dj*(u3 + w1)
2082 s21 = dj*(v1 + u2)
2083 s22 = dj*(v2 + v2)
2084 s23 = dj*(v3 + w2)
2085 s31 = dj*(w1 + u3)
2086 s32 = dj*(w2 + v3)
2087 s33 = dj*(w3 + w3)
2088
2089c Sum_j : (r_p/x_j) h1 J S_ij
2090
2091 ur(i,1,1)=rxm1(i,1,1,e)*s11+rym1(i,1,1,e)*s12+rzm1(i,1,1,e)*s13
2092 ur(i,2,1)=sxm1(i,1,1,e)*s11+sym1(i,1,1,e)*s12+szm1(i,1,1,e)*s13
2093 ur(i,3,1)=txm1(i,1,1,e)*s11+tym1(i,1,1,e)*s12+tzm1(i,1,1,e)*s13
2094
2095 ur(i,1,2)=rxm1(i,1,1,e)*s21+rym1(i,1,1,e)*s22+rzm1(i,1,1,e)*s23
2096 ur(i,2,2)=sxm1(i,1,1,e)*s21+sym1(i,1,1,e)*s22+szm1(i,1,1,e)*s23
2097 ur(i,3,2)=txm1(i,1,1,e)*s21+tym1(i,1,1,e)*s22+tzm1(i,1,1,e)*s23
2098
2099 ur(i,1,3)=rxm1(i,1,1,e)*s31+rym1(i,1,1,e)*s32+rzm1(i,1,1,e)*s33
2100 ur(i,2,3)=sxm1(i,1,1,e)*s31+sym1(i,1,1,e)*s32+szm1(i,1,1,e)*s33
2101 ur(i,3,3)=txm1(i,1,1,e)*s31+tym1(i,1,1,e)*s32+tzm1(i,1,1,e)*s33
2102
2103 enddo
2104 call local_grad3_t
2105 $ (au,ur(1,1,1),ur(1,2,1),ur(1,3,1),p,1,dxm1,dxtm1,av)
2106 call local_grad3_t
2107 $ (av,ur(1,1,2),ur(1,2,2),ur(1,3,2),p,1,dxm1,dxtm1,ur)
2108 call local_grad3_t
2109 $ (aw,ur(1,1,3),ur(1,2,3),ur(1,3,3),p,1,dxm1,dxtm1,ur)
2110
2111 do i=1,n
2112 au(i)=au(i) + h2(i)*bm1(i,1,1,e)*u(i)
2113 av(i)=av(i) + h2(i)*bm1(i,1,1,e)*v(i)
2114 aw(i)=aw(i) + h2(i)*bm1(i,1,1,e)*w(i)
2115 enddo
2116
2117 return
2118 end
2119c-----------------------------------------------------------------------
2120 subroutine axsf_e_2d(au,av,u,v,h1,h2,ur,e)
2121c
2122c du_i
2123c Compute the gradient tensor G_ij := ---- , for element e
2124c du_j
2125c
2126 include 'SIZE'
2127 include 'TOTAL'
2128
2129 real au(1),av(1),u(1),v(1),h1(1),h2(1)
2130
2131 parameter (l=lx1*ly1*lz1)
2132 real ur(l,ldim,ldim)
2133
2134 integer e,p
2135
2136 p = lx1-1 ! Polynomial degree
2137 n = lx1*ly1*lz1
2138
2139
2140 call local_grad2(ur(1,1,1),ur(1,2,1),u,p,1,dxm1,dxtm1)
2141 call local_grad2(ur(1,1,2),ur(1,2,2),v,p,1,dxm1,dxtm1)
2142
2143 do i=1,n
2144 dj = h1(i)*w3m1(i,1,1)*jacmi(i,e)
2145
2146 u1 = ur(i,1,1)*rxm1(i,1,1,e) + ur(i,2,1)*sxm1(i,1,1,e) !ux
2147 u2 = ur(i,1,1)*rym1(i,1,1,e) + ur(i,2,1)*sym1(i,1,1,e) !uy
2148 v1 = ur(i,1,2)*rxm1(i,1,1,e) + ur(i,2,2)*sxm1(i,1,1,e) !vx
2149 v2 = ur(i,1,2)*rym1(i,1,1,e) + ur(i,2,2)*sym1(i,1,1,e) !vy
2150
2151 s11 = dj*( u1 + u1 ) ! h1*rho*S_ij
2152 s12 = dj*( u2 + v1 )
2153 s21 = dj*( v1 + u2 )
2154 s22 = dj*( v2 + v2 )
2155
2156c Sum_j : (r_k/x_j) h1 J S_ij
2157
2158 ur(i,1,1)=rxm1(i,1,1,e)*s11+rym1(i,1,1,e)*s12 ! i=1,k=1
2159 ur(i,2,1)=sxm1(i,1,1,e)*s11+sym1(i,1,1,e)*s12 ! i=1,k=2
2160
2161 ur(i,1,2)=rxm1(i,1,1,e)*s21+rym1(i,1,1,e)*s22 ! i=2,k=1
2162 ur(i,2,2)=sxm1(i,1,1,e)*s21+sym1(i,1,1,e)*s22 ! i=2,k=2
2163
2164 enddo
2165
2166 call local_grad2_t(au,ur(1,1,1),ur(1,2,1),p,1,dxm1,dxtm1,av)
2167 call local_grad2_t(av,ur(1,1,2),ur(1,2,2),p,1,dxm1,dxtm1,ur)
2168
2169 do i=1,n
2170 au(i)=au(i) + h2(i)*bm1(i,1,1,e)*u(i)
2171 av(i)=av(i) + h2(i)*bm1(i,1,1,e)*v(i)
2172 enddo
2173
2174 return
2175 end
2176c-----------------------------------------------------------------------
2177 subroutine axsf_fast(au,av,aw,u,v,w,h1,h2,ifld)
2178 include 'SIZE'
2179 include 'TOTAL'
2180
2181 parameter (l=lx1*ly1*lz1)
2182 real au(l,1),av(l,1),aw(l,1),u(l,1),v(l,1),w(l,1),h1(l,1),h2(l,1)
2183
2184 common /btmp0/ ur(l,ldim,ldim)
2185
2186 integer e
2187
2188 nel = nelfld(ifld)
2189
2190 if (if3d) then
2191 do e=1,nel
2192 call axsf_e_3d(au(1,e),av(1,e),aw(1,e),u(1,e),v(1,e),w(1,e)
2193 $ ,h1(1,e),h2(1,e),ur,e)
2194 enddo
2195 else
2196 do e=1,nel
2197 call axsf_e_2d(au(1,e),av(1,e),u(1,e),v(1,e)
2198 $ ,h1(1,e),h2(1,e),ur,e)
2199 enddo
2200 endif
2201
2202 return
2203 end
2204c-----------------------------------------------------------------------
2205 subroutine ttxyz (ff,tx,ty,tz,nel)
2206C
2207 include 'SIZE'
2208 include 'DXYZ'
2209 include 'GEOM'
2210 include 'INPUT'
2211 include 'TSTEP'
2212 include 'WZ'
2213
2214 DIMENSION TX(LX1,LY1,LZ1,1)
2215 $ , TY(LX1,LY1,LZ1,1)
2216 $ , TZ(LX1,LY1,LZ1,1)
2217 $ , FF(LX1*LY1*LZ1,1)
2218
2219 common /scrsf/ fr(lx1*ly1*lz1,lelt)
2220 $ , fs(lx1*ly1*lz1,lelt)
2221 $ , ft(lx1*ly1*lz1,lelt)
2222 real wa(lx1,ly1,lz1,lelt)
2223 equivalence (wa,ft)
2224 real ys(lx1)
2225
2226 NXYZ1 = lx1*ly1*lz1
2227 NTOT1 = NXYZ1*NEL
2228
2229 CALL COL3 (FR,RXM1,TX,NTOT1)
2230 CALL ADDCOL3 (FR,RYM1,TY,NTOT1)
2231 CALL COL3 (FS,SXM1,TX,NTOT1)
2232 CALL ADDCOL3 (FS,SYM1,TY,NTOT1)
2233
2234 IF (ldim.EQ.3) THEN
2235 CALL ADDCOL3 (FR,RZM1,TZ,NTOT1)
2236 CALL ADDCOL3 (FS,SZM1,TZ,NTOT1)
2237 CALL COL3 (FT,TXM1,TX,NTOT1)
2238 CALL ADDCOL3 (FT,TYM1,TY,NTOT1)
2239 CALL ADDCOL3 (FT,TZM1,TZ,NTOT1)
2240 endif
2241C
2242 IF (IFAXIS) THEN
2243 DO 100 IEL=1,NEL
2244 IF ( IFRZER(IEL) ) THEN
2245 CALL MXM (YM1(1,1,1,IEL),lx1,DATM1,ly1,YS,1)
2246 DO 120 IX=1,lx1
2247 IY = 1
2248 WA(IX,IY,1,IEL)=YS(IX)*W2AM1(IX,IY)
2249 DO 120 IY=2,ly1
2250 DNR = 1.0 + ZAM1(IY)
2251 WA(IX,IY,1,IEL)=YM1(IX,IY,1,IEL)*W2AM1(IX,IY)/DNR
2252 120 CONTINUE
2253 ELSE
2254 CALL COL3 (WA(1,1,1,IEL),YM1(1,1,1,IEL),W2CM1,NXYZ1)
2255 endif
2256 100 CONTINUE
2257 CALL COL2 (FR,WA,NTOT1)
2258 CALL COL2 (FS,WA,NTOT1)
2259 else
2260 do 180 iel=1,nel
2261 call col2(fr(1,iel),w3m1,nxyz1)
2262 call col2(fs(1,iel),w3m1,nxyz1)
2263 call col2(ft(1,iel),w3m1,nxyz1)
2264 180 continue
2265 endif
2266
2267
2268 DO 200 IEL=1,NEL
2269 IF (IFAXIS) CALL SETAXDY ( IFRZER(IEL) )
2270 CALL TTRST (FF(1,IEL),FR(1,IEL),FS(1,IEL),
2271 $ FT(1,IEL),FR(1,IEL)) ! FR work array
2272 200 CONTINUE
2273C
2274 return
2275 end
2276c-----------------------------------------------------------------------
2277 subroutine ttrst (ff,fr,fs,ft,ta)
2278
2279 include 'SIZE'
2280 include 'DXYZ'
2281 include 'TSTEP'
2282
2283 DIMENSION FF(LX1,LY1,LZ1)
2284 $ , FR(LX1,LY1,LZ1)
2285 $ , FS(LX1,LY1,LZ1)
2286 $ , FT(LX1,LY1,LZ1)
2287 $ , TA(LX1,LY1,LZ1)
2288
2289 NXY1 = lx1*ly1
2290 NYZ1 = ly1*lz1
2291 NXYZ1 = NXY1*lz1
2292
2293 CALL MXM (DXTM1,lx1,FR,lx1,FF,NYZ1)
2294 IF (ldim.EQ.2) THEN
2295 CALL MXM (FS,lx1,DYM1,ly1,TA,ly1)
2296 CALL ADD2 (FF,TA,NXYZ1)
2297 ELSE
2298 DO 10 IZ=1,lz1
2299 CALL MXM (FS(1,1,IZ),lx1,DYM1,ly1,TA(1,1,IZ),ly1)
2300 10 CONTINUE
2301 CALL ADD2 (FF,TA,NXYZ1)
2302 CALL MXM (FT,NXY1,DZM1,lz1,TA,lz1)
2303 CALL ADD2 (FF,TA,NXYZ1)
2304 endif
2305
2306 return
2307 end
2308c-----------------------------------------------------------------------
2309 subroutine axitzz (vfy,tzz,nel)
2310C
2311 include 'SIZE'
2312 include 'DXYZ'
2313 include 'GEOM'
2314 include 'WZ'
2315 common /ctmp0/ phi(lx1,ly1)
2316C
2317 DIMENSION VFY(LX1,LY1,LZ1,1)
2318 $ , TZZ(LX1,LY1,LZ1,1)
2319C
2320 NXYZ1 = lx1*ly1*lz1
2321C
2322 DO 100 IEL=1,NEL
2323 CALL SETAXW1 ( IFRZER(IEL) )
2324 CALL COL4 (PHI,TZZ(1,1,1,IEL),JACM1(1,1,1,IEL),W3M1,NXYZ1)
2325 IF ( IFRZER(IEL) ) THEN
2326 DO 120 IX=1,lx1
2327 DO 120 IY=2,ly1
2328 DNR = PHI(IX,IY)/( 1.0 + ZAM1(IY) )
2329 DDS = WXM1(IX) * WAM1(1) * DATM1(IY,1) *
2330 $ JACM1(IX,1,1,IEL) * TZZ(IX,1,1,IEL)
2331 VFY(IX,IY,1,IEL)=VFY(IX,IY,1,IEL) + DNR + DDS
2332 120 CONTINUE
2333 ELSE
2334 CALL ADD2 (VFY(1,1,1,IEL),PHI,NXYZ1)
2335 endif
2336 100 CONTINUE
2337C
2338 return
2339 end
2340c-----------------------------------------------------------------------
2341 subroutine setaxdy (ifaxdy)
2342C
2343 include 'SIZE'
2344 include 'DXYZ'
2345C
2346 LOGICAL IFAXDY
2347C
2348 IF (IFAXDY) THEN
2349 CALL COPY (DYM1 ,DAM1 ,ly1*ly1)
2350 CALL COPY (DYTM1,DATM1,ly1*ly1)
2351 ELSE
2352 CALL COPY (DYM1 ,DCM1 ,ly1*ly1)
2353 CALL COPY (DYTM1,DCTM1,ly1*ly1)
2354 endif
2355C
2356 return
2357 end
2358c-----------------------------------------------------------------------
2359 function opnorm2w(v1,v2,v3,w)
2360 include 'SIZE'
2361 include 'TOTAL'
2362c
2363 real v1(1) , v2(1), v3(1), w(1)
2364
2365 n=lx1*ly1*lz1*nelv
2366 s=vlsc3(v1,w,v1,n)+vlsc3(v2,w,v2,n)
2367 if(if3d) s=s+vlsc3(v3,w,v3,n)
2368 s=glsum(s,1)
2369
2370 if (s.gt.0) s=sqrt(s/volvm1)
2371 opnorm2w = s
2372
2373 return
2374 end
2375c-----------------------------------------------------------------------
2376 subroutine strs_project_a(b1,b2,b3,h1,h2,wt,ifld,ierr,matmod)
2377
2378c Assumes if uservp is true and thus reorthogonalizes every step
2379
2380 include 'SIZE'
2381 include 'TOTAL'
2382 include 'ORTHOSTRS' ! Storage of approximation space
2383 include 'CTIMER'
2384
2385 real b1(1),b2(1),b3(1),h1(1),h2(1),wt(1)
2386 common /ctmp1/ w(lx1*ly1*lz1*lelt,ldim)
2387 real l2a,l2b
2388
2389 kmax = napproxstrs(1)
2390 k = napproxstrs(2)
2391 n = lx1*ly1*lz1*nelv
2392 m = n*ldim
2393
2394 if (k.eq.0.or.kmax.eq.0) return
2395
2396 etime0 = dnekclock()
2397
2398 l2b=opnorm2w(b1,b2,b3,binvm1)
2399
2400c Reorthogonalize basis
2401
2402 dh1max=difmax(bstrs(1 ),h1,n)
2403 call col3(bstrs,h2,bm1,n)
2404 dh2max=difmax(bstrs(1+n),bstrs,n)
2405
2406 if (dh1max.gt.0.or.dh2max.gt.0) then
2407 call strs_ortho_all(xstrs(1+m),bstrs(1+m),n,k,h1,h2,wt,ifld,w
2408 $ ,ierr,matmod)
2409 else
2410 call strs_ortho_one(xstrs(1+m),bstrs(1+m),n,k,h1,h2,wt,ifld,w
2411 $ ,ierr,matmod)
2412 endif
2413
2414 napproxstrs(2) = k
2415
2416 call opcopy(bstrs(1),bstrs(1+n),bstrs(1+2*n),b1,b2,b3)
2417 call opzero(xstrs(1),xstrs(1+n),xstrs(1+2*n))
2418
2419 do i=1,k
2420 i1 = 1 + 0*n + (i-1)*m + m
2421 i2 = 1 + 1*n + (i-1)*m + m
2422 i3 = 1 + 2*n + (i-1)*m + m
2423 alpha=op_glsc2_wt(bstrs(1),bstrs(1+n),bstrs(1+2*n)
2424 $ ,xstrs(i1),xstrs(i2),xstrs(i3),wt)
2425
2426 alphm=-alpha
2427 call opadds(bstrs(1),bstrs(1+n),bstrs(1+2*n)
2428 $ ,bstrs(i1),bstrs(i2),bstrs(i3),alphm,n,2)
2429
2430 call opadds(xstrs(1),xstrs(1+n),xstrs(1+2*n)
2431 $ ,xstrs(i1),xstrs(i2),xstrs(i3),alpha,n,2)
2432
2433
2434 enddo
2435
2436 call opcopy(b1,b2,b3,bstrs(1),bstrs(1+n),bstrs(1+2*n))
2437 l2a=opnorm2w(b1,b2,b3,binvm1)
2438
2439 call copy(bstrs(1 ),h1,n) ! Save h1 and h2 for comparison
2440 call col3(bstrs(1+n),bm1,h2,n)
2441
2442 ratio=l2b/l2a
2443 if (nio.eq.0) write(6,1) istep,' Project ' // 'HELM3 ',
2444 $ l2a,l2b,ratio,k,kmax
2445 1 format(i11,a,6x,1p3e13.4,i4,i4)
2446
2447c if (ierr.ne.0) call exitti(' h3proj quit$',ierr)
2448
2449 tproj = tproj + dnekclock() - etime0
2450
2451 return
2452 end
2453c-----------------------------------------------------------------------
2454 subroutine strs_project_b(x1,x2,x3,h1,h2,wt,ifld,ierr)
2455
2456c Reconstruct solution; don't bother to orthonomalize bases
2457
2458 include 'SIZE'
2459 include 'TOTAL'
2460 include 'ORTHOSTRS' ! Storage of approximation space
2461 include 'CTIMER'
2462
2463 real x1(1),x2(1),x3(1),h1(1),h2(1),wt(1)
2464 common /cptst/ xs(lx1*ly1*lz1*lelt*ldim)
2465
2466 kmax = napproxstrs(1)
2467 k = napproxstrs(2)
2468 n = lx1*ly1*lz1*nelv
2469 m = n*ldim
2470
2471 if (kmax.eq.0) return
2472
2473 etime0 = dnekclock()
2474
2475 if (k.eq.0) then ! _
2476 call opadd2(x1,x2,x3,xstrs(1),xstrs(1+n),xstrs(1+2*n)) ! x=dx+x
2477 k=1
2478 k1 = 1 + 0*n + (k-1)*m + m
2479 k2 = 1 + 1*n + (k-1)*m + m
2480 k3 = 1 + 2*n + (k-1)*m + m
2481 call opcopy(xstrs(k1),xstrs(k2),xstrs(k3),x1,x2,x3) ! x1=x^n
2482 elseif (k.eq.kmax) then ! _
2483 call opadd2(x1,x2,x3,xstrs(1),xstrs(1+n),xstrs(1+2*n)) ! x=dx+x
2484 k=1
2485 k1 = 1 + 0*n + (k-1)*m + m
2486 k2 = 1 + 1*n + (k-1)*m + m
2487 k3 = 1 + 2*n + (k-1)*m + m
2488 call opcopy(xstrs(k1),xstrs(k2),xstrs(k3),x1,x2,x3) ! x1=x^n
2489c k=2
2490c k1 = 1 + 0*n + (k-1)*m + m
2491c k2 = 1 + 1*n + (k-1)*m + m
2492c k3 = 1 + 2*n + (k-1)*m + m
2493c call opcopy(xstrs(k1),xstrs(k2),xstrs(k3),xs(1),xs(1+n),xs(1+2*n))
2494 else
2495 k=k+1
2496 k1 = 1 + 0*n + (k-1)*m + m
2497 k2 = 1 + 1*n + (k-1)*m + m
2498 k3 = 1 + 2*n + (k-1)*m + m
2499 call opcopy(xstrs(k1),xstrs(k2),xstrs(k3),x1,x2,x3) ! xk=dx _
2500 call opadd2(x1,x2,x3,xstrs(1),xstrs(1+n),xstrs(1+2*n)) ! x=dx + x
2501 endif
2502
2503c if (k.eq.kmax) call opcopy(xs(1),xs(1+n),xs(1+2*n),x1,x2,x3) ! presave
2504
2505 napproxstrs(2)=k
2506
2507 tproj = tproj + dnekclock() - etime0
2508
2509 return
2510 end
2511c-----------------------------------------------------------------------
2512 subroutine strs_orthok(x,b,n,k,h1,h2,wt,ifld,w,ierr,matmod)
2513
2514c Orthonormalize the kth element of X against x_j, j < k.
2515
2516 include 'SIZE'
2517 include 'TOTAL'
2518
2519 real x(n,ldim,k),b(n,ldim,k),h1(n),h2(n),wt(n),w(n,ldim)
2520 real al(mxprev),bt(mxprev)
2521
2522 m = n*ldim ! vector length
2523 nel = nelfld(ifld)
2524
2525 call axhmsf (b(1,1,k),b(1,2,k),b(1,3,k)
2526 $ ,x(1,1,k),x(1,2,k),x(1,3,k),h1,h2,matmod)
2527 call rmask (b(1,1,k),b(1,2,k),b(1,3,k),nel)
2528 call opdssum (b(1,1,k),b(1,2,k),b(1,3,k))
2529
2530 xax0 = op_glsc2_wt(b(1,1,k),b(1,2,k),b(1,3,k)
2531 $ ,x(1,1,k),x(1,2,k),x(1,3,k),wt)
2532 ierr = 3
2533
2534 if (xax0.le.0.and.nio.eq.0)write(6,*)istep,ierr,k,xax0,'Proj3 ERR'
2535 if (xax0.le.0) return
2536
2537 s = 0.
2538 do j=1,k-1! Modifed Gram-Schmidt
2539
2540 betaj = ( op_vlsc2_wt(b(1,1,j),b(1,2,j),b(1,3,j)
2541 $ ,x(1,1,k),x(1,2,k),x(1,3,k),wt)
2542 $ + op_vlsc2_wt(b(1,1,k),b(1,2,k),b(1,3,k)
2543 $ ,x(1,1,j),x(1,2,j),x(1,3,j),wt))/2.
2544 betam = -glsum(betaj,1)
2545 call add2s2 (x(1,1,k),x(1,1,j),betam,m) ! Full-vector-subtract
2546 call add2s2 (b(1,1,k),b(1,1,j),betam,m) !
2547
2548 s = s + betam**2
2549
2550 enddo
2551
2552 xax1 = xax0-s
2553 xax2 = op_glsc2_wt(b(1,1,k),b(1,2,k),b(1,3,k)
2554 $ ,x(1,1,k),x(1,2,k),x(1,3,k),wt)
2555 scale = xax2
2556
2557 eps = 1.e-8
2558 ierr = 0
2559 if (scale/xax0.lt.eps) ierr=1
2560
2561c if(nio.eq.0) write(6,*) istep,ierr,k,scale,xax0,' SCALE'
2562c if(nio.eq.0.and.(istep.lt.10.or.mod(istep,100).eq.0.or.ierr.gt.0))
2563c $write(6,3) istep,k,ierr,xax1/xax0,xax2/xax0
2564c 3 format(i9,2i3,1p2e12.4,' scale ortho')
2565
2566 if (scale.gt.0) then
2567 scale = 1./sqrt(scale)
2568 call cmult(x(1,1,k),scale,m)
2569 call cmult(b(1,1,k),scale,m)
2570 else
2571 ierr=2
2572 endif
2573
2574 return
2575 end
2576c-----------------------------------------------------------------------
2577 subroutine strs_ortho_one(x,b,n,k,h1,h2,wt,ifld,w,ierr,matmod)
2578
2579 include 'SIZE'
2580 include 'TOTAL'
2581
2582 real x(n,ldim,k),b(n,ldim,k),h1(n),h2(n),wt(n),w(n,ldim)
2583
2584 m = n*ldim
2585
2586 js=k
2587 do j=k,k
2588 if (js.lt.j) call copy(x(1,1,js),x(1,1,j),m)
2589 call strs_orthok(x,b,n,js,h1,h2,wt,ifld,w,ierr,matmod)
2590 if (ierr.eq.0) js=js+1
2591 enddo
2592 k = js-1
2593
2594 return
2595 end
2596c-----------------------------------------------------------------------
2597 subroutine strs_ortho_all(x,b,n,k,h1,h2,wt,ifld,w,ierr,matmod)
2598
2599 include 'SIZE'
2600 include 'TOTAL'
2601
2602 real x(n,ldim,k),b(n,ldim,k),h1(n),h2(n),wt(n),w(n,ldim)
2603
2604 m = n*ldim
2605
2606 js=1
2607 do j=1,k
2608 if (js.lt.j) call copy(x(1,1,js),x(1,1,j),m)
2609 call strs_orthok(x,b,n,js,h1,h2,wt,ifld,w,ierr,matmod)
2610 if (ierr.eq.0) js=js+1
2611 enddo
2612 k = js-1
2613
2614 return
2615 end
2616c-----------------------------------------------------------------------
2617 subroutine setprop
2618C------------------------------------------------------------------------
2619C
2620C Set variable property arrays
2621C
2622C------------------------------------------------------------------------
2623 include 'SIZE'
2624 include 'TOTAL'
2625 include 'CTIMER'
2626C
2627C Caution: 2nd and 3rd strainrate invariants residing in scratch
2628C common /SCREV/ are used in STNRINV and NEKASGN
2629C
2630 common /screv/ sii (lx1,ly1,lz1,lelt),siii(lx1,ly1,lz1,lelt)
2631
2632 if (nio.eq.0.and.loglevel.gt.2)
2633 $ write(6,*) 'setprop'
2634
2635#ifdef TIMER
2636 if (icalld.eq.0) tspro=0.0
2637 icalld=icalld+1
2638 nspro=icalld
2639 etime1=dnekclock()
2640#endif
2641
2642 NXYZ1 = lx1*ly1*lz1
2643 MFIELD=2
2644 IF (IFFLOW) MFIELD=1
2645 nfldt = nfield
2646 if (ifmhd) nfldt = nfield+1
2647
2648 ifld = ifield
2649
2650 do ifield=mfield,nfldt
2651 if (idpss(ifield-1).eq.-1) goto 100
2652 call vprops
2653 100 continue
2654 enddo
2655
2656 ifield = ifld
2657
2658#ifdef TIMER
2659 tspro=tspro+(dnekclock()-etime1)
2660#endif
2661
2662 return
2663 end
2664c-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.