Commit e31ecf15 authored by zhang's avatar zhang
Browse files

atm:ark4 loses mass

git-svn-id: https://geuz.org/svn/gmsh/trunk@20711 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent 9c616bb3
from dgpy import *
from math import *
import time, os, sys
import gmshPartition
import numpy as np
exportData=0
outFile=open('periodic_convection_2d_imexrk4_poly4_200s.dat','w')
#### Physical constants ###
gamma=1.4
Rd=287.0
p0=1.0e5
g=9.80616
Cv=Rd/(gamma-1.0)
Cp=Rd*(1.0+1.0/(gamma-1.0))
###########################
#### Begin and end time ###
Ti=0.0
Tf=200.0
###########################
#Period of the oscillations
period = 100
order=4
dimension=2
model = GModel()
name='line'
if Msg.GetCommSize()>1:
partStr='_part_%i' % Msg.GetCommSize()
else:
partStr=''
model.load(name + partStr + '_2d_XY.msh')
############# for convergence plot################
###################################################
errp=[]
erru=[]
errv=[]
errthetap=[]
errall=[]
stepsizes=[]
groups = dgGroupCollection(model, dimension, order)
groups.splitGroupsByPhysicalTag();
claw = dgEulerAtmLaw(dimension)
####Load reference solution from file###
reloaded =dgDofContainer(groups,claw.getNbFields())
reloaded.readMsh("refSolERK"+str(period)+"/subwindow4/subwindow4.idx")
solution = dgDofContainer(groups, claw.getNbFields())
XYZ = groups.getFunctionCoordinates();
xyzCircle = XYZ#transformToCircle(circleTransform)
integOrder = 2*order+1+2
xyzCircleInteg = functionPrecomputed(groups, integOrder, 3)
xyzCircleInteg.compute(xyzCircle)
def initialCondition(cmap, FCT, XYZ) :
FCT[:,:] = 0
def boundaryForcing(cmap, FCT, XYZ, time, rhoHs, thetaHs) :
thetac=5
xc=500.0
zc=0.0
rc=250.0
currentTime=time[0]
r=np.sqrt((XYZ[:,0]-xc)**2+(XYZ[:,1]-zc)**2)
thetaPert=np.empty(thetaHs.shape)
outR=(r>rc)
thetaPert[outR]=thetaHs[outR]
inR=np.invert(outR)
thetaPert[inR]=thetaHs[inR]+thetac/2.0*(1.0+np.cos(pi*r[inR]/rc))*np.sin(currentTime/period*2.0*pi)**2
FCT[:,0]=0
FCT[:,1]=0
FCT[:,2]=0
FCT[:,3]=rhoHs*thetaPert-rhoHs*thetaHs
def thetaHydrostatic(cmap, FCT) :
FCT[:]=300
def exnerHydrostatic(cmap, FCT, XYZ, thetaHs) :
FCT[:]=1.0-g/(Cp*thetaHs)*XYZ[:,1]
def rhoHydrostatic(cmap, FCT, thetaHs, exnerHs) :
FCT[:]=p0/(Rd*thetaHs)*exnerHs**(Cv/Rd)
def rhoThetaHydrostatic(cmap, FCT, rhoHs, thetaHs) :
FCT[:]=rhoHs*thetaHs
def getVelocity(cmap, FCT, sol, rhoHs) :
invrho=1/(rhoHs+sol[:,0])
FCT[:,0]=sol[:,1]*invrho
FCT[:,1]=sol[:,2]*invrho
FCT[:,2]=0
def getRhop(cmap, FCT, sol) :
FCT[:]=sol[:,0]
def getpp(cmap, FCT, sol, rhoHs, thetaHs) :
rho=rhoHs+sol[:,0]
rhoTheta = rhoHs*thetaHs+sol[:,3]
pHs=p0*(rhoHs*thetaHs*Rd/p0)**gamma
p=p0*(rhoTheta*Rd/p0)**gamma
FCT[:]=p-pHs
def getpHs(cmap, FCT, sol, rhoHs, thetaHs) :
pHs=p0*(rhoHs*thetaHs*Rd/p0)**gamma
FCT[:]=pHs
def getp(cmap, FCT, sol, rhoHs, thetaHs) :
rhoTheta = rhoHs*thetaHs+sol[:,3]
p=p0*(rhoTheta*Rd/p0)**gamma
FCT[:]=p
def getThetap(cmap, FCT, sol, rhoHs, thetaHs) :
rho=rhoHs+sol[:,0]
FCT[:]=(sol[:,3]+(rhoHs*thetaHs))/rho-thetaHs
thetaHsC = functionNumpy(1, thetaHydrostatic,[])
exnerHsC = functionNumpy(1, exnerHydrostatic, [XYZ, thetaHsC])
rhoHsC = functionNumpy(1, rhoHydrostatic, [thetaHsC, exnerHsC])
rhoThetaHsC = functionNumpy(1, rhoThetaHydrostatic, [rhoHsC, thetaHsC])
uv=functionNumpy(3, getVelocity, [solution.getFunction(), rhoHsC])
rhop=functionNumpy(1, getRhop, [solution.getFunction()])
pp=functionNumpy(1, getpp, [solution.getFunction(), rhoHsC, thetaHsC])
thetap=functionNumpy(1, getThetap, [solution.getFunction(), rhoHsC, thetaHsC])
initF=functionNumpy(claw.getNbFields(), initialCondition, [XYZ])
solution.interpolate(initF)
rhoHs = dgDofContainer(groups, 1)
rhoHs.interpolate(rhoHsC)
rhoThetaHs = dgDofContainer(groups, 1)
rhoThetaHs.interpolate(rhoThetaHsC)
claw.setHydrostaticState(rhoHs,rhoThetaHs)
claw.setPhysicalConstants(gamma,Rd,p0,g)
boundaryWall = claw.newBoundaryWall()
timeFunction = function.getTime()
boundaryForcingF=functionNumpy(4, boundaryForcing, [XYZ, timeFunction, rhoHsC, thetaHsC])
toReplace = VectorFunctorConst([function.getSolution()])
replaceBy = VectorFunctorConst([boundaryForcingF])
boundaryHeated = claw.newOutsideValueBoundaryGeneric("",toReplace,replaceBy)
claw.addBoundaryCondition('bottom_domain', boundaryWall)
claw.addBoundaryCondition('top_domain', boundaryHeated) #top is actually bottom due to the inverse extrusion
claw.addBoundaryCondition('left', boundaryWall)
claw.addBoundaryCondition('right', boundaryWall)
claw.setFilterMode(FILTER_LINEAR)
solution0 = dgDofContainer(groups, claw.getNbFields())
solution0.copy(solution)
#dt = 0.75
#dt = 1.0
#nbSteps = int(round((Tf-Ti)/dt))
ampf = 1.5
nbSteps = 400
dt=(Tf-Ti)/nbSteps
for iruns in range(0,9) :
petscIm = linearSystemPETScDouble()
petscIm.setParameter("petscOptions", "-pc_type lu -ksp_type preonly")
dofIm = dgDofManager.newDG(groups, claw, petscIm)
#timeIter = dgIMEXTSRK(claw, dofIm, 4) #timeorder
timeOrder=4
timeIter = dgIMEXRK(claw, dofIm, timeOrder) #timeorder
#timeIter = dgDIRK(claw, dofIm, 2) #timeorder
timeIter.getNewton().setVerb(1) #Verbosity
solution.copy(solution0)
#dt=claw.getMinOfTimeSteps(solution)*4.0/(2*timeOrder+1) #infinite because no initial velocity
if (iruns>0):
nbSteps = int(nbSteps*ampf)
dt=(Tf-Ti)/nbSteps
if (Msg.GetCommRank()==0):
print ("Time step:",dt,"No. of steps:",nbSteps)
# Compute relative error
def getRefSolSq(cmap, FCT, refSol) :
FCT[:,0]=refSol[:,0]**2
FCT[:,1]=refSol[:,1]**2
FCT[:,2]=refSol[:,2]**2
FCT[:,3]=refSol[:,3]**2
solRefSq=functionNumpy(4, getRefSolSq, [reloaded.getFunction()])
def getError(cmap, FCT, refSol, compSol) :
FCT[:,0]=(refSol[:,0]-compSol[:,0])**2
FCT[:,1]=(refSol[:,1]-compSol[:,1])**2
FCT[:,2]=(refSol[:,2]-compSol[:,2])**2
FCT[:,3]=(refSol[:,3]-compSol[:,3])**2
error = functionNumpy(4, getError, [reloaded.getFunction(), solution.getFunction()])
integratorError = dgFunctionIntegrator(groups, error)
intErr = fullMatrixDouble(4,1)
integratorRefSq = dgFunctionIntegrator(groups, solRefSq)
intRefSq = fullMatrixDouble(4,1)
#Export data
def getExp(cmap, FCT, uv, rhop, thetap, pp) :
FCT[:,0]=uv[:,0]
FCT[:,1]=uv[:,1]
FCT[:,2]=rhop[:]
FCT[:,3]=thetap[:]
FCT[:,4]=pp[:]
Exp=functionNumpy(5, getExp, [uv, rhop, thetap, pp, rhoHsC, thetaHsC])
nCompExp=[2,1,1,1,]
namesExp=["uv","rhop","thetap","pp"]
def densityOne(cmap, FCT, rhoHs, sol):
FCT[:,0] = rhoHs[:] + sol[:,0]
FCT[:,1] = 1
densityOneF=functionNumpy(2, densityOne, [rhoHs.getFunction(), solution.getFunction()])
massVolIntegrator=dgFunctionIntegrator(groups, densityOneF)
intMassVol = fullMatrixDouble(2,1)
massVolIntegrator.compute(intMassVol)
massInit = intMassVol(0,0)
prefixFileName = 'output_imexrk'+str(iruns)+'/export'
# sBV = 2.0/3*(order+1.0)
# muBV = 0.2
# filterBV = dgFilterBoydVandeven(groups, "", sBV, muBV, True, False)
# solutionBV = dgDofContainer(solution);
# minDiff = fullMatrixDouble(claw.getNbFields(),1)
# maxDiff = fullMatrixDouble(claw.getNbFields(),1)
t=Ti
#t=timeIter.starter(solution,dt,t,3)
n_export=0
timeStart=time.clock();
start = time.clock()
for i in range(0,nbSteps):
if (i%(250*2**iruns) == 0 and exportData):
solution.exportFunctionVtk(Exp,prefixFileName, t*1000, i,"solution",nCompExp,namesExp)
if (Msg.GetCommRank()==0):
print ('\nWriting output',n_export,'at time',t*1000,'ms and step',i,'over',nbSteps)
print("CFL: %.1e"%(dt/(claw.getMinOfTimeSteps(solution)*4.0)))
n_export=n_export+1
#We filter the dof manager and evaluate the max difference with unfiltered solution
# solutionBV.copy(solution)
# filterBV.apply(solutionBV)
# solutionBV.axpy(solution, -1)
# solutionBV.minmax(minDiff, maxDiff);
# if (Msg.GetCommRank()==0):
# #Print the max difference on density
# print("[%.1e]"%max(maxDiff(0,0),-minDiff(0,0)))
# print("CFL dt %.1e"%(claw.getMinOfTimeSteps(solution)*4.0/(2*timeOrder+1)))
##
timeIter.iterate(solution, dt, t)
t=t+dt
#if (Msg.GetCommRank()==0):
# sys.stdout.write('.')
# sys.stdout.flush()
if (Msg.GetCommRank()==0):
print 'dof',solution.getVector().size
print ('')
print 'Tf=',t,
elapsed = (time.clock() - start)
if (Msg.GetCommRank()==0):
print ('Time elapsed: ',elapsed,' s')
if(exportData) :
solution.exportFunctionVtk(Exp,prefixFileName, Tf*1000, nbSteps,"solution",nCompExp,namesExp)
integratorError.compute(intErr)
integratorRefSq.compute(intRefSq)
errp.append(sqrt(intErr(0,0)/intRefSq(0,0)))
erru.append(sqrt(intErr(1,0)/intRefSq(1,0)))
errv.append(sqrt(intErr(2,0)/intRefSq(2,0)))
errthetap.append(sqrt(intErr(3,0)/intRefSq(3,0)))
errall.append( sqrt((intErr(0,0)+intErr(1,0)+intErr(2,0)+intErr(3,0))/(intRefSq(0,0)+intRefSq(1,0)+intRefSq(2,0)+intRefSq(3,0))) )
stepsizes.append(dt)
massVolIntegrator.compute(intMassVol)
if (Msg.GetCommRank()==0):
print "\nError on rhop:",errp[-1],"rhou:",erru[-1],"rhov:", errv[-1],"rhothetap:",errthetap[-1],"all:",errall[-1]
if(iruns>0):
print "rhop order:",log(errp[-2]/errp[-1])/log(ampf)
print "rhou order:",log(erru[-2]/erru[-1])/log(ampf)
print "rhov order:",log(errv[-2]/errv[-1])/log(ampf)
print "rhothetap order:",log(errthetap[-2]/errthetap[-1])/log(ampf)
print "overall order:",log(errall[-2]/errall[-1])/log(ampf)
outFile.write('{} {} {} {} {} {} {}\n'.format(nbSteps,elapsed,errp[-1],erru[-1],errv[-1],errthetap[-1],errall[-1]))
print("------------------------\nMassInit %.25e \nMass %.25e \nDiff %e \nRelative diff %e"%(massInit, intMassVol(0,0), intMassVol(0,0)-massInit, (intMassVol(0,0)-massInit)/massInit))
outFile.close()
Msg.Exit(0)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment