-
Notifications
You must be signed in to change notification settings - Fork 12
/
m_les.f90
2396 lines (1953 loc) · 86.5 KB
/
m_les.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
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
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
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
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
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
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
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
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
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
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!================================================================================
! m_les - the module that contains all subroutines and variables that are
! needed to introduce Large Eddy Simulation (LES) into the code.
!
! The behaviour of the module is governed by the variable "les_mode" from the
! module m_parameters.f90
!
! Time-stamp: <2010-03-18 13:52:33 (chumakov)>
!================================================================================
module m_les
use m_openmpi
use m_io
use m_parameters
use m_fields
use m_work
use x_fftw
implicit none
!================================================================================
! Arrays and constants
!================================================================================
! indicator, whether we use LES at all
logical :: les = .false.
! model switch "les_mode" is contained by m_parameters
! integer(kind=4) :: les_model
! array for turbulent viscosity
real(kind=8), allocatable :: turb_visc(:,:,:)
! LES sources for velocities
real(kind=8), allocatable :: vel_source_les(:,:,:,:)
! separate array for the subgrid stress
real(kind=8), allocatable :: tauij(:,:,:,:)
! Smagorinsky constant
real(kind=8) :: c_smag = 0.18
! Scaling constant for the lag-model
real(kind=8) :: C_T = 1.d0
! Scaling constant for the mixed model
real(kind=8) :: C_mixed
! test filter width
real(kind=8) :: les_delta
! array with the test filter in it
real(kind=8), allocatable :: filter_g(:,:,:)
! model indicator for the output
character*3 :: les_model_name = ' '
! TEMP variables:
! - production
real(kind=8) :: energy, production, B, dissipation
!================================================================================
! SUBROUTINES
!================================================================================
contains
!================================================================================
! Allocation of LES arrays
!================================================================================
subroutine m_les_init
implicit none
integer :: n
real*8, allocatable :: sctmp(:)
! if les_model=0, do not initialize anything and return
if (les_model==0) return
write(out,*) 'Initializing LES...'
call flush(out)
! depending on the value of les_model, initialize different things
les = .true.
! initialize the filter width to be equal to the grid spaxing
les_delta = dx
write(out,*) 'LES_DELTA = ',les_delta
! initializeing stuff based on the model switch "les_model"
select case (les_model)
!--------------------------------------------------------------------------------
! Smagorinsky model
!--------------------------------------------------------------------------------
case(1)
write(out,*) ' - Smagorinsky model: initializing the eddy viscosity'
call flush(out)
allocate(turb_visc(nx,ny,nz),stat=ierr)
if (ierr/=0) then
write(out,*) 'Cannot allocate the turbulent viscosity array'
call flush(out)
call my_exit(-1)
end if
turb_visc = 0.0d0
n_les = 0
les_model_name = " SM"
!--------------------------------------------------------------------------------
! Dynamic localization model
!--------------------------------------------------------------------------------
case(2)
write(out,*) ' - DL model: initializing the eddy viscosity and adding extra transport equation'
call flush(out)
allocate(turb_visc(nx,ny,nz),stat=ierr)
if (ierr/=0) then
write(out,*) 'Cannot allocate the turbulent viscosity array'
call flush(out)
call my_exit(-1)
end if
turb_visc = 0.0d0
n_les = 1
les_model_name = "DLM"
!--------------------------------------------------------------------------------
! Dynamic localization model + lag model for the subgrid energy dissipation
!--------------------------------------------------------------------------------
case(3)
! Dynamic Localization model + lag model for the dissipation
write(out,*) ' - DL model + lag-model for dissipation'
call flush(out)
allocate(turb_visc(nx,ny,nz),stat=ierr)
if (ierr/=0) then
write(out,*) 'Cannot allocate the turbulent viscosity array'
call flush(out)
call my_exit(-1)
end if
turb_visc = 0.0d0
n_les = 3
les_model_name = "DLL"
!--------------------------------------------------------------------------------
! Dynamic Structure model + algebraic model for dissipation
!--------------------------------------------------------------------------------
case(4)
write(out,*) ' - DSt model + algebraic model for dissipation'
call flush(out)
allocate(turb_visc(nx,ny,nz), vel_source_les(nx+2,ny,nz,3), stat=ierr)
if (ierr/=0) then
write(out,*) 'Cannot allocate the turbulent viscosity array'
call flush(out)
call my_exit(-1)
end if
turb_visc = 0.0d0
n_les = 1
les_model_name = "DST"
! Fot this model we need a filter. The filter cannot be initialized
! without previous initialization of fields array. So initialization of
! the filter is done in m_les_begin
! Note that for this model we need a bigger wrk array
! this is taken care of in m_work.f90
!--------------------------------------------------------------------------------
! Dynamic Structure model + lag model for dissipation
!--------------------------------------------------------------------------------
case(5)
write(out,*) ' - DSt model + lag model for dissipation'
call flush(out)
allocate(turb_visc(nx,ny,nz), vel_source_les(nx+2,ny,nz,3), stat=ierr)
if (ierr/=0) then
write(out,*) 'Cannot allocate the turbulent viscosity array'
call flush(out)
call my_exit(-1)
end if
turb_visc = 0.0d0
n_les = 3
les_model_name = "STL"
! Fot this model we need a filter. The filter cannot be initialized
! without previous initialization of fields array. So initialization of
! the filter is done in m_les_begin
! Note that for this model we need a bigger wrk array
! this is taken care of in m_work.f90
!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> DEBUG+
inquire(file = 'c_t.in', exist = there)
if (.not.there) then
write(out,*) "Cannot find the file 'c_t.in', exiting"
call my_exit(-1)
end if
if (iammaster) then
open(900,file='c_t.in')
read(900,*) C_T
close(900)
end if
count = 1
call MPI_BCAST(C_T,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
write(out,*) "DSTM + lag model WITH C_T = ",C_T
call flush(out)
!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> DEBUG-
!--------------------------------------------------------------------------------
! Mixed model: Dynamic Structure + C-mixed * Dynamic Localization model
! Dissipation model is algebraic
!--------------------------------------------------------------------------------
case(6)
write(out,*) ' - MIXED MODEL: DSTM + DLM'
write(out,*) ' - Dissipation is algebraic'
call flush(out)
!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> DEBUG+
inquire(file = 'c_mixed.in', exist = there)
if (.not.there) then
write(out,*) "Cannot find the file 'c_mixed.in', exiting"
call my_exit(-1)
end if
if (iammaster) then
open(900,file='c_mixed.in')
read(900,*) C_mixed
close(900)
end if
count = 1
call MPI_BCAST(C_mixed,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
write(out,*) "MIXED MODEL WITH C_MIXED = ",C_mixed
call flush(out)
!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> DEBUG-
allocate(turb_visc(nx,ny,nz), vel_source_les(nx+2,ny,nz,3), stat=ierr)
if (ierr/=0) then
write(out,*) 'Cannot allocate the turbulent viscosity array'
call flush(out)
call my_exit(-1)
end if
turb_visc = 0.0d0
n_les = 1
les_model_name = "MMA"
! Fot this model we need a filter. The filter cannot be initialized
! without previous initialization of fields array. So initialization of
! the filter is done in m_les_begin
! Note that for this model we need a bigger wrk array
! this is taken care of in m_work.f90
!--------------------------------------------------------------------------------
! Mixed model: Dynamic Structure + some viscosity (about 15% of the usual)
! Dissipation model is lag-model
!--------------------------------------------------------------------------------
case(7)
write(out,*) ' - MIXED MODEL: DSTM + eddy viscosity'
write(out,*) ' - Lag-model for Dissipation'
call flush(out)
!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> DEBUG+
inquire(file = 'c_mixed.in', exist = there)
if (.not.there) then
write(out,*) "Cannot find the file 'c_mixed.in', exiting"
call my_exit(-1)
end if
if (iammaster) then
open(900,file='c_mixed.in')
read(900,*) C_mixed
close(900)
end if
count = 1
call MPI_BCAST(C_mixed,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
write(out,*) "MIXED MODEL WITH C_MIXED = ",C_mixed
call flush(out)
!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> DEBUG-
allocate(turb_visc(nx,ny,nz), vel_source_les(nx+2,ny,nz,3), stat=ierr)
if (ierr/=0) then
write(out,*) 'Cannot allocate the turbulent viscosity array'
call flush(out)
call my_exit(-1)
end if
turb_visc = 0.0d0
n_les = 3
les_model_name = "MML"
! Fot this model we need a filter. The filter cannot be initialized
! without previous initialization of fields array. So initialization of
! the filter is done in m_les_begin
! Note that for this model we need a bigger wrk array
! this is taken care of in m_work.f90
!--------------------------------------------------------------------------------
! LES MODEL # 10
!
! One-equation mixed model for the momentum closure (DSTM + C_mixed * SM)
! Lag model for the closure of the k-equation
! Harlow model for closure of the scalar tranport equation
!--------------------------------------------------------------------------------
case(10)
write(out,*) ' LES MODEL # 10'
write(out,*) ' - MIXED MODEL: DSTM + eddy viscosity'
write(out,*) ' Lag-model for Dissipation'
write(out,*) ' - SCALARS: Harlow model'
write(out,*) ' ============> file c_mixed.in is still required'
call flush(out)
!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> DEBUG+
inquire(file = 'c_mixed.in', exist = there)
if (.not.there) then
write(out,*) "Cannot find the file 'c_mixed.in', making default 0.5"
c_mixed = 0.5d0
else
if (iammaster) then
open(900,file='c_mixed.in')
read(900,*) C_mixed
close(900)
end if
count = 1
call MPI_BCAST(C_mixed,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
end if
write(out,*) "MIXED MODEL WITH C_MIXED = ",C_mixed
call flush(out)
!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> DEBUG-
allocate(turb_visc(nx,ny,nz), tauij(nx+2,ny,nz,6), stat=ierr)
if (ierr/=0) then
write(out,*) 'Cannot allocate the auxiliary arrays'
call flush(out)
call my_exit(-1)
end if
turb_visc = zip
tauij = zip
! number of extra LES quantities is 2: (BT), (eps * T)
n_les = 2
les_model_name = "#10"
! For this model we need a filter. The filter cannot be initialized
! without previous initialization of fields array. So initialization of
! the filter is done in m_les_begin
! Note that for this model we need a bigger wrk array
! this is taken care of in m_work.f90
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
case default
write(out,*) 'M_LES_INIT: invalid value of les_model:',les_model
call flush(out)
call my_exit(-1)
end select
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
write(out,*) "n_les = ", n_les
! If the number of LES quantities is non-zero, then for the sake of modularity
! we need to have SC and PE numbers defined for those quantities. For this we
! need to re-allocate the arrays SC and PE
additional_scalars: if (n_les .gt. 0) then
write(out,*) "Adding elements to arrays PE and SC for the LES-related scalars..."
call flush(out)
allocate(sctmp(1:n_scalars+n_les), stat=ierr)
sctmp(1:n_scalars) = sc(1:n_scalars)
sctmp(n_scalars+1:n_scalars+n_les) = one
if (allocated(sc)) deallocate(sc)
allocate(sc(1:n_scalars+n_les))
sc = sctmp
if (allocated(pe)) deallocate(pe)
allocate(pe(1:n_scalars+n_les))
pe = nu / sc
deallocate(sctmp)
write(out,*) " ...done."
call flush(out)
end if additional_scalars
write(out,*) 'initialized LES.'
call flush(out)
return
end subroutine m_les_init
!================================================================================
!================================================================================
! initialization of the LES arrays - part 2
! definition of the arrays
!
! called after the restart
!================================================================================
subroutine m_les_begin
use m_filter_xfftw
implicit none
write(out,*) "M_LES_BEGIN..."
call flush(out)
select case(les_model)
!--------------------------------------------------------------------------------
! if les_model=0, do not initialize anything and return
! also don't do anything if the model is Smagorinsky model - it does not
! need any additional initialization beside the array allocation which has
! been already done.
!--------------------------------------------------------------------------------
case(0)
return
case(1)
return
!--------------------------------------------------------------------------------
! Dynamic localization model
!--------------------------------------------------------------------------------
case(2)
write(out,*) "-- DLM model"
call flush(out)
if (itime.eq.0) then
write(out,*) "-- Initializing k_sgs"
call flush(out)
call m_les_dlm_k_init
end if
!--------------------------------------------------------------------------------
case(3)
write(out,*) "-- DLM model with lag model for dissipation"
call flush(out)
if (itime.gt.0) return
! call m_les_dlm_k_init
! COMMENT OUT THE FOLLOWING IN ORDER TO INITIALIZE B AND EPSILON AS ZERO
! making initial epsilon = k^(3/2)/Delta everywhere
! since k=const=0.5, just change one entry in epsilon array
! Note that the array itself contains not epsilon but (epsilon * T_epsilon)
! so the contents of the array is not presicely k^(3/2)/Delta. Some math is involved.
! (comment out to start from zero dissipation)
! if (iammaster) fields(1,1,1,3+n_scalars+3) = C_T * 0.5 * real(nxyz_all)
! Now initial conditions for B. We want B to be same as epsilon
! (kind of "starting from steady state"), but again the array contains B*T_B, not just B.
! Current implementation is T_B = 1/|S|. To get |S|, we call m_les_src_k_dlm, which
! gives us |S|^2 in wrk0.
! call m_les_k_src_dlm
! now wrk0 contains |S|^2 in x-space, and we can use it to get B
! fields(:,:,:,3+n_scalars+2) = 0.5**1.5d0 / les_delta / sqrt(wrk(:,:,:,0))
! call xFFT3d_fields(1,3+n_scalars+2)
! INITIALIZING K_SGS, B*T_B and eps*T_eps
! Initializing them so that k_sgs = B*T_b + eps*T_eps
! in fact this makes the equation for k_sgs unnecessary but we're keeping it for
! debug purposes and such
write(out,*) " Initializing k_sgs = 0.1, B*T_B = eps*T_eps = 0.05"
call flush(out)
! definition of k_sgs = 0.1 everywhere
if (iammaster) fields(1,1,1,3+n_scalars+1) = 0.1d0 * real(nxyz_all)
! definition of eps*T_eps = B*T_B = 0.5 k_sgs
if (iammaster) fields(1,1,1,3+n_scalars+2) = 0.5d0 * fields(1,1,1,3+n_scalars+1)
if (iammaster) fields(1,1,1,3+n_scalars+3) = 0.5d0 * fields(1,1,1,3+n_scalars+1)
!--------------------------------------------------------------------------------
case(4)
write(out,*) "-- Dynamic Structure model with algebraic model for dissipation"
call flush(out)
! initializing filtering arrays
! because filter_xfftw_init uses fields(1) as a temporary array, we need
! to store it before we initialize the filter, and the restore it to
! what it was.
wrk(:,:,:,LBOUND(wrk,4)) = fields(:,:,:,LBOUND(fields,4))
call filter_xfftw_init
fields(:,:,:,LBOUND(fields,4)) = wrk(:,:,:,LBOUND(wrk,4))
if (itime.eq.0) then
write(out,*) "-- Initializing k_sgs = 0.2"
call flush(out)
if (iammaster) fields(1,1,1,3+n_scalars+1) = 0.2d0 * real(nxyz_all)
end if
!--------------------------------------------------------------------------------
case(5)
write(out,*) "-- Dynamic Structure model with lag-model for dissipation"
call flush(out)
! initializing filtering arrays
! because filter_xfftw_init uses fields(1) as a temporary array, we need
! to store it before we initialize the filter, and the restore it to
! what it was.
wrk(:,:,:,LBOUND(wrk,4)) = fields(:,:,:,LBOUND(fields,4))
call filter_xfftw_init
fields(:,:,:,LBOUND(fields,4)) = wrk(:,:,:,LBOUND(wrk,4))
if (itime.eq.0) then
write(out,*) "-- Initializing k_sgs = 0.5"
call flush(out)
if (iammaster) fields(1,1,1,3+n_scalars+1) = 0.5d0 * real(nxyz_all)
! definition of eps*T_eps = 0, B*T_B = k_sgs
if (iammaster) fields(1,1,1,3+n_scalars+2) = fields(1,1,1,3+n_scalars+1)
if (iammaster) fields(1,1,1,3+n_scalars+3) = 0.d0*fields(1,1,1,3+n_scalars+1)
end if
!--------------------------------------------------------------------------------
case(6)
write(out,*) "-- MIXED MODEL (Dynamic Structure model + Smagorinsky)"
write(out,*) " Algebraic model for dissipation"
call flush(out)
! initializing filtering arrays
! because filter_xfftw_init uses fields(1) as a temporary array, we need
! to store it before we initialize the filter, and the restore it to
! what it was.
write(out,*) "Initializing filter"
call flush(out)
wrk(:,:,:,LBOUND(wrk,4)) = fields(:,:,:,LBOUND(fields,4))
call filter_xfftw_init
fields(:,:,:,LBOUND(fields,4)) = wrk(:,:,:,LBOUND(wrk,4))
if (itime.eq.0) then
write(out,*) "-- Initializing k_sgs = 0.5"
call flush(out)
if (iammaster) fields(1,1,1,3+n_scalars+1) = 0.5d0 * real(nxyz_all)
end if
!--------------------------------------------------------------------------------
case(7)
write(out,*) "-- MIXED MODEL (Dynamic Structure model + Smagorinsky)"
write(out,*) " Lag-model for dissipation"
call flush(out)
! initializing filtering arrays
! because filter_xfftw_init uses fields(1) as a temporary array, we need
! to store it before we initialize the filter, and the restore it to
! what it was.
write(out,*) " Initializing filter"
call flush(out)
wrk(:,:,:,LBOUND(wrk,4)) = fields(:,:,:,LBOUND(fields,4))
call filter_xfftw_init
fields(:,:,:,LBOUND(fields,4)) = wrk(:,:,:,LBOUND(wrk,4))
if (itime.eq.0) then
write(out,*) "-- Initializing k_sgs = 0.5, (BT)=k_s, (Eps T) = 0"
call flush(out)
! Initializing k
if (iammaster) fields(1,1,1,3+n_scalars+1) = 0.5d0 * real(nxyz_all)
! initializing (BT)
if (iammaster) fields(1,1,1,3+n_scalars+2) = fields(1,1,1,3+n_scalars+1)
! initializing (eps T)
if (iammaster) fields(1,1,1,3+n_scalars+3) = 0.d0*fields(1,1,1,3+n_scalars+1)
end if
!--------------------------------------------------------------------------------
case(10)
write(out,*) "-- MIXED MODEL (Dynamic Structure model + Smagorinsky)"
write(out,*) " Lag-model for dissipation"
call flush(out)
! initializing filtering arrays
! because filter_xfftw_init uses fields(1) as a temporary array, we need
! to store it before we initialize the filter, and the restore it to
! what it was.
write(out,*) " Initializing filter"
call flush(out)
wrk(:,:,:,LBOUND(wrk,4)) = fields(:,:,:,LBOUND(fields,4))
call filter_xfftw_init
fields(:,:,:,LBOUND(fields,4)) = wrk(:,:,:,LBOUND(wrk,4))
if (itime.eq.0) then
write(out,*) "-- (BT) = 0.1, (Eps T) = 0.1"
call flush(out)
! No need to initialize k_s, since k_s = (BT)+(eps T)
! if (iammaster) fields(1,1,1,3+n_scalars+1) = 0.5d0 * real(nxyz_all)
! initializing (BT)
if (iammaster) fields(1,1,1,3+n_scalars+1) = 0.1d0 * real(nxyz_all)
! initializing (eps T)
if (iammaster) fields(1,1,1,3+n_scalars+2) = 0.1d0 * real(nxyz_all)
end if
!--------------------------------------------------------------------------------
case default
write(out,*) "M_LES_BEGIN: invalid value of les_model: ",les_model
call flush(out)
call my_exit(-1)
end select
return
end subroutine m_les_begin
!================================================================================
!================================================================================
! Adding LES sources to the RHS of velocities
!================================================================================
subroutine les_rhs_velocity
use x_fftw
implicit none
integer :: i, j, k
select case (les_model)
case(1:3)
call les_rhsv_turb_visc
case(4:5)
! Dynamic Structure model
! add the velocity sources to RHS for velocitieis
wrk(:,:,:,1:3) = wrk(:,:,:,1:3) + vel_source_les(:,:,:,1:3)
case(6)
! Mixed model (DSTM + DLM)
! add the velocity sources to RHS for velocitieis
wrk(:,:,:,1:3) = wrk(:,:,:,1:3) + vel_source_les(:,:,:,1:3)
! also apply turbulent viscosity
call les_rhsv_turb_visc
case(7)
! Mixed model (DSTM + DLM) + lag-model for dissipation of k_sgs
! add the velocity sources to RHS for velocitieis
wrk(:,:,:,1:3) = wrk(:,:,:,1:3) + vel_source_les(:,:,:,1:3)
! also apply turbulent viscosity to velocities
call les_rhsv_turb_visc
case(10)
! Mixed model (DSTM + DLM) + lag-model for dissipation of k_sgs
! add the velocity sources to RHS for velocitieis
! The sources are directly computed from the array tauij
! so we do not need a separate array for them
call les_add_vel_source_from_tauij
case default
write(out,*) 'LES_RHS_VELOCITY: invalid value of les_model:',les_model
call flush(out)
call my_exit(-1)
end select
!--------------------------------------------------------------------------------
! Making sure that we are not getting any aliasing errors. That is done by
! making RHS zero for wavenumbers that are aliased.
!--------------------------------------------------------------------------------
do k = 1, nz
do j = 1, ny
do i = 1, nx + 2
if (ialias(i,j,k) .ne. 0) then
wrk(i,j,k,1:3) = zip
end if
end do
end do
end do
!--------------------------------------------------------------------------------
return
end subroutine les_rhs_velocity
!================================================================================
!================================================================================
! Adding LES sources to the RHS of scalars
!================================================================================
subroutine les_rhs_scalars
use x_fftw
implicit none
integer :: n, i, j, k
! note that the turbulent viscosity itself is computed in rhs_scalars.f90
! here we only modify the RHSs for scalars in case we're running LES
select case (les_model)
!--------------------------------------------------------------------------------
case(1)
call m_les_rhss_turb_visc
!--------------------------------------------------------------------------------
case(2)
! -- Dynamic Localization model with algebraic model for dissipation
call m_les_k_src_dlm
call m_les_k_diss_algebraic
call m_les_rhss_turb_visc
!--------------------------------------------------------------------------------
case(3)
! -- Dynamic Localization model with lag-model for dissipation
call m_les_k_src_dlm
call m_les_lag_model_sources
call m_les_rhss_turb_visc
!--------------------------------------------------------------------------------
case(4)
! -- Dynamic Structure Model with algebraic model for dissipation
! First taking care of the passive scalars (don't have them for now)
if (n_scalars .gt. 0) then
write(out,*) "*** Current version of the code cannot transport scalars"
write(out,*) "*** with the les_model=4 (Dynamic Structure Model)"
write(out,*) "Please specify a different LES model."
call flush(out)
call my_exit(-1)
end if
! now taking care of the LES-related scalars
! for the Dynamic Structure model, k-equation has turbulent viscosity
! so need to add that term.
call m_les_dstm_vel_k_sources
call m_les_k_diss_algebraic
call m_les_rhss_turb_visc
!--------------------------------------------------------------------------------
case(5)
! -- Dynamic Structure Model with lag-model for dissipation
! First taking care of the passive scalars (don't have them for now)
if (n_scalars .gt. 0) then
write(out,*) "*** Current version of the code cannot transport scalars"
write(out,*) "*** with the les_model=5 (Dynamic Structure Model)"
write(out,*) "Please specify a different LES model."
call flush(out)
call my_exit(-1)
end if
! now taking care of the LES-related scalars
! for the Dynamic Structure model, k-equation has turbulent viscosity
! so need to add that term.
! getting sources for velocities and production for k_s and B
call m_les_dstm_vel_k_sources
! getting |S|^2 and putting it into wrk0 (for the timescale for B)
call m_les_k_src_dlm
! getting the sources and sinks for (BT) and (eps T) and a sink for k_s
call m_les_lag_model_sources
! diffusing k_s, (BT) and (epsilon*T) with turbulent viscosity
call m_les_rhss_turb_visc
!--------------------------------------------------------------------------------
case(6)
! Mixed model (Dynamic Structure model + a fraction of Dynamic Localization model)
! The fraction is given by the constant C_mixed, which is read from the file
! This is taken care of in les_get_turb_visc
! First taking care of the passive scalars (don't have them for now)
if (n_scalars .gt. 0) then
write(out,*) "*** Current version of the code cannot transport scalars"
write(out,*) "*** with the les_model=6 (Mixed Model)"
write(out,*) "Please specify a different LES model."
call flush(out)
call my_exit(-1)
end if
! now taking care of the LES-related scalars
! for the mixed model, scalars and velocities have turbulent viscosity
! so need to add that term.
! getting DSTM sources for velocities and production for k_sgs
call m_les_dstm_vel_k_sources
! getting the DLM transfer term turb_visc*|S|^2 and adding to the RHS for k_sgs
call m_les_k_src_dlm
! diffusing k_sgs with turbulent viscosity
call m_les_rhss_turb_visc
! algebraic model for dissipation of k_sgs: k^{3/2}/Delta
call m_les_k_diss_algebraic
!--------------------------------------------------------------------------------
case(7)
! Mixed model (Dynamic Structure model + a fraction of Dynamic Localization model)
! The fraction is given by the constant C_mixed, which is read from the file
! This is taken care of in les_get_turb_visc
! Dissipation is via lag model
! First taking care of the passive scalars (don't have them for now)
if (n_scalars .gt. 0) then
write(out,*) "*** Current version of the code cannot transport scalars"
write(out,*) "*** with the les_model=6 (Mixed Model)"
write(out,*) "Please specify a different LES model."
call flush(out)
call my_exit(-1)
end if
! now taking care of the LES-related scalars
! for the Dynamic Structure model, k-equation has turbulent viscosity
! so need to add that term.
! getting sources for velocities and production for k_s and B
call m_les_dstm_vel_k_sources
! getting |S|^2 and putting it into wrk0 (for the timescale for B)
call m_les_k_src_dlm
! getting the sources and sinks for (BT) and (eps T) and a sink for k_s
call m_les_lag_model_sources
! diffusing k_s, (BT) and (epsilon*T) with turbulent viscosity
call m_les_rhss_turb_visc
!--------------------------------------------------------------------------------
case(10)
! Mixed model (Dynamic Structure model + a fraction of Smagorinsky model)
! The fraction is given by the constant C_mixed, which is read from the file
! c_mixed.in
! Dissipation term in the k-equation is closed via lag-model
! first we add the turbulent diffusion to (BT) and (eps T)
! (no turbulent diffusion for the passive scalars, since they are taken
! care of by the Harlow model later)
call m_les_rhss_turb_visc1(3+n_scalars+1)
call m_les_rhss_turb_visc1(3+n_scalars+2)
! now calculate the SGS stress tauij and calculate the source term Pi
! for (BT), the scalar with the number 3+n_scalars+1
! Add the source term for (BT) to the RHS for (BT)
! Store the SGS stress in the array tauij.
call m_les_get_tauij_dstm_source_BT
! get the SGS fluxes for all passive scalars, using Harlow model
if (n_scalars.gt.0) call m_les_rhss_sgs_flux_harlow
! get the lag-model sources for (BT) and (eps T) and add those to RHSs
! the "no_k" means that we calculate k_sgs as (BT)+(eps T), effectively
! reducing the number of extra equations.
call m_les_lag_model_sources_no_k
!--------------------------------------------------------------------------------
case default
write(out,*) 'LES_RHS_SCALARS: invalid value of les_model:',les_model
call flush(out)
call my_exit(-1)
end select
!--------------------------------------------------------------------------------
! Making sure that we are not getting any aliasing errors. That is done by
! making RHS zero for wavenumbers that produce aliasing.
!--------------------------------------------------------------------------------
do k = 1, nz
do j = 1, ny
do i = 1, nx + 2
if (ialias(i,j,k) .ne. 0) then
wrk(i,j,k,4:3+n_scalars+n_les) = zip
end if
end do
end do
end do
!--------------------------------------------------------------------------------
! outputting the LES quantities in the file les.gp
!--------------------------------------------------------------------------------
if(iammaster) then
if (mod(itime,iprint1)==0) then
open(999,file='les.gp', position='append')
if (n_les>0) then
if (les_model.lt.10) then
energy = fields(1,1,1,3+n_scalars+1) / real(nxyz_all)
else
energy = (fields(1,1,1,3+n_scalars+1) + fields(1,1,1,3+n_scalars+2))/ real(nxyz_all)
end if
write(999,"(i6,x,10e15.6)") itime, time, energy, production, B, dissipation
close(999)
end if
end if
B = zip
production = zip
dissipation = zip
end if
!--------------------------------------------------------------------------------
return
end subroutine les_rhs_scalars
!================================================================================
!================================================================================
! calculating LES sources for velocity field and adding them to the
! RHS's (wrk1...3)
!================================================================================
subroutine les_rhsv_turb_visc
use x_fftw
implicit none
integer :: i, j, k, n
real(kind=8) :: rtmp
! due to memory constraints we have only three work arrays wrk4..6,
! because the first three wrk1..3 contain already comptued velocity RHS's.
! Calculating S_11, S_12, S_13
do k = 1, nz
do j = 1, ny
do i = 1, nx + 1, 2
! S_11, du/dx
wrk(i ,j,k,4) = - akx(i+1) * fields(i+1,j,k,1)
wrk(i+1,j,k,4) = akx(i ) * fields(i ,j,k,1)
! S_12, 0.5 (du/dy + dv/dx)
wrk(i ,j,k,5) = -half * ( aky(k) * fields(i+1,j,k,1) + akx(i+1) * fields(i+1,j,k,2) )
wrk(i+1,j,k,5) = half * ( aky(k) * fields(i ,j,k,1) + akx(i ) * fields(i ,j,k,2) )
! S_13, 0.5 (du/dz + dw/dx)
wrk(i ,j,k,6) = -half * ( akz(j) * fields(i+1,j,k,1) + akx(i+1) * fields(i+1,j,k,3) )
wrk(i+1,j,k,6) = half * ( akz(j) * fields(i ,j,k,1) + akx(i ) * fields(i ,j,k,3) )
end do
end do
end do
! Multiplying them by -2 * turbulent viscosity to get tau_11, tau_12, tau_13
do n = 4,6
call xFFT3d(-1,n)
wrk(1:nx,1:ny,1:nz,n) = - two * turb_visc(1:nx,1:ny,1:nz) * wrk(1:nx,1:ny,1:nz,n)
call xFFT3d(1,n)
end do
! Taking d/dx tau_11, d/dy tau_12, d/dz tau_13 and subtracting from the current RHS
! note the sign reversal (-/+) because we subtract this from RHS
do k = 1, nz
do j = 1, ny
do i = 1, nx+1, 2
! Cutting off any wave modes that can introduce aliasing into the velocities
! This "dealiasing" is done in the most restrictive manner: only Fourier modes
! that have ialias=0 are added
if (ialias(i,j,k) .eq. 0) then
wrk(i ,j,k,1) = wrk(i ,j,k,1) + akx(i+1)*wrk(i+1,j,k,4) + aky(k)*wrk(i+1,j,k,5) + akz(j)*wrk(i+1,j,k,6)
wrk(i+1,j,k,1) = wrk(i+1,j,k,1) - akx(i )*wrk(i ,j,k,4) - aky(k)*wrk(i ,j,k,5) - akz(j)*wrk(i ,j,k,6)
wrk(i ,j,k,2) = wrk(i ,j,k,2) + akx(i+1)*wrk(i+1,j,k,5)
wrk(i+1,j,k,2) = wrk(i+1,j,k,2) - akx(i )*wrk(i ,j,k,5)
wrk(i ,j,k,3) = wrk(i ,j,k,3) + akx(i+1)*wrk(i+1,j,k,6)
wrk(i+1,j,k,3) = wrk(i+1,j,k,3) - akx(i )*wrk(i ,j,k,6)
end if
end do
end do
end do
! Calculating S_22, S_23, S_33
do k = 1, nz
do j = 1, ny
do i = 1, nx+1, 2
! S_22, dv/dy
wrk(i ,j,k,4) = - aky(k) * fields(i+1,j,k,2)
wrk(i+1,j,k,4) = aky(k) * fields(i ,j,k,2)
! S_23, 0.5 (dv/dz + dw/dy)
wrk(i ,j,k,5) = - half * ( akz(j) * fields(i+1,j,k,2) + aky(k) * fields(i+1,j,k,3) )
wrk(i+1,j,k,5) = half * ( akz(j) * fields(i ,j,k,2) + aky(k) * fields(i ,j,k,3) )
! S_33, de/dz
wrk(i ,j,k,6) = - akz(j) * fields(i+1,j,k,3)
wrk(i+1,j,k,6) = akz(j) * fields(i ,j,k,3)
end do
end do
end do
! Multiplying them by -2 * turbulent viscosity to get tau_22, tau_23, tau_33
do n = 4,6
call xFFT3d(-1,n)
wrk(1:nx,1:ny,1:nz,n) = - two * turb_visc(1:nx,1:ny,1:nz) * wrk(1:nx,1:ny,1:nz,n)
call xFFT3d(1,n)
end do
! Taking
! d/dy tau_22, d/dz tau_23
! d/dy tau_23, d/dz tau_33
! and subtracting from the current RHS for v and w
! note the sign reversal (-/+) because we subtract this from RHS
do k = 1, nz
do j = 1, ny
do i = 1, nx+1, 2