Commit d9d6c230 authored by zhang's avatar zhang
Browse files

atm 2d test for mass conservation

git-svn-id: https://geuz.org/svn/gmsh/trunk@20770 59c8cabf-10d6-4a28-b43f-3936b85956c4
parent ae3db8eb
......@@ -170,10 +170,7 @@ for iruns in range(0,9) :
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 = dgIMEXRK(claw, dofIm, IMEX_ARK4) #timeorder
timeIter.getNewton().setVerb(1) #Verbosity
solution.copy(solution0)
......@@ -235,7 +232,6 @@ for iruns in range(0,9) :
# maxDiff = fullMatrixDouble(claw.getNbFields(),1)
t=Ti
#t=timeIter.starter(solution,dt,t,3)
n_export=0
timeStart=time.clock();
......
nbx=15;
Point(1) = {0, 0, 0};
Point(2) = {1000, 0, 0};
Line(1) = {1, 2};
Transfinite Line {1} = nbx+1;
Physical Line("domain") = {1};
Physical Point("left") = {1};
Physical Point("right") = {2};
from partition1d import *
from dgpy import *
from extrude import *
from rotateZtoY import *
name='line'
nbPart = int(sys.argv[1])
nbLayers = 15
H = 1000.
res=H/nbLayers
def getLayersSigma (element, vertex) :
z=0
zTab=[z]
for i in range(nbLayers):
z=z-res
zTab.append(z)
return zTab
if nbPart>1:
partStr='_part_%i' % nbPart
else:
partStr=''
partition1d(name+'.msh', name + partStr + '.msh', nbPart)
extrude(mesh(name + partStr + '.msh'), getLayersSigma).write(name + partStr + '_2d.msh')
rotateZtoY(name + partStr + '_2d.msh',name + partStr + '_2d_XY.msh')
Msg.Exit(0)
from dgpy import *
from math import *
import time, os, sys
import gmshPartition
import numpy as np
exportData=0
outFile=open('periodic_convection_2d_imexark4_poly4_200s.dat','w')
massFile=open('periodic_convection_2d_imexark4_poly4_200s.mass','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();
def initialCondition(cmap, FCT, XYZ) :
FCT[:,:] = 0
def boundaryForcing(cmap, FCT, XYZ, time, rhoHs, thetaHs, sol, normals) :
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
vn = sol[:,1]*normals[:,0] + sol[:,2]*normals[:,1]
FCT[:,0]=sol[:,0]
FCT[:,1]=sol[:,1]-2*vn*normals[:,0]
FCT[:,2]=sol[:,2]-2*vn*normals[:,1]
FCT[:,3]=(rhoHs+sol[:,0])*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, function.getSolution(), function.getNormals()])
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 = dgIMEXRK(claw, dofIm, IMEX_ARK4) #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
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))
massFile.write('{0} {1:.15e} {2:.15e} {3:.15e} \n'.format(nbSteps,massInit,intMassVol(0,0),(intMassVol(0,0)-massInit)/massInit))
outFile.close()
massFile.close()
Msg.Exit(0)
from dgpy import *
from math import *
import time, os, sys
import gmshPartition
import numpy as np
#### 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=300.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')
groups = dgGroupCollection(model, dimension, order)
groups.splitGroupsByPhysicalTag();
claw = dgEulerAtmLaw(dimension)
solution = dgDofContainer(groups, claw.getNbFields())
XYZ = groups.getFunctionCoordinates();
def initialCondition(cmap, FCT, XYZ) :
FCT[:,:] = 0
def boundaryForcing(cmap, FCT, XYZ, time, rhoHs, thetaHs, sol, normals) :
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
vn = sol[:,1]*normals[:,0] + sol[:,2]*normals[:,1]
FCT[:,0]=sol[:,0]
FCT[:,1]=sol[:,1]-2*vn*normals[:,0]
FCT[:,2]=sol[:,2]-2*vn*normals[:,1]
FCT[:,3]=(rhoHs+sol[:,0])*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, function.getSolution(), function.getNormals()])
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)
timeIter = dgERK(claw,None,DG_ERK_44)
# claw.setFilterMode(FILTER_LINEAR)
# petscIm = linearSystemPETScBlockDouble()
# petscIm.setParameter("petscOptions","-ksp_type preonly -pc_type lu")
# #petscIm.setParameter("petscOptions", "-pc_type lu -ksp_type preonly -pc_factor_mat_solver_package mumps")
# dofIm = dgDofManager.newDGBlock(groups, claw.getNbFields(), petscIm)
# #timeIter = dgIMEXTSRK(claw, dofIm, 4) #timeorder
# timeOrder=3
# timeIter = dgIMEXRK(claw, dofIm, timeOrder) #timeorder
# #timeIter = dgDIRK(claw, dofIm, 2) #timeorder
# timeIter.getNewton().setVerb(1) #Verbosity
#dt=claw.getMinOfTimeSteps(solution)*4.0/(2*timeOrder+1) #infinite because no initial velocity
dt=Tf/60000
if (Msg.GetCommRank()==0):
print ("Time step:",dt)
nbSteps = int(round(Tf/dt))
#nbSteps=200
#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[:]
nCompExp=[2,1,1,1]
namesExp=["uv","rhop","thetap","pp"]
Exp=functionNumpy(sum(nCompExp), getExp, [uv, rhop, thetap, pp])
# 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
start = time.clock()
j=1
for i in range(0,nbSteps):
if (i%1000 == 0):
solution.exportFunctionVtk(Exp,'refSolERK100/export', 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)
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);</