| 1 |
# |
| 2 |
# upwinding test moving a Gaussian hill around |
| 3 |
# |
| 4 |
# we solve U_,t + v_i u_,i =0 |
| 5 |
# |
| 6 |
# the solution is given as u(x,t)=1/(4*pi*E*t)^{dim/2} * exp ( - |x-x_0(t)|^2/(4*E*t) ) |
| 7 |
# |
| 8 |
# where x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=2 and |
| 9 |
# |
| 10 |
# x_0(t) = [ cos(OMEGA0*T0)*0.5,-sin(OMEGA0*T0)*0.5 ] and v=[-y,x]*OMEGA0 for dim=3 |
| 11 |
# |
| 12 |
# the solution is started from some time T0>0. |
| 13 |
# |
| 14 |
# We are using five quality messurements for u_h |
| 15 |
# |
| 16 |
# - inf(u_h) > 0 |
| 17 |
# - sup(u_h)/sup(u(x,t)) = sup(u_h)*(4*pi*E*t)^{dim/2} ~ 1 |
| 18 |
# - integrate(u_h) ~ 1 |
| 19 |
# - | x_0h-x_0 | ~ 0 where x_0h = integrate(x*u_h) |
| 20 |
# - sigma_h/4*E*t ~ 1 where sigma_h=sqrt(integrate(length(x-x0h)**2 * u_h) * (DIM==3 ? sqrt(2./3.) :1 ) |
| 21 |
# |
| 22 |
# |
| 23 |
from esys.escript import * |
| 24 |
from esys.escript.linearPDEs import LinearSinglePDE, TransportPDE |
| 25 |
from esys.finley import Rectangle, Brick |
| 26 |
from math import pi, ceil |
| 27 |
NE=128 |
| 28 |
NE=4 |
| 29 |
DIM=2 |
| 30 |
THETA=0.5 |
| 31 |
OMEGA0=1. |
| 32 |
ALPHA=pi/4 |
| 33 |
T0=0 |
| 34 |
T_END=2.*pi |
| 35 |
dt=1e-3*10*10 |
| 36 |
E=1.e-3 |
| 37 |
TEST_SUPG=False or True |
| 38 |
|
| 39 |
|
| 40 |
if DIM==2: |
| 41 |
dom=Rectangle(NE,NE) |
| 42 |
else: |
| 43 |
dom=Brick(NE,NE,NE) |
| 44 |
u0=dom.getX()[0] |
| 45 |
# saveVTK("u.%s.xml"%0,u=u0) |
| 46 |
# print "XX"*80 |
| 47 |
dom.setX(2*dom.getX()-1) |
| 48 |
|
| 49 |
# set initial value |
| 50 |
x=dom.getX() |
| 51 |
r=sqrt(x[0]**2+(x[1]-1./3.)**2) |
| 52 |
# u0=whereNegative(r-1./3.)*wherePositive(wherePositive(abs(x[0])-0.05)+wherePositive(x[1]-0.5)) |
| 53 |
|
| 54 |
x=Function(dom).getX() |
| 55 |
if DIM == 2: |
| 56 |
V=OMEGA0*(x[0]*[0,-1]+x[1]*[1,0]) |
| 57 |
else: |
| 58 |
V=OMEGA0*(x[0]*[0,cos(ALPHA),0]+x[1]*[-cos(ALPHA),0,sin(ALPHA)]+x[2]*[0.,-sin(ALPHA),0.]) |
| 59 |
#=================== |
| 60 |
fc=TransportPDE(dom,num_equations=1,theta=THETA) |
| 61 |
x=Function(dom).getX() |
| 62 |
fc.setValue(M=Scalar(1.,Function(dom)),C=V) |
| 63 |
#============== |
| 64 |
if TEST_SUPG: |
| 65 |
supg=LinearSinglePDE(dom) |
| 66 |
supg.setValue(D=1.) |
| 67 |
supg.setSolverMethod(supg.LUMPING) |
| 68 |
dt_supg=inf(dom.getSize()/length(V)) |
| 69 |
u_supg=u0*1. |
| 70 |
|
| 71 |
c=0 |
| 72 |
# saveVTK("u.%s.xml"%c,u=u0) |
| 73 |
fc.setInitialSolution(u0) |
| 74 |
t=T0 |
| 75 |
print "QUALITY FCT: time = %s pi"%(t/pi),inf(u0),sup(u0),integrate(u0) |
| 76 |
while t<T_END: |
| 77 |
print "time step t=",t+dt |
| 78 |
u=fc.solve(dt, verbose=True) |
| 79 |
print "QUALITY FCT: time = %s pi"%(t+dt/pi),inf(u),sup(u),integrate(u) |
| 80 |
if TEST_SUPG: |
| 81 |
#========== supg tests ================ |
| 82 |
nn=max(ceil(dt/dt_supg),1.) |
| 83 |
dt2=dt/nn |
| 84 |
nnn=0 |
| 85 |
while nnn<nn : |
| 86 |
supg.setValue(Y=u_supg+dt2/2*inner(V,grad(u_supg))) |
| 87 |
u2=supg.getSolution() |
| 88 |
supg.setValue(Y=u_supg+dt2*inner(V,grad(u2))) |
| 89 |
u_supg=supg.getSolution() |
| 90 |
nnn+=1 |
| 91 |
c+=1 |
| 92 |
t+=dt |
| 93 |
if TEST_SUPG: |
| 94 |
print "QUALITY SUPG: time = %s pi"%(t/pi),inf(u_supg),sup(u_supg),integrate(u_supg) |
| 95 |
# saveVTK("u2.%s.xml"%c,u=u,u_supg=u_supg) |
| 96 |
else: |
| 97 |
# saveVTK("u.%s.xml"%c,u=u) |
| 98 |
pass |
| 99 |
# if c == 20: 1/0 |
Parent Directory
| 
Revision Log