-
Notifications
You must be signed in to change notification settings - Fork 6
/
Nonlinear-solver.edp
159 lines (135 loc) · 5.45 KB
/
Nonlinear-solver.edp
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
/*
Author(s): Johann Moulin <[email protected]>
Pierre Jolivet <[email protected]>
Olivier Marquet <[email protected]>
Date: 2019-03-22
Copyright (C) 2019- Office national d'études et de recherches aérospatiales
2019- Centre National de la Recherche Scientifique
This script is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
If you use this script, you are kindly asked to cite the following article:
"Augmented Lagrangian Preconditioning for Large-Scale Hydrodynamic Stability",
J. Moulin, P. Jolivet, and O. Marquet (2019).
*/
load "PETSc"
include "include/params.idp"
real tolNewton = getARGV("-Newton_tol", 1.0e-6);
if(mpirank == 0) cout << "Building mesh..." << endl;
real time = mpiWtime();
mesh3 th = buildlayers(square(1, 1), 1);
fespace Wh(th, Pk); // complete space [u, v, w, p]
fespace Qh(th, P1); // pressure space for Schur complement
int[int][int] intersection; // local-to-neighbors renumbering
real[int] D; // partition of unity
Wh [uc1, uc2, uc3, pc];
include "include/decomposition.idp"
time = mpiWtime() - time;
if(mpirank == 0) cout << " ### Building mesh done in " << time << "s" << endl;
include "include/varf.idp"
IFMACRO(original)
varf vRes([u1, u2, u3, p], [v1, v2, v3, q]) = int3d(th)(
nu * (grad(uc1)' * grad(v1) +
grad(uc2)' * grad(v2) +
grad(uc3)' * grad(v3))
+ UgradV(uc, uc)' * [v1, v2, v3]
- (pc - gamma * div(uc)) * div(v)
- div(uc) * q)
ENDIFMACRO
IFMACRO(opt)
varf vRes([u1, u2, u3, p], [v1, v2, v3, q]) = int3d(th, qforder = 3)(
+ dx(uc1) * ((nu + gamma) * dx(v1) + gamma * (dy(v2) + dz(v3)) - q)
+ dy(uc1) * (nu * dy(v1) )
+ dz(uc1) * (nu * dz(v1) )
+ dx(uc2) * (nu * dx(v2) )
+ dy(uc2) * ((nu + gamma) * dy(v2) + gamma * (dx(v1) + dz(v3)) - q)
+ dz(uc2) * (nu * dz(v2) )
+ dx(uc3) * (nu * dx(v3) )
+ dy(uc3) * (nu * dy(v3) )
+ dz(uc3) * ((nu + gamma) * dz(v3) + gamma * (dx(v1) + dy(v2)) - q)
- pc * (dx(v1) + dy(v2) + dz(v3))
)
+ int3d(th)(
dx(uc1) * (uc1 * v1)
+ dy(uc1) * (uc2 * v1)
+ dz(uc1) * (uc3 * v1)
+ dx(uc2) * (uc1 * v2)
+ dy(uc2) * (uc2 * v2)
+ dz(uc2) * (uc3 * v2)
+ dx(uc3) * (uc1 * v3)
+ dy(uc3) * (uc2 * v3)
+ dz(uc3) * (uc3 * v3)
)
ENDIFMACRO
// Dirichlet (inlet + plate)
+ on(1, u1 = uc1-1, u2 = uc2-0, u3 = uc3-0) + on(2, u1 = uc1-0, u2 = uc2-0, u3 = uc3-0)
// far field (zmax and ymax)
+ on(6, 7, u1 = uc1-1);
verbosity = 1;
Mat dJ(Wh.ndof, intersection, D);
verbosity = 0;
include "include/fieldsplit.idp"
/*# Schur #*/
matrix[int] S(1); // array with a single matrix
varf vSchur(p, q) = int3d(th, qforder = 3)
(-1.0/(gamma + 1.0/Re) * p * q); // %*\color{DarkGreen}{\cref{eq:approximatedshurcomplement}}*) with %*\color{DarkGreen}{$s=0$}*)
S[0] = vSchur(Qh, Qh); // matrix assembly
/*# EndSchur #*/
if(mpirank == 0) cout << "PETSc..." << endl;
time = mpiWtime();
include "include/paramsXYZ.idp"
/*# P #*/
string paramsP = "-prefix_push fieldsplit_p_ " +
"-ksp_type cg -ksp_max_it 5 -pc_type jacobi -prefix_pop";
/*# EndP #*/
/*# Krylov #*/
string paramsKrylov = "-ksp_type fgmres -ksp_monitor"
+ " -ksp_rtol 1.0e-1 -ksp_gmres_restart 200";
/*# EndKrylov #*/
/*# AllParams #*/
string params = paramsXYZ + " " + paramsP + " " + paramsKrylov +
" -pc_type fieldsplit -pc_fieldsplit_type multiplicative";
/*# EndAllParams #*/
set(dJ, sparams = params, fields = vX[], names = names, schurPreconditioner = S, schurList = listX[]);
func real[int] funcRes(real[int]& inPETSc) {
changeNumbering(dJ, uc1[], inPETSc, inverse = true, exchange = true);
real[int] out(Wh.ndof);
out = vRes(0, Wh, tgv = -1);
real[int] outPETSc;
changeNumbering(dJ, out, outPETSc);
return outPETSc;
}
func int funcJ(real[int]& inPETSc) {
changeNumbering(dJ, uc1[], inPETSc, inverse = true, exchange = true);
matrix J = vJ(Wh, Wh, tgv = -1);
dJ = J;
return 0;
}
real[int] xPETSc;
changeNumbering(dJ, uc1[], xPETSc);
SNESSolve(dJ, funcJ, funcRes, xPETSc, sparams = "-snes_monitor -snes_linesearch_order 1 -snes_atol " + tolNewton);
changeNumbering(dJ, uc1[], xPETSc, inverse = true, exchange = true);
time = mpiWtime() - time;
if(mpirank == 0) cout << " ### PETSc done in " << time << "s" << endl;
if(Re < 75) {
savemesh(th, SaveDirName + "/mesh_" + nf + "_" + mpirank + "-" + mpisize + ".mesh");
{
ofstream file(SaveDirName + "/setup_" + nf + "_" + mpirank + "-" + mpisize + ".dat");
file << D << endl;
file << intersection.n << endl;
for(int i = 0; i < intersection.n; ++i) {
file << intersection[i].n << endl;
file << intersection[i] << endl;
}
}
ofstream file(SaveDirName + "/base_" + nf + "_" + mpirank + "-" + mpisize + ".dat");
file.precision(16);
file.scientific;
file << uc1[] << endl;
}
else {
ofstream file(SaveDirName + "/sol_" + nf + "_" + mpirank + "-" + mpisize + ".dat");
file.precision(16);
file.scientific;
file << uc1[] << endl;
}