1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
load "BernardiRaugel"
real s0=clock();
// Parameters
real reylnods = 400.;
func g = (x)*(1-x)*4;
real nu = 1;
real dt = 0.1;
real penalty = 1.e-6;
int nbiter = 3;
real coefdt = 0.25^(1./nbiter);
real coefcut = 0.25^(1./nbiter);
real cut = 0.01;
real tol = 0.3;
real coeftol = 0.25^(1./nbiter);
// Mesh
mesh Th = square(2, 2);
// Fespace
fespace Vh2(Th, P2BR);
fespace Vh(Th, P0);
fespace Wh(Th, [P2BR, P0]);
Wh [u1, u2, p], [v1, v2, q], [up1, up2, pp];
// Macro
macro div(u1, u2) (dx(u1) + dy(u2)) //
// Problem
real alpha = 0;
int i = 0;
varf NS ([u1, u2, p], [v1, v2, q], init=i)
= int2d(Th)(
alpha*(u1*v1 + u2*v2)
+ nu*(
dx(u1)*dx(v1)
+ dy(u1)*dy(v1)
+ dx(u2)*dx(v2)
+ dy(u2)*dy(v2)
)
- penalty * p*q
- p*div(v1, v2)
- div(u1, u2)*q
)
- int2d(Th)(
- alpha*convect([up1, up2], -dt, up1)*v1
- alpha*convect([up1, up2], -dt, up2)*v2
)
+ on(3, u1=g, u2=0)
+ on(1, 2, 4, u1=0, u2=0)
;
// Initialization
[up1, up2, pp] = [0., 0., 0.];
// Solve (Stokes)
matrix A = NS(Wh, Wh, solver=UMFPACK);
real[int] b = NS(0, Wh);
u1[] = A^-1 * b;
// Plot
plot([u1,u2], p, coef=0.2, cmm="[u1, u2] and p", value=true, wait=true);
// Convergence loop
int iter = 0;
nu = 1./reylnods;
for (iter = 1; iter <= nbiter; iter++){
cout << "dt = " << dt << endl;
alpha = 1./dt;
A = NS(Wh, Wh, solver=UMFPACK);
// Time loop
for (i = 0; i <= 10; i++){
// Update
up1[] = u1[];
real[int] b = NS(0, Wh);
// Solve
u1[] = A^-1 * b;
// Plot
if (!(i % 10))
plot([u1,u2], p, coef=0.2, cmm="[u1,u2] and p");
// Display
cout << "CPU " << clock()-s0 << "s " << endl;
}
if (iter >= nbiter) break;
// Mesh adaptation
Th = adaptmesh(Th, [u1, u2], iso=0, abserror=0, cutoff=cut, err=tol, inquire=0, ratio=1.5, hmin=1./1000);
plot(Th);
// Update
dt = dt*coefdt;
tol = tol *coeftol;
cut = cut *coefcut;
// Re-interpolation
[u1, u2, p] = [u1, u2, p];
[up1, up2, pp] = [u1, u2, p];
}
// Display
cout << "CPU " << clock()-s0 << "s " << endl;
|
