diff --git a/examples/multi-layer/bowl-radial/setrun.py b/examples/multi-layer/bowl-radial/setrun.py index 20f4956f6..0f67603ab 100644 --- a/examples/multi-layer/bowl-radial/setrun.py +++ b/examples/multi-layer/bowl-radial/setrun.py @@ -131,7 +131,7 @@ def setrun(claw_pkg='geoclaw'): clawdata.output_t0 = True - clawdata.output_format = 'ascii' # 'ascii' or 'binary' + clawdata.output_format = 'ascii' # 'ascii', 'binary32' or 'binary64' clawdata.output_q_components = 'all' # could be list such as [True,True] clawdata.output_aux_components = 'none' # could be list @@ -404,11 +404,11 @@ def setgeo(rundata): # [fname] rundata.qinit_data.qinitfiles.append(['hump.xyz']) - # == setfixedgrids.data values == - fixedgrids = rundata.fixed_grid_data.fixedgrids - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] return rundata # end of function setgeo diff --git a/examples/multi-layer/plane_wave/setrun.py b/examples/multi-layer/plane_wave/setrun.py index 789e1f569..7fd9a7514 100644 --- a/examples/multi-layer/plane_wave/setrun.py +++ b/examples/multi-layer/plane_wave/setrun.py @@ -178,7 +178,7 @@ def setrun(claw_pkg='geoclaw'): clawdata.output_t0 = True - clawdata.output_format = 'ascii' # 'ascii' or 'binary' + clawdata.output_format = 'ascii' # 'ascii', 'binary32' or 'binary64' clawdata.output_q_components = 'all' # need all clawdata.output_aux_components = 'all' # eta=h+B is in q diff --git a/examples/storm-surge/ike/setrun.py b/examples/storm-surge/ike/setrun.py index 0921be20b..515b47258 100644 --- a/examples/storm-surge/ike/setrun.py +++ b/examples/storm-surge/ike/setrun.py @@ -386,11 +386,11 @@ def setgeo(rundata): topo_path = os.path.join(scratch_dir, 'gulf_caribbean.tt3') topo_data.topofiles.append([3, topo_path]) - # == setfixedgrids.data values == - rundata.fixed_grid_data.fixedgrids = [] - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] # ================ # Set Surge Data diff --git a/examples/storm-surge/isaac/setrun.py b/examples/storm-surge/isaac/setrun.py index 9df75b1b4..0112b09c1 100644 --- a/examples/storm-surge/isaac/setrun.py +++ b/examples/storm-surge/isaac/setrun.py @@ -385,11 +385,11 @@ def setgeo(rundata): topo_path = os.path.join(scratch_dir, 'gulf_caribbean.tt3') topo_data.topofiles.append([3, topo_path]) - # == setfixedgrids.data values == - rundata.fixed_grid_data.fixedgrids = [] - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] # ================ # Set Surge Data diff --git a/examples/tsunami/bowl-radial/setrun.py b/examples/tsunami/bowl-radial/setrun.py index 6f7a6d7ce..e6861f111 100644 --- a/examples/tsunami/bowl-radial/setrun.py +++ b/examples/tsunami/bowl-radial/setrun.py @@ -397,11 +397,11 @@ def setgeo(rundata): # [fname] rundata.qinit_data.qinitfiles.append(['hump.xyz']) - # == setfixedgrids.data values == - fixedgrids = rundata.fixed_grid_data.fixedgrids - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] return rundata # end of function setgeo diff --git a/examples/tsunami/bowl-slosh-netcdf/setrun.py b/examples/tsunami/bowl-slosh-netcdf/setrun.py index 7e9f9a23e..052e4c4dd 100644 --- a/examples/tsunami/bowl-slosh-netcdf/setrun.py +++ b/examples/tsunami/bowl-slosh-netcdf/setrun.py @@ -374,11 +374,11 @@ def setgeo(rundata): # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) # [fname] - # == setfixedgrids.data values == - fixedgrids = rundata.fixed_grid_data - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] return rundata # end of function setgeo diff --git a/examples/tsunami/bowl-slosh/setrun.py b/examples/tsunami/bowl-slosh/setrun.py index 594bfd306..b11753e01 100644 --- a/examples/tsunami/bowl-slosh/setrun.py +++ b/examples/tsunami/bowl-slosh/setrun.py @@ -374,11 +374,11 @@ def setgeo(rundata): # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) # [fname] - # == setfixedgrids.data values == - fixedgrids = rundata.fixed_grid_data - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] return rundata # end of function setgeo diff --git a/examples/tsunami/chile2010/setrun.py b/examples/tsunami/chile2010/setrun.py index 2c6418a03..20d66d9a2 100644 --- a/examples/tsunami/chile2010/setrun.py +++ b/examples/tsunami/chile2010/setrun.py @@ -410,11 +410,11 @@ def setgeo(rundata): # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) # [fname] - # == setfixedgrids.data values == - fixed_grids = rundata.fixed_grid_data - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] return rundata # end of function setgeo diff --git a/examples/tsunami/chile2010_adjoint/adjoint/setrun.py b/examples/tsunami/chile2010_adjoint/adjoint/setrun.py index e7e64d3ac..ff70a2e10 100644 --- a/examples/tsunami/chile2010_adjoint/adjoint/setrun.py +++ b/examples/tsunami/chile2010_adjoint/adjoint/setrun.py @@ -437,11 +437,11 @@ def setgeo(rundata): rundata.qinit_data.qinitfiles.append(['hump.xyz']) - # == setfixedgrids.data values == - fixed_grids = rundata.fixed_grid_data - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] return rundata # end of function setgeo diff --git a/examples/tsunami/chile2010_adjoint/setrun.py b/examples/tsunami/chile2010_adjoint/setrun.py index 7bebdb6bc..d40883c1a 100644 --- a/examples/tsunami/chile2010_adjoint/setrun.py +++ b/examples/tsunami/chile2010_adjoint/setrun.py @@ -487,11 +487,11 @@ def setgeo(rundata): # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) # [fname] - # == setfixedgrids.data values == - fixed_grids = rundata.fixed_grid_data - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] return rundata # end of function setgeo diff --git a/examples/tsunami/chile2010_fgmax-fgout/32412_notide.txt b/examples/tsunami/chile2010_fgmax-fgout/32412_notide.txt new file mode 100644 index 000000000..314d6058a --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/32412_notide.txt @@ -0,0 +1,1322 @@ +-1.361400000000000000e+05 7.166830903770460282e-03 +-1.352400000000000000e+05 5.192393528886896092e-03 +-1.343400000000000000e+05 3.049766923140850849e-03 +-1.334400000000000000e+05 3.167542919982224703e-03 +-1.325400000000000000e+05 3.978815151640446857e-03 +-1.316400000000000000e+05 1.914331188345386181e-03 +-1.307400000000000000e+05 3.956733680752222426e-04 +-1.298400000000000000e+05 -1.171423051346209832e-03 +-1.289400000000000000e+05 -3.403507979783171322e-03 +-1.280400000000000000e+05 -5.945507909927982837e-03 +-1.271400000000000000e+05 -5.476373226883879397e-03 +-1.262400000000000000e+05 -6.714191012179071549e-03 +-1.253400000000000000e+05 -5.420681160103413276e-03 +-1.244400000000000000e+05 -5.405003747000591829e-03 +-1.235400000000000000e+05 -4.526816206634975970e-03 +-1.226400000000000000e+05 -4.698531149188056588e-03 +-1.217400000000000000e+05 -3.886738296387193259e-03 +-1.208400000000000000e+05 -2.112767417202121578e-03 +-1.199400000000000000e+05 -2.452383099807775579e-03 +-1.190400000000000000e+05 -2.034615917182236444e-03 +-1.181400000000000000e+05 -2.039748725110257510e-03 +-1.172400000000000000e+05 -2.696490304515464231e-03 +-1.163400000000000000e+05 -6.278381737502058968e-03 +-1.154400000000000000e+05 -3.099493304944189731e-03 +-1.145400000000000000e+05 -2.509481320885242894e-03 +-1.136400000000000000e+05 -1.888084405436529778e-03 +-1.127400000000000000e+05 -6.391481965692946687e-04 +-1.118400000000000000e+05 2.815725126311008353e-03 +-1.109400000000000000e+05 2.043799870079965331e-03 +-1.100400000000000000e+05 4.607995247170038056e-03 +-1.091400000000000000e+05 5.073729485957301222e-03 +-1.082400000000000000e+05 5.015708634346083272e-03 +-1.073400000000000000e+05 4.024550333269871771e-03 +-1.064400000000000000e+05 5.713135005862568505e-03 +-1.055400000000000000e+05 5.722581552618066780e-03 +-1.046400000000000000e+05 6.727750677782751154e-03 +-1.037400000000000000e+05 7.442186356456659269e-03 +-1.028400000000000000e+05 5.622415405014180578e-03 +-1.019400000000000000e+05 7.071535187606059480e-03 +-1.010400000000000000e+05 6.642031114097335376e-03 +-1.001400000000000000e+05 7.237778303533559665e-03 +-9.924000000000000000e+04 6.815194795308343600e-03 +-9.834000000000000000e+04 5.383527478443284053e-03 +-9.744000000000000000e+04 4.004266485935659148e-03 +-9.654000000000000000e+04 2.789697306070593186e-03 +-9.564000000000000000e+04 1.900614667647460010e-03 +-9.474000000000000000e+04 1.543235537610598840e-03 +-9.384000000000000000e+04 -2.034638288932910655e-03 +-9.294000000000000000e+04 -2.548144789216166828e-03 +-9.204000000000000000e+04 -1.680499643953226041e-03 +-9.114000000000000000e+04 9.113371679632109590e-04 +-9.024000000000000000e+04 -4.094409205208648928e-04 +-8.934000000000000000e+04 -1.265924705876386724e-03 +-8.844000000000000000e+04 7.258440418809186667e-04 +-8.754000000000000000e+04 -2.049913673545233905e-03 +-8.664000000000000000e+04 7.844925448807771318e-04 +-8.574000000000000000e+04 -4.064658096467610449e-04 +-8.484000000000000000e+04 -1.278008501685690135e-03 +-8.394000000000000000e+04 -1.511195392595254816e-03 +-8.304000000000000000e+04 -3.818649332970380783e-03 +-8.214000000000000000e+04 -4.949757578287972137e-03 +-8.124000000000000000e+04 -3.695264224006677978e-03 +-8.034000000000000000e+04 -6.891176341923710424e-03 +-7.944000000000000000e+04 -8.421917232226405758e-03 +-7.854000000000000000e+04 -8.222672488045645878e-03 +-7.764000000000000000e+04 -6.280887874709151220e-03 +-7.674000000000000000e+04 -4.636891345398908015e-03 +-7.584000000000000000e+04 -2.383626302616903558e-03 +-7.494000000000000000e+04 -3.665497244583093561e-03 +-7.404000000000000000e+04 -3.676343832012207713e-03 +-7.314000000000000000e+04 -1.656573484069667757e-03 +-7.224000000000000000e+04 -8.894964821593021043e-04 +-7.134000000000000000e+04 3.030791685887379572e-04 +-7.044000000000000000e+04 -1.434076209989143535e-03 +-6.954000000000000000e+04 -4.839475877815857530e-04 +-6.864000000000000000e+04 -2.511090369807789102e-04 +-6.774000000000000000e+04 -1.551565173940616660e-04 +-6.684000000000000000e+04 -1.623845390895439778e-03 +-6.594000000000000000e+04 -1.086043585928564426e-03 +-6.504000000000000000e+04 1.035385454997594934e-03 +-6.414000000000000000e+04 4.330654251134546939e-03 +-6.324000000000000000e+04 3.409961786019266583e-03 +-6.234000000000000000e+04 5.910129978474287782e-03 +-6.144000000000000000e+04 4.500856836784805637e-03 +-6.054000000000000000e+04 8.904833739507012069e-04 +-5.964000000000000000e+04 -1.168821790088259149e-03 +-5.874000000000000000e+04 -1.876546377388876863e-03 +-5.784000000000000000e+04 -3.380063277290901169e-03 +-5.694000000000000000e+04 -3.771668473746103700e-03 +-5.604000000000000000e+04 -4.086479399120435119e-03 +-5.514000000000000000e+04 -7.301231650671979878e-03 +-5.424000000000000000e+04 -7.333997145906323567e-03 +-5.334000000000000000e+04 -7.044832545034296345e-03 +-5.244000000000000000e+04 -6.237351459276396781e-03 +-5.154000000000000000e+04 -6.661199165137077216e-03 +-5.064000000000000000e+04 -4.015394015368656255e-03 +-4.974000000000000000e+04 -9.524857805445208214e-04 +-4.884000000000000000e+04 9.165325009234948084e-04 +-4.794000000000000000e+04 2.016634205574519001e-03 +-4.704000000000000000e+04 1.802577170565200504e-03 +-4.614000000000000000e+04 7.522262303609750234e-04 +-4.524000000000000000e+04 3.359484570864879061e-03 +-4.434000000000000000e+04 5.126945161464391276e-03 +-4.344000000000000000e+04 3.558377699846460018e-03 +-4.254000000000000000e+04 4.151169969190959819e-03 +-4.164000000000000000e+04 2.388843199696566444e-03 +-4.074000000000000000e+04 3.733760093382443301e-03 +-3.984000000000000000e+04 3.620141766077722423e-03 +-3.894000000000000000e+04 1.447504466341342777e-03 +-3.804000000000000000e+04 2.574621161329559982e-03 +-3.714000000000000000e+04 1.314104478296940215e-03 +-3.624000000000000000e+04 3.927697733161039650e-03 +-3.534000000000000000e+04 3.622349888246390037e-03 +-3.444000000000000000e+04 8.547137715140706860e-03 +-3.354000000000000000e+04 8.791084773292823229e-03 +-3.264000000000000000e+04 8.381913150515174493e-03 +-3.174000000000000000e+04 8.285748347589105833e-03 +-3.084000000000000000e+04 7.407783040434878785e-03 +-2.994000000000000000e+04 7.593889942654641345e-03 +-2.904000000000000000e+04 5.633158491946232971e-03 +-2.814000000000000000e+04 6.261315938900224864e-03 +-2.724000000000000000e+04 7.164978484979656059e-03 +-2.634000000000000000e+04 8.986665346128575038e-03 +-2.544000000000000000e+04 7.330496126087382436e-03 +-2.454000000000000000e+04 5.768480594269931316e-03 +-2.364000000000000000e+04 5.847301020367012825e-03 +-2.274000000000000000e+04 2.095478522278426681e-03 +-2.184000000000000000e+04 2.030809244388365187e-03 +-2.094000000000000000e+04 1.167951368188369088e-03 +-2.004000000000000000e+04 1.026041643854114227e-03 +-1.914000000000000000e+04 -8.637805103717255406e-04 +-1.824000000000000000e+04 -9.510617401247145608e-04 +-1.734000000000000000e+04 -6.590536650037392974e-04 +-1.644000000000000000e+04 -1.378173702505591791e-03 +-1.554000000000000000e+04 1.539987219075555913e-03 +-1.464000000000000000e+04 7.880139683038578369e-04 +-1.374000000000000000e+04 -1.892712265544105321e-03 +-1.284000000000000000e+04 -7.081964513417915441e-04 +-1.194000000000000000e+04 -3.808876181210507639e-03 +-1.104000000000000000e+04 -4.287579728952550795e-03 +-1.014000000000000000e+04 -6.178427956001542043e-03 +-9.240000000000000000e+03 -2.456697719935618807e-03 +-8.340000000000000000e+03 -5.039648497586313169e-03 +-7.440000000000000000e+03 -6.788298409446724690e-03 +-6.540000000000000000e+03 -5.510120729923073668e-03 +-5.640000000000000000e+03 -5.962616538454312831e-03 +-5.640000000000000000e+03 -5.962616538454312831e-03 +-5.580000000000000000e+03 -4.683060635215952061e-03 +-5.520000000000000000e+03 -5.414283493337279651e-03 +-5.460000000000000000e+03 -6.156192179332720116e-03 +-5.400000000000000000e+03 -4.908692927529045846e-03 +-5.340000000000000000e+03 -4.671691250223375391e-03 +-5.280000000000000000e+03 -5.445091841465909965e-03 +-5.220000000000000000e+03 -4.228798661642940715e-03 +-5.160000000000000000e+03 -5.022714926781191025e-03 +-5.100000000000000000e+03 -4.826743041121517308e-03 +-5.040000000000000000e+03 -4.640784728508151602e-03 +-4.980000000000000000e+03 -4.464740964976954274e-03 +-4.920000000000000000e+03 -4.298511983506614342e-03 +-4.860000000000000000e+03 -6.141997314443869982e-03 +-4.800000000000000000e+03 -4.995095787307946011e-03 +-4.740000000000000000e+03 -5.857705488779174630e-03 +-4.740000000000000000e+03 -5.857705488779174630e-03 +-4.680000000000000000e+03 -5.729723853619361762e-03 +-4.620000000000000000e+03 -6.611047622754995245e-03 +-4.560000000000000000e+03 -6.501572832348756492e-03 +-4.500000000000000000e+03 -7.401194896374363452e-03 +-4.440000000000000000e+03 -6.309808503829117399e-03 +-4.380000000000000000e+03 -6.227307743756682612e-03 +-4.320000000000000000e+03 -6.153586042273673229e-03 +-4.260000000000000000e+03 -6.088536176321213134e-03 +-4.200000000000000000e+03 -7.032050317320681643e-03 +-4.140000000000000000e+03 -7.984020001458702609e-03 +-4.080000000000000000e+03 -6.944336142623797059e-03 +-4.020000000000000000e+03 -6.912889074556005653e-03 +-3.960000000000000000e+03 -6.889568535370926838e-03 +-3.900000000000000000e+03 -6.874263647659972776e-03 +-3.840000000000000000e+03 -7.866863009439839516e-03 +-3.840000000000000000e+03 -6.866863009236112703e-03 +-3.780000000000000000e+03 -7.867254576012783218e-03 +-3.720000000000000000e+03 -6.875325803775922395e-03 +-3.660000000000000000e+03 -7.890963590398314409e-03 +-3.600000000000000000e+03 -8.914054254091752227e-03 +-3.540000000000000000e+03 -7.944483624669373967e-03 +-3.480000000000000000e+03 -7.982136981809162535e-03 +-3.420000000000000000e+03 -8.026899079595750663e-03 +-3.360000000000000000e+03 -7.078654190991073847e-03 +-3.300000000000000000e+03 -6.137286090051929932e-03 +-3.240000000000000000e+03 -6.202678036061115563e-03 +-3.180000000000000000e+03 -5.274712822028959636e-03 +-3.120000000000000000e+03 -5.353272787942842115e-03 +-3.060000000000000000e+03 -5.438239775685360655e-03 +-3.000000000000000000e+03 -5.529495215341739822e-03 +-2.940000000000000000e+03 -5.626920054055517539e-03 +-2.940000000000000000e+03 -5.626920054055517539e-03 +-2.880000000000000000e+03 -4.730394849502772558e-03 +-2.820000000000000000e+03 -4.839799689762003254e-03 +-2.760000000000000000e+03 -3.955014265557110775e-03 +-2.700000000000000000e+03 -4.075917870068224147e-03 +-2.640000000000000000e+03 -5.202389403166307602e-03 +-2.580000000000000000e+03 -3.334307355544297025e-03 +-2.520000000000000000e+03 -4.471549847039568704e-03 +-2.460000000000000000e+03 -4.613994645296770614e-03 +-2.400000000000000000e+03 -4.761519135172420647e-03 +-2.340000000000000000e+03 -4.914000368444249034e-03 +-2.280000000000000000e+03 -4.071315071087155957e-03 +-2.220000000000000000e+03 -4.233339584061468486e-03 +-2.160000000000000000e+03 -4.399949987600848544e-03 +-2.100000000000000000e+03 -3.571022012692992575e-03 +-2.040000000000000000e+03 -3.746431107174430508e-03 +-2.040000000000000000e+03 -3.746431107174430508e-03 +-1.980000000000000000e+03 -3.926052432689175475e-03 +-1.920000000000000000e+03 -4.109760842766263522e-03 +-1.860000000000000000e+03 -4.297430929000256583e-03 +-1.800000000000000000e+03 -4.488937039241136517e-03 +-1.740000000000000000e+03 -4.684153235757548828e-03 +-1.680000000000000000e+03 -3.882953359607199673e-03 +-1.620000000000000000e+03 -3.085211033976520412e-03 +-1.560000000000000000e+03 -2.290799595357384533e-03 +-1.500000000000000000e+03 -3.499592238767945673e-03 +-1.440000000000000000e+03 -2.711461899707501289e-03 +-1.380000000000000000e+03 -1.926281353007652797e-03 +-1.320000000000000000e+03 -2.143923195035313256e-03 +-1.260000000000000000e+03 -1.364259792353550438e-03 +-1.200000000000000000e+03 -1.587163406838953961e-03 +-1.140000000000000000e+03 -8.125061240207287483e-04 +-1.080000000000000000e+03 -4.015986996819265187e-05 +-1.020000000000000000e+03 -2.699964516068575904e-04 +-9.600000000000000000e+02 -5.018875635869335383e-04 +-9.000000000000000000e+02 -7.357047315963427536e-04 +-8.400000000000000000e+02 -9.713194340292830020e-04 +-7.800000000000000000e+02 -2.086030362988822162e-04 +-7.200000000000000000e+02 -4.474267925616004504e-04 +-6.600000000000000000e+02 -6.876619290778762661e-04 +-6.000000000000000000e+02 -9.291795440731220879e-04 +-5.400000000000000000e+02 -1.718507255645818077e-04 +-4.800000000000000000e+02 -4.155465185249340720e-04 +-4.200000000000000000e+02 1.339862101303879172e-03 +-3.600000000000000000e+02 9.450415655010147020e-05 +-3.000000000000000000e+02 8.485086973450961523e-04 +-2.400000000000000000e+02 6.020047858328325674e-04 +-1.800000000000000000e+02 -6.448785325119388290e-04 +-1.200000000000000000e+02 1.079877729353029281e-04 +-6.000000000000000000e+01 8.607327317804447375e-04 +0.000000000000000000e+00 6.134853083494817838e-04 +6.000000000000000000e+01 1.366374439385253936e-03 +1.200000000000000000e+02 1.195290251416736282e-04 +1.800000000000000000e+02 1.873077864729566500e-03 +2.400000000000000000e+02 1.627149719752196688e-03 +3.000000000000000000e+02 3.818732366198673844e-04 +3.600000000000000000e+02 8.137376994454825763e-03 +4.200000000000000000e+02 -9.106210538448067382e-03 +4.800000000000000000e+02 1.165123898772435496e-02 +5.400000000000000000e+02 -1.590146212038234808e-03 +6.000000000000000000e+02 -2.883023802587558748e-02 +6.000000000000000000e+02 1.016976197388430592e-02 +6.000000000000000000e+02 2.616976197441545082e-02 +6.000000000000000000e+02 -1.468302380262684892e-01 +6.600000000000000000e+02 4.493109150280361064e-02 +6.600000000000000000e+02 2.889310915034002392e-01 +6.600000000000000000e+02 -2.330689084965342772e-01 +6.600000000000000000e+02 -2.280689084973346326e-01 +6.600000000000000000e+02 3.529310915027963347e-01 +7.200000000000000000e+02 6.269397015148570063e-02 +7.200000000000000000e+02 2.566939701519004302e-01 +7.200000000000000000e+02 -4.530602984868892236e-02 +7.200000000000000000e+02 -9.306029848630714696e-03 +7.200000000000000000e+02 4.969397015156573616e-02 +7.800000000000000000e+02 -1.354147441179520683e-02 +7.800000000000000000e+02 4.585255883284844458e-04 +7.800000000000000000e+02 2.945852558787009912e-02 +7.800000000000000000e+02 1.458525587622716557e-03 +7.800000000000000000e+02 -8.454147441170789534e-02 +8.400000000000000000e+02 -8.977511474950006232e-02 +8.400000000000000000e+02 -5.277511475014762254e-02 +9.000000000000000000e+02 4.199317643997346750e-02 +9.600000000000000000e+02 -4.236473771015880629e-03 +1.020000000000000000e+03 2.353606151245912770e-02 +1.080000000000000000e+03 -2.468909102117322618e-02 +1.140000000000000000e+03 2.408819507581938524e-02 +1.200000000000000000e+03 -1.813195399518008344e-02 +1.260000000000000000e+03 1.465058776830119314e-02 +1.320000000000000000e+03 3.043594610062427819e-02 +1.380000000000000000e+03 -3.775753549234650563e-03 +1.440000000000000000e+03 8.015614038413332310e-03 +1.500000000000000000e+03 -2.318982622909970814e-02 +1.560000000000000000e+03 -1.739194973742996808e-02 +1.620000000000000000e+03 4.093678744538920000e-04 +1.680000000000000000e+03 1.421425060289038811e-02 +1.740000000000000000e+03 -5.977177837849012576e-03 +1.800000000000000000e+03 1.083520591964770574e-02 +1.860000000000000000e+03 1.765152490224863868e-02 +1.920000000000000000e+03 -7.528098192778998055e-03 +1.980000000000000000e+03 -5.703541031834902242e-03 +2.040000000000000000e+03 -8.874681662746297661e-03 +2.100000000000000000e+03 1.495860151499073254e-02 +2.160000000000000000e+03 1.796429703063040506e-03 +2.220000000000000000e+03 -3.610763151300488971e-04 +2.280000000000000000e+03 -4.513796096034639049e-03 +2.340000000000000000e+03 2.733839031952811638e-02 +2.400000000000000000e+03 -8.804397477433667518e-03 +2.460000000000000000e+03 -4.942040362948318943e-03 +2.520000000000000000e+03 7.925580341179738753e-03 +2.580000000000000000e+03 6.798582913688733242e-03 +2.640000000000000000e+03 -3.229148787795566022e-04 +2.700000000000000000e+03 1.056120426892448450e-02 +2.760000000000000000e+03 -5.548942787754640449e-03 +2.820000000000000000e+03 6.346760295855347067e-03 +2.880000000000000000e+03 6.248429353945539333e-03 +2.940000000000000000e+03 -8.438202294200891629e-04 +3.000000000000000000e+03 -9.298736385972006246e-04 +3.060000000000000000e+03 6.990383433731039986e-03 +3.120000000000000000e+03 -2.082935216094483621e-03 +3.180000000000000000e+03 3.850283623251016252e-03 +3.240000000000000000e+03 7.790152677443984430e-03 +3.300000000000000000e+03 -2.263215939819929190e-03 +3.360000000000000000e+03 3.690289326186757535e-03 +3.420000000000000000e+03 5.650779483403312042e-03 +3.480000000000000000e+03 4.618364915586425923e-03 +3.540000000000000000e+03 -5.406844557910517324e-03 +3.600000000000000000e+03 2.575260280536895152e-03 +3.660000000000000000e+03 3.564788019502884708e-03 +3.720000000000000000e+03 -1.438153361959848553e-03 +3.780000000000000000e+03 5.566543494751385879e-03 +3.840000000000000000e+03 5.789852984889876097e-04 +3.900000000000000000e+03 8.599278125984710641e-03 +3.960000000000000000e+03 3.627527412390918471e-03 +4.020000000000000000e+03 2.663837894942844287e-03 +4.080000000000000000e+03 -1.291686317017592955e-03 +4.140000000000000000e+03 6.761058180018153507e-03 +4.200000000000000000e+03 8.221741190936882049e-04 +4.260000000000000000e+03 3.891763555657234974e-03 +4.320000000000000000e+03 3.969927804973849561e-03 +4.380000000000000000e+03 5.056767494352243375e-03 +4.440000000000000000e+03 -1.847617461862682831e-03 +4.500000000000000000e+03 7.256872136167658027e-03 +4.560000000000000000e+03 -3.629665278822358232e-03 +4.620000000000000000e+03 6.492868023087794427e-03 +4.680000000000000000e+03 4.624568991857813671e-03 +4.740000000000000000e+03 2.765533851743384730e-03 +4.800000000000000000e+03 9.158580332950805314e-04 +4.860000000000000000e+03 4.075636226843926124e-03 +4.920000000000000000e+03 1.244962330929411110e-03 +4.980000000000000000e+03 2.423929450742434710e-03 +5.040000000000000000e+03 6.126298994786338881e-04 +5.100000000000000000e+03 2.811155203744419850e-03 +5.160000000000000000e+03 2.019596083300712053e-03 +5.220000000000000000e+03 2.380424220973509364e-04 +5.280000000000000000e+03 1.466583304136293009e-03 +5.340000000000000000e+03 2.705306980715249665e-03 +5.400000000000000000e+03 1.954300832949229516e-03 +5.460000000000000000e+03 -1.786348530004033819e-03 +5.520000000000000000e+03 3.483444581434014253e-03 +5.580000000000000000e+03 1.763765049872745294e-03 +5.640000000000000000e+03 5.469686311698751524e-05 +5.700000000000000000e+03 3.356323141815664712e-03 +5.760000000000000000e+03 6.687261538900202140e-04 +5.820000000000000000e+03 9.919872700265841559e-04 +5.880000000000000000e+03 1.326186957157915458e-03 +5.940000000000000000e+03 1.671404813350818586e-03 +6.000000000000000000e+03 -2.972280484755174257e-03 +6.060000000000000000e+03 6.395208820322295651e-03 +6.120000000000000000e+03 -6.226050385521375574e-03 +6.180000000000000000e+03 1.164017818155116402e-03 +6.240000000000000000e+03 2.565488453910802491e-03 +6.300000000000000000e+03 4.978435593329777475e-03 +6.360000000000000000e+03 -2.597067645183415152e-03 +6.420000000000000000e+03 -1.160949032055214047e-03 +6.480000000000000000e+03 -4.713137342150730547e-03 +6.540000000000000000e+03 -2.253562298392353114e-03 +6.600000000000000000e+03 -2.782154560009075794e-03 +6.660000000000000000e+03 -2.298845778568647802e-03 +6.720000000000000000e+03 3.196431408468924928e-03 +6.780000000000000000e+03 -2.962565577036002651e-04 +6.840000000000000000e+03 -7.768442783344653435e-04 +6.900000000000000000e+03 -8.245267326856264845e-03 +6.960000000000000000e+03 4.298537744944042061e-03 +7.020000000000000000e+03 -2.145366650438518263e-03 +7.080000000000000000e+03 -2.576919053353776690e-03 +7.140000000000000000e+03 3.940942406188696623e-06 +7.200000000000000000e+03 -4.027272698294837028e-04 +7.260000000000000000e+03 2.031347121373983100e-04 +7.320000000000000000e+03 -1.784157648216933012e-04 +7.380000000000000000e+03 -1.547322365695436019e-03 +7.440000000000000000e+03 -4.903529785224236548e-03 +7.500000000000000000e+03 -2.246983776785782538e-03 +7.560000000000000000e+03 4.223688711135764606e-04 +7.620000000000000000e+03 -2.895419661399500910e-03 +7.680000000000000000e+03 -2.200298276875400916e-03 +7.740000000000000000e+03 -2.492216912287403829e-03 +7.800000000000000000e+03 2.288734449393814430e-04 +7.860000000000000000e+03 -4.036979295960918535e-03 +7.920000000000000000e+03 -1.289728247684251983e-03 +7.980000000000000000e+03 -1.529327641037525609e-03 +8.040000000000000000e+03 1.244267252332065254e-03 +8.100000000000000000e+03 -9.688999298305134289e-04 +8.160000000000000000e+03 -2.168786641050246544e-03 +8.220000000000000000e+03 -1.355351429992879275e-03 +8.280000000000000000e+03 -5.285539155011065304e-04 +8.340000000000000000e+03 -6.883548112455173396e-04 +8.400000000000000000e+03 2.165284043257997837e-03 +8.460000000000000000e+03 -2.967600261399638839e-03 +8.520000000000000000e+03 2.913028236434911378e-03 +8.580000000000000000e+03 1.807204415854357649e-03 +8.640000000000000000e+03 -1.285037988054682501e-03 +8.700000000000000000e+03 6.363336733556934632e-04 +8.760000000000000000e+03 -1.428649067747755907e-03 +8.820000000000000000e+03 1.520044182143465150e-03 +8.880000000000000000e+03 -5.175573060114402324e-04 +8.940000000000000000e+03 -5.414253655544598587e-04 +9.000000000000000000e+03 -1.551532963276258670e-03 +9.060000000000000000e+03 4.521457831287989393e-04 +9.120000000000000000e+03 4.696356527347234078e-04 +9.180000000000000000e+03 -1.499039731243101414e-03 +9.240000000000000000e+03 5.461421078507555649e-04 +9.300000000000000000e+03 1.605202551218098961e-03 +9.360000000000000000e+03 -3.218382062186719850e-04 +9.420000000000000000e+03 -1.234961082445806824e-03 +9.480000000000000000e+03 3.865851829687017016e-03 +9.540000000000000000e+03 9.806172965909354389e-04 +9.600000000000000000e+03 1.093509763450128958e-04 +9.660000000000000000e+03 1.252067327186523471e-03 +9.720000000000000000e+03 -5.912203314437647350e-04 +9.780000000000000000e+03 -1.420499824234866537e-03 +9.840000000000000000e+03 7.642398677489836700e-04 +9.900000000000000000e+03 -2.036991420027334243e-03 +9.960000000000000000e+03 1.758150492605636828e-04 +1.002000000000000000e+04 -1.597333203790185507e-03 +1.008000000000000000e+04 6.435701989175868221e-04 +1.014000000000000000e+04 -1.014695044432301074e-04 +1.020000000000000000e+04 1.675517150943051092e-04 +1.026000000000000000e+04 1.450636794288584497e-03 +1.032000000000000000e+04 7.477874460164457560e-04 +1.038000000000000000e+04 5.900423093407880515e-05 +1.044000000000000000e+04 1.384286575557780452e-03 +1.050000000000000000e+04 -2.763673010122147389e-04 +1.056000000000000000e+04 1.077039683877956122e-03 +1.062000000000000000e+04 -2.555496547756774817e-03 +1.068000000000000000e+04 -1.739812651067040861e-04 +1.074000000000000000e+04 -1.778420899427146651e-03 +1.080000000000000000e+04 -1.368823012853681576e-03 +1.086000000000000000e+04 -2.945196396467508748e-03 +1.092000000000000000e+04 -2.507550935661129188e-03 +1.098000000000000000e+04 -4.055897722537338268e-03 +1.104000000000000000e+04 -3.590249025364755653e-03 +1.110000000000000000e+04 -5.110618241815245710e-03 +1.116000000000000000e+04 -3.617019967350643128e-03 +1.122000000000000000e+04 -1.109469934817752801e-03 +1.128000000000000000e+04 5.412014954345067963e-03 +1.134000000000000000e+04 1.694741662504384294e-02 +1.140000000000000000e+04 3.649671583661984187e-02 +1.146000000000000000e+04 6.105989220850460697e-02 +1.152000000000000000e+04 9.763692419801373035e-02 +1.152000000000000000e+04 1.126369241983411484e-01 +1.158000000000000000e+04 1.412277890731274965e-01 +1.158000000000000000e+04 1.242277890733021195e-01 +1.158000000000000000e+04 1.352277890728146303e-01 +1.158000000000000000e+04 1.462277890732366359e-01 +1.158000000000000000e+04 1.582277890729528735e-01 +1.164000000000000000e+04 1.838324629907219787e-01 +1.164000000000000000e+04 1.678324629911003285e-01 +1.164000000000000000e+04 1.788324629906128393e-01 +1.164000000000000000e+04 1.898324629910348449e-01 +1.164000000000000000e+04 1.978324629908456700e-01 +1.170000000000000000e+04 2.184509209128009388e-01 +1.170000000000000000e+04 2.074509209123789333e-01 +1.170000000000000000e+04 2.164509209123934852e-01 +1.170000000000000000e+04 2.224509209127063514e-01 +1.170000000000000000e+04 2.274509209128154907e-01 +1.176000000000000000e+04 2.340831366736892960e-01 +1.176000000000000000e+04 2.340831366736892960e-01 +1.176000000000000000e+04 2.350831366738930228e-01 +1.176000000000000000e+04 2.340831366736892960e-01 +1.182000000000000000e+04 2.277290829279081663e-01 +1.188000000000000000e+04 2.053887312031292822e-01 +1.194000000000000000e+04 1.670620518771102070e-01 +1.200000000000000000e+04 1.247490141586240497e-01 +1.206000000000000000e+04 8.844958610643516295e-02 +1.212000000000000000e+04 6.016373466354707489e-02 +1.218000000000000000e+04 4.289142562447523233e-02 +1.224000000000000000e+04 3.963262360684893792e-02 +1.230000000000000000e+04 4.838729214316117577e-02 +1.236000000000000000e+04 6.315539358092792099e-02 +1.242000000000000000e+04 6.393688916250539478e-02 +1.248000000000000000e+04 5.273173900241090450e-02 +1.254000000000000000e+04 3.953990205445734318e-02 +1.260000000000000000e+04 4.936133618139137980e-02 +1.266000000000000000e+04 7.119599809902865672e-02 +1.272000000000000000e+04 8.704384339580428787e-02 +1.278000000000000000e+04 8.190482657028042013e-02 +1.284000000000000000e+04 7.377890098814532394e-02 +1.290000000000000000e+04 7.866601889145385940e-02 +1.296000000000000000e+04 7.756613143556023715e-02 +1.302000000000000000e+04 5.047918865693645785e-02 +1.308000000000000000e+04 2.240513946435385151e-02 +1.314000000000000000e+04 2.834393172361160396e-02 +1.320000000000000000e+04 4.229551215485116700e-02 +1.326000000000000000e+04 1.925982640750589781e-02 +1.332000000000000000e+04 -4.076318097031617071e-02 +1.338000000000000000e+04 -5.877356652763410239e-02 +1.344000000000000000e+04 -2.677138784474664135e-02 +1.350000000000000000e+04 3.243296360778913368e-03 +1.356000000000000000e+04 -2.729573684518982191e-03 +1.362000000000000000e+04 -6.690058808089816011e-03 +1.368000000000000000e+04 5.361779075428785291e-03 +1.374000000000000000e+04 4.258769822627073154e-04 +1.380000000000000000e+04 -1.849782908720953856e-02 +1.386000000000000000e+04 -1.540940423092251876e-02 +1.392000000000000000e+04 -1.430891459585836856e-02 +1.398000000000000000e+04 -5.219642733663931722e-02 +1.404000000000000000e+04 -9.707201070341398008e-02 +1.410000000000000000e+04 -9.493573392137477640e-02 +1.416000000000000000e+04 -6.078766729024209781e-02 +1.422000000000000000e+04 -5.562788214047031943e-02 +1.428000000000000000e+04 -7.745645079376117792e-02 +1.434000000000000000e+04 -9.227344661758252187e-02 +1.440000000000000000e+04 -9.207894397513882723e-02 +1.446000000000000000e+04 -6.987301823392044753e-02 +1.452000000000000000e+04 -3.065574578158702934e-02 +1.458000000000000000e+04 1.572795998981746379e-03 +1.464000000000000000e+04 6.812528757109248545e-03 +1.470000000000000000e+04 -7.936626887385500595e-03 +1.476000000000000000e+04 -2.867475123457552399e-02 +1.482000000000000000e+04 -5.340192561652656877e-02 +1.488000000000000000e+04 -4.111823231687594671e-02 +1.494000000000000000e+04 -7.823754569471930154e-03 +1.500000000000000000e+04 -2.651857659111556131e-02 +1.506000000000000000e+04 -6.820278353552566841e-02 +1.512000000000000000e+04 -4.687646149523061467e-02 +1.518000000000000000e+04 -7.539697528045508079e-03 +1.524000000000000000e+04 -3.019257961841503857e-02 +1.530000000000000000e+04 -8.383519667404470965e-02 +1.536000000000000000e+04 -8.346763851477589924e-02 +1.542000000000000000e+04 -3.208999590879102470e-02 +1.548000000000000000e+04 -1.670236050085804891e-02 +1.554000000000000000e+04 -2.730482487368135480e-02 +1.560000000000000000e+04 -1.689748247281386284e-02 +1.566000000000000000e+04 -4.748042766823346028e-02 +1.572000000000000000e+04 -8.205375568923045648e-02 +1.578000000000000000e+04 -5.361756264483119594e-02 +1.584000000000000000e+04 -5.717194553380977595e-02 +1.590000000000000000e+04 -8.471700221616629278e-02 +1.596000000000000000e+04 -4.525283139446401037e-02 +1.602000000000000000e+04 -2.877953262031951454e-02 +1.608000000000000000e+04 -5.029720631955569843e-02 +1.614000000000000000e+04 -3.680595371952222195e-02 +1.620000000000000000e+04 -3.630587690622633090e-02 +1.626000000000000000e+04 -3.579707876906468300e-02 +1.632000000000000000e+04 -2.427966301002015825e-02 +1.638000000000000000e+04 -1.875373416987713426e-02 +1.644000000000000000e+04 3.780602452025050297e-03 +1.650000000000000000e+04 2.432324073379277252e-02 +1.656000000000000000e+04 2.287407375752081862e-02 +1.662000000000000000e+04 3.243299386213038815e-02 +1.668000000000000000e+04 5.699989260210713837e-02 +1.674000000000000000e+04 5.257466074272088008e-02 +1.680000000000000000e+04 6.515718833361461293e-02 +1.686000000000000000e+04 8.174736465207388392e-02 +1.692000000000000000e+04 3.634507821607257938e-02 +1.698000000000000000e+04 1.995021687253029086e-02 +1.704000000000000000e+04 1.856266768299974501e-02 +1.710000000000000000e+04 1.118231702639604919e-02 +1.716000000000000000e+04 1.480905055905168410e-02 +1.722000000000000000e+04 -1.555724677564285230e-02 +1.728000000000000000e+04 -1.491669068036571844e-02 +1.734000000000000000e+04 -5.269397586744162254e-03 +1.740000000000000000e+04 -5.615484597001341172e-03 +1.746000000000000000e+04 1.704493054694466991e-02 +1.752000000000000000e+04 2.711729415750596672e-03 +1.758000000000000000e+04 1.638479292432748480e-02 +1.764000000000000000e+04 2.106400139109609881e-02 +1.770000000000000000e+04 5.749234454015095253e-03 +1.776000000000000000e+04 2.044037113500962732e-02 +1.782000000000000000e+04 2.413728985720808851e-02 +1.788000000000000000e+04 4.183986840871511959e-02 +1.794000000000000000e+04 4.354798401618609205e-02 +1.800000000000000000e+04 3.226151328817650210e-02 +1.806000000000000000e+04 4.798033223687525606e-02 +1.812000000000000000e+04 3.870431634459237102e-02 +1.818000000000000000e+04 4.143334049058466917e-02 +1.824000000000000000e+04 5.816727901947160717e-02 +1.830000000000000000e+04 4.190600573383562732e-02 +1.836000000000000000e+04 3.064939388423226774e-02 +1.842000000000000000e+04 1.239731618898076704e-02 +1.848000000000000000e+04 1.914964487514225766e-02 +1.854000000000000000e+04 3.790625161582283909e-02 +1.860000000000000000e+04 3.166700764359120512e-02 +1.866000000000000000e+04 2.443178365047060652e-02 +1.872000000000000000e+04 -6.799550151299627032e-03 +1.878000000000000000e+04 -5.027123981562908739e-03 +1.884000000000000000e+04 -1.125106858125946019e-02 +1.890000000000000000e+04 -2.447151510114053963e-02 +1.896000000000000000e+04 -4.168859514084033435e-02 +1.902000000000000000e+04 -6.390244075919326860e-02 +1.908000000000000000e+04 -6.211318443638447206e-02 +1.914000000000000000e+04 -7.232095904328161851e-02 +1.920000000000000000e+04 -4.952589787717442960e-02 +1.926000000000000000e+04 -4.072813463881175267e-02 +1.932000000000000000e+04 -2.592780339364253450e-02 +1.938000000000000000e+04 -5.125038614096411038e-03 +1.944000000000000000e+04 -7.319975109567167237e-03 +1.950000000000000000e+04 1.648725191898847697e-02 +1.956000000000000000e+04 1.129650694929296151e-02 +1.962000000000000000e+04 1.210765411542524816e-02 +1.968000000000000000e+04 -2.007944279012008337e-02 +1.974000000000000000e+04 -3.526492029777728021e-02 +1.980000000000000000e+04 -1.244891522219404578e-02 +1.986000000000000000e+04 -2.631564692819665652e-03 +1.992000000000000000e+04 9.186993857838388067e-03 +1.998000000000000000e+04 -6.993377250182675198e-03 +2.004000000000000000e+04 -8.172815996658755466e-03 +2.010000000000000000e+04 -1.935146062351122964e-02 +2.016000000000000000e+04 -2.452944957985891961e-02 +2.022000000000000000e+04 -3.770692160469479859e-02 +2.028000000000000000e+04 -1.088401562992658000e-02 +2.034000000000000000e+04 -2.606087083404418081e-02 +2.040000000000000000e+04 -3.223762659854401136e-02 +2.046000000000000000e+04 -2.741442249680403620e-02 +2.052000000000000000e+04 -4.459139827849867288e-02 +2.058000000000000000e+04 -3.776869391185755376e-02 +2.064000000000000000e+04 -2.794644949426583480e-02 +2.070000000000000000e+04 -4.112480530147877289e-02 +2.076000000000000000e+04 -3.730390177952358499e-02 +2.082000000000000000e+04 -2.248387944655405590e-02 +2.088000000000000000e+04 -1.066487902426160872e-02 +2.094000000000000000e+04 -8.470412976748775691e-04 +2.100000000000000000e+04 -1.403050718909071293e-02 +2.106000000000000000e+04 -1.621541772783530178e-02 +2.112000000000000000e+04 -2.640191397767921444e-02 +2.118000000000000000e+04 -3.159013713229796849e-02 +2.124000000000000000e+04 -7.780228442243242171e-03 +2.130000000000000000e+04 -4.972329188603907824e-03 +2.136000000000000000e+04 -9.166580726741813123e-03 +2.142000000000000000e+04 1.863687555305659771e-02 +2.148000000000000000e+04 1.143789827165164752e-02 +2.154000000000000000e+04 1.323634599066281226e-02 +2.160000000000000000e+04 1.203207728212873917e-02 +2.166000000000000000e+04 3.782495074301550630e-02 +2.172000000000000000e+04 2.061482493900257396e-02 +2.178000000000000000e+04 8.401558531659247819e-03 +2.184000000000000000e+04 9.185010166220308747e-03 +2.190000000000000000e+04 -1.603496146799443522e-02 +2.196000000000000000e+04 -5.258497583781718276e-03 +2.202000000000000000e+04 3.514260632982768584e-03 +2.208000000000000000e+04 2.831720994436182082e-04 +2.214000000000000000e+04 2.048095826467033476e-03 +2.220000000000000000e+04 8.088908980425912887e-04 +2.226000000000000000e+04 3.565416512174124364e-03 +2.232000000000000000e+04 -1.268246798463223968e-02 +2.238000000000000000e+04 -3.934903162189584691e-03 +2.244000000000000000e+04 -5.192029392674157862e-03 +2.250000000000000000e+04 -1.745398694583855104e-02 +2.256000000000000000e+04 -3.720915920894185547e-03 +2.262000000000000000e+04 1.600704378051887034e-02 +2.268000000000000000e+04 3.072975239592778962e-02 +2.274000000000000000e+04 1.444707038808701327e-02 +2.280000000000000000e+04 1.515885842218267499e-02 +2.286000000000000000e+04 2.386497735369630391e-02 +2.292000000000000000e+04 3.756528826579597080e-02 +2.298000000000000000e+04 3.325965249132423196e-02 +2.304000000000000000e+04 2.194793159924302017e-02 +2.310000000000000000e+04 1.162998739073373144e-02 +2.316000000000000000e+04 1.330568197317916201e-02 +2.322000000000000000e+04 5.974877678454504348e-03 +2.328000000000000000e+04 -2.362562831876857672e-03 +2.334000000000000000e+04 -1.370677664272079710e-02 +2.340000000000000000e+04 -1.505790050669020275e-02 +2.346000000000000000e+04 -1.141607085992291104e-02 +2.352000000000000000e+04 -9.781423831555002835e-03 +2.358000000000000000e+04 -1.154095220954332035e-03 +2.364000000000000000e+04 3.465779540420044214e-03 +2.370000000000000000e+04 -1.992193462774594082e-02 +2.376000000000000000e+04 -1.631737246407283237e-02 +2.382000000000000000e+04 -3.772066830151743488e-02 +2.388000000000000000e+04 -2.113195611127594020e-02 +2.394000000000000000e+04 -3.355136947357095778e-02 +2.400000000000000000e+04 -2.997904152107366826e-02 +2.406000000000000000e+04 -1.941510502911114600e-02 +2.412000000000000000e+04 3.140307697321986780e-03 +2.418000000000000000e+04 -1.631293523769272724e-02 +2.424000000000000000e+04 -2.377496529425116023e-02 +2.430000000000000000e+04 -2.224591344383952674e-02 +2.436000000000000000e+04 -1.772591021290281788e-02 +2.442000000000000000e+04 -1.021508566827833420e-02 +2.448000000000000000e+04 -1.371356936579104513e-02 +2.454000000000000000e+04 -2.221490385636570863e-03 +2.460000000000000000e+04 -1.173897732041950803e-02 +2.466000000000000000e+04 -1.826615819754806580e-02 +2.472000000000000000e+04 -2.380316057951858966e-02 +2.478000000000000000e+04 -2.535011146301258123e-02 +2.484000000000000000e+04 -2.790713731337746140e-02 +2.490000000000000000e+04 -2.047436406974156853e-02 +2.496000000000000000e+04 -2.805191705101606203e-02 +2.502000000000000000e+04 -3.363992106005753158e-02 +2.508000000000000000e+04 -3.223850031372421654e-02 +2.514000000000000000e+04 -2.884777841154573252e-02 +2.520000000000000000e+04 -2.546787838491582079e-02 +2.526000000000000000e+04 -2.009892265959933866e-02 +2.532000000000000000e+04 -1.874103300906426739e-02 +2.538000000000000000e+04 -1.339433061548334081e-02 +2.544000000000000000e+04 -1.605893602663854836e-02 +2.550000000000000000e+04 -3.473496912101836642e-02 +2.556000000000000000e+04 -4.042254916384990793e-02 +2.562000000000000000e+04 -3.712179470676346682e-02 +2.568000000000000000e+04 -1.483282366734783864e-02 +2.574000000000000000e+04 -1.855575329045677790e-02 +2.580000000000000000e+04 -2.829070010102441302e-02 +2.586000000000000000e+04 -1.903777995266864309e-02 +2.592000000000000000e+04 -8.797107981081353500e-03 +2.598000000000000000e+04 7.431201404870080296e-03 +2.604000000000000000e+04 6.647034491834347136e-03 +2.610000000000000000e+04 -4.149721689827856608e-03 +2.616000000000000000e+04 1.904082064902468119e-02 +2.622000000000000000e+04 3.921855000953655690e-02 +2.628000000000000000e+04 3.238335570040362654e-02 +2.634000000000000000e+04 1.853512776051502442e-02 +2.640000000000000000e+04 1.267375700354023138e-02 +2.646000000000000000e+04 4.799135052053316031e-03 +2.652000000000000000e+04 -1.208884571497037541e-02 +2.658000000000000000e+04 -1.599029208591673523e-02 +2.664000000000000000e+04 -1.490531005038064905e-02 +2.670000000000000000e+04 -8.340047934325411916e-04 +2.676000000000000000e+04 1.022351932624587789e-02 +2.682000000000000000e+04 1.126715879308903823e-02 +2.688000000000000000e+04 -5.703189061023294926e-03 +2.694000000000000000e+04 -2.068762609451368917e-02 +2.700000000000000000e+04 -1.568625325762695866e-02 +2.706000000000000000e+04 -1.699170679785311222e-03 +2.712000000000000000e+04 1.227352240766776958e-02 +2.718000000000000000e+04 1.123172763163893251e-02 +2.724000000000000000e+04 1.217534752959181787e-02 +2.730000000000000000e+04 1.810428552289522486e-02 +2.736000000000000000e+04 2.301844593512214487e-02 +2.742000000000000000e+04 2.091773398024088237e-02 +2.748000000000000000e+04 2.080205582024063915e-02 +2.754000000000000000e+04 1.967131853325554403e-02 +2.760000000000000000e+04 1.552543011439411202e-02 +2.766000000000000000e+04 1.836429951708851149e-02 +2.772000000000000000e+04 1.718783662636269582e-02 +2.778000000000000000e+04 1.999595226880046539e-02 +2.784000000000000000e+04 1.978855827383085852e-02 +2.790000000000000000e+04 2.056556739262305200e-02 +2.796000000000000000e+04 1.932689338173076976e-02 +2.802000000000000000e+04 1.107245095590769779e-02 +2.808000000000000000e+04 1.680215581836819183e-02 +2.814000000000000000e+04 1.751592468281160109e-02 +2.820000000000000000e+04 1.021367526664107572e-02 +2.826000000000000000e+04 -1.104673739973804913e-03 +2.832000000000000000e+04 1.556079741112625925e-02 +2.838000000000000000e+04 7.210009453046950512e-03 +2.844000000000000000e+04 9.842884147474251222e-03 +2.850000000000000000e+04 1.845934433094953420e-02 +2.856000000000000000e+04 1.105931382244307315e-02 +2.862000000000000000e+04 1.864271752947388450e-02 +2.868000000000000000e+04 1.920948138376843417e-02 +2.874000000000000000e+04 2.975953236455097795e-02 +2.880000000000000000e+04 1.729279856863286113e-02 +2.886000000000000000e+04 1.880920910480199382e-02 +2.892000000000000000e+04 2.430869417275971500e-02 +2.898000000000000000e+04 1.979118506915256148e-02 +2.904000000000000000e+04 1.725661415821377886e-02 +2.910000000000000000e+04 1.470491487998515368e-02 +2.916000000000000000e+04 4.136021831982361618e-03 +2.922000000000000000e+04 7.549870650109369308e-03 +2.928000000000000000e+04 1.094639809798536589e-02 +2.934000000000000000e+04 2.032554207198700169e-02 +2.940000000000000000e+04 2.768724154520896263e-02 +2.946000000000000000e+04 1.803143666893447516e-02 +2.952000000000000000e+04 1.635806867852807045e-02 +2.958000000000000000e+04 1.366707994384341873e-02 +2.964000000000000000e+04 1.695841400578501634e-02 +2.970000000000000000e+04 1.323201550439989660e-02 +2.976000000000000000e+04 9.487830260695773177e-03 +2.982000000000000000e+04 8.725805244466755539e-03 +2.988000000000000000e+04 1.794588855864276411e-02 +2.994000000000000000e+04 1.514802947622229112e-02 +3.000000000000000000e+04 1.933217846089974046e-02 +3.006000000000000000e+04 3.049828710481961025e-02 +3.012000000000000000e+04 2.364630819738522405e-02 +3.018000000000000000e+04 1.877619571860122960e-02 +3.024000000000000000e+04 1.888790479551971657e-02 +3.030000000000000000e+04 1.898139176955737639e-02 +3.036000000000000000e+04 1.305661414971837075e-02 +3.042000000000000000e+04 -1.886469341116026044e-03 +3.048000000000000000e+04 -9.847898781117692124e-03 +3.054000000000000000e+04 -1.482771307473740308e-02 +3.060000000000000000e+04 -1.182594990495999809e-02 +3.066000000000000000e+04 -7.842645744858600665e-03 +3.072000000000000000e+04 -7.877835847466485575e-03 +3.078000000000000000e+04 -1.093155428588943323e-02 +3.084000000000000000e+04 -9.003833871247479692e-03 +3.090000000000000000e+04 -8.094706232441239990e-03 +3.096000000000000000e+04 -1.720420174206083175e-02 +3.102000000000000000e+04 -2.533234957991226111e-02 +3.108000000000000000e+04 -3.347917766677710460e-02 +3.114000000000000000e+04 -3.964471269591740565e-02 +3.120000000000000000e+04 -3.382898013933299808e-02 +3.126000000000000000e+04 -4.203200423307862366e-02 +3.132000000000000000e+04 -3.825380795296950964e-02 +3.138000000000000000e+04 -2.949441303553612670e-02 +3.144000000000000000e+04 -2.975383999728364870e-02 +3.150000000000000000e+04 -3.703210807270806981e-02 +3.156000000000000000e+04 -4.932923527212551562e-02 +3.162000000000000000e+04 -4.864523836204170948e-02 +3.168000000000000000e+04 -4.398013280842860695e-02 +3.174000000000000000e+04 -4.133393287702347152e-02 +3.180000000000000000e+04 -4.470665153166919481e-02 +3.186000000000000000e+04 -4.509830051028984599e-02 +3.192000000000000000e+04 -4.150889028824167326e-02 +3.198000000000000000e+04 -4.193843005123198964e-02 +3.204000000000000000e+04 -4.738692773480579490e-02 +3.210000000000000000e+04 -4.985439003121427959e-02 +3.216000000000000000e+04 -5.434082232659420697e-02 +3.222000000000000000e+04 -5.884622877238143701e-02 +3.228000000000000000e+04 -4.437061225326033309e-02 +3.234000000000000000e+04 -3.891397433835663833e-02 +3.240000000000000000e+04 -3.347631538326822920e-02 +3.246000000000000000e+04 -3.405763442970055621e-02 +3.252000000000000000e+04 -3.465792927454458550e-02 +3.258000000000000000e+04 -3.627719644646276720e-02 +3.264000000000000000e+04 -2.691543115270178532e-02 +3.270000000000000000e+04 -2.457262737061682856e-02 +3.276000000000000000e+04 -7.248777810673345812e-03 +3.282000000000000000e+04 -2.943873865660862066e-03 +3.288000000000000000e+04 -1.065790569191449322e-02 +3.294000000000000000e+04 -8.390862165470025502e-03 +3.300000000000000000e+04 -1.714273084780870704e-02 +3.306000000000000000e+04 -2.291349807273945771e-02 +3.312000000000000000e+04 -2.070314888715074630e-02 +3.318000000000000000e+04 -2.051166704586648848e-02 +3.324000000000000000e+04 -7.339035063523624558e-03 +3.330000000000000000e+04 -1.018523412403737893e-02 +3.336000000000000000e+04 -7.050244191304955166e-03 +3.342000000000000000e+04 -4.934043945468147285e-03 +3.348000000000000000e+04 -3.836610764665238094e-03 +3.354000000000000000e+04 -9.757920794982055668e-03 +3.360000000000000000e+04 -2.697948890272527933e-03 +3.366000000000000000e+04 1.134333137633802835e-02 +3.372000000000000000e+04 4.365947673250047956e-03 +3.378000000000000000e+04 1.369928945678111631e-03 +3.384000000000000000e+04 -3.644694577815243974e-03 +3.390000000000000000e+04 -6.677891429717419669e-03 +3.396000000000000000e+04 -3.729628862856770866e-03 +3.402000000000000000e+04 -8.799872860436153132e-03 +3.408000000000000000e+04 -1.888588170913862996e-03 +3.414000000000000000e+04 -3.995738251433067489e-03 +3.420000000000000000e+04 -3.121285336419532541e-03 +3.426000000000000000e+04 -2.265190380967396777e-03 +3.432000000000000000e+04 -1.427413078090467025e-03 +3.438000000000000000e+04 2.392088101260014810e-03 +3.444000000000000000e+04 8.193355942239577416e-03 +3.450000000000000000e+04 4.976434497621085029e-03 +3.456000000000000000e+04 4.741369056318944786e-03 +3.462000000000000000e+04 6.488206130597973242e-03 +3.468000000000000000e+04 1.216993508023733739e-03 +3.474000000000000000e+04 -2.072219826914079022e-03 +3.480000000000000000e+04 4.620616387910558842e-03 +3.486000000000000000e+04 1.029555360219092108e-02 +3.492000000000000000e+04 9.526445292067364790e-04 +3.498000000000000000e+04 -7.408056894746550824e-03 +3.504000000000000000e+04 -1.786495526175713167e-03 +3.510000000000000000e+04 4.817385001842922065e-03 +3.516000000000000000e+04 7.403642284771194682e-03 +3.522000000000000000e+04 -2.027664866545819677e-03 +3.528000000000000000e+04 5.523523548617959023e-03 +3.534000000000000000e+04 5.057268765995104332e-03 +3.540000000000000000e+04 -3.426366790336032864e-03 +3.546000000000000000e+04 -7.927319519694719929e-03 +3.552000000000000000e+04 2.554475442593684420e-03 +3.558000000000000000e+04 4.019084080937318504e-03 +3.564000000000000000e+04 1.466573618927213829e-03 +3.570000000000000000e+04 3.897012439665559214e-03 +3.576000000000000000e+04 3.310470098767837044e-03 +3.582000000000000000e+04 1.370701736323098885e-02 +3.588000000000000000e+04 1.308672615505201975e-02 +3.594000000000000000e+04 1.144966954780102242e-02 +3.600000000000000000e+04 1.379592181274347240e-02 +3.606000000000000000e+04 1.412555835304374341e-02 +3.612000000000000000e+04 1.843865572391223395e-02 +3.618000000000000000e+04 1.073529167751985369e-02 +3.624000000000000000e+04 1.601554505123203853e-02 +3.630000000000000000e+04 1.727949587257171515e-02 +3.636000000000000000e+04 1.452722527119476581e-02 +3.642000000000000000e+04 1.375881551302882144e-02 +3.648000000000000000e+04 1.974350024283921812e-03 +3.654000000000000000e+04 -5.826086680826847441e-03 +3.660000000000000000e+04 -6.642408977313607465e-03 +3.666000000000000000e+04 -4.474530098377726972e-03 +3.672000000000000000e+04 -8.322362190483545419e-03 +3.678000000000000000e+04 -1.218581631837878376e-02 +3.684000000000000000e+04 -9.064802396096638404e-03 +3.690000000000000000e+04 -9.959229301784944255e-03 +3.696000000000000000e+04 -7.869004814892832655e-03 +3.702000000000000000e+04 -1.879403561906656250e-02 +3.708000000000000000e+04 -1.673422736348584294e-02 +3.714000000000000000e+04 -7.689484577895200346e-03 +3.720000000000000000e+04 -6.659710769781668205e-03 +3.726000000000000000e+04 -1.864480839776661014e-02 +3.732000000000000000e+04 -2.764467883753241040e-02 +3.738000000000000000e+04 -2.465922245937690604e-02 +3.744000000000000000e+04 -2.368833857235586038e-02 +3.750000000000000000e+04 -1.673192545604251791e-02 +3.756000000000000000e+04 -2.678988038860552479e-02 +3.762000000000000000e+04 -3.386209961718122941e-02 +3.768000000000000000e+04 -3.494847836736880708e-02 +3.774000000000000000e+04 -2.404891089190641651e-02 +3.780000000000000000e+04 -1.416329041876451811e-02 +3.786000000000000000e+04 -1.429150917829247192e-02 +3.792000000000000000e+04 -9.433458459170651622e-03 +3.798000000000000000e+04 -1.758902853089239215e-02 +3.804000000000000000e+04 -2.375810872217698488e-02 +3.810000000000000000e+04 -3.494058738579042256e-02 +3.816000000000000000e+04 -2.613635190755303483e-02 +3.822000000000000000e+04 -1.834528875133401016e-02 +3.828000000000000000e+04 -1.756728343116265023e-02 +3.834000000000000000e+04 -2.880222051771852421e-02 +3.840000000000000000e+04 -3.804998365194478538e-02 +3.846000000000000000e+04 -3.431045558318146504e-02 +3.852000000000000000e+04 -2.658351811260217801e-02 +3.858000000000000000e+04 -1.686905217320600059e-02 +3.864000000000000000e+04 -1.016693776728061493e-02 +3.870000000000000000e+04 -1.947705404381849803e-02 +3.876000000000000000e+04 -1.379927924790536053e-02 +3.882000000000000000e+04 -7.133490744308801368e-03 +3.888000000000000000e+04 -6.479565056906722020e-03 +3.894000000000000000e+04 -1.783737784808181459e-02 +3.900000000000000000e+04 -2.120680391726637026e-02 +3.906000000000000000e+04 -1.658771722213714384e-02 +3.912000000000000000e+04 -2.497999091337987920e-02 +3.918000000000000000e+04 -3.038349728376488201e-02 +3.924000000000000000e+04 -4.179810782716231188e-02 +3.930000000000000000e+04 -3.322369325178442523e-02 +3.936000000000000000e+04 -9.660123405410558917e-03 +3.942000000000000000e+04 2.892732593863911461e-03 +3.948000000000000000e+04 6.435006451283697970e-03 +3.954000000000000000e+04 7.966830619807296898e-03 +3.960000000000000000e+04 9.488338321716582868e-03 +3.966000000000000000e+04 1.999663553760910872e-03 +3.972000000000000000e+04 -8.499058943016279954e-03 +3.978000000000000000e+04 -1.000769371876231162e-02 +3.984000000000000000e+04 -5.526104539967491291e-03 +3.990000000000000000e+04 1.945845498084963765e-03 +3.996000000000000000e+04 1.240829402468079934e-02 +4.002000000000000000e+04 1.086137941183551447e-02 +4.008000000000000000e+04 1.430524067836813629e-02 +4.014000000000000000e+04 1.274001754245546181e-02 +4.020000000000000000e+04 1.616585041847429238e-02 +4.026000000000000000e+04 1.658288034377619624e-02 +4.032000000000000000e+04 1.899124908140947809e-02 +4.038000000000000000e+04 1.739109898244350916e-02 +4.044000000000000000e+04 2.178257305240549613e-02 +4.050000000000000000e+04 1.816581495040736627e-02 +4.056000000000000000e+04 2.454096892233792460e-02 +4.062000000000000000e+04 2.590817985492321895e-02 +4.068000000000000000e+04 2.826759322397265350e-02 +4.074000000000000000e+04 2.061935509482282214e-02 +4.080000000000000000e+04 9.963612108549568802e-03 +4.086000000000000000e+04 1.300511502449808177e-03 +4.092000000000000000e+04 5.630201040730753448e-03 +4.098000000000000000e+04 5.952829093075706623e-03 +4.104000000000000000e+04 1.026854452720726840e-02 +4.110000000000000000e+04 1.457749675410013879e-02 +4.116000000000000000e+04 1.387983572749362793e-02 +4.122000000000000000e+04 3.175711888616206124e-03 +4.128000000000000000e+04 3.465276201495726127e-03 +4.134000000000000000e+04 2.748680130935099442e-03 +4.140000000000000000e+04 -5.973924398858798668e-03 +4.146000000000000000e+04 -4.702384977463225368e-03 +4.152000000000000000e+04 -8.436548727331683040e-03 +4.158000000000000000e+04 -2.217626230867608683e-02 +4.164000000000000000e+04 -2.692137196845578728e-02 +4.170000000000000000e+04 -1.867172347738232929e-02 +4.176000000000000000e+04 -1.142716222693707095e-02 +4.182000000000000000e+04 -1.318753315808862681e-02 +4.188000000000000000e+04 -1.095268083736300468e-02 +4.194000000000000000e+04 -9.722449420223711058e-03 +4.200000000000000000e+04 -1.349668267812376143e-02 +4.206000000000000000e+04 -2.027522402840986615e-02 +4.212000000000000000e+04 -2.005791652391053503e-02 +4.218000000000000000e+04 -1.884460285236855270e-02 +4.224000000000000000e+04 -7.635125359229277819e-03 +4.230000000000000000e+04 -8.429326088844391052e-03 +4.236000000000000000e+04 -1.322704672566032968e-02 +4.242000000000000000e+04 -2.302812868401815649e-02 +4.248000000000000000e+04 -2.583241305364936125e-02 +4.254000000000000000e+04 -1.463974063426576322e-02 +4.260000000000000000e+04 -3.449951977017917670e-03 +4.266000000000000000e+04 -3.262887330492958426e-03 +4.272000000000000000e+04 -9.078386710825725459e-03 +4.278000000000000000e+04 2.103710105984646361e-03 +4.284000000000000000e+04 9.283563601456990000e-03 +4.290000000000000000e+04 1.461334478335629683e-03 +4.296000000000000000e+04 4.637183626982732676e-03 +4.302000000000000000e+04 3.811272173152246978e-03 +4.308000000000000000e+04 1.198376141201151768e-02 +4.314000000000000000e+04 1.815481281300890259e-02 +4.320000000000000000e+04 2.332458805722126272e-02 +4.326000000000000000e+04 1.749324893717130180e-02 +4.332000000000000000e+04 1.666095743075857172e-02 +4.338000000000000000e+04 1.682787562458543107e-02 +4.344000000000000000e+04 1.899416573451162549e-02 +4.350000000000000000e+04 2.315999013080727309e-02 +4.356000000000000000e+04 2.632551124952442478e-02 +4.362000000000000000e+04 1.749089162331074476e-02 +4.368000000000000000e+04 7.656293892068788409e-03 +4.374000000000000000e+04 7.821880748451803811e-03 +4.380000000000000000e+04 4.987814936612267047e-03 +4.386000000000000000e+04 5.154259300070407335e-03 +4.392000000000000000e+04 -2.678623328392859548e-03 +4.398000000000000000e+04 4.489329928219376598e-03 +4.404000000000000000e+04 4.658281976844591554e-03 +4.410000000000000000e+04 1.828395705160801299e-03 +4.416000000000000000e+04 -1.659527697484008968e-07 +4.422000000000000000e+04 1.727598701108945534e-04 +4.428000000000000000e+04 -3.652663966931868345e-03 +4.434000000000000000e+04 -3.476274626336817164e-03 +4.440000000000000000e+04 -1.029790934353513876e-02 +4.446000000000000000e+04 -1.011740539706806885e-02 +4.452000000000000000e+04 -2.934600120170216542e-03 +4.458000000000000000e+04 3.250669044064125046e-03 +4.464000000000000000e+04 6.438564566451532301e-03 +4.470000000000000000e+04 5.629248843433742877e-03 +4.476000000000000000e+04 8.822884115033957642e-03 +4.482000000000000000e+04 1.301963255082227988e-02 +4.488000000000000000e+04 1.821965614453802118e-02 +4.494000000000000000e+04 1.542311674438678892e-02 +4.500000000000000000e+04 1.363017607491201488e-02 +4.506000000000000000e+04 1.384099565530050313e-02 +4.512000000000000000e+04 7.055736845359206200e-03 +4.518000000000000000e+04 4.274560831618146040e-03 +4.524000000000000000e+04 4.976285717930295505e-04 +4.530000000000000000e+04 -6.274899185882532038e-03 +4.536000000000000000e+04 -1.704286189215054037e-02 +4.542000000000000000e+04 -1.080609923337760847e-02 +4.548000000000000000e+04 2.435548828543687705e-03 +4.554000000000000000e+04 2.682242124137701467e-03 +4.560000000000000000e+04 1.934140175762877334e-03 +4.566000000000000000e+04 1.119140227365278406e-02 +4.572000000000000000e+04 1.445418737785075791e-02 +4.578000000000000000e+04 1.172265418335882714e-02 +4.584000000000000000e+04 1.699696107607451268e-02 +4.590000000000000000e+04 1.527726610311219702e-02 +4.596000000000000000e+04 1.456372699067287613e-02 +4.602000000000000000e+04 5.856501134076097514e-03 +4.608000000000000000e+04 -1.844254414208990056e-03 +4.614000000000000000e+04 5.461617006403685082e-03 +4.620000000000000000e+04 5.774271720838441979e-03 +4.626000000000000000e+04 1.093865635994006880e-03 +4.632000000000000000e+04 -1.579445696734183002e-03 +4.638000000000000000e+04 3.754492854568525217e-03 +4.644000000000000000e+04 -2.904164003666664939e-03 +4.650000000000000000e+04 -8.555261966648686212e-03 +4.656000000000000000e+04 -1.198647174533107318e-03 +4.662000000000000000e+04 1.116583377461211057e-02 +4.668000000000000000e+04 -4.616661453837878071e-04 +4.674000000000000000e+04 -1.608099443637911463e-02 +4.680000000000000000e+04 -4.691999090027820785e-03 +4.686000000000000000e+04 -4.294528556783916429e-03 +4.692000000000000000e+04 -1.288843180554977152e-02 +4.698000000000000000e+04 -1.047355832361063221e-02 +4.704000000000000000e+04 -5.049758062341425102e-03 +4.710000000000000000e+04 -2.616881550238758791e-03 +4.716000000000000000e+04 -1.017477982259151759e-02 +4.722000000000000000e+04 -1.372330448793945834e-02 +4.728000000000000000e+04 -7.262307700329984073e-03 +4.734000000000000000e+04 5.208357848459854722e-03 +4.740000000000000000e+04 3.688838843117991928e-03 +4.746000000000000000e+04 -2.820718595103244297e-03 +4.752000000000000000e+04 -4.320168931371881627e-03 +4.758000000000000000e+04 -8.093672540780971758e-04 +4.764000000000000000e+04 1.711830717795237433e-03 +4.770000000000000000e+04 4.243568690071697347e-03 +4.776000000000000000e+04 7.785989691910799593e-03 +4.782000000000000000e+04 5.339236104191513732e-03 +4.788000000000000000e+04 1.903449684505176265e-03 +4.794000000000000000e+04 -2.521228512705420144e-03 +4.800000000000000000e+04 -5.934658059231878724e-03 +4.806000000000000000e+04 -1.133669925184221938e-02 +4.812000000000000000e+04 -7.727213076577754691e-03 +4.818000000000000000e+04 -1.060611884895479307e-04 +4.824000000000000000e+04 -1.473105996410595253e-03 +4.830000000000000000e+04 -3.828210607025539503e-03 +4.836000000000000000e+04 -4.171238845628977288e-03 +4.842000000000000000e+04 -1.050205530009407084e-02 +4.848000000000000000e+04 -1.382052531062072376e-02 +4.854000000000000000e+04 -8.126514946525276173e-03 +4.860000000000000000e+04 -1.041989108489360660e-02 +4.866000000000000000e+04 -8.700521320861298591e-03 +4.872000000000000000e+04 -1.296827408623357769e-02 +4.878000000000000000e+04 -2.022301859506114852e-02 +4.884000000000000000e+04 -8.464624836051370949e-03 +4.890000000000000000e+04 -6.692963652312755585e-03 +4.896000000000000000e+04 -2.907906674408877734e-03 +4.902000000000000000e+04 -3.109326362391584553e-03 +4.908000000000000000e+04 -6.297096034359128680e-03 +4.914000000000000000e+04 -1.247108984898659401e-02 +4.920000000000000000e+04 -7.631182806107972283e-03 +4.926000000000000000e+04 1.222749233420472592e-03 +4.932000000000000000e+04 3.090829511165793519e-03 +4.938000000000000000e+04 -2.681957812455948442e-05 +4.944000000000000000e+04 1.869923443337029312e-03 +4.950000000000000000e+04 1.781179184945358429e-03 +4.956000000000000000e+04 1.070706732298276620e-02 +4.962000000000000000e+04 1.164770664036041126e-02 +4.968000000000000000e+04 1.260321500467398437e-02 +4.974000000000000000e+04 1.557370933733182028e-02 +4.980000000000000000e+04 1.555930562517460203e-02 +4.986000000000000000e+04 1.256011892655806150e-02 +4.992000000000000000e+04 1.576263344759354368e-03 +4.998000000000000000e+04 7.607851998727710452e-03 +5.004000000000000000e+04 6.654997077021107543e-03 +5.010000000000000000e+04 7.717809770838357508e-03 +5.016000000000000000e+04 3.796400274040934164e-03 +5.022000000000000000e+04 -5.109122159410617314e-03 +5.028000000000000000e+04 1.350673301203642040e-06 +5.034000000000000000e+04 9.127926010478404351e-03 +5.040000000000000000e+04 1.127071005521429470e-02 +5.046000000000000000e+04 1.542980798603821313e-02 +5.052000000000000000e+04 1.260532397827773821e-02 +5.058000000000000000e+04 8.797361168035422452e-03 +5.064000000000000000e+04 6.006021634675562382e-03 +5.070000000000000000e+04 2.314064395250170492e-04 +5.076000000000000000e+04 1.473615561735641677e-03 +5.082000000000000000e+04 -2.672520849955617450e-04 +5.088000000000000000e+04 -2.991098605889419559e-03 +5.094000000000000000e+04 -6.978272340347757563e-04 +5.100000000000000000e+04 4.612657750840298831e-03 +5.106000000000000000e+04 1.094045096488116542e-02 +5.112000000000000000e+04 4.285645921299874317e-03 +5.118000000000000000e+04 -3.516649239827529527e-04 +5.124000000000000000e+04 -4.971390258106112015e-03 +5.130000000000000000e+04 4.265601182851241902e-04 +5.136000000000000000e+04 7.842275310395052657e-03 +5.142000000000000000e+04 1.127584328241937328e-02 +5.148000000000000000e+04 6.727350856635894161e-03 +5.154000000000000000e+04 3.196883767486724537e-03 +5.160000000000000000e+04 -2.315473439921333920e-03 +5.166000000000000000e+04 -1.180963735896511935e-02 +5.172000000000000000e+04 -1.928552571735053789e-02 +5.178000000000000000e+04 -2.074305743371951394e-02 +5.184000000000000000e+04 -9.182152538414811715e-03 +5.190000000000000000e+04 -4.602732282364740968e-03 +5.196000000000000000e+04 -8.004719073142041452e-03 +5.202000000000000000e+04 -1.738803648458997486e-02 +5.208000000000000000e+04 -1.675260930369404377e-02 +5.214000000000000000e+04 -1.509836348213866586e-02 +5.220000000000000000e+04 -2.242522616597852902e-02 +5.226000000000000000e+04 -2.573312572530994657e-02 +5.232000000000000000e+04 -2.402199172593100229e-02 +5.238000000000000000e+04 -1.829175492457579821e-02 +5.244000000000000000e+04 -1.354234732298209565e-02 +5.250000000000000000e+04 -1.377370209684158908e-02 +5.256000000000000000e+04 -1.198575369016907644e-02 +5.262000000000000000e+04 -1.117843776410154533e-02 +5.268000000000000000e+04 -5.351691186660900712e-03 +5.274000000000000000e+04 -8.505452098688692786e-03 +5.280000000000000000e+04 -1.263965985344839282e-02 +5.286000000000000000e+04 -9.754255049301718827e-03 +5.292000000000000000e+04 -2.849179548320535105e-03 +5.298000000000000000e+04 1.075623528777214233e-03 +5.304000000000000000e+04 -9.797901639103656635e-04 +5.310000000000000000e+04 -8.015366249310318381e-03 +5.316000000000000000e+04 9.689483731563086621e-04 +5.322000000000000000e+04 5.973205563350347802e-03 +5.328000000000000000e+04 7.997455881195492111e-03 +5.334000000000000000e+04 -1.958251347787154373e-03 +5.340000000000000000e+04 -9.893868116705561988e-03 +5.346000000000000000e+04 -1.080934767196595203e-02 +5.352000000000000000e+04 -1.270464454864850268e-02 +5.358000000000000000e+04 -1.257971457107487367e-02 +5.364000000000000000e+04 -1.343451486809499329e-02 +5.370000000000000000e+04 -1.026900382748863194e-02 +5.376000000000000000e+04 -7.083141152179450728e-03 +5.382000000000000000e+04 -6.876887844555312768e-03 +5.388000000000000000e+04 -4.650206195037753787e-03 +5.394000000000000000e+04 -5.403059804848453496e-03 +5.400000000000000000e+04 -3.135413596282887738e-03 +5.406000000000000000e+04 7.152766245781094767e-03 +5.412000000000000000e+04 4.461512179659621324e-03 +5.418000000000000000e+04 3.790855387705960311e-03 +5.424000000000000000e+04 5.140825714988750406e-03 +5.430000000000000000e+04 2.511451699319877662e-03 +5.436000000000000000e+04 1.290276055897265906e-02 +5.442000000000000000e+04 1.231477820965665160e-02 +5.448000000000000000e+04 1.747529214298992883e-03 +5.454000000000000000e+04 -2.798963152599753812e-03 +5.460000000000000000e+04 -4.324676961005025078e-03 +5.466000000000000000e+04 2.170408378333377186e-03 +5.472000000000000000e+04 7.686312186706345528e-03 +5.478000000000000000e+04 9.223052389643271454e-03 +5.484000000000000000e+04 7.780645634738903027e-03 +5.490000000000000000e+04 6.359107210300862789e-03 +5.496000000000000000e+04 8.958451065154804382e-03 +5.502000000000000000e+04 1.157868986592802685e-02 +5.508000000000000000e+04 1.021983490136335604e-02 +5.514000000000000000e+04 5.881896128812513780e-03 +5.520000000000000000e+04 8.564882206883339677e-03 +5.526000000000000000e+04 1.526880042365519330e-02 +5.532000000000000000e+04 1.099365673599095317e-02 +5.538000000000000000e+04 4.739455801427538972e-03 +5.544000000000000000e+04 -2.493799091098480858e-03 +5.550000000000000000e+04 -5.706105991521326359e-03 +5.646000000000000000e+04 -8.257657334979739971e-03 +5.736000000000000000e+04 -1.294846665223303717e-02 +5.826000000000000000e+04 -7.074896546328091063e-03 +5.916000000000000000e+04 4.193065290564845782e-03 +6.006000000000000000e+04 1.462318241283355746e-02 +6.096000000000000000e+04 2.924670918218907900e-03 +6.186000000000000000e+04 7.528095175075577572e-04 +6.276000000000000000e+04 1.071438542385294568e-02 +6.366000000000000000e+04 -1.262610964022314874e-02 +6.456000000000000000e+04 -5.739631480537354946e-03 +6.546000000000000000e+04 -9.125192001192772295e-03 +6.636000000000000000e+04 6.979586014494998381e-04 +6.726000000000000000e+04 1.019845200335112168e-02 +6.816000000000000000e+04 -2.158904887437529396e-03 +6.906000000000000000e+04 2.095205182740755845e-03 +6.996000000000000000e+04 -5.571460951614426449e-04 +7.086000000000000000e+04 -1.361311042728630127e-02 +7.176000000000000000e+04 1.145860972428636160e-02 +7.266000000000000000e+04 -1.077503794931544689e-02 +7.356000000000000000e+04 6.295159244473325089e-03 +7.446000000000000000e+04 2.326602399080002215e-03 +7.536000000000000000e+04 1.003001891240273835e-02 +7.626000000000000000e+04 1.517373321257764474e-02 +7.716000000000000000e+04 1.658700256120937411e-02 +7.806000000000000000e+04 -8.376347886951407418e-04 +7.896000000000000000e+04 1.685696384356560884e-02 +7.986000000000000000e+04 1.469282389734871686e-02 +8.076000000000000000e+04 7.560829544672742486e-04 +8.166000000000000000e+04 -3.480481571477866964e-02 +8.256000000000000000e+04 -7.819357142579974607e-04 +8.346000000000000000e+04 1.008837982499244390e-02 +8.436000000000000000e+04 6.120836686022812501e-03 +8.526000000000000000e+04 -2.324304773537733126e-03 +8.616000000000000000e+04 -8.474578135064803064e-04 +8.706000000000000000e+04 -1.016570184219744988e-03 +8.796000000000000000e+04 1.762549348859465681e-02 +8.886000000000000000e+04 8.553167517675319687e-03 +8.976000000000000000e+04 1.725012554379645735e-02 +9.066000000000000000e+04 2.201076116762124002e-03 +9.156000000000000000e+04 -4.116453088499838486e-03 +9.246000000000000000e+04 3.759790891308512073e-03 +9.336000000000000000e+04 -1.731053883304412011e-03 +9.426000000000000000e+04 -8.180423037629225291e-03 +9.516000000000000000e+04 -2.217281859884678852e-03 +9.606000000000000000e+04 -1.451446458759164670e-02 +9.696000000000000000e+04 -5.794293446342635434e-03 +9.786000000000000000e+04 7.166614861489506438e-03 +9.876000000000000000e+04 -3.466565482085570693e-03 +9.966000000000000000e+04 5.410172639130905736e-03 +1.005600000000000000e+05 2.837717569491360337e-03 +1.014600000000000000e+05 -9.207073198012949433e-03 +1.023600000000000000e+05 1.018877103342674673e-02 +1.032600000000000000e+05 8.755340377319953404e-04 +1.041600000000000000e+05 -6.356949212204199284e-03 +1.050600000000000000e+05 6.223949462764721829e-03 +1.059600000000000000e+05 2.978469774461700581e-04 +1.068600000000000000e+05 1.496431749728799332e-03 +1.077600000000000000e+05 6.409353354683844373e-03 +1.086600000000000000e+05 -3.409174560147221200e-03 +1.095600000000000000e+05 3.566836376194260083e-03 +1.104600000000000000e+05 -2.315714746509911492e-02 +1.113600000000000000e+05 -8.810956478555453941e-05 +1.122600000000000000e+05 -1.673717180165112950e-02 +1.131600000000000000e+05 3.888223836838733405e-04 +1.140600000000000000e+05 7.795694331434788182e-03 +1.149600000000000000e+05 -5.989849554680404253e-03 +1.158600000000000000e+05 -9.412296880327630788e-03 +1.167600000000000000e+05 -4.879836572399653960e-03 +1.176600000000000000e+05 2.124258357343933312e-02 +1.185600000000000000e+05 1.863953532028972404e-02 +1.194600000000000000e+05 -1.194937880245561246e-02 +1.203600000000000000e+05 1.427512742975522997e-02 +1.212600000000000000e+05 -1.682425932813202962e-02 +1.221600000000000000e+05 -1.331873392155102920e-02 +1.230600000000000000e+05 7.882993777457159013e-04 +1.239600000000000000e+05 -3.438177172938594595e-03 +1.248600000000000000e+05 1.347117413388332352e-04 +1.257600000000000000e+05 1.706155614556337241e-03 +1.266600000000000000e+05 1.539015512207697611e-03 +1.275600000000000000e+05 1.395615599540178664e-02 +1.284600000000000000e+05 -7.664193280106701422e-03 +1.293600000000000000e+05 1.010602402675431222e-02 +1.302600000000000000e+05 -5.261662406155664939e-03 +1.311600000000000000e+05 -4.259319705852249172e-03 +1.320600000000000000e+05 -2.350289103560498916e-03 +1.329600000000000000e+05 -7.977350707733421586e-03 +1.338600000000000000e+05 1.428751766070490703e-03 +1.347600000000000000e+05 8.440559807240788359e-03 +1.356600000000000000e+05 3.625059850492107216e-03 +1.365600000000000000e+05 5.534892904506705236e-03 +1.374600000000000000e+05 -3.300247534752998035e-03 +1.383600000000000000e+05 8.618106423455174081e-03 +1.392600000000000000e+05 -7.250626406857918482e-03 +1.401600000000000000e+05 -1.349332137942838017e-02 +1.410600000000000000e+05 -6.749700095497246366e-03 +1.419600000000000000e+05 9.281914337407215498e-03 +1.428600000000000000e+05 -1.416026971492101438e-02 +1.437600000000000000e+05 1.009458604130486492e-02 +1.446600000000000000e+05 -9.852923979451588821e-03 +1.455600000000000000e+05 1.002575051006715512e-02 +1.464600000000000000e+05 -3.135820761599461548e-04 +1.473600000000000000e+05 -8.987438848635065369e-03 +1.482600000000000000e+05 -1.183138265332672745e-03 +1.491600000000000000e+05 8.438082559223403223e-04 +1.500600000000000000e+05 -2.266519841214176267e-04 +1.509600000000000000e+05 1.225524485562345944e-03 +1.518600000000000000e+05 7.766027524667151738e-03 +1.527600000000000000e+05 1.491256866029289085e-02 +1.536600000000000000e+05 6.141980730717477854e-03 +1.545600000000000000e+05 -1.210207033273036359e-02 +1.554600000000000000e+05 1.059911638094490627e-02 +1.563600000000000000e+05 -4.352161000497289933e-03 +1.572600000000000000e+05 -7.561220061688800342e-03 +1.581600000000000000e+05 7.367878652075887658e-03 +1.590600000000000000e+05 2.841149690539168660e-03 +1.599600000000000000e+05 1.228330634148733225e-02 +1.608600000000000000e+05 9.146053481345006730e-03 +1.617600000000000000e+05 -1.408408164570573717e-02 +1.626600000000000000e+05 -7.878504583459289279e-03 +1.635600000000000000e+05 -6.660009203187655658e-03 diff --git a/examples/tsunami/chile2010_fgmax-fgout/Makefile b/examples/tsunami/chile2010_fgmax-fgout/Makefile new file mode 100644 index 000000000..b548cfba7 --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/Makefile @@ -0,0 +1,72 @@ +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = geoclaw # Clawpack package to use +EXE = xgeoclaw # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow +include $(GEOLIB)/Makefile.geoclaw + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + +# ---------------------------------------- +# List of custom sources for this program: +# ---------------------------------------- + + +MODULES = \ + +SOURCES = \ + $(CLAW)/riemann/src/rpn2_geoclaw.f \ + $(CLAW)/riemann/src/rpt2_geoclaw.f \ + $(CLAW)/riemann/src/geoclaw_riemann_utils.f \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +# Construct the topography data +.PHONY: topo all fgout_plots +topo: + python maketopo.py + +fgout_plots: + $(MAKE) plots SETPLOT_FILE=setplot_fgout.py PLOTDIR=_plots_fgout + +all: + $(MAKE) topo + $(MAKE) .plots + $(MAKE) fgout_plots + $(MAKE) .htmls + diff --git a/examples/tsunami/chile2010_fgmax-fgout/README.rst b/examples/tsunami/chile2010_fgmax-fgout/README.rst new file mode 100644 index 000000000..c627cd0e2 --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/README.rst @@ -0,0 +1,149 @@ + +.. _geoclaw_examples_tsunami_chile2010_fgmax-fgout: + +Chile 2010 test case for fgmax and fgout routines +================================================= + +This example illustrates the use of an fgmax grid and an fgout grid, +which can be used independently of one another. + +For details about this test problem, see +`$CLAW/geoclaw/examples/tsunami/chile2010 +`__. + +fgmax grid +---------- + +See https://www.clawpack.org/fgmax.html for general documentation of +fgmax grids. In this example, an fgmax grid is used to monitor the +maximum amplitude of the wave at each point in the domain and the +arrival times, in order to make a plot displaying these over the +computational domain, a portion of the south Pacific. + +The fgmax grid is specified in `setrun.py` and doing `make data` +(or `make .output`) leads to the creation of a file `fgmax_grids.data` +that is read into GeoClaw. + +To test:: + + make topo + make .output + make plots # to make frame plots and _PlotIndex.html + +Or simply:: + + make all + +This should produce +`_plots/amplitude_times.png <./_plots/amplitude_times.png>`_, +a color map of maximum amplitudes along with contours of arrival +time. This is generated by the code in `plot_fgmax.py` and +a link to this plot should show up in `_plots/_PlotIndex.html` +along with the usual time frame plots. + +*Note:* + +- See http://www.clawpack.org/fgmax.html for more information about + specifying fgmax parameters. + +- The time `fg.tstart_max` in `setrun.py` is set to 10 seconds so that the + topography in the source region has been finalized following the + earthquake before we start monitoring the maxima. (Since the topo on the + fixed grid must also be stored for later postprocessing.) + +- The refinement parameters and regions are set so that the maximum + amplitude we wish to capture always appears on a level 3 grid and + `fg.min_level_check = 3` is set in `setrun.py`. Other choices of these + parameters may give misleading or bizarre results. The fgmax capabilities + were designed with the assumption that the region of interest will always + be refined to the maximum level allowed. + +- The code in `plot_fgmax.py` is used to plot the fgmax results. Also the file + `setplot.py` includes the lines:: + + otherfigure = plotdata.new_otherfigure(name='max amplitude and arrival times', + fname='amplitude_times.png') + + def makefig(plotdata): + plot_fgmax.plot_fgmax_grid(plotdata.outdir, plotdata.plotdir) + + otherfigure.makefig = makefig + + This results in `plot_fgmax.plot_fgmax_grid` being run and + the link to the resulting figure appearing in `_plots/_PlotIndex.html`. + +fgout grid +---------- + +See https://www.clawpack.org/fgout.html for general documentation of +fgout grids. Here a single fgout grid is specified in `setrun.py` +that covers most of the computational domain at a fixed resolution. +The solution interpolated to this grid is output every 15 minutes, +more frequently than the usual GeoClaw frames, which in this example +are output only every 2 hours. + + +To test:: + + make topo + make .output + make fgout_plots # creates _plots_fgout directory + python make_fgout_animation.py # creates fgout_animation.mp4 + +Note that `make fgout_plots` is defined in the `Makefile` as shorthand for:: + + make plots SETPLOT_FILE=setplot_fgout.py PLOTDIR=_plots_fgout + +This illustrates one approach to plotting fgout grid results: A setplot +function is specified (in this case by `setplot_fgout.py`) that has the +same form as a setplot function for plotting standard GeoClaw/Clawpack +output frames, but setting :: + + plotdata.file_prefix = 'fgout0001' # for fgout grid fgno==1 + +to indicate that instead of the usual output files with names like +`fort.t*` and `fort.q*` (and also `fort.b*` in the case of binary output), +as described at https://www.clawpack.org/output_styles.html, +the data is in files named `fgout0001.t*`, etc. but with the same format. + +This creates the usual sort of `_plots` directory displaying all the +resulting fgout plots. In this example we have called this directory +`_plots_fgout <./_plots_fgout/_PlotIndex.html>`_ +to differentiate it from the directory +`_plots <./_plots/_PlotIndex.html>`_ +which contains the usual plots from output times. + +Loading and plotting fgout results directly, and making an mp4 animation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Alternatively, since every fgout frame consists of only a single +uniform grid of data, it is much easier to manipulate or plot +directly than general AMR data. The `clawpack.geoclaw.fgout_tools` +module described at https://www.clawpack.org/fgout_tools_module.html +provides tools for reading frames and producing +arrays that can then be worked with directly. + + +An example of how this might be done is provided in `plot_fgout.py`, +where a single frame is plotted. The sample code in +`make_fgout_animation.py` reads in all the frames of fgout data +and produces an animation as a stand-alone mp4 file. The use of +fgout grids provides a way to produce frequent outputs on a fixed grid +resolution, as often desired for making smooth animations of a portion of +the computational domain. + +Saving a sequence of fgout frames to a single netCDF file +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The script `make_netcdf.py` illustrates how to combine multiple fgout +frames of data into a single netCDF file using `fgout_tools.write_netcdf`. +This is easily done since all the fgout results are on the same uniform +grid. You can also select which quantities of interest to store and use +32-bit floats to store them, so that the resulting file takes much less +space than the original collection of fgout files. + +The script `make_netcdf.py` also illustrates how to read the arrays back in +from the netCDF file. + +This example requires the Python module `netCDF4 +`__. diff --git a/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py new file mode 100644 index 000000000..f6324689d --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py @@ -0,0 +1,133 @@ +""" +Make an mp4 animation of fgout grid results. +This is done in a way that makes the animation quickly and with minimum +storage required, by making one plot and then defining an update function +that only changes the parts of the plot that change in each frame. +The tuple update_artists contains the list of Artists that must be changed +in update. Modify this as needed. + +""" + +from pylab import * +import os +from clawpack.visclaw import plottools, geoplot +from clawpack.visclaw import animation_tools +from matplotlib import animation, colors +from datetime import timedelta + +from clawpack.geoclaw import fgout_tools + +fgno = 1 # which fgout grid + +outdir = '_output' +format = 'binary' # format of fgout grid output + +fgframes = range(1,26) # frames of fgout solution to use in animation + +figsize = (10,8) + +# Instantiate object for reading fgout frames: +fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format) + + +# Plot one frame of fgout data and define the Artists that will need to +# be updated in subsequent frames: + +fgout1 = fgout_grid.read_frame(fgframes[0]) + +plot_extent = fgout1.extent_edges + +ylat = fgout1.Y.mean() # for aspect ratio of plots + +fig,ax = subplots(figsize=figsize) + +ax.set_xlim(plot_extent[:2]) +ax.set_ylim(plot_extent[2:]) + +plottools.pcolorcells(fgout1.X,fgout1.Y,fgout1.B, cmap=geoplot.land_colors) +clim(0,100) + +eta = ma.masked_where(fgout1.h<0.001, fgout1.eta) + +eta_plot = plottools.pcolorcells(fgout1.X,fgout1.Y,eta, + cmap=geoplot.tsunami_colormap) +clim(-0.3,0.3) +cb = colorbar(eta_plot, extend='both', shrink=0.7) +cb.set_label('meters') +title_text = title('Surface at time %s' % timedelta(seconds=fgout1.t)) + +ax.set_aspect(1./cos(ylat*pi/180.)) +ticklabel_format(useOffset=False) +xticks(rotation=20) +ax.set_xlim(plot_extent[:2]) +ax.set_ylim(plot_extent[2:]) + + +# The artists that will be updated for subsequent frames: +update_artists = (eta_plot, title_text) + + +def update(fgframeno, *update_artists): + """ + Update an exisiting plot with solution from fgout frame fgframeno. + The artists in update_artists must have new data assigned. + """ + + fgout = fgout_grid.read_frame(fgframeno) + print('Updating plot at time %s' % timedelta(seconds=fgout.t)) + + # unpack update_artists (must agree with definition above): + eta_plot, title_text = update_artists + + # reset title to current time: + title_text.set_text('Surface at time %s' % timedelta(seconds=fgout.t)) + + # reset surface eta to current state: + eta = ma.masked_where(fgout.h<0.001, fgout.eta) + eta_plot.set_array(eta.T.flatten()) + + update_artists = (eta_plot, title_text) + return update_artists + +def plot_fgframe(fgframeno): + """ + Convenience function for plotting one frame. + But if you use this function in IPython and then try to make the animation, + it may get into an infinite loop (not sure why). Close the figure to abort. + """ + update(fgframeno, *update_artists) + + +def make_anim(): + print('Making anim...') + anim = animation.FuncAnimation(fig, update, + frames=fgframes, + fargs=update_artists, + interval=200, blit=True) + return anim + +if __name__ == '__main__': + + anim = make_anim() + + # Output files: + name = 'fgout_animation' + + fname_mp4 = name + '.mp4' + fname_html = None # name + '.html' + + if fname_mp4: + fps = 5 + print('Making mp4...') + writer = animation.writers['ffmpeg'](fps=fps) + anim.save(fname_mp4, writer=writer) + print("Created %s" % fname_mp4) + + if fname_html: + # html version: + fname_html = name + '.html' + animation_tools.make_html(anim, file_name=fname_html, title=name) + + + + diff --git a/examples/tsunami/chile2010_fgmax-fgout/make_netcdf.py b/examples/tsunami/chile2010_fgmax-fgout/make_netcdf.py new file mode 100644 index 000000000..048a21c56 --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/make_netcdf.py @@ -0,0 +1,67 @@ + +""" +Test script to read in a sequence of fgout frames, and write a netCDF file +containing only eta from each fgout frame, along with the final topography B. + +Then print info about it, and read it back in. +With datatype = 'f8' there should be no difference from the original. +With datatype = 'f4' (32-bit floats) the netCDF file is only half as large. + +This example should be turned into a notebook with more discussin... +""" + +from clawpack.geoclaw import fgout_tools + +try: + import netCDF4 +except: + 'You must install netCDF4 in order to write a netCDF file' + raise + +fgno = 1 # which fgout grid + +outdir = '_output' +format = 'binary' # format of fgout grid output + +fgframenos = range(1,5) # frames of fgout solution to write to netCDF file + +# Instantiate object for reading fgout frames: +fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format) + +fgout_frames = [] +for fgframeno in fgframenos: + fgout = fgout_grid.read_frame(fgframeno) + fgout_frames.append(fgout) + +datatype = 'f4' # 'f4' or 'f8' + +fgout_tools.write_netcdf(fgout_frames, + fname_nc='fgout_frames.nc', + qois=['eta'], + datatype=datatype, + include_B0=True, + include_Bfinal=False, + description='Chile 2010 test problem', + verbose=True) + +# ================================================ +# Print out info about this file: +fgout_tools.print_netcdf_info('fgout_frames.nc') + +# ================================================ +# Read back the saved fgout frames as arrays: + +x,y,t,qoi_arrays = fgout_tools.read_netcdf_arrays('fgout_frames.nc', + ['eta','B0'], verbose=True) + +B0 = qoi_arrays['B0'] +eta = qoi_arrays['eta'] + +# Compare h in first frame with the original data: +h0_original = fgout_frames[0].h +h0 = eta[0,:,:] - B0 + +dh0 = abs(h0 - h0_original).max() +print("With datatype='%s', abs(h0 - h0_original).max() = %.3e" % (datatype,dh0)) + + \ No newline at end of file diff --git a/examples/tsunami/chile2010_fgmax-fgout/maketopo.py b/examples/tsunami/chile2010_fgmax-fgout/maketopo.py new file mode 100644 index 000000000..c77ec4dfa --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/maketopo.py @@ -0,0 +1,113 @@ +""" +Create topo and dtopo files needed for this example: + etopo10min120W60W60S0S.asc download from GeoClaw topo repository + dtopo_usgs100227.tt3 create using Okada model +Prior to Clawpack 5.2.1, the fault parameters we specified in a .cfg file, +but now they are explicit below. + +Call functions with makeplots==True to create plots of topo, slip, and dtopo. +""" + +from __future__ import absolute_import +from __future__ import print_function +import os + +import clawpack.clawutil.data + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW enviornment variable") + +# Scratch directory for storing topo and dtopo files: +scratch_dir = os.path.join(CLAW, 'geoclaw', 'scratch') + +def get_topo(makeplots=False): + """ + Retrieve the topo file from the GeoClaw repository. + """ + from clawpack.geoclaw import topotools + topo_fname = 'etopo10min120W60W60S0S.asc' + url = 'http://depts.washington.edu/clawpack/geoclaw/topo/etopo/' + topo_fname + clawpack.clawutil.data.get_remote_file(url, output_dir=scratch_dir, + file_name=topo_fname, verbose=True) + + if makeplots: + from matplotlib import pyplot as plt + topo = topotools.Topography(os.path.join(scratch_dir,topo_fname), topo_type=2) + topo.plot() + fname = os.path.splitext(topo_fname)[0] + '.png' + plt.savefig(fname) + print("Created ",fname) + + + +def make_dtopo(makeplots=False): + """ + Create dtopo data file for deformation of sea floor due to earthquake. + Uses the Okada model with fault parameters and mesh specified below. + """ + from clawpack.geoclaw import dtopotools + import numpy + + dtopo_fname = os.path.join(scratch_dir, "dtopo_usgs100227.tt3") + + # Specify subfault parameters for this simple fault model consisting + # of a single subfault: + + usgs_subfault = dtopotools.SubFault() + usgs_subfault.strike = 16. + usgs_subfault.length = 450.e3 + usgs_subfault.width = 100.e3 + usgs_subfault.depth = 35.e3 + usgs_subfault.slip = 15. + usgs_subfault.rake = 104. + usgs_subfault.dip = 14. + usgs_subfault.longitude = -72.668 + usgs_subfault.latitude = -35.826 + usgs_subfault.coordinate_specification = "top center" + + fault = dtopotools.Fault() + fault.subfaults = [usgs_subfault] + + print("Mw = ",fault.Mw()) + + if os.path.exists(dtopo_fname): + print("*** Not regenerating dtopo file (already exists): %s" \ + % dtopo_fname) + else: + print("Using Okada model to create dtopo file") + + x = numpy.linspace(-77, -67, 100) + y = numpy.linspace(-40, -30, 100) + times = [1.] + + fault.create_dtopography(x,y,times) + dtopo = fault.dtopo + dtopo.write(dtopo_fname, dtopo_type=3) + + + if makeplots: + from matplotlib import pyplot as plt + if fault.dtopo is None: + # read in the pre-existing file: + print("Reading in dtopo file...") + dtopo = dtopotools.DTopography() + dtopo.read(dtopo_fname, dtopo_type=3) + x = dtopo.x + y = dtopo.y + plt.figure(figsize=(12,7)) + ax1 = plt.subplot(121) + ax2 = plt.subplot(122) + fault.plot_subfaults(axes=ax1,slip_color=True) + ax1.set_xlim(x.min(),x.max()) + ax1.set_ylim(y.min(),y.max()) + dtopo.plot_dZ_colors(1.,axes=ax2) + fname = os.path.splitext(os.path.split(dtopo_fname)[-1])[0] + '.png' + plt.savefig(fname) + print("Created ",fname) + + +if __name__=='__main__': + get_topo(False) + make_dtopo(False) diff --git a/examples/tsunami/chile2010_fgmax-fgout/plot_fgmax.py b/examples/tsunami/chile2010_fgmax-fgout/plot_fgmax.py new file mode 100644 index 000000000..b6fb84368 --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/plot_fgmax.py @@ -0,0 +1,50 @@ +""" +Plot fgmax output from GeoClaw run. + +""" + +import matplotlib.pyplot as plt +import os +import numpy +from clawpack.geoclaw import fgmax_tools +from clawpack.visclaw import geoplot + +def plot_fgmax_grid(outdir,plotdir): + + fg = fgmax_tools.FGmaxGrid() + fg.outdir = outdir + data_file = os.path.join(outdir, 'fgmax_grids.data') + fg.read_fgmax_grids_data(fgno=1, data_file=data_file) + fg.read_output() + + clines_zeta = [0.01] + list(numpy.linspace(0.05,0.3,6)) + [0.5,1.0,10.0] + colors = geoplot.discrete_cmap_1(clines_zeta) + plt.figure(1) + plt.clf() + zeta = numpy.where(fg.B>0, fg.h, fg.h+fg.B) # surface elevation in ocean + plt.contourf(fg.X,fg.Y,zeta,clines_zeta,colors=colors) + plt.colorbar() + plt.contour(fg.X,fg.Y,fg.B,[0.],colors='k') # coastline + + # plot arrival time contours and label: + arrival_t = fg.arrival_time/3600. # arrival time in hours + clines_t = numpy.linspace(0,8,17) # hours + clines_t_label = clines_t[2::2] # which ones to label + clines_t_colors = ([.5,.5,.5],) + con_t = plt.contour(fg.X,fg.Y,arrival_t, clines_t,colors=clines_t_colors) + plt.clabel(con_t, clines_t_label) + + # fix axes: + plt.ticklabel_format(style='plain',useOffset=False) + plt.xticks(rotation=20) + plt.gca().set_aspect(1./numpy.cos(fg.Y.mean()*numpy.pi/180.)) + plt.title("Maximum amplitude / arrival times") + + if not os.path.isdir(plotdir): + os.mkdir(plotdir) + fname = os.path.join(plotdir, "amplitude_times.png") + plt.savefig(fname) + print("Created ",fname) + +if __name__=="__main__": + plot_fgmax_grid(outdir='_output', plotdir='.') diff --git a/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py b/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py new file mode 100644 index 000000000..8b61666d3 --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py @@ -0,0 +1,40 @@ + +from pylab import * +import os +from clawpack.visclaw import plottools, geoplot +from clawpack.visclaw import animation_tools +from matplotlib import animation, colors +from clawpack.geoclaw import fgout_tools +from datetime import timedelta + +fgno = 1 +outdir = '_output' +output_format = 'binary' + +# Instantiate object for reading fgout frames: +fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, output_format) + +# Plot one frame of fgout data +fgframe = 20 +fgout = fgout_grid.read_frame(2) + +fgout = fgout_grid.read_frame(fgframe) + +figure(1, figsize=(8,8)) +imshow(flipud(fgout.B.T), extent=fgout.extent_edges, + cmap=geoplot.land_colors) + +clim(0,100) + +eta_water = where(fgout.h > 0, fgout.eta, nan) +imshow(flipud(eta_water.T), extent=fgout.extent_edges, + cmap=geoplot.tsunami_colormap) + +clim(-0.2, 0.2) + +title('Surface at time %s' % timedelta(seconds=fgout.t)) + +fname = 'fgout_frame%s.png' % str(fgframe).zfill(4) +savefig(fname) +print('Created ',fname) + diff --git a/examples/tsunami/chile2010_fgmax-fgout/setplot.py b/examples/tsunami/chile2010_fgmax-fgout/setplot.py new file mode 100644 index 000000000..74a2d7133 --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/setplot.py @@ -0,0 +1,184 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + +import numpy as np +import matplotlib.pyplot as plt + +from clawpack.geoclaw import topotools + +import plot_fgmax + +try: + TG32412 = np.loadtxt('32412_notide.txt') +except: + print("*** Could not load DART data file") + +#-------------------------- +def setplot(plotdata): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + from clawpack.visclaw import colormaps, geoplot + from numpy import linspace + + plotdata.clearfigures() # clear any old figures,axes,items data + + + # To plot gauge locations on pcolor or contour plot, use this as + # an afteraxis function: + + def addgauges(current_data): + from clawpack.visclaw import gaugetools + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos='all', format_string='ko', add_labels=True) + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.title = 'Surface' + plotaxes.scaled = True + + def fixup(current_data): + from pylab import title + from datetime import timedelta + addgauges(current_data) + t = current_data.t + title('Surface at time %s' % timedelta(seconds=t)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = -0.2 + plotitem.pcolor_cmax = 0.2 + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 1 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [1,0,0] + plotitem.patchedges_show = 1 + plotaxes.xlimits = [-120,-60] + plotaxes.ylimits = [-60,0] + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface at gauges', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Surface' + + # Plot surface as blue curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 3 + plotitem.plotstyle = 'b-' + + # Plot topo as green curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.show = False + + def gaugetopo(current_data): + q = current_data.q + h = q[0,:] + eta = q[3,:] + topo = eta - h + return topo + + plotitem.plot_var = gaugetopo + plotitem.plotstyle = 'g-' + + def add_zeroline(current_data): + from pylab import plot, legend, xticks, floor, axis, xlabel + t = current_data.t + gaugeno = current_data.gaugeno + + if gaugeno == 32412: + try: + plot(TG32412[:,0], TG32412[:,1], 'r') + legend(['GeoClaw','Obs'],loc='lower right') + except: pass + axis((0,t.max(),-0.3,0.3)) + + plot(t, 0*t, 'k') + n = int(floor(t.max()/3600.) + 2) + xticks([3600*i for i in range(n)], ['%i' % i for i in range(n)]) + xlabel('time (hours)') + + plotaxes.afteraxes = add_zeroline + + + #----------------------------------------- + # Figures for fgmax - max values on fixed grids + #----------------------------------------- + otherfigure = plotdata.new_otherfigure(name='max amplitude and arrival times', + fname='amplitude_times.png') + + def makefig(plotdata): + plot_fgmax.plot_fgmax_grid(plotdata.outdir, plotdata.plotdir) + + otherfigure.makefig = makefig + + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via pyclaw.plotters.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_gaugenos = 'all' # list of gauges to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + + return plotdata + diff --git a/examples/tsunami/chile2010_fgmax-fgout/setplot_fgout.py b/examples/tsunami/chile2010_fgmax-fgout/setplot_fgout.py new file mode 100644 index 000000000..d84e327a6 --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/setplot_fgout.py @@ -0,0 +1,181 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + +import numpy as np +import matplotlib.pyplot as plt + +from clawpack.geoclaw import topotools + +try: + TG32412 = np.loadtxt('32412_notide.txt') +except: + print("*** Could not load DART data file") + +#-------------------------- +def setplot(plotdata): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + from clawpack.visclaw import colormaps, geoplot + from numpy import linspace + + plotdata.clearfigures() # clear any old figures,axes,items data + + # To plot fgout files rather than usual frames of AMR solution: + plotdata.file_prefix = 'fgout0001' # for fgout grid fgno==1 + plotdata.format = 'binary' # 'ascii' or 'binary' to match fgout + + # To plot gauge locations on pcolor or contour plot, use this as + # an afteraxis function: + + def addgauges(current_data): + from clawpack.visclaw import gaugetools + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos='all', format_string='ko', add_labels=True) + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.title = 'Surface' + plotaxes.scaled = True + + def fixup(current_data): + from pylab import title + from datetime import timedelta + addgauges(current_data) + t = current_data.t + title('Surface at time %s' % timedelta(seconds=t)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = -0.2 + plotitem.pcolor_cmax = 0.2 + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 1 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [1,0,0] + plotitem.patchedges_show = 1 + plotaxes.xlimits = [-120,-60] + plotaxes.ylimits = [-60,0] + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface at gauges', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Surface' + + # Plot surface as blue curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 3 + plotitem.plotstyle = 'b-' + + # Plot topo as green curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.show = False + + def gaugetopo(current_data): + q = current_data.q + h = q[0,:] + eta = q[3,:] + topo = eta - h + return topo + + plotitem.plot_var = gaugetopo + plotitem.plotstyle = 'g-' + + def add_zeroline(current_data): + from pylab import plot, legend, xticks, floor, axis, xlabel + t = current_data.t + gaugeno = current_data.gaugeno + + if gaugeno == 32412: + try: + plot(TG32412[:,0], TG32412[:,1], 'r') + legend(['GeoClaw','Obs'],loc='lower right') + except: pass + axis((0,t.max(),-0.3,0.3)) + + plot(t, 0*t, 'k') + n = int(floor(t.max()/3600.) + 2) + xticks([3600*i for i in range(n)], ['%i' % i for i in range(n)]) + xlabel('time (hours)') + + plotaxes.afteraxes = add_zeroline + + + #----------------------------------------- + # Figures for fgmax - max values on fixed grids + #----------------------------------------- + otherfigure = plotdata.new_otherfigure(name='max amplitude and arrival times', + fname='amplitude_times.png') + + + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via pyclaw.plotters.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' #range(25) # list of frames to print + plotdata.print_gaugenos = 'all' # list of gauges to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + + return plotdata + diff --git a/examples/tsunami/chile2010_fgmax-fgout/setrun.py b/examples/tsunami/chile2010_fgmax-fgout/setrun.py new file mode 100644 index 000000000..f4e54cab5 --- /dev/null +++ b/examples/tsunami/chile2010_fgmax-fgout/setrun.py @@ -0,0 +1,465 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +import os +import numpy as np +from clawpack.geoclaw import fgmax_tools +from clawpack.geoclaw import fgout_tools + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW enviornment variable") + +# Scratch directory for storing topo and dtopo files: +scratch_dir = os.path.join(CLAW, 'geoclaw', 'scratch') + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + #probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + + + #------------------------------------------------------------------ + # GeoClaw specific parameters are set later + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + clawdata.lower[0] = -120.0 # west longitude + clawdata.upper[0] = -60.0 # east longitude + + clawdata.lower[1] = -60.0 # south latitude + clawdata.upper[1] = 0.0 # north latitude + + # Number of grid cells: Coarsest grid + clawdata.num_cells[0] = 30 + clawdata.num_cells[1] = 30 + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 2 + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # Note: If restarting, you must also change the Makefile to set: + # RESTART = True + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = 'fort.chk00036' # File to use for restart data + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 5 + clawdata.tfinal = 10*3600.0 + clawdata.output_t0 = False # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [30., 60., 300., 600.] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 3 + clawdata.output_t0 = True + + + clawdata.output_format = 'binary32' # 'ascii', 'binary32', 'binary64' + + clawdata.output_q_components = 'all' # need all + clawdata.output_aux_components = 'none' # eta=h+B is in q + clawdata.output_aux_onlyonce = False # output aux arrays each frame + + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # If dt_variable==0 then dt=dt_initial for all steps: + clawdata.dt_initial = 0.2 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + clawdata.cfl_desired = 0.75 + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 5000 + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = ['mc', 'mc', 'mc'] + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' + clawdata.bc_upper[0] = 'extrap' + + clawdata.bc_lower[1] = 'extrap' + clawdata.bc_upper[1] = 'extrap' + + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif np.abs(clawdata.checkpt_style) == 1: + # Checkpoint only at tfinal. + pass + + elif np.abs(clawdata.checkpt_style) == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = [0.1,0.15] + + elif np.abs(clawdata.checkpt_style) == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + # maximum size of patches in each direction (matters in parallel): + amrdata.max1d = 60 + + # max number of refinement levels: + amrdata.amr_levels_max = 3 + + # List of refinement ratios at each level (length at least mxnest-1) + amrdata.refinement_ratios_x = [4,3] + amrdata.refinement_ratios_y = [4,3] + amrdata.refinement_ratios_t = [4,3] + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = ['center','capacity','yleft'] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag2refine = True + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 0 + + + # --------------- + # Regions: + # --------------- + rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + #rundata.regiondata.regions.append([1, 2, 0., 1e9, -360,360,-90,90]) + rundata.regiondata.regions.append([3, 3, 0., 600., -76,-72,-38,-33]) + + # --------------- + # Gauges: + # --------------- + rundata.gaugedata.gauges = [] + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + rundata.gaugedata.gauges.append([32412, -86.392, -17.975, 0., 1.e10]) + + # -------------------- + # GeoClaw parameters: + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 2 + geo_data.earth_radius = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 0.0 + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = True + geo_data.manning_coefficient =.025 + geo_data.friction_depth = 1e6 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = True + refinement_data.wave_tolerance = 0.01 + + # == settopo.data values == + topo_data = rundata.topo_data + # for topography, append lines of the form + # [topotype, fname] + topo_path = os.path.join(scratch_dir, 'etopo10min120W60W60S0S.asc') + topo_data.topofiles.append([2, topo_path]) + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + # for moving topography, append lines of the form : (<= 1 allowed for now!) + # [topotype, fname] + dtopo_path = os.path.join(scratch_dir, 'dtopo_usgs100227.tt3') + dtopo_data.dtopofiles.append([3,dtopo_path]) + dtopo_data.dt_max_dtopo = 0.2 + + + # == setqinit.data values == + rundata.qinit_data.qinit_type = 0 + rundata.qinit_data.qinitfiles = [] + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # == setfixedgrids.data values == + # rundata.fixed_grid_data removed in v5.9.0; instead use fgmax and/or fgout + + + # == fgmax_grids.data values == + # NEW STYLE STARTING IN v5.7.0 + + # set num_fgmax_val = 1 to save only max depth, + # 2 to also save max speed, + # 5 to also save max hs,hss,hmin + rundata.fgmax_data.num_fgmax_val = 2 # Save depth and speed + + fgmax_grids = rundata.fgmax_data.fgmax_grids # empty list to start + + # Now append to this list objects of class fgmax_tools.FGmaxGrid + # specifying any fgmax grids. + + # Points on a uniform 2d grid: + dx_fine = 2./(3.*4.) # grid resolution at finest level + + fg = fgmax_tools.FGmaxGrid() + fg.point_style = 2 # uniform rectangular x-y grid + fg.x1 = -120. + dx_fine/2. # specify pts to align with FV cell centers + fg.x2 = -60. - dx_fine/2. + fg.y1 = -60. + dx_fine/2. + fg.y2 = 0. - dx_fine/2. + fg.dx = dx_fine + fg.tstart_max = 10. # when to start monitoring max values + fg.tend_max = 1.e10 # when to stop monitoring max values + fg.dt_check = 60. # target time (sec) increment between updating + # max values + fg.min_level_check = 3 # which levels to monitor max on + fg.arrival_tol = 1.e-2 # tolerance for flagging arrival + + fg.interp_method = 0 # 0 ==> pw const in cells, recommended + fgmax_grids.append(fg) # written to fgmax_grids.data + + + # == fgout_grids.data values == + # NEW IN v5.9.0 + # Set rundata.fgout_data.fgout_grids to be a list of + # objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + fgout_grids = rundata.fgout_data.fgout_grids # empty list initially + + fgout = fgout_tools.FGoutGrid() + fgout.fgno = 1 + fgout.point_style = 2 # will specify a 2d grid of points + fgout.output_format = 'binary32' # 4-byte, float32 + fgout.nx = 200 + fgout.ny = 250 + fgout.x1 = -115. # specify edges (fgout pts will be cell centers) + fgout.x2 = -70. + fgout.y1 = -55. + fgout.y2 = -10. + fgout.tstart = 0. + fgout.tend = 6.*3600 + fgout.nout = 25 + fgout_grids.append(fgout) # written to fgout_grids.data + + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = True # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + # More AMR parameters can be set -- see the defaults in pyclaw/data.py + + return rundata + + # end of function setrun + # ---------------------- + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() + diff --git a/examples/tsunami/island-particles/setrun.py b/examples/tsunami/island-particles/setrun.py index dc231ab8f..b6fb66c1f 100644 --- a/examples/tsunami/island-particles/setrun.py +++ b/examples/tsunami/island-particles/setrun.py @@ -396,11 +396,11 @@ def setgeo(rundata): # [fname] rundata.qinit_data.qinitfiles.append(['qinit.xyz']) - # == setfixedgrids.data values == - fixedgrids = rundata.fixed_grid_data.fixedgrids - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] + # == fgout grids == + # new style as of v5.9.0 (old rundata.fixed_grid_data is deprecated) + # set rundata.fgout_data.fgout_grids to be a + # list of objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + #rundata.fgout_data.fgout_grids = [] return rundata # end of function setgeo diff --git a/src/2d/shallow/Makefile.geoclaw b/src/2d/shallow/Makefile.geoclaw index 53785a983..f1f79c6d6 100644 --- a/src/2d/shallow/Makefile.geoclaw +++ b/src/2d/shallow/Makefile.geoclaw @@ -16,7 +16,7 @@ COMMON_MODULES += \ $(GEOLIB)/topo_module.f90 \ $(GEOLIB)/qinit_module.f90 \ $(GEOLIB)/refinement_module.f90 \ - $(GEOLIB)/fixedgrids_module.f90 \ + $(GEOLIB)/fgout_module.f90 \ $(GEOLIB)/fgmax_module.f90 \ $(GEOLIB)/surge/model_storm_module.f90 \ $(GEOLIB)/surge/data_storm_module.f90 \ diff --git a/src/2d/shallow/advanc.f b/src/2d/shallow/advanc.f index aa8ccaca6..96098565b 100644 --- a/src/2d/shallow/advanc.f +++ b/src/2d/shallow/advanc.f @@ -4,7 +4,7 @@ subroutine advanc (level,nvar,dtlevnew,vtime,naux) c use amr_module - use fixedgrids_module +c use fgout_module, only: ? use topo_module, only: topo_finalized implicit double precision (a-h,o-z) diff --git a/src/2d/shallow/amr2.f90 b/src/2d/shallow/amr2.f90 index 09b90daeb..139e01f1e 100644 --- a/src/2d/shallow/amr2.f90 +++ b/src/2d/shallow/amr2.f90 @@ -97,12 +97,12 @@ program amr2 use geoclaw_module, only: set_geo use topo_module, only: read_topo_settings, read_dtopo_settings use qinit_module, only: set_qinit - use fixedgrids_module, only: set_fixed_grids use refinement_module, only: set_refinement 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 @@ -491,7 +491,7 @@ program amr2 call read_dtopo_settings() ! specifies file with dtopo from earthquake call read_topo_settings(rest) ! specifies topography (bathymetry) files call set_qinit() ! specifies file with dh if this used instead - call set_fixed_grids() ! Fixed grid settings + call set_fgout(rest) ! Fixed grid settings call setup_variable_friction() ! Variable friction parameter !call set_multilayer() ! Set multilayer SWE parameters call set_storm() ! Set storm parameters @@ -526,7 +526,7 @@ program amr2 call read_dtopo_settings() ! specifies file with dtopo from earthquake call read_topo_settings(rest) ! specifies topography (bathymetry) files call set_qinit() ! specifies file with dh if this used instead - call set_fixed_grids() ! Fixed grid settings + call set_fgout(rest) ! Fixed grid settings call setup_variable_friction() ! Variable friction parameter call set_multilayer() ! Set multilayer SWE parameters call set_storm() ! Set storm parameters diff --git a/src/2d/shallow/fgout_module.f90 b/src/2d/shallow/fgout_module.f90 new file mode 100644 index 000000000..f36544f56 --- /dev/null +++ b/src/2d/shallow/fgout_module.f90 @@ -0,0 +1,626 @@ +module fgout_module + + implicit none + save + + ! Container for fixed grid data, geometry and output settings + type fgout_grid + ! Grid data + real(kind=8), pointer :: early(:,:,:) + real(kind=8), pointer :: late(:,:,:) + + ! Geometry + integer :: num_vars(2),mx,my,point_style,fgno,output_format + real(kind=8) :: dx,dy,x_low,x_hi,y_low,y_hi + + ! Time Tracking and output types + integer :: num_output,next_output_index + real(kind=8) :: start_time,end_time,dt + + integer, allocatable :: output_frames(:) + real(kind=8), allocatable :: output_times(:) + end type fgout_grid + + + logical, private :: module_setup = .false. + + ! Fixed grid arrays and sizes + integer :: FGOUT_num_grids + type(fgout_grid), target, allocatable :: FGOUT_fgrids(:) + real(kind=8) :: FGOUT_tcfmax + real(kind=8), parameter :: FGOUT_ttol = 1.d-6 ! tolerance for times + + +contains + + + ! Setup routine that reads in the fixed grids data file and sets up the + ! appropriate data structures + + subroutine set_fgout(rest,fname) + + use amr_module, only: parmunit, tstart_thisrun + + implicit none + + ! Subroutine arguments + logical :: rest ! restart? + character(len=*), optional, intent(in) :: fname + + ! Local storage + integer, parameter :: unit = 7 + integer :: i,k + type(fgout_grid), pointer :: fg + real(kind=8) :: ts + + + if (.not.module_setup) then + + write(parmunit,*) ' ' + write(parmunit,*) '--------------------------------------------' + write(parmunit,*) 'SETFGOUT:' + write(parmunit,*) '-----------' + + ! Open data file + if (present(fname)) then + call opendatafile(unit,fname) + else + call opendatafile(unit,'fgout_grids.data') + endif + + ! Read in data + read(unit,'(i2)') FGOUT_num_grids + write(parmunit,*) ' mfgrids = ',FGOUT_num_grids + if (FGOUT_num_grids == 0) then + write(parmunit,*) ' No fixed grids specified for output' + return + endif + + ! Allocate fixed grids (not the data yet though) + allocate(FGOUT_fgrids(FGOUT_num_grids)) + + ! Read in data for each fixed grid + do i=1,FGOUT_num_grids + fg => FGOUT_fgrids(i) + ! Read in this grid's data + read(unit,*) fg%fgno + read(unit,*) fg%start_time + read(unit,*) fg%end_time + read(unit,*) fg%num_output + read(unit,*) fg%point_style + read(unit,*) fg%output_format + read(unit,*) fg%mx, fg%my + read(unit,*) fg%x_low, fg%y_low + read(unit,*) fg%x_hi, fg%y_hi + + allocate(fg%output_times(fg%num_output)) + allocate(fg%output_frames(fg%num_output)) + + ! Initialize next_output_index + ! (might be reset below in case of a restart) + fg%next_output_index = 1 + + if (fg%point_style .ne. 2) then + print *, 'set_fgout: ERROR, unrecognized point_style = ',\ + fg%point_style + endif + + ! Setup data for this grid + ! Set dtfg (the timestep length between outputs) for each grid + if (fg%end_time <= fg%start_time) then + if (fg%num_output > 1) then + print *,'set_fgout: ERROR for fixed grid', i + print *,'start_time <= end_time yet num_output > 1' + print *,'set end_time > start_time or set num_output = 1' + stop + else + fg%dt = 0.d0 + endif + else + if (fg%num_output < 2) then + print *,'set_fgout: ERROR for fixed grid', i + print *,'end_time > start_time, yet num_output = 1' + print *,'set num_output > 2' + stop + else + fg%dt = (fg%end_time - fg%start_time) & + / (fg%num_output - 1) + do k=1,fg%num_output + fg%output_times(k) = fg%start_time + (k-1)*fg%dt + if (rest) then + ! don't write initial time or earlier + ts = tstart_thisrun+FGOUT_ttol + else + ! do write initial time + ts = tstart_thisrun-FGOUT_ttol + endif + + if (fg%output_times(k) < ts) then + ! will not output this time in this run + ! (might have already be done when restarting) + fg%output_frames(k) = -2 + fg%next_output_index = k+1 + else + ! will be reset to frameno when this is written + fg%output_frames(k) = -1 + endif + enddo + endif + endif + + + + ! Set spatial intervals dx and dy on each grid + if (fg%mx > 1) then + !fg%dx = (fg%x_hi - fg%x_low) / (fg%mx - 1) ! points + fg%dx = (fg%x_hi - fg%x_low) / fg%mx ! cells + else if (fg%mx == 1) then + fg%dx = 0.d0 + else + print *,'set_fgout: ERROR for fixed grid', i + print *,'x grid points mx <= 0, set mx >= 1' + endif + + if (fg%my > 1) then + !fg%dy = (fg%y_hi - fg%y_low) / (fg%my - 1) ! points + fg%dy = (fg%y_hi - fg%y_low) / fg%my ! cells + else if (fg%my == 1) then + fg%dy = 0.d0 + else + print *,'set_fgout: ERROR for fixed grid', i + print *,'y grid points my <= 0, set my >= 1' + endif + + ! set the number of variables stored for each grid + ! this should be (the number of variables you want to write out + 1) + fg%num_vars(1) = 6 + + ! Allocate new fixed grid data array + allocate(fg%early(fg%num_vars(1), fg%mx,fg%my)) + fg%early = nan() + + allocate(fg%late(fg%num_vars(1), fg%mx,fg%my)) + fg%late = nan() + + enddo + close(unit) + + FGOUT_tcfmax=-1.d16 + + module_setup = .true. + end if + + end subroutine set_fgout + + + !=====================FGOUT_INTERP======================================= + ! This routine interpolates q and aux on a computational grid + ! to an fgout grid not necessarily aligned with the computational grid + ! using bilinear interpolation defined on computational grid + !======================================================================= + subroutine fgout_interp(fgrid_type,fgrid, & + t,q,meqn,mxc,myc,mbc,dxc,dyc,xlowc,ylowc, & + maux,aux) + + use geoclaw_module, only: dry_tolerance + implicit none + + ! Subroutine arguments + integer, intent(in) :: fgrid_type + type(fgout_grid), intent(inout) :: fgrid + integer, intent(in) :: meqn,mxc,myc,mbc,maux + real(kind=8), intent(in) :: t,dxc,dyc,xlowc,ylowc + real(kind=8), intent(in) :: q(meqn,1-mbc:mxc+mbc,1-mbc:myc+mbc) + real(kind=8), intent(in) :: aux(maux,1-mbc:mxc+mbc,1-mbc:myc+mbc) + + integer, parameter :: method = 0 ! interpolate in space? + + ! Indices + integer :: ifg,jfg,m,ic1,ic2,jc1,jc2 + integer :: bathy_index,eta_index + + ! Tolerances + real(kind=8) :: total_depth,depth_indicator,nan_check + + ! Geometry + real(kind=8) :: xfg,yfg,xc1,xc2,yc1,yc2,xhic,yhic + real(kind=8) :: geometry(4) + + real(kind=8) :: points(2,2), eta_tmp + + ! Work arrays for eta interpolation + real(kind=8) :: eta(2,2),h(2,2) + + + ! Alias to data in fixed grid + integer :: num_vars + real(kind=8), pointer :: fg_data(:,:,:) + + + ! Setup aliases for specific fixed grid + if (fgrid_type == 1) then + num_vars = fgrid%num_vars(1) + fg_data => fgrid%early + else if (fgrid_type == 2) then + num_vars = fgrid%num_vars(1) + fg_data => fgrid%late + else + write(6,*) '*** Unexpected fgrid_type = ', fgrid_type + stop + ! fgrid_type==3 is deprecated, use fgmax grids instead + endif + + xhic = xlowc + dxc*mxc + yhic = ylowc + dyc*myc + + ! Find indices of various quantities in the fgrid arrays + bathy_index = meqn + 1 + eta_index = meqn + 2 + + !write(59,*) '+++ ifg,jfg,eta,geometry at t = ',t + + ! Primary interpolation loops + do ifg=1,fgrid%mx + xfg=fgrid%x_low + (ifg-0.5d0)*fgrid%dx ! cell centers + do jfg=1,fgrid%my + yfg=fgrid%y_low + (jfg-0.5d0)*fgrid%dy ! cell centers + + ! Check to see if this coordinate is inside of this grid + if (.not.((xfg < xlowc.or.xfg > xhic) & + .or.(yfg < ylowc.or.yfg > yhic))) then + + ! find where xfg,yfg is in the computational grid and + ! compute the indices + ! (Note: may be subject to rounding error if fgout point + ! is right on a cell edge!) + ic1 = int((xfg - xlowc + dxc)/dxc) + jc1 = int((yfg - ylowc + dyc)/dyc) + + if (method == 0) then + + ! piecewise constant: take values from cell (ic1,jc1): + + forall (m=1:meqn) + fg_data(m,ifg,jfg) = q(m,ic1,jc1) + end forall + + fg_data(bathy_index,ifg,jfg) = aux(1,ic1,jc1) + + ! for pw constant we take B, h, eta from same cell, + ! so setting eta = h+B should be fine even near shore: + fg_data(eta_index,ifg,jfg) = fg_data(1,ifg,jfg) & + + fg_data(bathy_index,ifg,jfg) + + + else if (method == 1) then + + ! bilinear used to interpolate to xfg,yfg + ! (not recommended) + + ! define constant parts of bilinear + if (ic1.eq.mxc) ic1=mxc-1 + if (jc1.eq.myc) jc1=myc-1 + ic2 = ic1 + 1 + jc2 = jc1 + 1 + + xc1 = xlowc + dxc * (ic1 - 0.5d0) + yc1 = ylowc + dyc * (jc1 - 0.5d0) + xc2 = xlowc + dxc * (ic2 - 0.5d0) + yc2 = ylowc + dyc * (jc2 - 0.5d0) + + geometry = [(xfg - xc1) / dxc, & + (yfg - yc1) / dyc, & + (xfg - xc1) * (yfg - yc1) / (dxc*dyc), & + 1.d0] + + + ! Interpolate all conserved quantities and bathymetry + forall (m=1:meqn) + fg_data(m,ifg,jfg) = & + interpolate(q(m,ic1:ic2,jc1:jc2), geometry) + end forall + + fg_data(bathy_index,ifg,jfg) = & + interpolate(aux(1,ic1:ic2,jc1:jc2),geometry) + + + ! surface eta = h + B: + + ! Note that for pw bilinear interp there may + ! be a problem interpolating each separately since + ! interpolated h + interpolated B may be much larger + ! than eta should be offshore. + eta = q(1,ic1:ic2,jc1:jc2) + aux(1,ic1:ic2,jc1:jc2) + fg_data(eta_index,ifg,jfg) = interpolate(eta,geometry) + ! NEED TO FIX + endif + + + ! save the time this fgout point was computed: + fg_data(num_vars,ifg,jfg) = t + + + endif ! if fgout point is on this grid + enddo ! fgout y-coordinate loop + enddo ! fgout x-coordinte loop + + end subroutine fgout_interp + + + !================ fgout_write ========================================== + ! This routine interpolates in time and then outputs a grid at + ! time=out_time + ! + ! files now have the same format as frames produced by outval + !======================================================================= + subroutine fgout_write(fgrid,out_time,out_index) + + use geoclaw_module, only: dry_tolerance + implicit none + + ! Subroutine arguments + type(fgout_grid), intent(inout) :: fgrid + real(kind=8), intent(in) :: out_time + integer, intent(in) :: out_index + + ! I/O + integer, parameter :: unit = 87 + character(len=15) :: fg_filename + character(len=4) :: cfgno, cframeno + character(len=8) :: file_format + integer :: grid_number,ipos,idigit,out_number,columns + integer :: ifg,ifg1, iframe,iframe1 + + integer, parameter :: method = 0 ! interpolate in time? + + ! Output format strings + ! These are now the same as in outval for frame data, for compatibility + ! For fgout grids there is only a single grid (ngrids=1) + ! and we set AMR_level=0, naux=0, nghost=0 (so no extra cells in binary) + + character(len=*), parameter :: header_format = & + "(i6,' grid_number',/," // & + "i6,' AMR_level',/," // & + "i6,' mx',/," // & + "i6,' my',/" // & + "e26.16,' xlow', /, " // & + "e26.16,' ylow', /," // & + "e26.16,' dx', /," // & + "e26.16,' dy',/)" + + character(len=*), parameter :: t_file_format = "(e18.8,' time', /," // & + "i6,' meqn'/," // & + "i6,' ngrids'/," // & + "i6,' naux'/," // & + "i6,' ndim'/," // & + "i6,' nghost'/," // & + "a10,' format'/,/)" + + ! Other locals + integer :: i,j,m + real(kind=8) :: t0,tf,tau, qaug(6) + real(kind=8), allocatable :: qeta(:,:,:) + real(kind=4), allocatable :: qeta4(:,:,:) + real(kind=8) :: h_early,h_late,topo_early,topo_late + + allocate(qeta(4, fgrid%mx, fgrid%my)) ! to store h,hu,hv,eta + + + ! Interpolate the grid in time, to the output time, using + ! the solution in fgrid1 and fgrid2, which represent the + ! solution on the fixed grid at the two nearest computational times + do j=1,fgrid%my + do i=1,fgrid%mx + + ! Check for small numbers + forall(m=1:fgrid%num_vars(1)-1,abs(fgrid%early(m,i,j)) < 1d-90) + fgrid%early(m,i,j) = 0.d0 + end forall + + if (method == 0) then + + ! no interpolation in time, use solution from full step: + qaug = fgrid%early(:,i,j) + + ! note that CFL condition ==> waves can't move more than 1 + ! cell per time step on each level, so solution from nearest + ! full step should be correct to within a cell width + ! Better to use early than late since for refinement tracking + ! wave moving out into still water. + + else if (method == 1) then + + ! interpolate in time. May have problems near shore? + + ! Fetch times for interpolation, this is done per grid point + ! since each grid point may come from a different source + t0 = fgrid%early(fgrid%num_vars(1),i,j) + tf = fgrid%late(fgrid%num_vars(1),i,j) + tau = (out_time - t0) / (tf - t0) + + ! check for small values: + forall(m=1:fgrid%num_vars(1)-1,abs(fgrid%late(m,i,j)) < 1d-90) + fgrid%late(m,i,j) = 0.d0 + end forall + + ! linear interpolation: + qaug = (1.d0-tau)*fgrid%early(:,i,j) + tau*fgrid%late(:,i,j) + + ! If resolution changed between early and late time, may be + ! problems near shore when interpolating B, h, eta + ! separately (at least in case when B changed and point + ! was dry at one time and wet the other). + ! Switch back to fgrid%early values, only in this case. + ! This is implemented below but not extensively tested. + + if (qaug(1) > 0.d0) then + topo_early = fgrid%early(4,i,j) + topo_late = fgrid%late(4,i,j) + if (topo_early .ne. topo_late) then + ! AMR resolution changed between times + h_early = fgrid%early(1,i,j) + h_late = fgrid%late(1,i,j) + if (((h_early < dry_tolerance) & + .and. (h_late >= dry_tolerance)) & + .or. ((h_late < dry_tolerance) & + .and. (h_early >= dry_tolerance))) then + ! point changed between wet and dry + qaug = fgrid%early(:,i,j) ! don't interpolate + endif + endif + endif + endif + + ! Output the conserved quantities and eta value + qeta(1,i,j) = qaug(1) ! h + qeta(2,i,j) = qaug(2) ! hu + qeta(3,i,j) = qaug(3) ! hv + qeta(4,i,j) = qaug(5) ! eta + + enddo + enddo + + + ! Make the file names and open output files + cfgno = '0000' + ifg = fgrid%fgno + ifg1 = ifg + do ipos=4,1,-1 + idigit = mod(ifg1,10) + cfgno(ipos:ipos) = char(ichar('0') + idigit) + ifg1 = ifg1/10 + enddo + + cframeno = '0000' + iframe = out_index + iframe1 = iframe + do ipos=4,1,-1 + idigit = mod(iframe1,10) + cframeno(ipos:ipos) = char(ichar('0') + idigit) + iframe1 = iframe1/10 + enddo + + fg_filename = 'fgout' // cfgno // '.q' // cframeno + + open(unit,file=fg_filename,status='unknown',form='formatted') + + ! Determine number of columns that will be written out + columns = fgrid%num_vars(1) - 1 + if (fgrid%num_vars(2) > 1) then + columns = columns + 2 + endif + + !write(6,*) '+++ fgout out_time = ',out_time + !write(6,*) '+++ fgrid%num_vars: ',fgrid%num_vars(1),fgrid%num_vars(2) + + ! Write out header in .q file: + !write(unit,header_format) out_time,fgrid%mx,fgrid%my, & + ! fgrid%x_low,fgrid%y_low, fgrid%x_hi,fgrid%y_hi, columns + + write(unit,header_format) fgrid%fgno, 0, fgrid%mx,fgrid%my, & + fgrid%x_low,fgrid%y_low, fgrid%dx, fgrid%dy + + if (fgrid%output_format == 1) then + ! ascii output added to .q file: + do j=1,fgrid%my + do i=1,fgrid%mx + write(unit, "(50e26.16)") qeta(1,i,j),qeta(2,i,j), & + qeta(3,i,j),qeta(4,i,j) + enddo + write(unit,*) ' ' ! blank line required between rows + enddo + endif + + close(unit) + + if (fgrid%output_format == 3) then + ! binary output goes in .b file as full 8-byte (float64): + fg_filename = 'fgout' // cfgno // '.b' // cframeno + open(unit=unit, file=fg_filename, status="unknown", & + access='stream') + write(unit) qeta + close(unit) + else if (fgrid%output_format == 2) then + ! binary output goes in .b file as 4-byte (float32): + fg_filename = 'fgout' // cfgno // '.b' // cframeno + open(unit=unit, file=fg_filename, status="unknown", & + access='stream') + allocate(qeta4(4, fgrid%mx, fgrid%my)) ! for 4-byte binary output + qeta4 = real(qeta, kind=4) + write(unit) qeta4 + deallocate(qeta4) + close(unit) + endif + + deallocate(qeta) + + ! time info .t file: + + + if (fgrid%output_format == 1) then + file_format = 'ascii' + else if (fgrid%output_format == 2) then + file_format = 'binary32' + else if (fgrid%output_format == 3) then + file_format = 'binary64' + else + write(6,*) '*** unrecognized fgrid%output_format = ', & + fgrid%output_format + write(6,*) '*** should be ascii, binary32, or binary64' + endif + + fg_filename = 'fgout' // cfgno // '.t' // cframeno + open(unit=unit, file=fg_filename, status='unknown', form='formatted') + ! time, num_eqn+1, num_grids, num_aux, num_dim, num_ghost: + write(unit, t_file_format) out_time, 4, 1, 0, 2, 0, file_format + close(unit) + + print "(a,i4,a,i4,a,e18.8)",'Writing fgout grid #',fgrid%fgno, & + ' frame ',out_index,' at time =',out_time + + ! Index into qeta for binary output + ! Note that this implicitly assumes that we are outputting only h, hu, hv + ! and will not output more (change num_eqn parameter above) + + end subroutine fgout_write + + + ! ========================================================================= + ! Utility functions for this module + ! Returns back a NaN + + real(kind=8) function nan() + real(kind=8) dnan + integer inan(2) + equivalence (dnan,inan) + inan(1)=2147483647 + inan(2)=2147483647 + nan=dnan + end function nan + + ! Interpolation function (in space) + ! Given 4 points (points) and geometry from x,y,and cross terms + + real(kind=8) pure function interpolate(points,geometry) result(interpolant) + + implicit none + + ! Function signature + real(kind=8), intent(in) :: points(2,2) + real(kind=8), intent(in) :: geometry(4) + integer :: icell, jcell + + ! pw bilinear + ! This is set up as a dot product between the approrpriate terms in + ! the input data. This routine could be vectorized or a BLAS routine + ! used instead of the intrinsics to ensure that the fastest routine + ! possible is being used + interpolant = sum([points(2,1)-points(1,1), & + points(1,2)-points(1,1), & + points(1,1) + points(2,2) - (points(2,1) + points(1,2)), & + points(1,1)] * geometry) + + end function interpolate + + +end module fgout_module diff --git a/src/2d/shallow/fixedgrids_module.f90 b/src/2d/shallow/fixedgrids_module.f90 deleted file mode 100644 index 600b9e51b..000000000 --- a/src/2d/shallow/fixedgrids_module.f90 +++ /dev/null @@ -1,560 +0,0 @@ -module fixedgrids_module - - implicit none - save - - ! Container for fixed grid data, geometry and output settings - type fixedgrid_type - ! Grid data - real(kind=8), pointer :: early(:,:,:) - real(kind=8), pointer :: late(:,:,:) - real(kind=8), pointer :: often(:,:,:) - - ! Geometry - integer :: num_vars(2),mx,my - real(kind=8) :: dx,dy,x_low,x_hi,y_low,y_hi - - ! Time Tracking and output types - integer :: num_output,last_output_index - integer :: output_arrival_times,output_surface_max - real(kind=8) :: last_output_time,start_time,end_time,dt - end type fixedgrid_type - - - logical, private :: module_setup = .false. - - ! Fixed grid arrays and sizes - integer :: num_fixed_grids - type(fixedgrid_type), allocatable :: fgrids(:) - real(kind=8) :: tcfmax - -contains - - ! Setup routine that reads in the fixed grids data file and sets up the - ! appropriate data structures - subroutine set_fixed_grids(fname) - - use amr_module, only: parmunit - - implicit none - - ! Subroutine arguments - character(len=*), optional, intent(in) :: fname - - ! Local storage - integer, parameter :: unit = 7 - integer :: i - - if (.not.module_setup) then - - write(parmunit,*) ' ' - write(parmunit,*) '--------------------------------------------' - write(parmunit,*) 'SETFIXEDGRIDS:' - write(parmunit,*) '-----------' - - ! Open data file - if (present(fname)) then - call opendatafile(unit,fname) - else - call opendatafile(unit,'fixed_grids.data') - endif - - ! Read in data - read(unit,'(i2)') num_fixed_grids - write(parmunit,*) ' mfgrids = ',num_fixed_grids - if (num_fixed_grids == 0) then - write(parmunit,*) ' No fixed grids specified for output' - return - endif - - ! Allocate fixed grids (not the data yet though) - allocate(fgrids(num_fixed_grids)) - - ! Read in data for each fixed grid - do i=1,num_fixed_grids - ! Read in this grid's data - read(unit,*) fgrids(i)%start_time, & - fgrids(i)%end_time, & - fgrids(i)%num_output, & - fgrids(i)%x_low, & - fgrids(i)%x_hi , & - fgrids(i)%y_low, & - fgrids(i)%y_hi , & - fgrids(i)%mx , & - fgrids(i)%my , & - fgrids(i)%output_arrival_times, & - fgrids(i)%output_surface_max - - ! Setup data for this grid - ! Set dtfg (the timestep length between outputs) for each grid - if (fgrids(i)%end_time <= fgrids(i)%start_time) then - if (fgrids(i)%num_output > 1) then - print *,'SETFIXEDGRIDS: ERROR for fixed grid', i - print *,'start_time <= end_time yet num_output > 1' - print *,'set end_time > start_time or set num_output = 1' - stop - else - fgrids(i)%dt = 0.d0 - endif - else - if (fgrids(i)%num_output < 2) then - print *,'SETFIXEDGRIDS: ERROR for fixed grid', i - print *,'end_time > start_time, yet num_output = 1' - print *,'set num_output > 2' - stop - else - fgrids(i)%dt = (fgrids(i)%end_time - fgrids(i)%start_time) & - / (fgrids(i)%num_output - 1) - endif - endif - - ! Initialize last_output_time and index - fgrids(i)%last_output_time = fgrids(i)%start_time - fgrids(i)%dt - fgrids(i)%last_output_index = 0 - - ! Set spatial intervals dx and dy on each grid - if (fgrids(i)%mx > 1) then - fgrids(i)%dx = (fgrids(i)%x_hi - fgrids(i)%x_low) / (fgrids(i)%mx - 1) - else if (fgrids(i)%mx == 1) then - fgrids(i)%dx = 0.d0 - else - print *,'SETFIXEDGRIDS: ERROR for fixed grid', i - print *,'x grid points mx <= 0, set mx >= 1' - endif - - if (fgrids(i)%my > 1) then - fgrids(i)%dy = (fgrids(i)%y_hi - fgrids(i)%y_low) / (fgrids(i)%my - 1) - else if (fgrids(i)%my == 1) then - fgrids(i)%dy = 0.d0 - else - print *,'SETFIXEDGRIDS: ERROR for fixed grid', i - print *,'y grid points my <= 0, set my >= 1' - endif - - ! set the number of variables stored for each grid - ! this should be (the number of variables you want to write out + 1) - fgrids(i)%num_vars(1) = 6 - fgrids(i)%num_vars(2) = 3*fgrids(i)%output_surface_max & - + fgrids(i)%output_arrival_times - - ! Allocate new fixed grid data array - allocate(fgrids(i)%early(fgrids(i)%num_vars(1),fgrids(i)%mx,fgrids(i)%my)) - fgrids(i)%early = nan() - allocate(fgrids(i)%late(fgrids(i)%num_vars(1),fgrids(i)%mx,fgrids(i)%my)) - fgrids(i)%late = nan() - allocate(fgrids(i)%often(fgrids(i)%num_vars(2),fgrids(i)%mx,fgrids(i)%my)) - fgrids(i)%often = nan() - enddo - close(unit) - - tcfmax=-1.d16 - - module_setup = .true. - end if - - end subroutine set_fixed_grids - - - !=====================FGRIDINTERP======================================= - ! This routine interpolates q and aux on a computational grid - ! to a fgrid not necessarily aligned with the computational grid - ! using bilinear interpolation defined on computation grid - !======================================================================= - subroutine fgrid_interp(fgrid_type,fgrid, & - t,q,meqn,mxc,myc,mbc,dxc,dyc,xlowc,ylowc, & - maux,aux,maxcheck) - - use geoclaw_module, only: dry_tolerance - implicit none - - ! Subroutine arguments - integer, intent(in) :: fgrid_type - type(fixedgrid_type), intent(inout) :: fgrid - integer, intent(in) :: meqn,mxc,myc,mbc,maux,maxcheck - real(kind=8), intent(in) :: t,dxc,dyc,xlowc,ylowc - real(kind=8), intent(in) :: q(meqn,1-mbc:mxc+mbc,1-mbc:myc+mbc) - real(kind=8), intent(in) :: aux(maux,1-mbc:mxc+mbc,1-mbc:myc+mbc) - - ! Indices - integer :: ifg,jfg,m,ic1,ic2,jc1,jc2 - integer :: bathy_index,eta_index,arrival_index - integer :: eta_min_index,eta_max_index,eta_now_index - - ! Tolerances - real(kind=8), parameter :: arrival_tolerance = 1.d-2 - real(kind=8) :: total_depth,depth_indicator,nan_check - - ! Geometry - real(kind=8) :: xfg,yfg,xc1,xc2,yc1,yc2,xhic,yhic - real(kind=8) :: geometry(4) - - ! Work arrays for eta interpolation - real(kind=8) :: eta(2,2),h(2,2) - - - ! Alias to data in fixed grid - integer :: num_vars - real(kind=8), pointer :: fg_data(:,:,:) - - ! Setup aliases for specific fixed grid - if (fgrid_type == 1) then - num_vars = fgrid%num_vars(1) - fg_data => fgrid%early - else if (fgrid_type == 2) then - num_vars = fgrid%num_vars(1) - fg_data => fgrid%late - else - num_vars = fgrid%num_vars(2) - fg_data => fgrid%often - endif - - xhic = xlowc + dxc*mxc - yhic = ylowc + dyc*myc - - ! Find indices of various quantities in the fgrid arrays - bathy_index = meqn + 1 - eta_index = meqn + 2 - - if (maxcheck > 0) then - arrival_index = 0 - eta_now_index = 0 - eta_min_index = 0 - eta_max_index = 0 - - if (fgrid%output_arrival_times > 0) then - arrival_index = 1 - endif - - if (fgrid%output_surface_max > 0) then - eta_now_index = arrival_index + 1 - eta_min_index = arrival_index + 2 - eta_max_index = arrival_index + 3 - endif - endif - - ! Primary interpolation loops - do ifg=1,fgrid%mx - xfg=fgrid%x_low + (ifg-1)*fgrid%dx - do jfg=1,fgrid%my - yfg=fgrid%y_low + (jfg-1)*fgrid%dy - - ! Check to see if this coordinate is inside of this grid - if (.not.((xfg < xlowc.or.xfg > xhic).or.(yfg < ylowc.or.yfg > yhic))) then - - ! find where xfg,yfg is in the computational grid and compute the indices - ! and relevant coordinates of each corner - ic1 = int((xfg-(xlowc+0.5d0*dxc))/(dxc))+1 - jc1 = int((yfg-(ylowc+0.5d0*dyc))/(dyc))+1 - if (ic1.eq.mxc) ic1=mxc-1 - if (jc1.eq.myc) jc1=myc-1 - ic2 = ic1 + 1 - jc2 = jc1 + 1 - - xc1 = xlowc + dxc * (ic1 - 0.5d0) - yc1 = ylowc + dyc * (jc1 - 0.5d0) - xc2 = xlowc + dxc * (ic2 - 0.5d0) - yc2 = ylowc + dyc * (jc2 - 0.5d0) - - ! Calculate geometry of interpolant - ! interpolate bilinear used to interpolate to xfg,yfg - ! define constant parts of bilinear - geometry = [(xfg - xc1) / dxc, & - (yfg - yc1) / dyc, & - (xfg - xc1) * (yfg - yc1) / (dxc*dyc), & - 1.d0] - - ! Interpolate for all conserved quantities and bathymetry - if (maxcheck == 0) then - forall (m=1:meqn) - fg_data(m,ifg,jfg) = interpolate([[q(m,ic1,jc1),q(m,ic1,jc2)], & - [q(m,ic2,jc1),q(m,ic2,jc2)]], geometry) - end forall - fg_data(bathy_index,ifg,jfg) = interpolate([[aux(1,ic1,jc1),aux(1,ic1,jc2)], & - [aux(1,ic2,jc1),aux(1,ic2,jc2)]], geometry) - endif - - ! If eta max/min are saved on this grid initialize if necessary - if (fgrid%output_surface_max > 0 .and. maxcheck == 2) then - if (.not.(fg_data(eta_min_index,ifg,jfg) == fg_data(eta_min_index,ifg,jfg))) then - fg_data(eta_min_index,ifg,jfg) = 0.d0 - endif - if (.not.(fg_data(eta_max_index,ifg,jfg) == fg_data(eta_max_index,ifg,jfg))) then - fg_data(eta_max_index,ifg,jfg) = 0.d0 - endif - endif - - ! Interpolate surface eta, only use wet eta points near the shoreline - eta(1,:) = [aux(1,ic1,jc1) + q(1,ic1,jc1), aux(1,ic1,jc2) + q(1,ic1,jc2)] - eta(2,:) = [aux(1,ic2,jc1) + q(1,ic2,jc1), aux(1,ic2,jc2) + q(1,ic2,jc2)] - h(1,:) = [q(1,ic1,jc1),q(1,ic1,jc2)] - h(2,:) = [q(1,ic2,jc1),q(1,ic2,jc2)] - - depth_indicator= min(h(1,1),h(1,2),h(2,1),h(2,2)) - total_depth = sum(h) - - ! We are near shoreline - if (depth_indicator < dry_tolerance .and. & - total_depth > 4.d0 * dry_tolerance) then - ! Check to see if each cell around fixed grid point is - ! wet, if not re-balance - if (h(1,1) < dry_tolerance) then - eta(1,1) = (h(1,2)*eta(1,2) & - + h(2,1)*eta(2,1) & - + h(2,2)*eta(2,2)) / total_depth - endif - if (h(1,2) < dry_tolerance) then - eta(1,2) = (h(1,1)*eta(1,1) & - + h(2,1)*eta(2,1) & - + h(2,2)*eta(2,2)) / total_depth - endif - if (h(2,1) < dry_tolerance) then - eta(2,1) = (h(1,1)*eta(1,1) & - + h(1,2)*eta(1,2) & - + h(2,2)*eta(2,2)) / total_depth - endif - if (h(2,2) < dry_tolerance) then - eta(2,2)= (h(1,2)*eta(1,2) & - + h(2,1)*eta(2,1) & - + h(1,1)*eta(1,1)) / total_depth - endif - endif - if (total_depth <= 4.d0*dry_tolerance) then - eta(2,2) = nan() - endif - - ! Check which task to perform and evaluate the interpolant - ! or evaluate the eta min and max functions - if (maxcheck == 0) then - fg_data(eta_index,ifg,jfg) = interpolate(eta,geometry) - fg_data(num_vars,ifg,jfg) = t - else if (maxcheck.eq.1.and.fgrid%output_surface_max.gt.0) then - fg_data(eta_now_index,ifg,jfg) = interpolate(eta,geometry) - else if (maxcheck.eq.2.and.fgrid%output_surface_max.gt.0) then - fg_data(eta_min_index,ifg,jfg) = min(fg_data(eta_min_index,ifg,jfg), & - fg_data(eta_now_index,ifg,jfg)) - fg_data(eta_max_index,ifg,jfg) = max(fg_data(eta_max_index,ifg,jfg), & - fg_data(eta_now_index,ifg,jfg)) - endif - - ! If arrival times are saved on this grid - if (maxcheck == 1 .and. fgrid%output_arrival_times > 0) then - ! TODO: It would be nice to replace this with an - ! intrinsic such as ieee_is_nan but this is not widely - ! implemented yet - nan_check = fg_data(arrival_index,ifg,jfg) - ! if nan_check = NaN: Waves haven't arrived previously - if (.not.(nan_check == nan_check)) then - if (abs(fg_data(eta_index,ifg,jfg)) > arrival_tolerance) then - fg_data(arrival_index,ifg,jfg)= t - endif - endif - endif - - endif ! Enclosing if statement to see if fixed grid point is - ! in this computational grid - enddo ! Fixed grid y-coordinate loop - enddo ! Fixed grid x-coordinte loop - - end subroutine fgrid_interp - - - !=====================FGRIDOUT========================================== - ! This routine interpolates in time and then outputs a grid at - ! time=toutfg - ! - ! files have a header, followed by columns of data - !======================================================================= - subroutine fgrid_out(grid_index,fgrid,out_time,out_index,out_flag) - - implicit none - - ! Subroutine arguments - type(fixedgrid_type), intent(inout) :: fgrid - real(kind=8), intent(in) :: out_time - integer, intent(in) :: grid_index,out_index, out_flag - - ! I/O - integer, parameter :: unit = 95 - character(len=30) :: fg_filename - integer :: grid_number,pos,digit,out_number,columns - - ! Out format strings - character(len=*), parameter :: header_format = "(e18.8,' time', /," // & - "i5,' mx', /," // & - "i5,' my', /," // & - "e18.8,' xlow',/" // & - "e18.8,' ylow',/" // & - "e18.8,' xhi',/," // & - "e18.8,' yhi',/," // & - "i5,' columns',/)" - character(len=*), parameter :: data_format = "(8e26.16)" - character(len=*), parameter :: arrival_header_format = & - "(i5,' mx', /," // & - "i5,' my', /," // & - "e18.8,' xlow',/" // & - "e18.8,' ylow',/" // & - "e18.8,' xhi',/," // & - "e18.8,' yhi',/," // & - "i5,' columns',/)" - - ! Other locals - integer :: i,j,m - integer :: eta_min_index,eta_max_index - real(kind=8) :: t0,tf,tau - - - ! Make the file names and open output files - fg_filename = 'fort.fgnn_xxxx' - grid_number= grid_index - do pos = 9, 8, -1 - digit = mod(grid_number,10) - fg_filename(pos:pos) = char(ichar('0') + digit) - grid_number = grid_number/ 10 - enddo - - out_number = out_index - do pos = 14, 11, -1 - digit = mod(out_number,10) - fg_filename(pos:pos) = char(ichar('0') + digit) - out_number = out_number / 10 - enddo - - open(unit,file=fg_filename,status='unknown',form='formatted') - - ! Determine number of columns that will be written out - columns = fgrid%num_vars(1) - 1 - if (fgrid%num_vars(2) > 1) then - columns = columns + 2 - endif - - ! Write out header - write(unit,header_format) out_time,fgrid%mx,fgrid%my,fgrid%x_low,fgrid%y_low,fgrid%x_hi,fgrid%y_hi,columns - - ! Surface max/min index - eta_min_index = fgrid%output_arrival_times + 2 - eta_max_index = fgrid%output_arrival_times + 3 - - ! Interpolate the grid in time, to the output time, using - ! the solution in fgrid1 and fgrid2, which represent the - ! solution on the fixed grid at the two nearest computational times - do j=1,fgrid%my - do i=1,fgrid%mx - ! Fetch times for interpolation, this is done per grid point - ! since each grid point may come from a different source - t0 = fgrid%early(fgrid%num_vars(1),i,j) - tf = fgrid%late(fgrid%num_vars(1),i,j) - tau = (out_time - t0) / (tf - t0) - - ! Check for small numbers - forall(m=1:fgrid%num_vars(1)-1,abs(fgrid%early(m,i,j)) < 1d-90) - fgrid%early(m,i,j) = 0.d0 - end forall - forall(m=1:fgrid%num_vars(1)-1,abs(fgrid%late(m,i,j)) < 1d-90) - fgrid%late(m,i,j) = 0.d0 - end forall - - ! Check which output form we are doing - if (columns == fgrid%num_vars(1) - 1) then - ! Output only the conserved quantities - write(unit,data_format) interpolate_time(fgrid%num_vars(1), & - fgrid%early(:,i,j), & - fgrid%late(:,i,j), & - tau) - else - ! Output min/max of eta as well as the conserved quantities - if (abs(fgrid%often(eta_min_index,i,j)) < 1d-90) then - fgrid%often(eta_min_index,i,j) = 0.d0 - endif - if (abs(fgrid%often(eta_max_index,i,j)) < 1d-90) then - fgrid%often(eta_max_index,i,j) = 0.d0 - endif - write(unit,data_format) interpolate_time(fgrid%num_vars(1), & - fgrid%early(:,i,j), & - fgrid%late(:,i,j), & - tau), & - fgrid%often(eta_min_index,i,j), & - fgrid%often(eta_max_index,i,j) - endif - enddo - enddo - - close(unit) - print "(a,i2,a,i2,a,e18.8)",' FGRIDOUT: Fixed Grid ',grid_index, ' frame ',out_index,' at time =',out_time - - ! ==================== Output for arrival times============ - if (out_flag == 1) then - ! Make the file name and open output file for arrival times - fg_filename = 'fort.fgnn_arrivaltimes' - grid_number= grid_index - do pos = 9, 8, -1 - digit= mod(grid_number,10) - fg_filename(pos:pos) = char(ichar('0') + digit) - grid_number = grid_number/ 10 - enddo - open(unit,file=fg_filename,status='unknown',form='formatted') - - write(95,arrival_header_format) fgrid%mx,fgrid%my,fgrid%x_low,fgrid%y_low,fgrid%x_hi,fgrid%y_hi - - do j=1,fgrid%my - do i=1,fgrid%mx - write(unit,"(1e26.16)") fgrid%often(1,i,j) - enddo - enddo - close(unit) - - print "(a,i2,a)", ' FGRIDOUT: Fixed Grid ', grid_index, ' arrival times output' - endif - - end subroutine fgrid_out - - ! ========================================================================= - ! Utility functions for this module - ! Returns back a NaN - real(kind=8) function nan() - real(kind=8) dnan - integer inan(2) - equivalence (dnan,inan) - inan(1)=2147483647 - inan(2)=2147483647 - nan=dnan - end function nan - - ! Interpolation function (in space) - ! Given 4 points (points) and geometry from x,y,and cross terms - real(kind=8) pure function interpolate(points,geometry) result(interpolant) - - implicit none - - ! Function signature - real(kind=8), intent(in) :: points(2,2) - real(kind=8), intent(in) :: geometry(4) - - ! This is set up as a dot product between the approrpriate terms in - ! the input data. This routine could be vectorized or a BLAS routine - ! used instead of the intrinsics to ensure that the fastest routine - ! possible is being used - interpolant = sum([points(2,1)-points(1,1), & - points(1,2)-points(1,1), & - points(1,1) + points(2,2) - (points(2,1) + points(1,2)), & - points(1,1)] * geometry) - - end function interpolate - - ! Interpolation function in time - pure function interpolate_time(num_vars,early,late,tau) result(interpolant) - - implicit none - - ! Input arguments - integer, intent(in) :: num_vars - real(kind=8), intent(in) :: early(num_vars),late(num_vars),tau - - ! Return value - real(kind=8) :: interpolant(num_vars) - - interpolant = (1.d0 - tau) * early(:) + tau * late(:) - - end function interpolate_time - - -end module fixedgrids_module diff --git a/src/2d/shallow/multilayer/Makefile.multilayer b/src/2d/shallow/multilayer/Makefile.multilayer index bd399b709..ca23bbb72 100644 --- a/src/2d/shallow/multilayer/Makefile.multilayer +++ b/src/2d/shallow/multilayer/Makefile.multilayer @@ -50,7 +50,7 @@ COMMON_MODULES += \ $(GEOLIB)/topo_module.f90 \ $(GEOLIB)/qinit_module.f90 \ $(GEOLIB)/refinement_module.f90 \ - $(GEOLIB)/fixedgrids_module.f90 \ + $(GEOLIB)/fgout_module.f90 \ $(GEOLIB)/fgmax_module.f90 \ $(GEOLIB)/friction_module.f90 \ $(GEOLIB)/adjointsup_module.f90 diff --git a/src/2d/shallow/multilayer/stepgrid.f b/src/2d/shallow/multilayer/stepgrid.f index 81d6a1c00..c0ac2c009 100644 --- a/src/2d/shallow/multilayer/stepgrid.f +++ b/src/2d/shallow/multilayer/stepgrid.f @@ -2,7 +2,7 @@ c ------------------------------------------------------------- c subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, - & nvar,xlow,ylow,time,mptr,maux,aux) + & nvar,xlow,ylow,time,mptr,maux,aux,actualstep) c c c ::::::::::::::::::: STEPGRID :::::::::::::::::::::::::::::::::::: @@ -23,23 +23,21 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, c c c This version of stepgrid, stepgrid_geo.f allows output on -c fixed grids specified in setfixedgrids.data +c fgout grids specified in fgout_grids.data c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - use amr_module - use geoclaw_module, only: grav, rho + use amr_module + use fgout_module, only: FGOUT_num_grids, FGOUT_fgrids, + & FGOUT_tcfmax, fgout_interp, fgout_grid, + & FGOUT_ttol use multilayer_module, only: num_layers, dry_tolerance - - use fixedgrids_module implicit double precision (a-h,o-z) external rpn2,rpt2 -c parameter (msize=max1d+4) -c parameter (mwork=msize*(maxvar*maxvar + 13*maxvar + 3*maxaux +2)) dimension q(nvar,mitot,mjtot) dimension fp(nvar,mitot,mjtot),gp(nvar,mitot,mjtot) @@ -49,17 +47,22 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, logical :: debug = .false. logical :: dump = .false. + logical, intent (in) :: actualstep + type(fgout_grid), pointer :: fgout + logical, allocatable :: fgout_interp_needed(:) + real(kind=8) :: fgout_tnext c - tcfmax = -rinfinity + FGOUT_tcfmax = -rinfinity level = node(nestlevel,mptr) + if (dump) then write(outunit,*)" at start of stepgrid: dumping grid ",mptr do i = 1, mitot do j = 1, mjtot - write(outunit,545) i,j,(q(ivar,i,j),ivar=1,nvar) -c write(*,545) i,j,(q(ivar,i,j),ivar=1,nvar) - 545 format(2i4,4e15.7) + write(outunit,545) i,j,(q(ivar,i,j),ivar=1,nvar), + . (aux(iaux,i,j),iaux=1,maux) + 545 format(2i4,4e15.7,/,8x,4e15.7) end do end do endif @@ -110,104 +113,61 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, C mwork1 = mwork - mused !# remaining space (passed to step2) c - -c::::::::::::::::::::::::Fixed Grid Output::::::::::::::::::::::::::::::::: tc0=time !# start of computational step tcf=tc0+dt !# end of computational step + +c Check if fgout interpolation needed before and after step: + + allocate(fgout_interp_needed(FGOUT_num_grids)) !$OMP CRITICAL (FixedGrids) -c # see if any f-grids should be written out - do ng=1,num_fixed_grids - if (tc0 > fgrids(ng)%start_time .and. - & fgrids(ng)%last_output_index < fgrids(ng)%num_output) then -c # fgrid ng may need to be written out -c # find the first output number that has not been written out and -c # find the first output number on a fixed grid that is >= tc0 -c # which will not be written out - if (fgrids(ng)%dt > 0.d0) then - ioutfgend= 1+max(0,nint((tc0 - fgrids(ng)%start_time) - & / fgrids(ng)%dt)) - else - ioutfgend=1 - endif - ioutfgend = min(ioutfgend,fgrids(ng)%num_output) - ioutfgstart = fgrids(ng)%last_output_index + 1 -c # write-out fgrid times that are less than tc0, and have not been written yet -c # these should be the most accurate values at any given point in the fgrid -c # since tc0> output time - do ioutfg=ioutfgstart,ioutfgend - toutfg=fgrids(ng)%start_time+(ioutfg-1)*fgrids(ng)%dt - if (toutfg < tc0) then -c # write out the solution for fixed grid ng -c # test if arrival times should be output - ioutflag = fgrids(ng)%output_arrival_times* - & (fgrids(ng)%num_output- - & fgrids(ng)%last_output_index) - call fgrid_out(ng,fgrids(ng),toutfg,ioutfg,ioutflag) - - fgrids(ng)%last_output_time = toutfg - fgrids(ng)%last_output_index = - & fgrids(ng)%last_output_index + 1 - endif - enddo - - endif + + do ng=1,FGOUT_num_grids + fgout => FGOUT_fgrids(ng) + if (fgout%next_output_index > fgout%num_output) then + fgout_interp_needed(ng) = .false. + else + fgout_tnext = fgout%output_times(fgout%next_output_index) + fgout_interp_needed(ng) = + & ((fgout%x_low < xlowmbc + mx * dx) .and. + & (fgout%x_hi > xlowmbc) .and. + & (fgout%y_low < ylowmbc + my * dy) .and. + & (fgout%y_hi > ylowmbc) .and. + & (fgout_tnext >= tc0 - FGOUT_ttol) .and. + & (fgout_tnext <= tcf + FGOUT_ttol)) + endif +c write(6,*) '+++ level, before- tc0,tcf,needed: ', +c & level,tc0,tcf,fgout_interp_needed(ng) +c write(6,*) '+++ next index: ',fgout%next_output_index +c write(6,*) '+++ fgout_tnext, tcf + FGOUT_ttol: ', +c & fgout_tnext, tcf + FGOUT_ttol enddo !$OMP END CRITICAL (FixedGrids) + + +c::::::::::::::::::::::::fgout Output::::::::::::::::::::::::::::::::: +c This has been moved to tick.f, after advancing all patches on +c finest level. No need to check on each patch separately. c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: call b4step2(mbc,mx,my,nvar,q, - & xlowmbc,ylowmbc,dx,dy,time,dt,maux,aux) + & xlowmbc,ylowmbc,dx,dy,time,dt,maux,aux,actualstep) -c::::::::::::::::::::::::FIXED GRID DATA before step::::::::::::::::::::::: -c # fill in values at fixed grid points effected at time tc0 -!$OMP CRITICAL (FixedGrids) - do ng=1,num_fixed_grids - - if ( (fgrids(ng)%x_low < xlowmbc + mx*dx) .and. - & (fgrids(ng)%x_hi > xlowmbc) .and. - & (fgrids(ng)%y_low < ylowmbc + my*dy) .and. - & (fgrids(ng)%y_hi > ylowmbc) .and. - & (fgrids(ng)%last_output_index < fgrids(ng)%num_output) .and. - & (tcf >= fgrids(ng)%start_time) ) then - - if (fgrids(ng)%last_output_time + fgrids(ng)%dt >= tc0 .and. - & fgrids(ng)%last_output_time + fgrids(ng)%dt <= tcf) then - -c # fixedgrid ng has an output time within [tc0,tcf] interval -c # and it overlaps this computational grid spatially - call fgrid_interp(1,fgrids(ng),tc0,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,0) - -c # routine to spatially interpolate computational solution -c # at tc0 to the fixed grid spatial points, -c #saving solution, variables and tc0 at every grid point - endif - -c # set maxima or minima if this is a new coarse step -c if (tc0.ge.tcfmax) then - -c # RJL: rewrote to set min/max every time a grid at level 1 -c # is about to be taken. The previous code failed if there was more than one grid -c # at level 1. Note that all grids are up to date at start of step on level 1. -c # New feature added at end of this routine to check more frequently if -c # levelcheck > 0. - if (level .eq. 1) then - if (fgrids(ng)%output_surface_max - & + fgrids(ng)%output_arrival_times > 0) then - - call fgrid_interp(3,fgrids(ng),tc0,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,2) - +c::::::::::::::::::::::::FGOUT DATA before step::::::::::::::::::::::: +c # fill in values at fgout points affected at time tc0 + do ng=1,FGOUT_num_grids + if (fgout_interp_needed(ng)) then +c # fgout grid ng has an output time within [tc0,tcf] interval +c # and it overlaps this computational grid spatially +!$OMP CRITICAL (FixedGrids) + fgout => FGOUT_fgrids(ng) +c write(6,*) '+++ fout_interp(1), tc0, level: ',tc0,level + call fgout_interp(1,fgout,tc0,q,nvar,mx,my,mbc, + & dx,dy,xlowmbc,ylowmbc,maux,aux) + +!$OMP END CRITICAL (FixedGrids) endif - endif - - endif enddo - tcfmax=max(tcfmax,tcf) - -!$OMP END CRITICAL (FixedGrids) - c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c New fixed grid stuff: Update fixed grid info from this patch... @@ -226,8 +186,8 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, c mptr_level = node(nestlevel,mptr) - write(outunit,811) mptr, mptr_level, cflgrid - 811 format(" Courant # of grid ",i5," level",i3," is ",d12.4) +c write(outunit,811) mptr, mptr_level, cflgrid +c811 format(" Courant # of grid ",i5," level",i3," is ",d12.4) c !$OMP CRITICAL (cflm) @@ -271,12 +231,14 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, c 50 continue c c # Copied here from b4step2 since need to do before saving to qc1d: +c # (This is the only place there's a difference for multilayer case) forall(i=1:mitot, j=1:mjtot, k=1:num_layers, & q(3*(k-1)+1,i,j) / rho(k) < dry_tolerance(k)) q(3*(k-1)+1,i,j) = max(q(3*(k-1)+1,i,j), 0.d0) q(3*(k-1)+2,i,j) = 0.d0 q(3*(k-1)+3,i,j) = 0.d0 end forall + c if (method(5).eq.1) then c # with source term: use Godunov splitting @@ -284,66 +246,25 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, & q,maux,aux,time,dt) endif -!$OMP CRITICAL (FixedGrids) -c ::::::::::::::::::::::::Fixed Grid data afterstep::::::::::::::::::::::: -c # fill in values at fixed grid points effected at time tcf - do ng=1,num_fixed_grids - if ((fgrids(ng)%x_low < xlowmbc + mx * dx) .and. - & (fgrids(ng)%x_hi > xlowmbc) .and. - & (fgrids(ng)%y_low < ylowmbc + my * dy) .and. - & (fgrids(ng)%y_hi > ylowmbc) .and. - & (fgrids(ng)%last_output_index < fgrids(ng)%num_output) .and. - & (tcf >= fgrids(ng)%start_time)) then - - if (fgrids(ng)%last_output_time + fgrids(ng)%dt >= tc0 .and. - & fgrids(ng)%last_output_time + fgrids(ng)%dt <= tcf) then - -c # fixedgrid ng has an output time within [tc0,tcf] interval -c # and it overlaps this computational grid spatially -C i0=i0fg(ng) !# index into the ng grid in the work array - - call fgrid_interp(2,fgrids(ng),tcf,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,0) - -c # routine to interpolate solution -c # at tcf to the fixed grid storage array, -c #saving solution and tcf at every grid point - - endif - -c # fill in values for eta if they need to be saved for later checking max/mins -c # check for arrival times - if (fgrids(ng)%output_surface_max - & + fgrids(ng)%output_arrival_times > 0) then - - call fgrid_interp(3,fgrids(ng),tc0,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,1) - - endif - -c # RJL: Modified 8/20/11 -c # If levelcheck > 0 then update max/mins at end of step on this grid. -c # Note that if there are finer grids then fgridoften will not have been updated -c # properly yet by those grids. This modification allows checking max/min more -c # frequently than the original code (equivalent to levelcheck==0) when you know -c # what level is most relevant for this fixed grid. Note also that if there are no -c # grids at levelcheck overlapping a portion of the fixed grid then the max/min values -c # will be updated only at start of next level 1 step. - - levelcheck = 0 - if (level == levelcheck) then - if (fgrids(ng)%output_arrival_times - & + fgrids(ng)%output_surface_max > 0) then - - call fgrid_interp(3,fgrids(ng),tc0,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,2) - endif +c ::::::::::::::::::::::::fgout data afterstep::::::::::::::::::::::: +c # fill in values at fgout points affected at time tcf + do ng=1,FGOUT_num_grids + if (fgout_interp_needed(ng)) then +c # fgout grid ng has an output time within [tc0,tcf] interval +c # and it overlaps this computational grid spatially +!$OMP CRITICAL (FixedGrids) + fgout => FGOUT_fgrids(ng) +c write(6,*) '+++ fout_interp(2), tcf, level: ',tcf,level + call fgout_interp(2,fgout,tcf,q,nvar,mx,my,mbc, + & dx,dy,xlowmbc,ylowmbc,maux,aux) + +!$OMP END CRITICAL (FixedGrids) endif - - - endif +c write(6,*) '+++ level,after- tc0,tcf,needed: ', +c & level,tc0,tcf,fgout_interp_needed(ng) +c write(6,*) '+++ next index: ',fgout%next_output_index +c write(6,*) '+++ fgout_tnext: ',fgout_tnext enddo -!$OMP END CRITICAL (FixedGrids) c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -399,3 +320,4 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, return end + diff --git a/src/2d/shallow/multilayer/valout.f90 b/src/2d/shallow/multilayer/valout.f90 index 063a569f5..241b5b24a 100644 --- a/src/2d/shallow/multilayer/valout.f90 +++ b/src/2d/shallow/multilayer/valout.f90 @@ -12,7 +12,7 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) use amr_module, only: node, rnode, ndilo, ndihi, ndjlo, ndjhi use amr_module, only: cornxlo, cornylo, levelptr, mxnest use amr_module, only: timeValout, timeValoutCPU, tvoll, tvollCPU, rvoll - use amr_module, only: timeTick, tick_clock_start, t0 + use amr_module, only: timeTick, tick_clock_start, t0, timeTickCPU use storm_module, only: storm_specification_type, output_storm_location use storm_module, only: output_storm_location @@ -34,6 +34,7 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) ! Locals logical :: timing_file_exists integer, parameter :: out_unit = 50 + integer, parameter :: binary_unit = 51 integer :: i, j, k, m, level, output_aux_num, num_stop, digit integer :: index, grid_ptr, num_cells(2), num_grids, q_loc, aux_loc real(kind=8) :: lower_corner(2), delta(2) @@ -43,6 +44,8 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) real(kind=8) :: h(num_layers), hu(num_layers), hv(num_layers) real(kind=8) :: eta(num_layers) real(kind=8), allocatable :: qeta(:) + real(kind=4), allocatable :: qeta4(:), aux4(:) + integer :: lenaux4 #ifdef HDF5 ! HDF writing @@ -52,12 +55,13 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) #endif ! Timing - integer :: clock_start, clock_finish, clock_rate - integer tick_clock_finish, tick_clock_rate, timeTick_int + integer(kind=8) :: clock_start, clock_finish, clock_rate + integer(kind=8) :: tick_clock_finish, tick_clock_rate, timeTick_int real(kind=8) :: cpu_start, cpu_finish, t_CPU_overall, timeTick_overall character(len=128) :: console_format character(len=512) :: timing_line, timing_substr character(len=*), parameter :: timing_file_name = "timing.csv" + character(len=8) :: file_format character(len=*), parameter :: header_format = & "(i6,' grid_number',/," // & @@ -73,7 +77,9 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) "i6,' ngrids'/," // & "i6,' naux'/," // & "i6,' ndim'/," // & - "i6,' nghost'/,/)" + "i6,' nghost'/," // & + "a10,' format'/,/)" + ! Output timing call system_clock(clock_start,clock_rate) @@ -120,8 +126,8 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) ! Write out fort.q file (and fort.bXXXX and clawxxxx.h5 files if necessary) ! Here we let fort.q be out_unit and the the other two be out_unit + 1 open(unit=out_unit, file=file_name(1), status='unknown', form='formatted') - if (output_format == 3) then - open(unit=out_unit + 1, file=file_name(4), status="unknown", & + if (output_format == 2 .or. output_format == 3) then + open(unit=binary_unit, file=file_name(4), status="unknown", & access='stream') else if (output_format == 4) then #ifdef HDF5 @@ -158,18 +164,19 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) delta(1), delta(2) ! Output grids - select case(output_format) - ! ASCII output - case(1) - ! Round off if nearly zero - forall (m = 1:num_eqn, & - i=num_ghost + 1:num_cells(1) + num_ghost, & - j=num_ghost + 1:num_cells(2) + num_ghost, & - abs(alloc(iadd(m, i, j))) < 1d-90) - - alloc(iadd(m, i, j)) = 0.d0 - end forall - + + ! Round off if nearly zero + ! (Do this for all output_format's) + forall (m = 1:num_eqn, & + i=num_ghost + 1:num_cells(1) + num_ghost, & + j=num_ghost + 1:num_cells(2) + num_ghost, & + abs(alloc(iadd(m, i, j))) < 1d-90) + + alloc(iadd(m, i, j)) = 0.d0 + end forall + + if (output_format == 1) then + ! ascii output do j = num_ghost + 1, num_cells(2) + num_ghost do i = num_ghost + 1, num_cells(1) + num_ghost @@ -198,52 +205,69 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) write(out_unit, *) ' ' end do - ! What is case 2? - case(2) - stop "Unknown format." - - ! Binary output - case(3) - ! Need to add eta to the output data - allocate(qeta((num_eqn + num_layers) & - * (num_cells(1) + 2 * num_ghost) & - * (num_cells(2) + 2 * num_ghost))) - do j = 1, num_cells(2) + 2 * num_ghost - do i = 1, num_cells(1) + 2 * num_ghost - - ! Extract depth and momenta - do k=1,num_layers - index = 3 * (k - 1) - h(k) = alloc(iadd(index + 1,i,j)) / rho(k) - hu(k)= alloc(iadd(index + 2,i,j)) / rho(k) - hv(k) = alloc(iadd(index + 3,i,j)) / rho(k) - end do - - ! Calculate surfaces - eta(num_layers) = h(num_layers) & - + alloc(iaddaux(1,i,j)) - do k=num_layers-1,1,-1 - eta(k) = h(k) + eta(k+1) - enddo - - do m=1,num_layers - index = 3*(m - 1) - qeta(iaddqeta(index+1,i,j)) = h(m) - qeta(iaddqeta(index+2,i,j)) = hu(m) - qeta(iaddqeta(index+3,i,j)) = hv(m) - qeta(iaddqeta(3*num_layers+m,i,j)) = eta(m) - end do - + else if (output_format==2 .or. output_format==3) then + ! binary32 or binary64 + + ! Note: We are writing out ghost cell data also, + ! so need to update this + call bound(time,num_eqn,num_ghost,alloc(q_loc), & + num_cells(1) + 2*num_ghost, & + num_cells(2) + 2*num_ghost, & + grid_ptr,alloc(aux_loc),num_aux) + + + ! Need to add eta to the output data + allocate(qeta((num_eqn + num_layers) & + * (num_cells(1) + 2 * num_ghost) & + * (num_cells(2) + 2 * num_ghost))) + do j = 1, num_cells(2) + 2 * num_ghost + do i = 1, num_cells(1) + 2 * num_ghost + + ! Extract depth and momenta + do k=1,num_layers + index = 3 * (k - 1) + h(k) = alloc(iadd(index + 1,i,j)) / rho(k) + hu(k)= alloc(iadd(index + 2,i,j)) / rho(k) + hv(k) = alloc(iadd(index + 3,i,j)) / rho(k) end do - end do - ! Note: We are writing out ghost cell data also - write(out_unit + 1) qeta + ! Calculate surfaces + eta(num_layers) = h(num_layers) & + + alloc(iaddaux(1,i,j)) + do k=num_layers-1,1,-1 + eta(k) = h(k) + eta(k+1) + enddo + + do m=1,num_layers + index = 3*(m - 1) + qeta(iaddqeta(index+1,i,j)) = h(m) + qeta(iaddqeta(index+2,i,j)) = hu(m) + qeta(iaddqeta(index+3,i,j)) = hv(m) + qeta(iaddqeta(3*num_layers+m,i,j)) = eta(m) + end do - deallocate(qeta) + end do + end do + + if (output_format==2) then + ! binary32 (shorten to 4-byte) + allocate(qeta4((num_eqn + num_layers) & + * (num_cells(1) + 2 * num_ghost) & + * (num_cells(2) + 2 * num_ghost))) + qeta4 = real(qeta, kind=4) + write(binary_unit) qeta4 + deallocate(qeta4) + + else if (output_format==3) then + ! binary64 (full 8-byte) + write(binary_unit) qeta + endif + + deallocate(qeta) + + else if (output_format == 4) then ! HDF5 output - case(4) #ifdef HDF5 ! Create data space - handles dimensions of the corresponding ! data set - annoyingling need to stick grid size into other @@ -264,15 +288,20 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) call h5dclose_f(data_set, hdf_error) call h5sclose_f(data_space, hdf_error) #endif - case default - print *, "Unsupported output format", output_format,"." - stop + else + print *, "Unsupported output format", output_format,"." + stop - end select + endif grid_ptr = node(levelptr, grid_ptr) end do end do close(out_unit) + + if (output_format==2 .or. output_format==3) then + close(binary_unit) + endif + #ifdef HDF5 if (output_format == 4) then call h5gclose_f(q_group, hdf_error) @@ -285,13 +314,13 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) if (output_format == 1) then open(unit=out_unit, file=file_name(3), status='unknown', & form='formatted') - else if (output_format == 3) then + else if (output_format==2 .or. output_format==3) then open(unit=out_unit, file=file_name(3), status='unknown', & access='stream') - else if (output_format == 4) then #ifdef HDF5 - ! Create group for aux - call h5gcreate_f(hdf_file, "/aux", aux_group, hdf_error) + else if (output_format == 4) then + ! Create group for aux + call h5gcreate_f(hdf_file, "/aux", aux_group, hdf_error) #endif end if @@ -339,12 +368,19 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) write(out_unit, *) ' ' end do - ! What is case 2? case(2) - stop "Unknown format." - - ! Binary output + ! binary32 + ! Note: We are writing out ghost cell data also + i = (iaddaux(num_aux, num_cells(1) + 2 * num_ghost, & + num_cells(2) + 2 * num_ghost)) + lenaux4 = i - iaddaux(1, 1, 1) + 1 + allocate(aux4(lenaux4)) + aux4 = real(alloc(iaddaux(1, 1, 1):i), kind=4) + write(out_unit) aux4 + deallocate(aux4) + case(3) + ! binary64 ! Note: We are writing out ghost cell data also i = (iaddaux(num_aux, num_cells(1) + 2 * num_ghost, & num_cells(2) + 2 * num_ghost)) @@ -384,6 +420,7 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) end do end do end if + close(out_unit) #ifdef HDF5 if (out_aux) then call h5gclose_f(aux_group, hdf_error) @@ -391,21 +428,32 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) call h5fclose_f(hdf_file, hdf_error) #endif - ! ========================================================================== ! Write fort.t file open(unit=out_unit, file=file_name(2), status='unknown', form='formatted') ! Note: We need to print out num_ghost too in order to strip ghost cells ! from q array when reading in pyclaw.io.binary + + if (output_format == 1) then + file_format = 'ascii' + else if (output_format == 2) then + file_format = 'binary32' + else if (output_format == 3) then + file_format = 'binary64' + else if (output_format == 4) then + file_format = 'hdf' + endif + write(out_unit, t_file_format) time, num_eqn + num_layers, num_grids, & - num_aux, 2, num_ghost + num_aux, 2, num_ghost, file_format close(out_unit) ! ========================================================================== ! Write out timing stats open(unit=out_unit, file=timing_file_name, form='formatted', & - status='old', action='write', position='append') + status='unknown', action='write', position='append') + !status='old', action='write', position='append') timing_line = "(e16.6, ', ', e16.6, ', ', e16.6," do level=1, mxnest @@ -414,15 +462,18 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) end do timing_line = trim(timing_line) // ")" - if (time == t0) then + if (abs(time - t0) < 1d-15) then t_CPU_overall = 0.d0 timeTick_overall = 0.d0 - else + else call cpu_time(t_CPU_overall) + ! if this is a restart, need to adjust add in time from previous run: + t_CPU_overall = t_CPU_overall + timeTickCPU + call system_clock(tick_clock_finish,tick_clock_rate) timeTick_int = timeTick + tick_clock_finish - tick_clock_start timeTick_overall = real(timeTick_int, kind=8)/real(clock_rate,kind=8) - endif + endif write(out_unit, timing_line) time, timeTick_overall, t_CPU_overall, & (real(tvoll(i), kind=8) / real(clock_rate, kind=8), & @@ -434,11 +485,11 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) ! Print output info if (display_landfall_time) then ! Convert time to days relative to landfall - console_format = "('AMRCLAW: Frame ',i4,' output files done at " // & + console_format = "('GEOCLAW: Frame ',i4,' output files done at " // & "time t = ', f5.2,/)" - print console_format, frame, time / (3.3d3 * 24.d0) + print console_format, frame, time / (24.d0 * 60d0**2) else - console_format = "('AMRCLAW: Frame ',i4,' output files done at " // & + console_format = "('GEOCLAW: Frame ',i4,' output files done at " // & "time t = ', d13.6,/)" print console_format, frame, time end if diff --git a/src/2d/shallow/stepgrid.f b/src/2d/shallow/stepgrid.f index ed4c30af2..990319871 100644 --- a/src/2d/shallow/stepgrid.f +++ b/src/2d/shallow/stepgrid.f @@ -23,12 +23,14 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, c c c This version of stepgrid, stepgrid_geo.f allows output on -c fixed grids specified in setfixedgrids.data +c fgout grids specified in fgout_grids.data c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: use geoclaw_module use amr_module - use fixedgrids_module + use fgout_module, only: FGOUT_num_grids, FGOUT_fgrids, + & FGOUT_tcfmax, fgout_interp, fgout_grid, + & FGOUT_ttol implicit double precision (a-h,o-z) external rpn2,rpt2 @@ -44,9 +46,13 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, logical :: debug = .false. logical :: dump = .false. logical, intent (in) :: actualstep + type(fgout_grid), pointer :: fgout + logical, allocatable :: fgout_interp_needed(:) + real(kind=8) :: fgout_tnext c - tcfmax = -rinfinity + FGOUT_tcfmax = -rinfinity level = node(nestlevel,mptr) + if (dump) then write(outunit,*)" at start of stepgrid: dumping grid ",mptr @@ -105,104 +111,61 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, C mwork1 = mwork - mused !# remaining space (passed to step2) c - -c::::::::::::::::::::::::Fixed Grid Output::::::::::::::::::::::::::::::::: tc0=time !# start of computational step tcf=tc0+dt !# end of computational step + +c Check if fgout interpolation needed before and after step: + + allocate(fgout_interp_needed(FGOUT_num_grids)) !$OMP CRITICAL (FixedGrids) -c # see if any f-grids should be written out - do ng=1,num_fixed_grids - if (tc0 > fgrids(ng)%start_time .and. - & fgrids(ng)%last_output_index < fgrids(ng)%num_output) then -c # fgrid ng may need to be written out -c # find the first output number that has not been written out and -c # find the first output number on a fixed grid that is >= tc0 -c # which will not be written out - if (fgrids(ng)%dt > 0.d0) then - ioutfgend= 1+max(0,nint((tc0 - fgrids(ng)%start_time) - & / fgrids(ng)%dt)) - else - ioutfgend=1 - endif - ioutfgend = min(ioutfgend,fgrids(ng)%num_output) - ioutfgstart = fgrids(ng)%last_output_index + 1 -c # write-out fgrid times that are less than tc0, and have not been written yet -c # these should be the most accurate values at any given point in the fgrid -c # since tc0> output time - do ioutfg=ioutfgstart,ioutfgend - toutfg=fgrids(ng)%start_time+(ioutfg-1)*fgrids(ng)%dt - if (toutfg < tc0) then -c # write out the solution for fixed grid ng -c # test if arrival times should be output - ioutflag = fgrids(ng)%output_arrival_times* - & (fgrids(ng)%num_output- - & fgrids(ng)%last_output_index) - call fgrid_out(ng,fgrids(ng),toutfg,ioutfg,ioutflag) - - fgrids(ng)%last_output_time = toutfg - fgrids(ng)%last_output_index = - & fgrids(ng)%last_output_index + 1 - endif - enddo - - endif + + do ng=1,FGOUT_num_grids + fgout => FGOUT_fgrids(ng) + if (fgout%next_output_index > fgout%num_output) then + fgout_interp_needed(ng) = .false. + else + fgout_tnext = fgout%output_times(fgout%next_output_index) + fgout_interp_needed(ng) = + & ((fgout%x_low < xlowmbc + mx * dx) .and. + & (fgout%x_hi > xlowmbc) .and. + & (fgout%y_low < ylowmbc + my * dy) .and. + & (fgout%y_hi > ylowmbc) .and. + & (fgout_tnext >= tc0 - FGOUT_ttol) .and. + & (fgout_tnext <= tcf + FGOUT_ttol)) + endif +c write(6,*) '+++ level, before- tc0,tcf,needed: ', +c & level,tc0,tcf,fgout_interp_needed(ng) +c write(6,*) '+++ next index: ',fgout%next_output_index +c write(6,*) '+++ fgout_tnext, tcf + FGOUT_ttol: ', +c & fgout_tnext, tcf + FGOUT_ttol enddo !$OMP END CRITICAL (FixedGrids) + + +c::::::::::::::::::::::::fgout Output::::::::::::::::::::::::::::::::: +c This has been moved to tick.f, after advancing all patches on +c finest level. No need to check on each patch separately. c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: call b4step2(mbc,mx,my,nvar,q, & xlowmbc,ylowmbc,dx,dy,time,dt,maux,aux,actualstep) -c::::::::::::::::::::::::FIXED GRID DATA before step::::::::::::::::::::::: -c # fill in values at fixed grid points effected at time tc0 -!$OMP CRITICAL (FixedGrids) - do ng=1,num_fixed_grids - - if ( (fgrids(ng)%x_low < xlowmbc + mx*dx) .and. - & (fgrids(ng)%x_hi > xlowmbc) .and. - & (fgrids(ng)%y_low < ylowmbc + my*dy) .and. - & (fgrids(ng)%y_hi > ylowmbc) .and. - & (fgrids(ng)%last_output_index < fgrids(ng)%num_output) .and. - & (tcf >= fgrids(ng)%start_time) ) then - - if (fgrids(ng)%last_output_time + fgrids(ng)%dt >= tc0 .and. - & fgrids(ng)%last_output_time + fgrids(ng)%dt <= tcf) then - -c # fixedgrid ng has an output time within [tc0,tcf] interval -c # and it overlaps this computational grid spatially - call fgrid_interp(1,fgrids(ng),tc0,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,0) - -c # routine to spatially interpolate computational solution -c # at tc0 to the fixed grid spatial points, -c #saving solution, variables and tc0 at every grid point - endif - -c # set maxima or minima if this is a new coarse step -c if (tc0.ge.tcfmax) then - -c # RJL: rewrote to set min/max every time a grid at level 1 -c # is about to be taken. The previous code failed if there was more than one grid -c # at level 1. Note that all grids are up to date at start of step on level 1. -c # New feature added at end of this routine to check more frequently if -c # levelcheck > 0. - if (level .eq. 1) then - if (fgrids(ng)%output_surface_max - & + fgrids(ng)%output_arrival_times > 0) then - - call fgrid_interp(3,fgrids(ng),tc0,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,2) - - endif +c::::::::::::::::::::::::FGOUT DATA before step::::::::::::::::::::::: +c # fill in values at fgout points affected at time tc0 + do ng=1,FGOUT_num_grids + if (fgout_interp_needed(ng)) then +c # fgout grid ng has an output time within [tc0,tcf] interval +c # and it overlaps this computational grid spatially +!$OMP CRITICAL (FixedGrids) + fgout => FGOUT_fgrids(ng) +c write(6,*) '+++ fout_interp(1), tc0, level: ',tc0,level + call fgout_interp(1,fgout,tc0,q,nvar,mx,my,mbc, + & dx,dy,xlowmbc,ylowmbc,maux,aux) + +!$OMP END CRITICAL (FixedGrids) endif - - endif enddo - tcfmax=max(tcfmax,tcf) - -!$OMP END CRITICAL (FixedGrids) - c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c New fixed grid stuff: Update fixed grid info from this patch... @@ -277,66 +240,25 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, & q,maux,aux,time,dt) endif -!$OMP CRITICAL (FixedGrids) -c ::::::::::::::::::::::::Fixed Grid data afterstep::::::::::::::::::::::: -c # fill in values at fixed grid points effected at time tcf - do ng=1,num_fixed_grids - if ((fgrids(ng)%x_low < xlowmbc + mx * dx) .and. - & (fgrids(ng)%x_hi > xlowmbc) .and. - & (fgrids(ng)%y_low < ylowmbc + my * dy) .and. - & (fgrids(ng)%y_hi > ylowmbc) .and. - & (fgrids(ng)%last_output_index < fgrids(ng)%num_output) .and. - & (tcf >= fgrids(ng)%start_time)) then - - if (fgrids(ng)%last_output_time + fgrids(ng)%dt >= tc0 .and. - & fgrids(ng)%last_output_time + fgrids(ng)%dt <= tcf) then - -c # fixedgrid ng has an output time within [tc0,tcf] interval -c # and it overlaps this computational grid spatially -C i0=i0fg(ng) !# index into the ng grid in the work array - - call fgrid_interp(2,fgrids(ng),tcf,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,0) - -c # routine to interpolate solution -c # at tcf to the fixed grid storage array, -c #saving solution and tcf at every grid point - - endif - -c # fill in values for eta if they need to be saved for later checking max/mins -c # check for arrival times - if (fgrids(ng)%output_surface_max - & + fgrids(ng)%output_arrival_times > 0) then - - call fgrid_interp(3,fgrids(ng),tc0,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,1) - - endif - -c # RJL: Modified 8/20/11 -c # If levelcheck > 0 then update max/mins at end of step on this grid. -c # Note that if there are finer grids then fgridoften will not have been updated -c # properly yet by those grids. This modification allows checking max/min more -c # frequently than the original code (equivalent to levelcheck==0) when you know -c # what level is most relevant for this fixed grid. Note also that if there are no -c # grids at levelcheck overlapping a portion of the fixed grid then the max/min values -c # will be updated only at start of next level 1 step. - - levelcheck = 0 - if (level == levelcheck) then - if (fgrids(ng)%output_arrival_times - & + fgrids(ng)%output_surface_max > 0) then - - call fgrid_interp(3,fgrids(ng),tc0,q,nvar,mx,my,mbc,dx,dy, - & xlowmbc,ylowmbc,maux,aux,2) - endif +c ::::::::::::::::::::::::fgout data afterstep::::::::::::::::::::::: +c # fill in values at fgout points affected at time tcf + do ng=1,FGOUT_num_grids + if (fgout_interp_needed(ng)) then +c # fgout grid ng has an output time within [tc0,tcf] interval +c # and it overlaps this computational grid spatially +!$OMP CRITICAL (FixedGrids) + fgout => FGOUT_fgrids(ng) +c write(6,*) '+++ fout_interp(2), tcf, level: ',tcf,level + call fgout_interp(2,fgout,tcf,q,nvar,mx,my,mbc, + & dx,dy,xlowmbc,ylowmbc,maux,aux) + +!$OMP END CRITICAL (FixedGrids) endif - - - endif +c write(6,*) '+++ level,after- tc0,tcf,needed: ', +c & level,tc0,tcf,fgout_interp_needed(ng) +c write(6,*) '+++ next index: ',fgout%next_output_index +c write(6,*) '+++ fgout_tnext: ',fgout_tnext enddo -!$OMP END CRITICAL (FixedGrids) c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: diff --git a/src/2d/shallow/tick.f b/src/2d/shallow/tick.f index 9a176d440..fc1e3912e 100644 --- a/src/2d/shallow/tick.f +++ b/src/2d/shallow/tick.f @@ -14,8 +14,8 @@ subroutine tick(nvar,cut,nstart,vtime,time,naux,start_time, use storm_module, only: landfall, display_landfall_time use fgmax_module, only: FG_num_fgrids, FG_fgrids, fgrid - - + use fgout_module, only: FGOUT_num_grids, FGOUT_fgrids, + & fgout_write,fgout_grid, FGOUT_ttol implicit double precision (a-h,o-z) @@ -27,6 +27,7 @@ subroutine tick(nvar,cut,nstart,vtime,time,naux,start_time, character(len=128) :: time_format real(kind=8) cpu_start,cpu_finish type(fgrid), pointer :: fg + type(fgout_grid), pointer :: fgout c c :::::::::::::::::::::::::::: TICK ::::::::::::::::::::::::::::: @@ -343,6 +344,39 @@ subroutine tick(nvar,cut,nstart,vtime,time,naux,start_time, possk(level) = possk(level-1)/ntogo(level) go to 60 endif + +c write(6,*) '+++ in tick, done with level',level, +c & ' tlevel = ',tlevel(1:level) + +c When we reach here, we are done with grid patches at +c the finest level with grids present at this time + +c ---- fgout output ---- +c This used to be done in stepgrid.f, but only needs to be done +c after all patches at finest level have been advanced. + + tc0 = tlevel(level) ! current time on finest level present + !write(6,*) '+++ tick: tc0 = ',tc0 + + do ng=1,FGOUT_num_grids + fgout => FGOUT_fgrids(ng) + + do ioutfg=1,fgout%num_output + if (fgout%output_frames(ioutfg) == -1) then + ! this time not yet written out + if (fgout%output_times(ioutfg) < + & tc0+FGOUT_ttol) then + toutfg = fgout%output_times(ioutfg) +c write(6,*) '+++ tick call fgrid_out, frame, t: ', +c & ioutfg,toutfg + call fgout_write(fgout,toutfg,ioutfg) + fgout%output_frames(ioutfg) = ioutfg + fgout%next_output_index = ioutfg+1 + endif + endif + enddo + enddo +c ---- end fgout output ---- c 105 if (level .eq. 1) go to 110 if (ntogo(level) .gt. 0) then diff --git a/src/2d/shallow/valout.f90 b/src/2d/shallow/valout.f90 index f778cee3e..c7b58b463 100644 --- a/src/2d/shallow/valout.f90 +++ b/src/2d/shallow/valout.f90 @@ -40,6 +40,8 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) real(kind=8) :: h, hu, hv, eta real(kind=8), allocatable :: qeta(:) + real(kind=4), allocatable :: qeta4(:), aux4(:) + integer :: lenaux4 #ifdef HDF5 ! HDF writing @@ -55,6 +57,7 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) character(len=128) :: console_format character(len=512) :: timing_line, timing_substr character(len=*), parameter :: timing_file_name = "timing.csv" + character(len=8) :: file_format character(len=*), parameter :: header_format = & "(i6,' grid_number',/," // & @@ -70,7 +73,8 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) "i6,' ngrids'/," // & "i6,' naux'/," // & "i6,' ndim'/," // & - "i6,' nghost'/,/)" + "i6,' nghost'/," // & + "a10,' format'/,/)" ! Output timing @@ -118,7 +122,7 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) ! Write out fort.q file (and fort.bXXXX and clawxxxx.h5 files if necessary) ! Here we let fort.q be out_unit and the the other two be out_unit + 1 open(unit=out_unit, file=file_name(1), status='unknown', form='formatted') - if (output_format == 3) then + if (output_format == 2 .or. output_format == 3) then open(unit=binary_unit, file=file_name(4), status="unknown", & access='stream') #ifdef HDF5 @@ -157,18 +161,19 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) delta(1), delta(2) ! Output grids - select case(output_format) - ! ASCII output - case(1) - ! Round off if nearly zero - forall (m = 1:num_eqn, & - i=num_ghost + 1:num_cells(1) + num_ghost, & - j=num_ghost + 1:num_cells(2) + num_ghost, & - abs(alloc(iadd(m, i, j))) < 1d-90) - - alloc(iadd(m, i, j)) = 0.d0 - end forall - + + ! Round off if nearly zero + ! (Do this for all output_format's) + forall (m = 1:num_eqn, & + i=num_ghost + 1:num_cells(1) + num_ghost, & + j=num_ghost + 1:num_cells(2) + num_ghost, & + abs(alloc(iadd(m, i, j))) < 1d-90) + + alloc(iadd(m, i, j)) = 0.d0 + end forall + + if (output_format == 1) then + ! ascii output do j = num_ghost + 1, num_cells(2) + num_ghost do i = num_ghost + 1, num_cells(1) + num_ghost @@ -188,40 +193,49 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) write(out_unit, *) ' ' end do - ! What is case 2? - case(2) - stop "Unknown format." - - ! Binary output - case(3) - - ! Updating ghost cell data - call bound(time,num_eqn,num_ghost,alloc(q_loc), & - num_cells(1) + 2*num_ghost, & - num_cells(2) + 2*num_ghost, & - grid_ptr,alloc(aux_loc),num_aux) - - ! Need to add eta to the output data - allocate(qeta((num_eqn + 1) & - * (num_cells(1) + 2 * num_ghost) & - * (num_cells(2) + 2 * num_ghost))) - do j = 1, num_cells(2) + 2 * num_ghost - do i = 1, num_cells(1) + 2 * num_ghost - do m = 1, num_eqn - qeta(iaddqeta(m, i, j)) = alloc(iadd(m, i, j)) - end do - eta = alloc(iadd(1, i, j)) + alloc(iaddaux(1, i ,j)) - qeta(iaddqeta(num_eqn + 1, i, j)) = eta + else if (output_format==2 .or. output_format==3) then + ! binary32 or binary64 + + ! Note: We are writing out ghost cell data also, + ! so need to update this + call bound(time,num_eqn,num_ghost,alloc(q_loc), & + num_cells(1) + 2*num_ghost, & + num_cells(2) + 2*num_ghost, & + grid_ptr,alloc(aux_loc),num_aux) + + + ! Need to add eta to the output data + allocate(qeta((num_eqn + 1) & + * (num_cells(1) + 2 * num_ghost) & + * (num_cells(2) + 2 * num_ghost))) + do j = 1, num_cells(2) + 2 * num_ghost + do i = 1, num_cells(1) + 2 * num_ghost + do m = 1, num_eqn + qeta(iaddqeta(m, i, j)) = alloc(iadd(m, i, j)) end do + eta = alloc(iadd(1, i, j)) + alloc(iaddaux(1, i ,j)) + qeta(iaddqeta(num_eqn + 1, i, j)) = eta end do + end do + + if (output_format==2) then + ! binary32 (shorten to 4-byte) + allocate(qeta4((num_eqn + 1) & + * (num_cells(1) + 2 * num_ghost) & + * (num_cells(2) + 2 * num_ghost))) + qeta4 = real(qeta, kind=4) + write(binary_unit) qeta4 + deallocate(qeta4) - ! Note: We are writing out ghost cell data also + else if (output_format==3) then + ! binary64 (full 8-byte) write(binary_unit) qeta + endif + + deallocate(qeta) - deallocate(qeta) - + else if (output_format == 4) then ! HDF5 output - case(4) #ifdef HDF5 ! Create data space - handles dimensions of the corresponding ! data set - annoyingling need to stick grid size into other @@ -242,17 +256,17 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) call h5dclose_f(data_set, hdf_error) call h5sclose_f(data_space, hdf_error) #endif - case default - print *, "Unsupported output format", output_format,"." - stop + else + print *, "Unsupported output format", output_format,"." + stop - end select + endif grid_ptr = node(levelptr, grid_ptr) end do end do close(out_unit) - if (output_format == 3) then + if (output_format==2 .or. output_format==3) then close(binary_unit) #ifdef HDF5 else if (output_format == 4) then @@ -266,7 +280,7 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) if (output_format == 1) then open(unit=out_unit, file=file_name(3), status='unknown', & form='formatted') - else if (output_format == 3) then + else if (output_format==2 .or. output_format==3) then open(unit=out_unit, file=file_name(3), status='unknown', & access='stream') #ifdef HDF5 @@ -320,12 +334,19 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) write(out_unit, *) ' ' end do - ! What is case 2? case(2) - stop "Unknown format." - - ! Binary output + ! binary32 + ! Note: We are writing out ghost cell data also + i = (iaddaux(num_aux, num_cells(1) + 2 * num_ghost, & + num_cells(2) + 2 * num_ghost)) + lenaux4 = i - iaddaux(1, 1, 1) + 1 + allocate(aux4(lenaux4)) + aux4 = real(alloc(iaddaux(1, 1, 1):i), kind=4) + write(out_unit) aux4 + deallocate(aux4) + case(3) + ! binary64 ! Note: We are writing out ghost cell data also i = (iaddaux(num_aux, num_cells(1) + 2 * num_ghost, & num_cells(2) + 2 * num_ghost)) @@ -380,8 +401,19 @@ subroutine valout(level_begin, level_end, time, num_eqn, num_aux) ! Note: We need to print out num_ghost too in order to strip ghost cells ! from q array when reading in pyclaw.io.binary + + if (output_format == 1) then + file_format = 'ascii' + else if (output_format == 2) then + file_format = 'binary32' + else if (output_format == 3) then + file_format = 'binary64' + else if (output_format == 4) then + file_format = 'hdf' + endif + write(out_unit, t_file_format) time, num_eqn + 1, num_grids, num_aux, & - 2, num_ghost + 2, num_ghost, file_format close(out_unit) ! ========================================================================== diff --git a/src/python/geoclaw/data.py b/src/python/geoclaw/data.py index fdb4b59c9..12739f31a 100755 --- a/src/python/geoclaw/data.py +++ b/src/python/geoclaw/data.py @@ -212,9 +212,12 @@ def write(self,data_source='setrun.py', out_file='topo.data'): self.close_data_file() - class FixedGridData(clawpack.clawutil.data.ClawData): + """ + Deprecated, starting in 5.9.0 use FGoutData instead. + """ + def __init__(self): super(FixedGridData,self).__init__() @@ -225,12 +228,49 @@ def __init__(self): def write(self,data_source='setrun.py', out_file='fixed_grids.data'): # Fixed grid settings + msg = 'rundata.fixed_grid_data is deprecated starting in v5.9.0,' \ + + ' use rundata.fgout_data instead' + #warnings.warn(msg) + if len(self.fixedgrids) > 0: + raise AttributeError(msg) + + +class FGoutData(clawpack.clawutil.data.ClawData): + + def __init__(self): + + super(FGoutData,self).__init__() + + # File name for fgout points and parameters: + self.add_attribute('fgout_grids',[]) + + + def write(self,data_source='setrun.py', out_file='fgout_grids.data'): self.open_data_file(out_file, data_source) - nfixedgrids = len(self.fixedgrids) - self.data_write(value=nfixedgrids,alt_name='nfixedgrids') + num_fgout_grids = len(self.fgout_grids) + self.data_write(value=num_fgout_grids,alt_name='num_fgout_grids') self.data_write() - for fixedgrid in self.fixedgrids: - self._out_file.write(11*"%g " % tuple(fixedgrid) +"\n") + + fgno_unset = 0 # to use if fg.fgno not set by user + fgno_list = [] # to check for uniqueness of fgno's + + for fg in self.fgout_grids: + # if path is relative in setrun, assume it's relative to the + # same directory that out_file comes from + + if fg.fgno is None: + # not set by user in setrun + fgno_unset += 1 + fg.fgno = fgno_unset + + if fg.fgno in fgno_list: + msg = 'Trying to set fgout grid number to fgno = %i' % fg.fgno \ + + '\n but this fgno was already used' \ + + '\n Set unique fgno for each fgout grid' + raise ValueError(msg) + + fgno_list.append(fg.fgno) + fg.write_to_fgout_data(self._out_file) self.close_data_file() class FGmaxData(clawpack.clawutil.data.ClawData): diff --git a/src/python/geoclaw/fgmax_tools.py b/src/python/geoclaw/fgmax_tools.py index 77c264835..c98b55230 100644 --- a/src/python/geoclaw/fgmax_tools.py +++ b/src/python/geoclaw/fgmax_tools.py @@ -179,6 +179,10 @@ def read_fgmax_grids_data(self, fgno, data_file='fgmax_grids.data'): def write_to_fgmax_data(self, fid): + """ + Write the fgmax grid data to the file specified by `fid`, normally + the `fgmax_grids.data` file that is read in by the GeoClaw Fortran code. + """ print("\n---------------------------------------------- ") point_style = self.point_style @@ -507,6 +511,9 @@ def set_q_time(ind_q, ind_q_time): def bounding_box(self): + """ + Return the bounding box of the grid as a list [x1,x2,y1,y2] + """ x1 = self.X.min() x2 = self.X.max() y1 = self.Y.min() @@ -598,6 +605,11 @@ def interp_dz(self,dtopo_path,dtopo_type): def adjust_fgmax_1d(x1_desired, x2_desired, x1_domain, dx): """ + Adjust the upper and lower limits of a grid so that equally spaced + grid points with spacing `dx` lie exactly at cell centers, so that + no interpolation is needed for fgmax values. Note that parameter + names refer to `x` limits, but works equally well for `y` values. + :Input: - x1_desired, x2_desired: approximate desired limits of fgmax grid - x1_domain: lower edge of computational domain @@ -605,7 +617,8 @@ def adjust_fgmax_1d(x1_desired, x2_desired, x1_domain, dx): :Output: - x1_new, x2_new: limits to set so (x2-x1) is integer multiple of dx and points are at cell centers of computational grid - - npoints: number of points + - npoints: number of points to specify, so that + `linspace(x1_new, x2_new, npoints) gives points with spacing `dx`. """ i1 = numpy.floor((x1_desired-x1_domain - 0.5*dx)/dx) diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py new file mode 100644 index 000000000..f5df57c90 --- /dev/null +++ b/src/python/geoclaw/fgout_tools.py @@ -0,0 +1,869 @@ +r""" +fgout_tools module: $CLAW/geoclaw/src/python/geoclaw/fgout_tools.py + +Tools to specify and work with fgout grids, used to output GeoClaw solutions +on a fixed grid at a sequence of times, regardless of the AMR structure. + +Includes: + +- class FGoutFrame: used to hold a single frame of fgout output data +- class FGoutGrid: used to specify and store info about an fgout grid, with + methods to read and write info to fgout_grids.data +- function make_fgout_fcn_xy: Takes an FGoutFrame object and produces an + interpolating function that can be evaluated for any (x,y). +- function make_fgout_fcn_xyt: Takes 2 FGoutFrame objects and produces an + interpolating function that can be evaluated for any (x,y,t) + at intermediate times. +- function write_netcdf: Write a specified set of qoi's from a list of + fgout frames, as a single netCDF file +- function read_netcdf: Read a netCDF file and return a list of fgout frames, + assuming the file contains all the qoi's needed to reconstruct q. +- function read_netcdf_arrays: Read a netCDF file and extract the + requested quantities of interest as numpy arrays. +- print_netcdf_info: Print info about the contents of a netCDF file containing + some fgout frames. +""" + +import os +from numpy import sqrt, ma, mod +import numpy +from six.moves import range + + +class FGoutFrame(object): + + """ + Class to hold a single frame of fgout data at one output time. + Several attributes are defined as properties that can be evaluated + and stored only when needed by the user. + """ + + def __init__(self, fgout_grid, frameno=None): + self.fgout_grid = fgout_grid + self.frameno = frameno + self.t = None + + # private attributes for those that are only created if + # needed by the user: + self._h = None + self._hu = None + self._hv = None + self._eta = None + self._B = None + self._u = None + self._v = None + self._s = None + self._hss = None + + # Define shortcuts to attributes of self.fgout_grid that are the same + # for all frames (e.g. X,Y) to avoid storing grid for every frame. + + @property + def x(self): + return self.fgout_grid.x + + @property + def y(self): + return self.fgout_grid.y + + @property + def X(self): + return self.fgout_grid.X + + @property + def Y(self): + return self.fgout_grid.Y + + @property + def delta(self): + return self.fgout_grid.delta + + @property + def extent_centers(self): + return self.fgout_grid.extent_centers + + @property + def extent_edges(self): + return self.fgout_grid.extent_edges + + @property + def drytol(self): + return self.fgout_grid.drytol + + # Define attributes such as h as @properties with lazy evaluation: + # the corresponding array is created and stored only when first + # accessed by the user. Those not needed are not created. + + @property + def h(self): + """depth""" + if self._h is None: + #print('+++ setting _h...') + self._h = self.q[0,:,:] + #print('+++ getting _h...') + return self._h + + @property + def hu(self): + """momentum h*u""" + if self._hu is None: + self._hu = self.q[1,:,:] + return self._hu + + @property + def u(self): + """speed u, computed as hu/h or set to 0 if hself.drytol)) + return self._u + + @property + def hv(self): + """momentum h*v""" + if self._hv is None: + self._hv = self.q[2,:,:] + return self._hv + + @property + def v(self): + """speed v, computed as hv/h or set to 0 if hself.drytol)) + return self._v + + @property + def eta(self): + """surface eta = h+B""" + if self._eta is None: + self._eta = self.q[-1,:,:] + return self._eta + + @property + def B(self): + """topography""" + if self._B is None: + self._B = self.q[-1,:,:] - self.q[0,:,:] + return self._B + + @property + def s(self): + """speed s = sqrt(u**2 + v**2)""" + if self._s is None: + self._s = numpy.sqrt(self.u**2 + self.v**2) + return self._s + + @property + def hss(self): + """momentum flux h*s**2""" + if self._hss is None: + self._hss = self.h * self.s**2 + return self._hss + + +class FGoutGrid(object): + + """ + New class introduced in 5.9.0 to keep store information both about the + fgout input data and the output generated by a GeoClaw run. + """ + + def __init__(self,fgno=None,outdir=None,output_format=None): + + + # GeoClaw input values: + self.id = '' # identifier, optional + self.point_style = 2 # only option currently supported + self.npts = None + self.nx = None + self.ny = None + self.tstart = None + self.tend = None + self.nout = None + self.fgno = fgno + self.outdir = outdir + self.output_format = output_format + + self.drytol = 1e-3 # used for computing u,v from hu,hv + + # private attributes for those that are only created if + # needed by the user: + + self._X = None + self._Y = None + self._x = None + self._y = None + self._delta = None # (dx,dy) + self._extent_centers = None # defined by fgout points + self._extent_edges = None # extended so points are cell centers + + self._plotdata = None + + @property + def x(self): + """1D array x of longitudes (cell centers)""" + if self._x is None: + dx = (self.x2 - self.x1)/self.nx + self._x = numpy.linspace(self.x1+dx/2, self.x2-dx/2, self.nx) + return self._x + + @property + def y(self): + """1D array y of latitudes (cell centers)""" + if self._y is None: + dy = (self.y2 - self.y1)/self.ny + self._y = numpy.linspace(self.y1+dy/2, self.y2-dy/2, self.ny) + return self._y + + @property + def X(self): + """2D array X of longitudes (cell centers)""" + if self._X is None: + self._X = numpy.meshgrid(self.x, self.y)[0].T + return self._X + + @property + def Y(self): + """2D array Y of latitudes (cell centers)""" + if self._Y is None: + self._Y = numpy.meshgrid(self.x, self.y)[1].T + return self._Y + + @property + def delta(self): + if self._delta is None: + dx = (self.x2 - self.x1)/self.nx + dy = (self.y2 - self.y1)/self.ny + self._delta = (dx,dy) + return self._delta + + @property + def extent_centers(self): + """Extent of cell centers [xmin,xmax,ymin,ymax]""" + if self._extent_centers is None: + self._extent_centers = [self.x.min(), self.x.max(),\ + self.y.min(), self.y.max()] + return self._extent_centers + + @property + def extent_edges(self): + """Extent of cell edges [xmin,xmax,ymin,ymax]""" + if self._extent_edges is None: + dx,dy = self.delta + self._extent_edges = [self.x.min()-dx/2, self.x.max()+dx/2,\ + self.y.min()-dy/2, self.y.max()+dy/2] + return self._extent_edges + + + # Create plotdata of class clawpack.visclaw.ClawPlotData + # only when needed for reading GeoClaw output, + # since this must be done after fgno, outdir, output_format + # have been specified: + + @property + def plotdata(self): + if self._plotdata is None: + self._plotdata = self.set_plotdata() + return self._plotdata + + def set_plotdata(self): + """ + Create a plotdata, assuming attributes fgno, outdir, output_format + have all been set. + """ + from clawpack.visclaw.data import ClawPlotData + + assert self.fgno is not None, '*** fgno must be set' + assert self.outdir is not None, '*** outdir must be set' + assert self.output_format is not None, '*** output_format must be set' + if 0: + print('+++ creating plotdata for fgno=%i, outdir=%s, format=%s' \ + % (self.fgno,self.outdir,self.output_format)) + + plotdata = ClawPlotData() + plotdata.outdir = self.outdir + plotdata.format = self.output_format + plotdata.file_prefix = 'fgout%s' % str(self.fgno).zfill(4) + if self.output_format[:6]=='binary': + # could be 'binary', 'binary32' or 'binary64' + file_prefix_str = plotdata.file_prefix + '.b' + else: + file_prefix_str = plotdata.file_prefix + '.q' + + # Contrary to usual default, set save_frames to False so that all + # the fgout frames read in are not saved in plotdata.framesoln_dict. + # Otherwise this might be a memory hog when making an animation with + # many frames on a fine fgout grid. + plotdata.save_frames = False + + return plotdata + + + def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'): + """ + Read input info for fgout grid number fgno from the data file + fgout_grids.data, which should have been created by setrun.py. + This file now contains info about all fgout grids. + """ + if fgno is not None: + self.fgno = fgno + assert self.fgno is not None, '*** fgno must be set' + + with open(data_file) as filep: + lines = filep.readlines() + fgout_input = None + for lineno,line in enumerate(lines): + if 'fgno' in line: + if int(line.split()[0]) == self.fgno: + fgout_input = lines[lineno+1:] + #print('Found line %i: %s' % (lineno,line)) + break + + if fgout_input is None: + raise ValueError('fgout grid fgno = %i not found in %s' \ + % (fgno, data_file)) + + self.tstart = float(fgout_input[0].split()[0]) + self.tend = float(fgout_input[1].split()[0]) + self.nout = int(fgout_input[2].split()[0]) + self.point_style = point_style = int(fgout_input[3].split()[0]) + output_format = int(fgout_input[4].split()[0]) + if output_format == 1: + self.output_format = 'ascii' + elif output_format == 3: + self.output_format = 'binary' + print('Reading input for fgno=%i, point_style = %i ' \ + % (self.fgno, self.point_style)) + if point_style == 2: + self.nx = nx = int(fgout_input[5].split()[0]) + self.ny = ny = int(fgout_input[5].split()[1]) + self.x1 = float(fgout_input[6].split()[0]) + self.y1 = float(fgout_input[6].split()[1]) + self.x2 = float(fgout_input[7].split()[0]) + self.y2 = float(fgout_input[7].split()[1]) + else: + raise NotImplementedError("fgout not implemented for point_style %i" \ + % point_style) + + + def write_to_fgout_data(self, fid): + """ + Convert fgout data specified in setrun.py to file `fgout_grids.data` + read in by GeoClaw fortran code. + """ + + print("\n---------------------------------------------- ") + assert self.point_style is not None, 'Need to set point_style' + + point_style = self.point_style + if point_style not in [2]: + errmsg = 'point_style %s is not implement, only point_style==2' \ + % point_style + raise NotImplementedError(errmsg) + + if self.output_format == 'ascii': + output_format = 1 + elif self.output_format == 'binary32': + output_format = 2 + elif self.output_format in ['binary','binary64']: + output_format = 3 + else: + errmsg = "fgout output_format must be ascii, binary32, or binary63" + raise NotImplementedError(errmsg) + + assert self.tstart is not None, 'Need to set tstart' + assert self.tend is not None, 'Need to set tend' + assert self.nout is not None, 'Need to set nout' + assert self.point_style is not None, 'Need to set point_style' + + # write header, independent of point_style: + #fid = open(self.input_file_name,'w') + fid.write("\n") + fid.write("%i # fgno\n" % self.fgno) + fid.write("%16.10e # tstart\n" % self.tstart) + fid.write("%16.10e # tend\n" % self.tend) + fid.write("%i %s # nout\n" % (self.nout, 11*" ")) + fid.write("%i %s # point_style\n" \ + % (self.point_style,12*" ")) + fid.write("%i %s # output_format\n" \ + % (output_format,12*" ")) + + print('fgout grid %i has point_style = %i' % (self.fgno, point_style)) + + + if point_style == 2: + # 2d grid of points + x1,x2 = self.x1, self.x2 + y1,y2 = self.y1, self.y2 + nx,ny = self.nx, self.ny + + dx = (x2-x1)/nx # x1,x2 are cell edges + dy = (y2-y1)/ny # y1,y2 are cell edges + + npts = nx*ny + + fid.write("%i %i %s # nx,ny\n" \ + % (nx,ny,10*" ")) + fid.write("%16.10e %20.10e # x1, y1\n" % (x1,y1)) + fid.write("%16.10e %20.10e # x2, y2\n" % (x2,y2)) + + print(" specifying fgout grid with shape %i by %i, with %i points" \ + % (nx,ny,npts)) + print(" lower left = (%15.10f,%15.10f)" % (x1,y1)) + print(" upper right = (%15.10f,%15.10f)" % (x2,y2)) + print(" dx = %15.10e, dy = %15.10e" % (dx,dy)) + + + + def read_frame(self, frameno): + """ + Read a single frame of fgout data. + """ + + from datetime import timedelta + + try: + fr = self.plotdata.getframe(frameno) + except: + print('*** Could not read fgout grid %i frame %i from %s' \ + % (frameno,self.fgno,self.plotdata.outdir)) + raise + state = fr.states[0] # only 1 AMR grid + patch = state.patch + + fgout_frame = FGoutFrame(self, frameno) + fgout_frame.fgout_grid = self + + fgout_frame.q = state.q + + fgout_frame.t = state.t + + fgout_frame.frameno = frameno + + X,Y = patch.grid.p_centers[:2] + + try: + fgoutX = self.X + fgoutY = self.Y + except: + # presumably x,y, etc. were not set for this fgout_grid + # (e.g. by read_fgout_grids_data) so infer from this frame and set + + # reset most attributes including private _x etc to None: + self.__init__(fgno=self.fgno,outdir=self.outdir, + output_format=self.output_format) + + print(' Setting grid attributes based on Frame %i:' % frameno) + x = X[:,0] + y = Y[0,:] + dx = x[1] - x[0] + dy = y[1] - y[0] + self.x1 = x[0] - dx/2 + self.x2 = x[-1] + dx/2 + self.y1 = y[0] - dy/2 + self.y2 = y[-1] + dy/2 + self.nx = len(x) + self.ny = len(y) + fgoutX = self.X # will be computed from info above + fgoutY = self.Y # will be computed from info above + print(' Grid edges: ',self.extent_edges) + print(' with %i cells in x and %i cells in y' \ + % (self.nx,self.ny)) + + + if not numpy.allclose(X, fgoutX): + errmsg = '*** X read from output does not match fgout_grid.X' + raise ValueError(errmsg) + + if not numpy.allclose(Y, fgoutY): + errmsg = '*** Y read from output does not match fgout_grid.Y' + raise ValueError(errmsg) + + return fgout_frame + + +# ======================== +# Functions for interpolating from fgout grid to arbitrary points, +# useful for example if using velocity field to model particle/debris motion + +def make_fgout_fcn_xy(fgout, qoi, method='nearest', + bounds_error=False, fill_value=numpy.nan): + """ + Create a function that can be called at (x,y) and return the qoi + interpolated in space from the fgout array. + + qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to + an attribute of fgout. + + The function returned takes arguments x,y that can be floats or + (equal length) 1D arrays of values that lie within the spatial + extent of fgout. + + bounds_error and fill_value determine the behavior if (x,y) is not in + the bounds of the data, as in scipy.interpolate.RegularGridInterpolator. + """ + + from scipy.interpolate import RegularGridInterpolator + + try: + q = getattr(fgout,qoi) + except: + print('*** fgout missing attribute qoi = %s?' % qoi) + + err_msg = '*** q must have same shape as fgout.X\n' \ + + 'fgout.X.shape = %s, q.shape = %s' % (fgout.X.shape,q.shape) + assert fgout.X.shape == q.shape, err_msg + + x1 = fgout.X[:,0] + y1 = fgout.Y[0,:] + fgout_fcn1 = RegularGridInterpolator((x1,y1), q, method=method, + bounds_error=bounds_error, fill_value=fill_value) + + def fgout_fcn(x,y): + """ + Function that can be evaluated at single point or arrays (x,y). + """ + from numpy import array, vstack + xa = array(x) + ya = array(y) + xyout = vstack((xa,ya)).T + qout = fgout_fcn1(xyout) + if len(qout) == 1: + qout = qout[0] # return scalar + return qout + + return fgout_fcn + + +def make_fgout_fcn_xyt(fgout1, fgout2, qoi, method_xy='nearest', + method_t='linear', bounds_error=False, + fill_value=numpy.nan): + """ + Create a function that can be called at (x,y,t) and return the qoi + interpolated in space and time between the two frames fgout1 and fgout2. + + qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to + an attribute of fgout. + + method_xy is the method used in creating the spatial interpolator, + and is passed to make_fgout_fcn_xy. + + method_t is the method used for interpolation in time, currently only + 'linear' is supported, which linearly interpolates. + + bounds_error and fill_value determine the behavior if (x,y,t) is not in + the bounds of the data. + + The function returned takes arguments x,y (floats or equal-length 1D arrays) + of values that lie within the spatial extent of fgout1, fgout2 + (which are assumed to cover the same uniform grid at different times) + and t should be a float that lies between fgout1.t and fgout2.t. + """ + + assert numpy.allclose(fgout1.X, fgout2.X), \ + '*** fgout1 and fgout2 must have same X' + assert numpy.allclose(fgout1.Y, fgout2.Y), \ + '*** fgout1 and fgout2 must have same Y' + + t1 = fgout1.t + t2 = fgout2.t + #assert t1 < t2, '*** expected fgout1.t < fgout2.t' + + fgout1_fcn_xy = make_fgout_fcn_xy(fgout1, qoi, method=method_xy, + bounds_error=bounds_error, fill_value=fill_value) + fgout2_fcn_xy = make_fgout_fcn_xy(fgout2, qoi, method=method_xy, + bounds_error=bounds_error, fill_value=fill_value) + + def fgout_fcn(x,y,t): + """ + Function that can be evaluated at single point or arrays (x,y) + at a single time t. + """ + from numpy import array, ones + xa = array(x) + ya = array(y) + tol = 1e-6 # to make sure it works ok when called with t=t1 or t=t2 + if t1-tol <= t <= t2+tol: + alpha = (t-t1)/(t2-t1) + elif bounds_error: + errmsg = '*** argument t=%g should be between t1=%g and t2=%g' \ + % (t,t1,t2) + raise InputError(errmsg) + else: + qout = fill_value * ones(xa.shape) + return qout + + qout1 = fgout1_fcn_xy(x,y) + qout2 = fgout2_fcn_xy(x,y) + + if method_t == 'linear': + if t1 <= t <= t2: + alpha = (t-t1)/(t2-t1) + qout = (1-alpha)*qout1 + alpha*qout2 + else: + raise NotImplementedError('method_t = %s not supported' % method_t) + + return qout + + return fgout_fcn + +# =============================== +# Functions for writing a set of fgout frames as a netCDF file, and +# reading such a file: + +def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc', + qois = ['h','hu','hv','eta'], datatype='f4', + include_B0=False, include_Bfinal=False, + description='', verbose=True): + """ + Write a list of fgout frames (at different times on the same rectangular + grid) to a single netCDF file, with some metadata and the topography, + if desired. + + fgout_frames should be a list of FGoutFrame objects, all of the same size + and at increasing times. + + fname_nc is the name of the file to write. + + qois is a list of strings, the quantities of interest to include in file. + This could include any of: + + 'h', 'eta', 'hu', 'hv', 'u', 'v', 's', 'hss', 'B'. + + All other quantities can be computed from h, hu, hv, eta, + the original fgout variables from GeoClaw, but for some applications + you might only want to save 'h' and 's', for example. + + datatype should be 'f4' [default] or 'f8', specifying bytes per qoi value. + 'f8' has full precision of the original data, but the file will be + twice as large and may not be needed for downstream applications. + + Note that the topography B = eta - h, so it is not necessary to store all + three of these. Also, B is often the same for all frames, so rather than + storing B at each frame as a qoi, two other options are also provided + (and then storing eta or h for all frames allows calculating the other): + + include_Bfinal: If True, include the topography B array from the final frame + as the Bfinal array. + + include_B0: If True, include the topography B array from the first frame + as the B0 array. This is only useful if, e.g., the first frame is initial + topography before co-seismic deformation, and at later times the topography + is always equal to Bfinal. + + `description` is a string that will be added as metadata. + A metadata field `history` will also be added, which includes the + time the file was created and the path to the directory where it was made. + """ + + import netCDF4 + from datetime import datetime, timedelta + from cftime import num2date, date2num + import time + timestr = time.ctime(time.time()) # current time for metadata + + fg_times = numpy.array([fg.t for fg in fgout_frames]) + + if verbose: + print('Creating %s with fgout frames at times: ' % fname_nc) + print(fg_times) + + fg0 = fgout_frames[0] + x = fg0.x + y = fg0.y + + xs = numpy.array([fg.x for fg in fgout_frames]) + ys = numpy.array([fg.y for fg in fgout_frames]) + # assert same for all times + + units = {'h':'meters', 'eta':'meters', 'hu':'m^2/s', 'hv':'m^2/s', + 'u':'m/s', 'v':'m/s', 's':'m/s', 'hss':'m^3/s^2', 'B':'meters'} + + with netCDF4.Dataset(fname_nc, 'w') as rootgrp: + + rootgrp.description = description + rootgrp.history = "Created " + timestr + rootgrp.history += " in %s; " % os.getcwd() + + lon = rootgrp.createDimension('lon', len(x)) + longitudes = rootgrp.createVariable('lon','f8',('lon',)) + longitudes[:] = x + longitudes.units = 'degrees_east' + + lat = rootgrp.createDimension('lat', len(y)) + latitudes = rootgrp.createVariable('lat','f8',('lat',)) + latitudes[:] = y + latitudes.units = 'degrees_north' + + time = rootgrp.createDimension('time', len(fg_times)) + times = rootgrp.createVariable('time','f8',('time',)) + times[:] = fg_times + times.units = 'seconds' + + if 0: + # Could make times be datetimes relative to some event time, e.g.: + times.units = 'seconds since 1700-01-26 21:00:00.0' + times.calendar = 'gregorian' + dates = [datetime(1700,1,26,21) + timedelta(seconds=ss) \ + for ss in fg_times] + times[:] = date2num(dates,units=times.units,calendar=times.calendar) + + if include_B0: + B0 = rootgrp.createVariable('B0',datatype,('lon','lat',)) + B0[:,:] = fg0.B + B0.units = 'meters' + + if include_Bfinal: + fg_final = fgout_frames[-1] + Bfinal = rootgrp.createVariable('Bfinal',datatype,('lon','lat',)) + Bfinal[:,:] = fg_final.B + Bfinal.units = 'meters' + + for qoi in qois: + qoi_frames = [getattr(fgout,qoi) for fgout in fgout_frames] + qoi_var = rootgrp.createVariable(qoi,datatype,('time','lon','lat',)) + qoi_var[:,:,:] = qoi_frames + qoi_var.units = units[qoi] + +def get_as_array(var, rootgrp, verbose=True): + """ + Utility function to retrieve variable from netCDF file and convert to + numpy array. + """ + a = rootgrp.variables.get(var, None) + if a is not None: + if verbose: print(' Loaded %s with shape %s' % (var,repr(a.shape))) + return numpy.array(a) + else: + if verbose: print(' Did not find %s' % var) + return None + + +def read_netcdf_arrays(fname_nc, qois, verbose=True): + """ + Read a netCDF file and extract the quantities of interest denoted by + strings in the list qois, which can include: + + 'h', 'eta', 'hu', 'hv', 'u', 'v', 's', 'hss', 'B'. + + qois can also include the time-independent 'B0' and/or 'Bfinal' + + Returns + + x, y, t, qoi_arrays + + where x,y define the longitude, latitudes, t is the times of the frames, + and qoi_arrays is a dictionary indexed by the strings from qois. + + Each dict element is an array with shape (len(t), len(x), len(y)) + for time-dependent qoi's, or (len(x), len(y)) for B0 or Bfinal, + or None if that qoi was not found in the netCDF file. + """ + + import netCDF4 + + with netCDF4.Dataset(fname_nc, 'r') as rootgrp: + if verbose: + print('Reading data to fgout frames from nc file',fname_nc) + print(' nc file description: ', rootgrp.description) + print('History: ', rootgrp.history) + + x = get_as_array('lon', rootgrp, verbose) + y = get_as_array('lat', rootgrp, verbose) + t = get_as_array('time', rootgrp, verbose) + + vars = list(rootgrp.variables) + + qoi_arrays = {} + for qoi in qois: + qoi_array = get_as_array(qoi, rootgrp, verbose) + qoi_arrays[qoi] = qoi_array + + return x,y,t,qoi_arrays + + + +def read_netcdf(fname_nc, fgout_grid=None, verbose=True): + """ + Read a netCDF file and return a list of FGoutFrame instances. + This will only be possible if the netCDF file contains at least + the qoi's 'h','hu','hv','eta' required to reconstruct the q array + as output by GeoClaw. + """ + + import netCDF4 + + with netCDF4.Dataset(fname_nc, 'r') as rootgrp: + if 1: + print('Reading data to fgout frames from nc file',fname_nc) + print(' nc file description: ', rootgrp.description) + print('History: ', rootgrp.history) + + x = get_as_array('lon', rootgrp, verbose) + y = get_as_array('lat', rootgrp, verbose) + t = get_as_array('time', rootgrp, verbose) + + vars = list(rootgrp.variables) + + + for qoi in ['h','hu','hv','eta']: + errmsg = '*** Cannot reconstruct fgout frame without %s' % qoi + assert qoi in vars, errmsg + h = get_as_array('h', rootgrp, verbose) + hu = get_as_array('hu', rootgrp, verbose) + hv = get_as_array('hv', rootgrp, verbose) + eta = get_as_array('eta', rootgrp, verbose) + + if (x is None) or (y is None): + print('*** Could not create grid') + else: + X,Y = numpy.meshgrid(x,y) + + fgout_frames = [] + + for k in range(eta.shape[0]): + fgout = FGoutFrame(fgout_grid=fgout_grid, frameno=k) + fgout.x = x + fgout.y = y + fgout.t = t[k] + fgout.q = numpy.empty((4,eta.shape[1],eta.shape[2])) + fgout.q[0,:,:] = h[k,:,:] + fgout.q[1,:,:] = hu[k,:,:] + fgout.q[2,:,:] = hv[k,:,:] + fgout.q[3,:,:] = eta[k,:,:] + fgout.X = X + fgout.Y = Y + fgout_frames.append(fgout) + + print('Created fgout_frames as list of length %i' % len(fgout_frames)) + + return fgout_frames + +def print_netcdf_info(fname_nc): + """ + Print out info about the contents of a netCDF file contining fgout frames, + written using write_netcdf. + """ + import netCDF4 + + with netCDF4.Dataset(fname_nc, 'r') as rootgrp: + x = get_as_array('lon', rootgrp, verbose=False) + y = get_as_array('lat', rootgrp, verbose=False) + t = get_as_array('time', rootgrp, verbose=False) + + vars = list(rootgrp.variables) + + print('===================================================') + print('netCDF file %s contains:' % fname_nc) + print('description: \n', rootgrp.description) + print('history: \n', rootgrp.history) + print('%i longitudes from %.6f to %.6f' % (len(x),x[0],x[-1])) + print('%i latitudes from %.6f to %.6f' % (len(y),y[0],y[-1])) + print('%i times from %.3f to %.3f' % (len(t),t[0],t[-1])) + print('variables: ',vars) + print('===================================================') + + diff --git a/src/python/geoclaw/kmltools.py b/src/python/geoclaw/kmltools.py index 27171a597..6d4c69f49 100644 --- a/src/python/geoclaw/kmltools.py +++ b/src/python/geoclaw/kmltools.py @@ -21,6 +21,7 @@ - topo2kml - create a kml outline for each topo grid specified in setrun - dtopo2kml - create a kml outline for each dtopo grid specified in setrun - fgmax2kml - create a kml outline for each fgmax grid specified in setrun + - fgout2kml - create a kml outline for each fgout grid specified in setrun - make_input_data_kmls - make kml files for many things specified in setrun - pcolorcells_for_kml - version of pcolormesh with appropriate dpi and size - png2kml - create kml file wrapping a png figure to be viewed on GE @@ -962,6 +963,53 @@ def fgmax2kml(rundata=None,fname='fgmax_grids.kml',verbose=True,combined=False): % fg.point_style) +def fgout2kml(rundata=None,fname='fgout_grids.kml',verbose=True,combined=False): + + """ + Create a KML box for each fgout grid specified for a GeoClaw run. + + :Inputs: + + - *rundata* - an object of class *ClawRunData* or None + + If *rundata==None*, try to create based on executing function *setrun* + from the `setrun.py` file in the current directory. + + - *fname* (str) - resulting kml file. + + - *verbose* (bool) - If *True*, print out info about each region found + + - *combined* (bool) - If *True*, combine into single kml file with + name given by *fname*. NOT YET IMPLEMENTED. + If False, *fname* is ignored and individual files are created for + each fgout grid. + + """ + + from numpy import cos,pi,floor + + if rundata is None: + try: + import setrun + reload(setrun) + rundata = setrun.setrun() + except: + raise IOError("*** cannot execute setrun file") + + if combined: + fname_combined = 'fgout_grids.kml' + print('*** combined fgout kml files not yet supported') + print(' making a kml file for each fgout grid') + + fgout_grids = rundata.fgout_data.fgout_grids + + for fg in fgout_grids: + fname_root = 'fgout%s' % str(fg.fgno).zfill(4) + kml_file = fname_root + '.kml' + xy = ([fg.x1,fg.x2], [fg.y1,fg.y2]) + box2kml(xy, kml_file, fname_root, color='8888FF') + + def make_input_data_kmls(rundata=None, combined=False): """ Produce kml files for the computational domain, all gauges and regions @@ -993,6 +1041,7 @@ def make_input_data_kmls(rundata=None, combined=False): regions2kml(rundata, combined=combined) gauges2kml(rundata) fgmax2kml(rundata) + fgout2kml(rundata) topofiles = rundata.topo_data.topofiles for f in topofiles: diff --git a/tests/bowl_slosh/setrun.py b/tests/bowl_slosh/setrun.py index 13864be4d..1d2120913 100644 --- a/tests/bowl_slosh/setrun.py +++ b/tests/bowl_slosh/setrun.py @@ -381,12 +381,12 @@ def setgeo(rundata): # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) # [fname] - # == setfixedgrids.data values == - fixedgrids = rundata.fixed_grid_data - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] - + # == fgout_grids.data values == + # NEW IN v5.9.0 + # Set rundata.fgout_data.fgout_grids to be a list of + # objects of class clawpack.geoclaw.fgout_tools.FGoutGrid: + fgout_grids = rundata.fgout_data.fgout_grids # empty list initially + # == fgmax.data values == rundata.fgmax_data.num_fgmax_val = 2 fgmax_grids = rundata.fgmax_data.fgmax_grids # empty list to start