integration_example.usr 16.9 KB
Newer Older
Cody Permann's avatar
Cody Permann committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
c-----------------------------------------------------------------------
C
C  USER SPECIFIED ROUTINES:
C
C     - boundary conditions
C     - initial conditions
C     - variable properties
C     - local acceleration for fluid (a)
C     - forcing function for passive scalar (q)
C     - general purpose routine for checking errors etc.
C
c-----------------------------------------------------------------------
      subroutine uservp (ix,iy,iz,eg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
      integer e,f,eg
      udiff =0.
      utrans=0.
      return
      end
c-----------------------------------------------------------------------
23
24
25
! Sets the force term in the momentum equation. Use this to set body
! forces such as gravity. Note that ffx, ffy, and ffz will later be
! multiplied by the density.
Cody Permann's avatar
Cody Permann committed
26
27
28
29
30
31
      subroutine userf  (ix,iy,iz,eg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
      ffx=0.0
      ffy=0.0
32
      ffz=0.0
Cody Permann's avatar
Cody Permann committed
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
      return
      end
c-----------------------------------------------------------------------
      subroutine userq  (ix,iy,iz,eg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer e,f,eg
c     e = gllel(eg)

      qvol   = 0.0

      return
      end
c-----------------------------------------------------------------------
      subroutine userchk
      include 'SIZE'
      include 'TOTAL'
52
53
54
55
56
57
58
      parameter (nl_max=100)
      parameter (nf_max=100)
      common/expansion_tdata/n_legendre, m_fourier
      common/expansion_tcoef/coeff_tij(nl_max,nf_max)
      common/expansion_fcoef/coeff_fij(nl_max,nf_max)
      common/expansion_recon/flux_recon(lx1,ly1,lz1,lelt)
      common/test_passing/ flux_moose, temp_nek 
Cody Permann's avatar
Cody Permann committed
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
      real sint, sint1, sarea, sarea1, wtmp
      integer e, f 

      n1=nelt*lx1*ly1*lz1
      n2=nelt*lx2*ly2*lz2

      pmin=glmin(pr,n2)
      pmax=glmax(pr,n2)
      wmin=glmin(vz,n1)
      wmax=glmax(vz,n1)
      tmin=glmin(t,n1)
      tmax=glmax(t,n1)

      ifflow=.false. 

74
75
76
c      if (istep.eq.0) then
c         call rzero(t,n1) 
c      endif
77

Cody Permann's avatar
Cody Permann committed
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
      sint1=0.0
      sarea1=0.0

      do e=1,lelt
        do f=1,6
          call surface_int(sint,sarea,t,e,f)
          if (cbc(f,e,1).eq.'W  ') then
           sint1=sint1+sint  
           sarea1=sarea1+sarea
          endif          
         enddo 
      enddo

      call  gop(sint1,wtmp,'+  ',1)
      call  gop(sarea1,wtmp,'+  ',1)

      temp_nek=sint1/sarea1
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112

      sint1=0.0
      sarea1=0.0

      do e=1,lelt
        do f=1,6
          call surface_int(sint,sarea,flux_recon,e,f)
          if (cbc(f,e,1).eq.'W  ') then
           sint1=sint1+sint
           sarea1=sarea1+sarea
          endif
         enddo
      enddo

      call  gop(sint1,wtmp,'+  ',1)
      call  gop(sarea1,wtmp,'+  ',1)

      flux_moose=sint1/sarea1
113
114

      call heat_balance(flux_moose)
Cody Permann's avatar
Cody Permann committed
115
116
117
      
      if (nid.eq.0) then
         write(6,*)"*** Temperature: ",tmin," - ",tmax
118
119
         write(6,*)"*** Average Temperature: ",temp_nek,coeff_tij(1,1)
         write(6,*)"*** Average Flux: ",flux_moose,coeff_fij(1,1)
Cody Permann's avatar
Cody Permann committed
120
121
      endif

122
c     Will this be overwritten -----------------
123
124
125
c      write(6,*)"*** Setting polynomial orders ..." 
      n_legendre=10
      m_fourier=5
126
127
128
129
130
131
132
133
134
135
136
c-----------------------------------------------

c This is for testing -----------
c      call flux_reconstruction()
c      call nek_expansion()
c      call nek_diag()
c      if (nid.eq.0) write(6,*)" *** Expansion done"
c      call nek_testp() 
c      call exitt()
c--------------------------------

137
138
139
140
141
142
143
144
c This is to test coefficients with MOOSE --
c      do i0=1,n_legendre
c        do j0=1,m_fourier
c        if (nid.eq.0) write(6,*)"C(",i0,",",j0,")=",
c     &      coeff_fij(i0,j0)
c        enddo 
c      enddo
c------------------------------------------
Cody Permann's avatar
Cody Permann committed
145
146
147
148
149
150
151
152
153
      return
      end
c-----------------------------------------------------------------------
      subroutine userbc (ix,iy,iz,iside,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'
      
      common/test_passing/ flux_moose, temp_nek
154
      common/expansion_recon/flux_recon(lx1,ly1,lz1,lelt)
Cody Permann's avatar
Cody Permann committed
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169

      integer e,ieg
      real ucx, ucy, ucz, ucy_e, yy 

      e=gllel(ieg)

      zz  = zm1(ix,iy,iz,e)
      xx  = xm1(ix,iy,iz,e)
      yy  = ym1(ix,iy,iz,e)
      rr  = sqrt(xx**2 + yy**2)

      ux = 0.0
      uy = 0.0
      uz = 1.0
      temp = 0.0
170
      flux = flux_recon(ix,iy,iz,e) !flux_moose
Cody Permann's avatar
Cody Permann committed
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247

      return
      end
c-----------------------------------------------------------------------
      subroutine useric (ix,iy,iz,ieg)
      include 'SIZE'
      include 'TOTAL'
      include 'NEKUSE'

      integer idum, e
      save    idum
      data    idum / 0 /

      if (idum.eq.0) idum = 99 + nid
      eps = .35

      uz=0.0
      uy=0.0
      ux=0.0
      temp=0

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat
      include 'SIZE'
      include 'TOTAL'

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat2
      include 'SIZE'
      include 'TOTAL'

      param(66) = 4.   ! These give the std nek binary i/o and are 
      param(67) = 4.   ! good default values
      ifuservp=.false.

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat3
      include 'SIZE'
      include 'TOTAL'
c
      return
      end

C=======================================================================
      subroutine nek_init_step()
      include 'SIZE'
      include 'TSTEP'
      include 'INPUT'
      include 'CTIMER'
      real*4 papi_mflops
      integer*8 papi_flops
      integer icall, kstep, i, pstep
      common /cht_coupler/ pstep
      save icall
 
      if (icall.eq.0) then
      call nekgsync()
      if (instep.eq.0) then
        if(nid.eq.0) write(6,'(/,A,/,A,/)')
     &     ' nsteps=0 -> skip time loop',
     &     ' running solver in post processing mode'
      else
        if(nio.eq.0) write(6,'(/,A,/)') 'Starting time loop ...'
      endif
      isyc  = 0
      itime = 0
      if(ifsync) isyc=1
      itime = 1
      call nek_comm_settings(isyc,itime)
      call nek_comm_startstat()
      istep  = 0
248
      icall  = 1
Cody Permann's avatar
Cody Permann committed
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
      endif

      istep=istep+1

      if (lastep .eq. 1) then
        pstep=2
      else
        call nek_advance
        pstep=2
      endif

      return
      end
C=======================================================================
      subroutine nek_step()

      include 'SIZE'
      include 'TSTEP'
      include 'INPUT'
      include 'CTIMER'
      common /cht_coupler/ pstep
      integer pstep

      pstep=pstep+1
273
      call heat(pstep)
Cody Permann's avatar
Cody Permann committed
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300

      return
      end
C=======================================================================
      subroutine nek_finalize_step()
      include 'SIZE'
      include 'TSTEP'
      include 'INPUT'
      include 'CTIMER'
      common /cht_coupler/ pstep
      integer pstep 

      real*4 papi_mflops
      integer*8 papi_flops
      integer icall, kstep, knstep, i

      if (param(103).gt.0)   call q_filter      (param(103))
      call setup_convect (pstep)  ! Save convective velocity _after_ filter

      call userchk
      call prepost (.false.,'his')
      call in_situ_check()

      if (mod(istep,nsteps).eq.0) lastep=1
      call nek_comm_settings(isyc,0)
      call comment

301
302
303
304
305
306
307
308
309
      return
      end
C=======================================================================
      subroutine nek_expansion()

      parameter (nl_max=100)
      parameter (nf_max=100)

      include 'SIZE'
310
      include 'TOTAL'
311
      common/expansion_tdata/n_legendre, m_fourier
312
      common/expansion_tcoef/coeff_tij(nl_max,nf_max)
313
314
      integer e,f

315
316
317
      real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
      real*8 sint, sint1, sarea, sarea1
      real*8 pi
318

319
      pi=4.0*atan(1.0)
320
321
322
      ntot=nx1*ny1*nz1*nelt

      call rzero(fmode,ntot)
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
      do i0=1,n_legendre
        do j0=1,m_fourier
          call nek_mode(fmode,i0,j0)
          sarea1=0.0
          sint1= 0.0
          do e=1,nelt
            do f=1,6
              sint=0.0
              sarea=0.0
              if (cbc(f,e,1).eq.'W  ') then
                call col3(cache,fmode,t,ntot)
                call surface_int(sint,sarea,cache,e,f)
                sint1=sint1+sint
                sarea1=sarea1+sarea
              endif
            enddo
          enddo
340
          call  gop(sint1,wtmp,'+  ',1)
341
342
343
344
345
          call  gop(sarea1,wtmp,'+  ',1)
!
          coeff_tij(i0,j0)=sint1*4.0*pi/sarea1
!
        enddo
346
347
      enddo

348
c     For Testing
349
      call nek_testp()
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
      return
      end
C=======================================================================
      subroutine nek_diag()
      parameter (nl_max=100)
      parameter (nf_max=100)
      include 'SIZE'
      include 'TOTAL'
      common/expansion_tdata/n_legendre, m_fourier
      common/expansion_tcoef/coeff_tij(nl_max,nf_max)
      common/diags_coeff/diag_c(nl_max,nf_max)
      integer e,f

      real*8 fmode(lx1,ly1,lz1,lelt)
      real*8 cache(lx1,ly1,lz1,lelt)
      real*8 sint, sint1

      ntot=nx1*ny1*nz1*nelt

      zmin=glmin(zm1,ntot)
      zmax=glmax(zm1,ntot)

      call rzero(fmode,ntot)
         do i0=1,n_legendre
           do j0=1,m_fourier
             call nek_mode(fmode,i0,j0)  
             sint1=0.0
             do e=1,nelt
               do f=1,6
                 sint=0.0
                 if (cbc(f,e,1).eq.'W  ') then
                   call col3(cache,fmode,fmode,ntot)
                   call surface_int(sint,sarea,cache,e,f)
                   sint1=sint1+sint
                 endif
               enddo
             enddo
          call  gop(sint1,wtmp,'+  ',1)
          diag_c(i0,j0)=sint1*2.0/(0.5*(zmax-zmin))
          if (nid.eq.0) write(6,*)i0,j0,diag_c(i0,j0)
        enddo
      enddo

      return
      end
C=======================================================================
      subroutine nek_testp()
      parameter (nl_max=100)
      parameter (nf_max=100)
      include 'SIZE'
      include 'TOTAL'
      common/expansion_tdata/n_legendre, m_fourier
      common/expansion_tcoef/coeff_tij(nl_max,nf_max)
      integer e,f
      real*8 fmode(lx1,ly1,lz1,lelt), cache(lx1,ly1,lz1,lelt)
      real*8 fun(lx1,ly1,lz1,lelt) 
      ntot=nx1*ny1*nz1*nelt
      call rzero(fmode,ntot)
      call rzero(fun,ntot)   
         do i0=1,n_legendre
           do j0=1,m_fourier
             call nek_mode(fmode,i0,j0)
             do i=1,ntot
             fun(i,1,1,1)=fun(i,1,1,1)+fmode(i,1,1,1)*coeff_tij(i0,j0) 
             enddo  
          enddo
         enddo
      call rzero(cache,ntot) 
      call sub3(cache,fun,t,ntot)
419
420
421
422
      do i=1,ntot
        c1=cache(i,1,1,1)**2
        cache(i,1,1,1)=c1 
      enddo
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
 
      sint1=0.0
      sarea1=0.0
      do e=1,lelt
        do f=1,6
          call surface_int(sint,sarea,cache,e,f)
          if (cbc(f,e,1).eq.'W  ') then
           sint1=sint1+sint
           sarea1=sarea1+sarea
          endif
         enddo
      enddo
      call  gop(sint1,wtmp,'+  ',1)
      call  gop(sarea1,wtmp,'+  ',1)

438
      er_avg=sqrt(sint1/sarea1)
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
      if (nid.eq.0) write(6,*)"Error: ",er_avg
      return
      end
C=======================================================================
      subroutine nek_mode(fmode,im,jm)
      include 'SIZE'
      include 'TOTAL'
      integer e,f
      real*8 fmode(lx1,ly1,lz1,lelt) 
      ntot=nx1*ny1*nz1*nelt
      zmin=glmin(zm1,ntot)
      zmax=glmax(zm1,ntot)
             do e=1,nelt
               do f=1,6
                 if (cbc(f,e,1).eq.'W  ') then
                   call dsset(nx1,ny1,nz1)
                   iface  = eface1(f)
                   js1    = skpdat(1,iface)
                   jf1    = skpdat(2,iface)
                   jskip1 = skpdat(3,iface)
                   js2    = skpdat(4,iface)
                   jf2    = skpdat(5,iface)
                   jskip2 = skpdat(6,iface)
                   do j2=js2,jf2,jskip2
                     do j1=js1,jf1,jskip1
                       x=xm1(j1,j2,1,e)
                       y=ym1(j1,j2,1,e)
                       z=zm1(j1,j2,1,e)
                       z_leg=2*((z-zmin)/zmax)-1
                       theta=atan2(y,x)
                       fmode(j1,j2,1,e)=
     &                 pl_leg(z_leg,im-1)*fl_four(theta,jm-1)
                     enddo
                   enddo
                 endif
              enddo
            enddo
      return
      end 
c-----------------------------------------------------------------------
      subroutine flux_reconstruction()
      include 'SIZE'
      include 'TOTAL'
      parameter (nl_max=100)
      parameter (nf_max=100) 
      common/expansion_tdata/n_legendre, m_fourier
      common/expansion_tcoef/coeff_tij(nl_max,nf_max)
      common/expansion_fcoef/coeff_fij(nl_max,nf_max)
      common/expansion_recon/flux_recon(lx1,ly1,lz1,lelt)
      integer e,f
489
      real*8 coeff_base, flux_base, flux_moose 
490
491
492
      real*8 fmode(lx1,ly1,lz1,lelt)
      ntot=nx1*ny1*nz1*nelt

493
494
495
496
c     --------------------------------------
c     The flux from MOOSE must have proper sign
c     --------------------------------------
      coeff_base=-1.0 
497
      flux_base=0.25 ! from energy conservation in the problem
498
499
c     --------------------------------------

500
501
502
      call rzero(flux_recon,ntot)
         do i0=1,n_legendre
           do j0=1,m_fourier
503
             call rzero(fmode,ntot)
504
505
506
             call nek_mode(fmode,i0,j0)
             do i=1,ntot
             flux_recon(i,1,1,1)= flux_recon(i,1,1,1)
507
     &                  + coeff_base*fmode(i,1,1,1)*coeff_fij(i0,j0)
508
509
510
511
512
513
             enddo
          enddo
         enddo

c---- Below is for testing
c             do i=1,ntot
514
c             flux_recon(i,1,1,1)= 0.0
515
516
c             enddo
c--------------------------
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539

c---- Renormalization
      sint1=0.0
      sarea1=0.0

      do e=1,lelt
        do f=1,6
          call surface_int(sint,sarea,flux_recon,e,f)
          if (cbc(f,e,1).eq.'W  ') then
           sint1=sint1+sint
           sarea1=sarea1+sarea
          endif
         enddo
      enddo

      call  gop(sint1,wtmp,'+  ',1)
      call  gop(sarea1,wtmp,'+  ',1)

      flux_moose=sint1/sarea1
      do i=1,ntot
         flux_recon(i,1,1,1)= flux_recon(i,1,1,1)*flux_base/flux_moose
      enddo
c-----------------------
Cody Permann's avatar
Cody Permann committed
540
541
542
      return
      end
c-----------------------------------------------------------------------
543
544
! calculates Legendre polynomials using a recurrence relationship. If
! n > the maximum Legendre order, the function returns 0.0.
545
      function pl_leg(x,n)
546
      parameter (nl_max=100)
547
548
549
550
      real*8 pl,pl_leg
      real*8 x
      real*8 pln(0:n)
      integer n, k
551

552
553
      pln(0) = 1.0
      pln(1) = x
554

555
556
      if (n.le.1) then
        pl = pln(n)
557
558
559
560
561
562
563
      else if (n.le.nl_max) then
        do k=1,n-1
          pln(k+1) = ((2.0*k+1.0)*x*pln(k)-dble(k)*pln(k-1))/(dble(k+1))
        end do
        pl = pln(n)
      else
        pl = 0.0
564
      end if
565
566
567
568

      pl_leg = pl*sqrt(dble(2*n+1)/2.0)

      return
569
      end
570
571
c-----------------------------------------------------------------------
! calculates Fourier polynomials Fn(x)
572
573
574
575
      function fl_four(x,n)
      real*8 fl_four,pi
      real*8 x, A
      integer n, k
576

577
      pi=4.0*atan(1.0)
578
579
580
581
582
583
      A=1.0/sqrt(pi)

      if (n.eq.0) A=1.0/sqrt(2*pi)
      fl_four=A*cos(n*x)

      return
584
      end
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
c-----------------------------------------------------------------------
      subroutine heat_balance(fflux)
      include 'SIZE'
      include 'TOTAL'
      real sint, sint1, sarea, sarea1, wtmp, fflux
      real cache(lx1,ly1,lz1,lelt),dtdz(lx1,ly1,lz1,lelt),
     & dtdy(lx1,ly1,lz1,lelt),
     & dtdx(lx1,ly1,lz1,lelt)
      integer e, f

      n1=nelt*lx1*ly1*lz1
      n2=nelt*lx2*ly2*lz2
      nxyz=lx1*ly1*lz1
      call gradm1(dtdx,dtdy,dtdz,t(1,1,1,1,1))

      sint1=0.0
      sarea1=0.0
      do e=1,lelt
        do i=1,nxyz
          cache(i,1,1,e)=t(i,1,1,e,1)*vz(i,1,1,e)
        enddo
        do f=1,6
          call surface_int(sint,sarea,cache,e,f)
          if (cbc(f,e,1).eq.'O  ') then
           sint1=sint1+sint
           sarea1=sarea1+sarea
          endif
         enddo
      enddo
      call  gop(sint1,wtmp,'+  ',1)
      test_flux=sint1
      
      sint1=0.0
      sarea1=0.0
      do e=1,lelt
        do f=1,6
          call surface_int(sint,sarea,dtdz,e,f)
          if (cbc(f,e,1).eq.'v  ') then
           sint1=sint1+sint
           sarea1=sarea1+sarea
          endif
         enddo
      enddo
      call  gop(sint1,wtmp,'+  ',1)
      call  gop(sarea1,wtmp,'+  ',1)
      flux_inlet=sint1
631

632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
      sint1=0.0
      sarea1=0.0
      do e=1,lelt
        do f=1,6
          call surface_int(sint,sarea,vz,e,f)
          if (cbc(f,e,1).eq.'W  ') then
           sarea1=sarea1+sarea
          endif
         enddo
      enddo
      call  gop(sarea1,wtmp,'+  ',1)

      if (nid.eq.0) then
         rr_loss=
     &   flux_inlet*param(8)/(test_flux+flux_inlet*param(8))
         write(6,*)"*** Heat balance: ",
     &   sarea1*fflux,test_flux+flux_inlet*param(8)
         write(6,*)"*** Heat loss: ",rr_loss*100.0," %"
      endif

      return
      end
c-----------------------------------------------------------------------