/
Helmholtz-2d-FEM-BEM-coupling-PETSc-composite.edp
141 lines (101 loc) · 4.28 KB
/
Helmholtz-2d-FEM-BEM-coupling-PETSc-composite.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
//ff-mpirun -np 4 Helmholtz-2d-FEM-BEM-coupling-PETSc-composite.edp -wg
/* example of wave guiding with gradient-index lenses */
load "bem"
load "msh3"
load "PETSc-complex"
complex k = 10; // wavenumber
real lambda = 2*pi/real(k);
real nloc = 10./lambda;
real n = nloc*2*pi;
int[int] nsl(10); // number of lenses
nsl = n;
real theta = 10*2*pi/360; // angular shift between lenses
real[int] nsx(nsl.n), nsy(nsl.n);
nsx[0] = 3.2;
nsy[0] = 0;
for (int i=1; i< nsl.n; i++) {
nsx[i] = nsx[i-1] + 2.01*cos(i*theta);
nsy[i] = nsy[i-1] + 2.01*sin(i*theta);
}
int interface = 1;
int waveguide = 2;
border circle(t=0, 2*pi; i){x=cos(t)+nsx[i]; y=sin(t)+nsy[i]; z=0; label=interface;}
real L = 2;
real l = 0.3;
real dd = 0.1;
func finc = exp(-100*((x+L-0.1)^2+y^2+z^2)); // source for waveguide excitation
// interface around the waveguide :
border a(t=-L-dd, L+dd){x=t; y=-l-dd; z=0; label=interface;}
border b(t=-l-dd, l+dd){x=L+dd; y=t; z=0; label=interface;}
border c(t=L+dd, -L-dd){x=t; y=l+dd; z=0; label=interface;}
border d(t=l+dd, -l-dd){x=-L-dd; y=t; z=0; label=interface;}
// waveguide :
border ga(t=-L, L){x=t; y=-l; z=0; label=waveguide;}
border gc(t=L, -L){x=t; y=l; z=0; label=waveguide;}
border gd(t=l, -l){x=-L; y=t; z=0; label=waveguide;}
// Fem mesh :
mesh Th = buildmesh(a(nloc*(2*L+2*dd))+b(nloc*(2*l+2*dd))+c(nloc*(2*L+2*dd))+d(nloc*(2*l+2*dd))
+ga(nloc*2*L)+gc(nloc*2*L)+gd(nloc*2*l)+circle(nsl));
plot(Th, wait=1, dim=2, cmm="FEM mesh Th");
int[int] lab = [interface];
meshL ThL = extract(Th, label=lab); // BEM mesh
ThL = OrientNormal(ThL,unbounded=1); // BEM mesh
plot(ThL, wait=1, dim=2, cmm="BEM mesh ThL");
// change region labels in the fem mesh to define gradient index in lenses
fespace Ph(Th,P0);
Ph reg = region;
int[int] regs(2*nsl.n+2);
for (int i=0; i< nsl.n; i++){
int regi = reg(nsx[i], nsy[i]); // the lenses
regs[2*i] = regi;
regs[2*i+1] = i+1;
}
regs[2*nsl.n] = reg(0,0); // the waveguide
regs[2*nsl.n+1] = -1;
Th = change(Th,region=regs);
reg = region;
func ind = reg == -1 ? 1 : 2./(1+((x-nsx[max(0.,reg-1)])^2+(y-nsy[max(0.,reg-1)])^2)); // gradient index in lenses
fespace Uh(Th,P1);
fespace UhL(ThL,P1);
fespace Ch=Uh*UhL;
macro Grad(u) [dx(u),dy(u)] // EOM
/* the coupled problem we want to solve is :
( F TDL ) (ufem) = (Frhs)
( mass -SL ) (ubem) ( 0 ) */
Uh<complex> ufem,v1;
UhL<complex> ubem,v2;
// problem formulation
varf vfLenses(<[ufem],[ubem]>,<[v1],[v2]>) =
int2d(Th)(-ind*k^2*ufem*v1 + Grad(ufem)'*Grad(v1)) // F
+ int1dx1d(ThL)(ThL)(BEM(BemKernel("TDL",k=k),ubem,v1)) + int1d(ThL)(0.5*ubem*v1) // TDL
+ int1d(ThL)(ufem*v2) // mass
+ int1dx1d(ThL)(ThL)(BEM(-1*BemKernel("SL",k=k),ubem,v2)) // -SL
+ int2d(Th)(finc*v1) + on(waveguide,ufem=0) ; // RHS
// the linear system is assembled in parallel ; the distributed matrix is then passed to PETSc as a nested Mat (MatNest)
Mat<complex> HC = vfLenses(Ch,Ch);
complex[int] rhs = vfLenses(0,Ch);
// example of fieldsplit preconditioner ; the fields are automatically deduced from the composite FE space Ch
set(HC,sparams="-ksp_view -ksp_monitor -ksp_type fgmres -ksp_view_final_residual -ksp_gmres_restart 200 -pc_type fieldsplit -fieldsplit_0_pc_type lu -fieldsplit_0_pc_mat_solver_type mumps "+"-fieldsplit_1_ksp_type gmres -fieldsplit_1_ksp_max_it 20");
complex[int] u = HC^-1*rhs;
[ufem[],ubem[]] = u; // dispatch solution
plot(ufem, fill=1, value=1, wait=1, dim=2, cmm="FEM solution");
plot(ubem, fill=1, value=1, wait=1, dim=2, cmm="BEM ansatz on ThL");
// output mesh for visualization
int np = 200/2;
real R = 20;
real rr = 20;
border b1(t=-rr, R){x=t; y=-rr;}
border b2(t=-rr, rr){x=R; y=t;}
border b3(t=R, -rr){x=t; y=rr;}
border b4(t=rr, -rr){x=-rr; y=t;}
nsl = -nsl;
// exterior mesh
mesh ThOut = buildmesh(b1(np*R/rr)+b2(np)+b3(np*R/rr)+b4(np)+circle(nsl)
+a(-nloc*(2*L+2*dd))+b(-nloc*(2*l+2*dd))+c(-nloc*(2*L+2*dd))+d(-nloc*(2*l+2*dd)));
fespace UhOut(ThOut,P1);
varf vp(u,v)=int1d(ThL)(POT(BemPotential("SL",k=k),u,v));
HMatrix<complex> B = vp(UhL,UhOut);
if (mpirank == 0) cout << B.infos << endl;
UhOut<complex> uext;
uext[] = B*ubem[];
plot(ufem, uext, dim=2, fill=1, value=1, nbiso=40, cmm="u");