-
Notifications
You must be signed in to change notification settings - Fork 1
/
respcompute.edp
118 lines (114 loc) · 4 KB
/
respcompute.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
//
// respcompute.edp
// Chris Douglas
//
// EXAMPLE USAGE:
// Solve linear forced response problem, store solution:
// ff-mpirun -np 4 respcompute.edp -eps_target 0+1.0i -forcing 1 -fi <FILEIN> -fo <VECFILE>
//
// Solve linear forced response problem over a sequence of shifts, store solutions:
// ff-mpirun -np 4 respcompute.edp -eps_target 0+0.2i -ntarget 10 -targetf 0+2.0i -forcing 1 -fi <FILEIN> -fo <VECFILE>
//
// NOTE: This file should not be changed unless you know what you're doing.
//
load "iovtk"
load "PETSc-complex"
include "settings.idp"
include "macros_bifbox.idp"
// arguments
string meshin = getARGV("-mi", ""); // input meshfile
string filein = getARGV("-fi", "");
string fileout = getARGV("-fo", "");
string statout = getARGV("-so", "");
string symstr = getARGV("-sym", "0");
real omega = getARGV("-omega", 1.0);
int nomega = getARGV("-nomega", 1);
real omegaf = getARGV("-omegaf", 1.0);
// Load mesh, make FE basis
string fileroot, fileext = parsefilename(filein, fileroot); //extract file name and extension
parsefilename(statout, statout); // trim extension from output file, if given
parsefilename(fileout, fileout); // trim extension from output file, if given
if(fileext == "mode" || fileext == "resp" || fileext == "rslv" || fileext == "tdls"){
filein = readbasename(workdir + filein);
fileext = parsefilename(filein, fileroot);
}
if(filein != "" && meshin == "") meshin = readmeshname(workdir + filein); // get mesh file
Th = readmeshN(workdir + meshin);
Thg = Th;
buildDmesh(Th);
restu = restrict(XMh, XMhg, n2o);
XMh defu(ub), defu(um2), defu(um3);
if(fileext == "base") {
ub[] = loadbase(fileroot, meshin);
}
else if(fileext == "fold") {
real[string] alpha;
real beta;
real[int] qm, qma;
ub[] = loadfold(fileroot, meshin, qm, qma, alpha, beta);
}
else if(fileext == "hopf") {
real omega;
complex[string] alpha;
complex beta;
complex[int] qm, qma;
ub[] = loadhopf(fileroot, meshin, qm, qma, sym, omega, alpha, beta);
}
else if(fileext == "foho") {
real omega;
complex[string] alpha1;
real[string] alpha2;
complex beta1, gamma12, gamma13;
real beta22, beta23, gamma22, gamma23;
complex[int] q1m, q1ma;
real[int] q2m, q2ma;
ub[] = loadfoho(fileroot, meshin, q1m, q1ma, q2m, q2ma, sym, omega, alpha1, alpha2, beta1, beta22, beta23, gamma12, gamma13, gamma22, gamma23);
}
else if(fileext == "hoho") {
real[int] sym1(sym.n), sym2(sym.n);
real omega1, omega2;
complex[string] alpha1, alpha2;
complex beta1, beta2, gamma1, gamma2, gamma12, gamma13, gamma22, gamma23;
complex[int] q1m, q1ma, q2m, q2ma;
ub[] = loadhoho(fileroot, meshin, q1m, q1ma, q2m, q2ma, sym1, sym2, omega1, omega2, alpha1, alpha2, beta1, beta2, gamma1, gamma2, gamma12, gamma13, gamma22, gamma23);
}
else if(fileext == "tdns") {
real time;
ub[] = loadtdns(fileroot, meshin, time);
}
else if(fileext == "porb") {
int Nh=1;
real omega;
complex[int, int] qh(ub[].n, Nh);
ub[] = loadporb(fileroot, meshin, qh, sym, omega, Nh);
}
XMh<complex> defu(um);
Mat<complex> J;
createMatu(Th, J, Pk);
sym = parsesymstr(symstr);
complex[int] ik(sym.n), ik2(sym.n), ik3(sym.n);
ik.im = sym;
complex iomega, iomega2 = 0.0, iomega3 = 0.0; // Let PETSc/SLEPc do the shift
include "eqns.idp"
real omegas = omega;
for (int n = 0; n < nomega; ++n){
if (nomega > 1) omega = omegas + (omegaf - omegas)*real(n)/real(nomega - 1);
iomega = 1i*omega;
complex[int] RHS = vMq(0, XMh, tgv = -1);
J = vJ(XMh, XMh, tgv = -1);
// set up Mat parameters
IFMACRO(Jprecon) Jprecon(iomega); ENDIFMACRO
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO sparams = KSPparams + " -options_left no");
um[] = J^-1*RHS;
real[int] ubsave = ub[];
complex[int] ubcplx(ub[].n), temp;
ubcplx.re = ub[];
ChangeNumbering(J, ubcplx, temp);
ChangeNumbering(J, ubcplx, temp, inverse = true);
ub[] = ubcplx.re;
ChangeNumbering(J, um[], temp);
ChangeNumbering(J, um[], temp, inverse = true);
saveresp(fileout + ((nomega > 1) ? ("_" + n) : ""), statout, filein, meshin, sym, omega, (fileout != ""), true);
ub[] = ubsave;
}