doc/examples/mount.py

from esys.escript import *
from esys.escript.models import Mountains
from esys.finley import Brick,Rectangle
from math import pi, ceil

NE=16
DIM=3
H=1.
L=2*H
TOL=1.e-4
OMEGA=10
EPS=0.01
t=0
T_END=(2*pi)/OMEGA
n=0
if DIM==2:
  mydomain=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=1, useFullElementOrder=True,optimize=True)
else:
  mydomain=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=1, useFullElementOrder=True,optimize=True)
x=mydomain.getX()
v = Vector(0.0, Solution(mydomain))
if DIM==2:
    a0=1
    n0=1
    n1=0.5
    a1=-(a0*n0)/n1
    v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
    v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
else:
    a0=1
    a1=1
    n0=2
    n1=2
    n2=0.5
    a2=-(a0*n0+a1*n1)/n2
    v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
    v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
    v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])


H_t=Scalar(0.0, Solution(mydomain))
mts=Mountains(mydomain,v,eps=EPS,z=H)
dt=0.
while t<T_END:
    print "STEP ", t
    u=v*cos(OMEGA*t)
    u,Z=mts.update(u=u,H_t=H_t,verbose=False)
    
    saveVTK("state.%d.vtu"%n,sol=Z)
    print "Integral(Z)=",integrate(Z),Lsup(u[DIM-1])
    n+=1
    H_t=Z
    dt=mts.getDt()
    t+=dt

1	from esys.escript import *
2	from esys.escript.models import Mountains
3	from esys.finley import Brick,Rectangle
4	from math import pi, ceil
5
6	NE=16
7	DIM=3
8	H=1.
9	L=2*H
10	TOL=1.e-4
11	OMEGA=10
12	EPS=0.01
13	t=0
14	T_END=(2*pi)/OMEGA
15	n=0
16	if DIM==2:
17	mydomain=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=1, useFullElementOrder=True,optimize=True)
18	else:
19	mydomain=Brick(int(ceil(LNE/H)),int(ceil(LNE/H)),NE,l0=L,l1=L,l2=H,order=1, useFullElementOrder=True,optimize=True)
20	x=mydomain.getX()
21	v = Vector(0.0, Solution(mydomain))
22	if DIM==2:
23	a0=1
24	n0=1
25	n1=0.5
26	a1=-(a0*n0)/n1
27	v[0]=a0sin(pin0x[0]) cos(pin1x[1])
28	v[1]=a1cos(pin0x[0]) sin(pin1x[1])
29	else:
30	a0=1
31	a1=1
32	n0=2
33	n1=2
34	n2=0.5
35	a2=-(a0n0+a1n1)/n2
36	v[0]=a0sin(pin0x[0]) cos(pin1x[1])* cos(pin2x[2])
37	v[1]=a1cos(pin0x[0]) sin(pin1x[1])* cos(pin2x[2])
38	v[2]=a2cos(pin0x[0]) cos(pin1x[1])* sin(pin2x[2])
39
40
41	H_t=Scalar(0.0, Solution(mydomain))
42	mts=Mountains(mydomain,v,eps=EPS,z=H)
43	dt=0.
44	while t<T_END:
45	print "STEP ", t
46	u=vcos(OMEGAt)
47	u,Z=mts.update(u=u,H_t=H_t,verbose=False)
48
49	saveVTK("state.%d.vtu"%n,sol=Z)
50	print "Integral(Z)=",integrate(Z),Lsup(u[DIM-1])
51	n+=1
52	H_t=Z
53	dt=mts.getDt()
54	t+=dt
55