SinglerateScalingTsunami.py 2.26 KB
Newer Older
1
2
3
4
5
6
7
8
9
from dgpy import *
from mpi4py import MPI
from SRInputTsunami import *
from math import *
import time, os, sys

model = GModel()
model.load(filename + '_sr_partitioned%d.msh'%(Msg.GetCommSize()))

bseny's avatar
   
bseny committed
10
groups = dgGroupCollection(model, dimension, order)
11
12
13
14
15
16
XYZ = groups.getFunctionCoordinates();
bath = dgDofContainer(groups, 1);
bathname=outputDir+"/"+filename+"_bath_smooth"+"/"+filename+"_bath_smooth.idx"
Msg.Barrier()
if(os.path.exists(bathname)):
  Msg.Info('.... Import Default Smoothed Bathymetry')
17
  bath.importIdx(bathname)
18
19
20
21
22
23
24
25
26
27
28
29
else: 
  Msg.Fatal("No correct bathymetry file,  launch DiffuseBathymetryTsunami.py")

claw = dgConservationLawShallowWater2d()
solution = dgDofContainer(groups, claw.getNbFields())
Finit = functionC(tlib, "initial_condition_okada_honshu",3, [XYZ])
solution.interpolate(Finit)
FCoriolis= functionC(tlib, "coriolis", 1, [XYZ])
claw.setBathymetry(bath.getFunction())
claw.setBathymetryGradient(bath.getFunctionGradient())
claw.setCoriolisFactor(FCoriolis)
claw.setIsSpherical(6371220.)
bseny's avatar
   
bseny committed
30
claw.setIsLinear(True)
31
32
wall=claw.newBoundaryWall()
claw.addBoundaryCondition("Coast", wall)
bseny's avatar
   
bseny committed
33
34
35
36
#bottomDrag = functionC(tlib, "bottomDrag", 1, [function.getSolution(), bath.getFunction()])
#claw.setQuadraticDissipation(bottomDrag)
bottomDragLinear = functionC(tlib, "bottomDragLinear", 1, [function.getSolution(), bath.getFunction()])
claw.setLinearDissipation(bottomDragLinear)
37
38
39
40
41
42
43

rk=dgMultirateERK(groups, claw, RK_TYPE)
dt=rk.splitGroupsForMultirate(mL, solution, [solution, bath])
rk.printMultirateInfo(0)

t=Ti
if(Msg.GetCommRank() == 0):
44
  print 'SU=', rk.speedUp()
45
46
47
48
49
50
51
52
53
54
55
56
Msg.Barrier()
rk.initialize()
Msg.Barrier()
timer=dgTimer.root()
timer.reset()
Msg.Barrier()
timer.start()
for i in range (0,nbSteps) :
  rk.iterate(solution, dt, t)
  t = t + dt
timer.stop()
Msg.Barrier()
57
58
minRatioProc=rk.nbInterfaceTerms2()/float(rk.nbVolumeTerms());
maxRatioProc=rk.nbInterfaceTerms2()/float(rk.nbVolumeTerms());
59
60
61
62
63
64
minRatio = MPI.COMM_WORLD.allreduce(minRatioProc,   op=MPI.MIN)
maxRatio = MPI.COMM_WORLD.allreduce(maxRatioProc,   op=MPI.MAX)
maxTime = MPI.COMM_WORLD.allreduce(timer.get(),   op=MPI.MAX)
minTime = MPI.COMM_WORLD.allreduce(timer.get(),   op=MPI.MIN)
norm=solution.norm()
if (Msg.GetCommRank()==0):
bseny's avatar
   
bseny committed
65
    print "%.8f %.8f %.8f %.8f %.8f %.8f"%(t, norm,  minRatio, maxRatio, minTime,  maxTime)
66
Msg.Exit(0)