forked from geoflows/dclaw
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathamr2.f90
853 lines (706 loc) · 29.4 KB
/
amr2.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
!
! Use adaptive mesh refinement to solve the hyperbolic 2-d equation:
!
! u + f(u) + g(u) = 0
! t x y
!
! or the more general non-conservation law form:
! u + A u + B u = 0
! t x y
!
! using the wave propagation method as in CLAWPACK in combination
! with the locally uniform embedded grids of AMR.
!
! Estimate error with Richardson extrap. (in errest.f)
! + gradient checking (in errsp.f). Initial conditions set
! in (qinit.f), b.c.'s in (physbd.f).
!
! Specify rectangular domain from
! (xlower,ylower) to (xupper,yupper).
!
!
! =========================================================================
!
! This software is made available for research and instructional use only.
! You may copy and use this software without charge for these non-commercial
! purposes, provided that the copyright notice and associated text is
! reproduced on all copies. For all other uses (including distribution of
! modified versions), please contact the author at the address given below.
!
! *** This software is made available "as is" without any assurance that it
! *** will work for your purposes. The software may in fact have defects, so
! *** use the software at your own risk.
!
! --------------------------------------
! AMRCLAW Version 5.0, 2012
! Homepage: http://www.clawpack.org
! --------------------------------------
!
! Authors:
!
! Marsha J. Berger
! Courant Institute of Mathematical Sciences
! New York University
! 251 Mercer St.
! New York, NY 10012
!
! Randall J. LeVeque
! Applied Mathematics
! Box 352420
! University of Washington,
! Seattle, WA 98195-2420
!
! =========================================================================
program amr2
use amr_module, only: dbugunit, parmunit, outunit, inunit, matunit
use amr_module, only: mxnest, rinfinity, iinfinity
use amr_module, only: xupper, yupper, xlower, ylower
use amr_module, only: hxposs, hyposs, intratx, intraty, kratio
use amr_module, only: cfl, cflv1, cflmax, evol
use amr_module, only: checkpt_style, checkpt_interval, tchk, nchkpt
use amr_module, only: rstfile
use amr_module, only: max1d, maxvar, maxlv
use amr_module, only: method, mthlim, use_fwaves, numgrids
use amr_module, only: nghost, mwaves, mcapa, auxtype
use amr_module, only: tol, tolsp, flag_richardson, flag_gradient
use amr_module, only: nghost, mthbc
use amr_module, only: xperdom, yperdom, spheredom
use amr_module, only: nstop, nout, iout, tfinal, tout, output_style
use amr_module, only: output_format, printout, verbosity_regrid
use amr_module, only: output_q_components, output_aux_components
use amr_module, only: output_aux_onlyonce, matlabu
use amr_module, only: lfine, lentot, iregridcount, avenumgrids
use amr_module, only: tvoll, tvollCPU, rvoll, rvol, mstart, possk, ibuff
use amr_module, only: timeRegridding,timeUpdating, timeValout
use amr_module, only: timeBound,timeStepgrid, timeFlagger,timeBufnst
use amr_module, only: timeBoundCPU,timeStepGridCPU,timeRegriddingCPU
use amr_module, only: timeValoutCPU,timeTick,timeTickCPU
use amr_module, only: kcheck, iorder, lendim, lenmax, memsize
use amr_module, only: dprint, eprint, edebug, gprint, nprint, pprint
use amr_module, only: rprint, sprint, tprint, uprint
use amr_module, only: t0, tstart_thisrun
use amr_module, only: tick_clock_start, tick_cpu_start
! Data modules
use geoclaw_module, only: set_geo
use topo_module, only: read_topo_settings, read_dtopo_settings
use qinit_module, only: set_qinit, set_qinit_dig
use refinement_module, only: set_refinement, set_flow_grades
use storm_module, only: set_storm
use friction_module, only: setup_variable_friction
use gauges_module, only: set_gauges, num_gauges
use regions_module, only: set_regions
use fgout_module, only: set_fgout
use fgmax_module, only: set_fgmax, FG_num_fgrids
use multilayer_module, only: set_multilayer
use adjoint_module, only: read_adjoint_data
use auxinit_module, only: set_auxinit
use digclaw_module, only: set_dig, set_pinit
implicit none
! Local variables
integer :: i, iaux, mw, level
integer :: ndim, nvar, naux, mcapa1, mindim, dimensional_split
integer :: nstart, nsteps, nv1, nx, ny, lentotsave, num_gauge_SAVE
integer :: omp_get_max_threads, maxthreads
real(kind=8) :: time, ratmet, cut, dtinit, dt_max
logical :: vtime, rest, output_t0
! Timing variables
integer(kind=8) :: clock_start, clock_finish, clock_rate, ttotal, count_max
real(kind=8) :: ttotalcpu
integer(kind=8) :: tick_clock_finish
integer, parameter :: timing_unit = 48
character(len=512) :: timing_line, timing_substr
character(len=*), parameter :: timing_base_name = "timing."
character(len=*), parameter :: timing_header_format = &
"(' wall time (', i2,')," // &
" CPU time (', i2,'), " // &
"cells updated (', i2,'),')"
! Common block variables
real(kind=8) :: dxmin, dymin
common /comfine/ dxmin,dymin
character(len=364) :: format_string
character(len=*), parameter :: clawfile = 'claw.data'
character(len=*), parameter :: amrfile = 'amr.data'
character(len=*), parameter :: outfile = 'fort.amr'
character(len=*), parameter :: dbugfile = 'fort.debug'
character(len=*), parameter :: matfile = 'fort.nplot'
character(len=*), parameter :: parmfile = 'fort.parameters'
! Open parameter and debug files
open(dbugunit,file=dbugfile,status='unknown',form='formatted')
open(parmunit,file=parmfile,status='unknown',form='formatted')
maxthreads = 1 !! default, if no openmp
! Open AMRClaw primary parameter file
call opendatafile(inunit,clawfile)
! Number of space dimensions, not really a parameter but we read it in and
! check to make sure everyone is on the same page.
read(inunit,"(i1)") ndim
if (ndim /= 2) then
print *,'Error *** ndim = 2 is required, ndim = ',ndim
print *,'*** Are you sure input has been converted'
print *,'*** to Clawpack 5.x form?'
stop
endif
! Domain variables
read(inunit,*) xlower, ylower
read(inunit,*) xupper, yupper
read(inunit,*) nx, ny
read(inunit,*) nvar ! meqn
read(inunit,*) mwaves
read(inunit,*) naux
read(inunit,*) t0
! ==========================================================================
! Output Options
! Output style
read(inunit,*) output_style
if (output_style == 1) then
read(inunit,*) nout
read(inunit,*) tfinal
read(inunit,*) output_t0
iout = 0
else if (output_style == 2) then
read(inunit,*) nout
allocate(tout(nout))
read(inunit,*) (tout(i), i=1,nout)
output_t0 = (abs(tout(1) - t0) < 1e-15)
! Move output times down one index
if (output_t0) then
nout = nout - 1
do i=1,nout
tout(i) = tout(i+1)
enddo
endif
iout = 0
tfinal = tout(nout)
else if (output_style == 3) then
read(inunit,*) iout
read(inunit,*) nstop
read(inunit,*) output_t0
nout = 0
tfinal = rinfinity
else
stop "Error *** Invalid output style."
endif
! Error checking
if ((output_style == 1) .and. (nout > 0)) then
allocate(tout(nout))
do i=1,nout
tout(i) = t0 + i * (tfinal - t0) / real(nout,kind=8)
enddo
endif
! What and how to output
read(inunit,*) output_format
allocate(output_q_components(nvar))
read(inunit,*) (output_q_components(i),i=1,nvar)
if (naux > 0) then
allocate(output_aux_components(naux))
read(inunit,*) (output_aux_components(i),i=1,naux)
read(inunit,*) output_aux_onlyonce
endif
! ==========================================================================
! ==========================================================================
! Algorithm parameters
read(inunit,*) possk(1) ! dt_initial
read(inunit,*) dt_max ! largest allowable dt
read(inunit,*) cflv1 ! cfl_max
read(inunit,*) cfl ! clf_desired
read(inunit,*) nv1 ! steps_max
if (output_style /= 3) then
nstop = nv1
endif
read(inunit,*) vtime ! dt_variable
if (vtime) then
method(1) = 2
else
method(1) = 1
endif
read(inunit,*) method(2) ! order
iorder = method(2)
read(inunit,*) method(3) ! order_trans
read(inunit,*) dimensional_split
if (dimensional_split > 0) then
print *, '*** ERROR *** dimensional_split = ', dimensional_split
print *, ' dimensional splitting not supported in amrclaw'
stop
endif
read(inunit,*) method(4) ! verbosity
read(inunit,*) method(5) ! src_split
read(inunit,*) mcapa1
read(inunit,*) use_fwaves
allocate(mthlim(mwaves))
read(inunit,*) (mthlim(mw), mw=1,mwaves)
! Boundary conditions
read(inunit,*) nghost
read(inunit,*) mthbc(1),mthbc(3)
read(inunit,*) mthbc(2),mthbc(4)
! 1 = left, 2 = right 3 = bottom 4 = top boundary
xperdom = (mthbc(1) == 2 .and. mthbc(2) == 2)
yperdom = (mthbc(3) == 2 .and. mthbc(4) == 2)
spheredom = (mthbc(3) == 5 .and. mthbc(4) == 5)
if ((mthbc(1).eq.2 .and. mthbc(2).ne.2) .or. &
(mthbc(2).eq.2 .and. mthbc(1).ne.2)) then
print *, '*** ERROR *** periodic boundary conditions: '
print *, ' mthbc(1) and mthbc(2) must BOTH be set to 2'
stop
endif
if ((mthbc(3).eq.2 .and. mthbc(4).ne.2) .or. &
(mthbc(4).eq.2 .and. mthbc(3).ne.2)) then
print *, '*** ERROR *** periodic boundary conditions: '
print *, ' mthbc(3) and mthbc(4) must BOTH be set to 2'
stop
endif
if ((mthbc(3).eq.5 .and. mthbc(4).ne.5) .or. &
(mthbc(4).eq.5 .and. mthbc(3).ne.5)) then
print *, '*** ERROR *** sphere bcs at top and bottom: '
print *, ' mthbc(3) and mthbc(4) must BOTH be set to 5'
stop
endif
if (spheredom .and. .not. xperdom) then
print *,'*** ERROR *** sphere bcs at top and bottom: '
print *,'must go with periodic bcs at left and right '
stop
endif
! ==========================================================================
! Restart and Checkpointing
read(inunit,*) rest
read(inunit,*) rstfile
read(inunit,*) checkpt_style
! default value unless set by checkpt_style==3:
checkpt_interval = iinfinity
if (abs(checkpt_style) == 2) then
read(inunit,*) nchkpt
allocate(tchk(nchkpt))
read(inunit,*) (tchk(i), i=1,nchkpt)
else if (abs(checkpt_style) == 3) then
! Checkpoint every checkpt_interval steps on coarse grid
read(inunit,*) checkpt_interval
endif
close(inunit)
! ==========================================================================
! Refinement Control
call opendatafile(inunit, amrfile)
read(inunit,*) memsize ! initial size of alloc array
read(inunit,*) max1d ! max size of each grid patch
read(inunit,*) mxnest
if (mxnest <= 0) then
stop 'Error *** mxnest (amrlevels_max) <= 0 not allowed'
endif
if (mxnest > maxlv) then
stop 'Error *** mxnest > max. allowable levels (maxlv) in common'
endif
! Anisotropic refinement always allowed in 5.x:
read(inunit,*) (intratx(i),i=1,max(1,mxnest-1))
read(inunit,*) (intraty(i),i=1,max(1,mxnest-1))
read(inunit,*) (kratio(i), i=1,max(1,mxnest-1))
read(inunit,*)
do i=1,mxnest-1
if ((intratx(i) > max1d) .or. (intraty(i) > max1d)) then
print *, ""
format_string = "(' *** Error: Refinement ratios must be no " // &
"larger than max1d = ',i5,/,' (set max1d" // &
" in amr_module.f90)')"
print format_string, max1d
stop
endif
enddo
if (naux > 0) then
allocate(auxtype(naux))
read(inunit,*) (auxtype(iaux), iaux=1,naux)
endif
read(inunit,*)
read(inunit,*) flag_richardson
read(inunit,*) tol ! for richardson
read(inunit,*) flag_gradient
read(inunit,*) tolsp ! for gradient
read(inunit,*) kcheck
read(inunit,*) ibuff
read(inunit,*) cut
read(inunit,*) verbosity_regrid
! read verbose/debugging flags
read(inunit,*) dprint
read(inunit,*) eprint
read(inunit,*) edebug
read(inunit,*) gprint
read(inunit,*) nprint
read(inunit,*) pprint
read(inunit,*) rprint
read(inunit,*) sprint
read(inunit,*) tprint
read(inunit,*) uprint
close(inunit)
! Finished with reading in parameters
! ==========================================================================
! Look for capacity function via auxtypes:
mcapa = 0
do iaux = 1, naux
if (auxtype(iaux) == "capacity") then
if (mcapa /= 0) then
print *, " only 1 capacity allowed"
stop
else
mcapa = iaux
endif
endif
! Change to Version 4.1 terminology:
if (auxtype(iaux) == "leftface") auxtype(iaux) = "xleft"
if (auxtype(iaux) == "bottomface") auxtype(iaux) = "yleft"
if (.not. (auxtype(iaux) .eq."xleft" .or. &
auxtype(iaux) .eq. "yleft".or. &
auxtype(iaux) .eq. "capacity".or. &
auxtype(iaux) .eq. "center")) then
print *," unknown type for auxiliary variables"
print *," use xleft/yleft/center/capacity"
stop
endif
enddo
! Error checking of input data
if (mcapa /= mcapa1) then
stop 'Error *** mcapa does not agree with auxtype'
endif
if (nvar > maxvar) then
stop 'Error *** nvar > maxvar in common'
endif
if (2*nghost > min(nx,ny) .and. ny /= 1) then
mindim = 2 * nghost
print *, 'Error *** need finer domain >', mindim, ' cells'
stop
endif
if (mcapa > naux) then
stop 'Error *** mcapa > naux in input file'
endif
if (.not. vtime .and. nout /= 0) then
print *, ' cannot specify output times with fixed dt'
stop
endif
! Write out parameters
write(parmunit,*) ' '
write(parmunit,*) 'Running amrclaw with parameter values:'
write(parmunit,*) ' '
print *, ' '
print *, 'Running amrclaw ... '
print *, ' '
hxposs(1) = (xupper - xlower) / nx
hyposs(1) = (yupper - ylower) / ny
! initialize frame number for output.
! Note: might be reset in restrt if this is a restart
if (output_t0) then
matlabu = 0
else
matlabu = 1
endif
if (rest) then
open(outunit, file=outfile, status='unknown', position='append', &
form='formatted')
! moved upt before restrt or won't properly initialize
call set_fgmax()
! need these before set_gauges so num_out_vars set right for gauges
call set_geo() ! sets basic parameters g and coord system
! Set gauge output, note that restrt might reset x,y for lagrangian:
call set_gauges(rest, nvar, naux)
call restrt(nsteps,time,nvar,naux)
nstart = nsteps
tstart_thisrun = time
print *, ' '
print *, 'Restarting from previous run'
print *, ' at time = ',time
print *, ' '
! Call user routine to set up problem parameters:
call setprob()
! Non-user defined setup routine
call set_refinement() ! sets refinement control parameters
call read_dtopo_settings() ! specifies file with dtopo from earthquake
call read_topo_settings(rest) ! specifies topography (bathymetry) files
call set_auxinit() ! specifies file(s) for auxiliary variables
call set_fgout(rest,nvar) ! Fixed grid settings
call set_regions() ! Set refinement regions
call set_dig(naux)
call set_pinit()
call set_qinit_dig()
call set_flow_grades()
else
! Create new timing file
open(unit=timing_unit, file=timing_base_name//"csv", &
form='formatted', status='unknown', action='write')
! Construct header string
timing_line = 'output_time,total_wall_time,total_cpu_time,'
timing_substr = ""
do level=1, mxnest
write(timing_substr, timing_header_format) level, level, level
timing_line = trim(timing_line) // trim(timing_substr)
end do
write(timing_unit, "(a)") timing_line
close(timing_unit)
open(outunit, file=outfile, status='unknown', form='formatted')
tstart_thisrun = t0
! Call user routine to set up problem parameters:
call setprob()
! Non-user defined setup routine
call set_geo() ! sets basic parameters g and coord system
call set_refinement() ! sets refinement control parameters
call read_dtopo_settings() ! specifies file with dtopo from earthquake
call read_topo_settings(rest) ! specifies topography (bathymetry) files
call set_auxinit() ! specifies file(s) for auxiliary variables
call set_fgout(rest,nvar) ! Fixed grid settings
call set_regions() ! Set refinement regions
call set_gauges(rest, nvar, naux) ! Set gauge output
call set_fgmax()
call set_dig(naux)
call set_pinit()
call set_qinit_dig()
call set_flow_grades()
cflmax = 0.d0 ! otherwise use previously checkpointed val
lentot = 0
lenmax = 0
lendim = 0
rvol = 0.0d0
do i = 1, mxnest
rvoll(i) = 0.0d0
enddo
evol = 0.0d0
call stst1()
! changed 4/24/09: store dxmin,dymin for setaux before
! grids are made, in order to average up from finest grid.
dxmin = hxposs(mxnest)
dymin = hyposs(mxnest)
call domain(nvar,vtime,nx,ny,naux,t0)
! Hold off on gauges until grids are set.
! The fake call to advance at the very first timestep
! looks at the gauge array but it is not yet built
num_gauge_SAVE = num_gauges
num_gauges = 0
call setgrd(nvar,cut,naux,dtinit,t0)
num_gauges = num_gauge_SAVE
! commented out to match 4-x version
!!$ if (possk(1) .gt. dtinit*cflv1/cfl .and. vtime) then
!!$ ! initial time step was too large. reset to dt from setgrd
!!$ print *, "*** Initial time step reset for desired cfl"
!!$ possk(1) = dtinit
!!$ do i = 2, mxnest-1
!!$ possk(i) = possk(i-1)*kratio(i-1)
!!$ end do
!!$ endif
time = t0
nstart = 0
endif
write(parmunit,*) ' '
write(parmunit,*) '--------------------------------------------'
write(parmunit,*) ' '
write(parmunit,*) ' rest = ', rest, ' (restart?)'
write(parmunit,*) ' start time = ',time
write(parmunit,*) ' '
!$ maxthreads = omp_get_max_threads()
write(outunit,*)" max threads set to ",maxthreads
print *," max threads set to ",maxthreads
!
! print out program parameters for this run
!
format_string = "(/' amrclaw parameters:',//," // &
"' error tol ',e12.5,/," // &
"' spatial error tol ',e12.5,/," // &
"' order of integrator ',i9,/," // &
"' error checking interval ',i9,/," // &
"' buffer zone size ',i9,/," // &
"' nghost ',i9,/," // &
"' volume ratio cutoff ',e12.5,/," // &
"' max. refinement level ',i9,/," // &
"' user sub. calling times ',i9,/," // &
"' cfl # (if var. delt) ',e12.5,/)"
write(outunit,format_string) tol,tolsp,iorder,kcheck,ibuff,nghost,cut, &
mxnest,checkpt_interval,cfl
format_string = "(' xupper(upper corner) ',e12.5,/," // &
"' yupper(upper corner) ',e12.5,/," // &
"' xlower(lower corner) ',e12.5,/," // &
"' ylower(lower corner) ',e12.5,/," // &
"' nx = no. cells in x dir.',i9,/," // &
"' ny = no. cells in y dir.',i9,/," // &
"' refinement ratios ',6i5,/,/)"
write(outunit,format_string) xupper,yupper,xlower,ylower,nx,ny
write(outunit,"(' refinement ratios: ',6i5,/)" ) &
(intratx(i),i=1,mxnest)
write(outunit,"(' refinement ratios: ',6i5,/)" ) &
(intraty(i),i=1,mxnest)
write(outunit,"(' no. auxiliary vars. ',i9)") naux
write(outunit,"(' var ',i5,' of type ', a10)") &
(iaux,auxtype(iaux),iaux=1,naux)
if (mcapa > 0) write(outunit,"(' capacity fn. is aux. var',i9)") mcapa
print *, ' '
print *, 'Done reading data, starting computation ... '
print *, ' '
! initialize timers before calling valout:
call system_clock(tick_clock_start, clock_rate)
call cpu_time(tick_cpu_start)
call outtre (mstart,printout,nvar,naux)
write(outunit,*) " original total mass ..."
call conck(1,nvar,naux,time,rest)
if (output_t0) then
call valout(1,lfine,time,nvar,naux)
endif
close(parmunit)
! --------------------------------------------------------
! Tick is the main routine which drives the computation:
! --------------------------------------------------------
call tick(nvar,cut,nstart,vtime,time,naux,t0,rest,dt_max)
! --------------------------------------------------------
! Done with computation, cleanup:
! Print out the fgmax files
if (FG_num_fgrids > 0) call fgmax_finalize()
write(*,"('See fort.amr for more info on this run and memory usage')")
! call system_clock to get clock_finish and count_max for debug output:
call system_clock(clock_finish,clock_rate,count_max)
!output timing data
open(timing_unit, file=timing_base_name//"txt", status='unknown', &
form='formatted')
write(*,*)
write(timing_unit,*)
format_string="('============================== Timing Data ==============================')"
write(timing_unit,format_string)
write(*,format_string)
write(*,*)
write(timing_unit,*)
!Integration time
format_string="('Integration Time (stepgrid + BC + overhead)')"
write(timing_unit,format_string)
write(*,format_string)
!Advanc time
format_string="('Level Wall Time (seconds) CPU Time (seconds) Total Cell Updates')"
write(timing_unit,format_string)
write(*,format_string)
! level counters are cumulative after restart, so initialize sums to zero
! even after restart to sum up time over all levels
ttotalcpu=0.d0
ttotal=0
do level=1,mxnest
format_string="(i3,' ',1f15.3,' ',1f15.3,' ', e17.3)"
write(timing_unit,format_string) level, &
real(tvoll(level),kind=8) / real(clock_rate,kind=8), tvollCPU(level), rvoll(level)
write(*,format_string) level, &
real(tvoll(level),kind=8) / real(clock_rate,kind=8), tvollCPU(level), rvoll(level)
ttotalcpu=ttotalcpu+tvollCPU(level)
ttotal=ttotal+tvoll(level)
end do
format_string="('total ',1f15.3,' ',1f15.3,' ', e17.3)"
write(timing_unit,format_string) &
real(ttotal,kind=8) / real(clock_rate,kind=8), ttotalCPU, rvol
write(*,format_string) &
real(ttotal,kind=8) / real(clock_rate,kind=8), ttotalCPU, rvol
write(*,*)
write(timing_unit,*)
format_string="('All levels:')"
write(*,format_string)
write(timing_unit,format_string)
!stepgrid
format_string="('stepgrid ',1f15.3,' ',1f15.3,' ',e17.3)"
write(timing_unit,format_string) &
real(timeStepgrid,kind=8) / real(clock_rate,kind=8), timeStepgridCPU
write(*,format_string) &
real(timeStepgrid,kind=8) / real(clock_rate,kind=8), timeStepgridCPU
!bound
format_string="('BC/ghost cells',1f15.3,' ',1f15.3)"
write(timing_unit,format_string) &
real(timeBound,kind=8) / real(clock_rate,kind=8), timeBoundCPU
write(*,format_string) &
real(timeBound,kind=8) / real(clock_rate,kind=8), timeBoundCPU
!regridding time
format_string="('Regridding ',1f15.3,' ',1f15.3,' ')"
write(timing_unit,format_string) &
real(timeRegridding,kind=8) / real(clock_rate,kind=8), timeRegriddingCPU
write(*,format_string) &
real(timeRegridding,kind=8) / real(clock_rate,kind=8), timeRegriddingCPU
!output time
format_string="('Output (valout)',1f14.3,' ',1f15.3,' ')"
write(timing_unit,format_string) &
real(timeValout,kind=8) / real(clock_rate,kind=8), timeValoutCPU
write(*,format_string) &
real(timeValout,kind=8) / real(clock_rate,kind=8), timeValoutCPU
write(*,*)
write(timing_unit,*)
!Total Time
format_string="('Total time: ',1f15.3,' ',1f15.3,' ')"
write(*,format_string) real(timeTick,kind=8)/real(clock_rate,kind=8), &
timeTickCPU
write(timing_unit,format_string) real(timeTick,kind=8)/real(clock_rate,kind=8), &
timeTickCPU
format_string="('Using',i3,' thread(s)')"
write(timing_unit,format_string) maxthreads
write(*,format_string) maxthreads
write(*,*)
write(timing_unit,*)
write(*,"('Note: The CPU times are summed over all threads.')")
write(timing_unit,"('Note: The CPU times are summed over all threads.')")
write(*,"(' Total time includes more than the subroutines listed above')")
write(timing_unit,"(' Total time includes more than the subroutines listed above')")
if (rest) then
write(*,"(' Times for restart runs are cumulative')")
write(timing_unit,"(' Times for restart runs are cumulative')")
endif
write(*, "('Note: timings are also recorded for each output step')")
write(*, "(' in the file timing.csv.')")
write(timing_unit, "('Note: timings are also recorded for each output step')")
write(timing_unit, "(' in the file timing.csv.')")
write(*,*)
write(timing_unit,*)
! output clock_rate etc. useful in debugging or if negative time reported
format_string="('clock_rate = ',i10, ' per second, count_max = ',i23)"
!write(*,format_string) clock_rate, count_max
write(timing_unit,format_string) clock_rate, count_max
format_string="('clock_start = ',i20, ', clock_finish = ',i20)"
!write(*,format_string) tick_clock_start, clock_finish
write(timing_unit,format_string) tick_clock_start, clock_finish
format_string="('=========================================================================')"
write(timing_unit,format_string)
write(*,format_string)
write(*,*)
write(timing_unit,*)
close(timing_unit)
! Done with computation, cleanup:
lentotsave = lentot
call cleanup(nvar,naux)
if (lentot /= 0) then
write(outunit,*) lentot," words not accounted for in memory cleanup"
print *, lentot," words not accounted for in memory cleanup"
endif
!
! report on statistics
!
open(matunit,file=matfile,status='unknown',form='formatted')
write(matunit,*) matlabu-1
write(matunit,*) mxnest
close(matunit)
write(outunit,*)
write(outunit,*)
do i = 1, mxnest
if (iregridcount(i) > 0) then
write(outunit,801) i,avenumgrids(i)/iregridcount(i),iregridcount(i)
801 format("for level ",i3, " average num. grids = ",f10.2," over ",i10, &
" regridding steps")
write(outunit,802) i,numgrids(i)
802 format("for level ",i3," current num. grids = ",i7)
endif
end do
write(outunit,*)
write(outunit,*)
write(outunit,"('alloc array statistics:')")
write(outunit,"(' current alloc usage = ',i12)") lentotsave
write(outunit,"(' maximum alloc usage = ',i12)") lenmax
write(outunit,"('required alloc memsize >= ',i12,/)") lendim
write(outunit,"('number of cells advanced for time integration = ',f20.6)")&
rvol
do level = 1,mxnest
write(outunit,"(3x,'# cells advanced on level ',i4,' = ',f20.2)") &
level, rvoll(level)
enddo
write(outunit,"('number of cells advanced for error estimation = ',f20.6,/)") &
evol
if (evol + rvol > 0.d0) then
ratmet = rvol / (evol + rvol) * 100.0d0
else
ratmet = 0.0d0
endif
write(outunit,"(' percentage of cells advanced in time = ', f10.2)") ratmet
write(outunit,"(' maximum Courant number seen = ', f10.5)") cflmax
write(*,"(' maximum Courant number seen = ', f10.5)") cflmax
write(outunit,"(//,' ------ end of AMRCLAW integration -------- ')")
! Close output and debug files.
close(outunit)
close(dbugunit)
end program amr2