1 /*
2 * Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include "operator_ext_upml.h"
19 #include "FDTD/operator_cylindermultigrid.h"
20 #include "engine_ext_upml.h"
21 #include "tools/array_ops.h"
22 #include "fparser.hh"
23
24 using namespace std;
25
Operator_Ext_UPML(Operator * op)26 Operator_Ext_UPML::Operator_Ext_UPML(Operator* op) : Operator_Extension(op)
27 {
28 m_GradingFunction = new FunctionParser();
29 //default grading function
30 SetGradingFunction(" -log(1e-6)*log(2.5)/(2*dl*Z*(pow(2.5,W/dl)-1)) * pow(2.5, D/dl) ");
31
32 for (int n=0; n<6; ++n)
33 {
34 m_BC[n]=0;
35 m_Size[n]=0;
36 }
37 for (int n=0; n<3; ++n)
38 {
39 m_StartPos[n]=0;
40 m_numLines[n]=0;
41 }
42
43 vv = NULL;
44 vvfo = NULL;
45 vvfn = NULL;
46 ii = NULL;
47 iifo = NULL;
48 iifn = NULL;
49 }
50
~Operator_Ext_UPML()51 Operator_Ext_UPML::~Operator_Ext_UPML()
52 {
53 delete m_GradingFunction;
54 m_GradingFunction = NULL;
55 DeleteOp();
56 }
57
SetBoundaryCondition(const int * BCs,const unsigned int size[6])58 void Operator_Ext_UPML::SetBoundaryCondition(const int* BCs, const unsigned int size[6])
59 {
60 for (int n=0; n<6; ++n)
61 {
62 m_BC[n]=BCs[n];
63 m_Size[n]=size[n];
64 }
65 }
66
SetRange(const unsigned int start[3],const unsigned int stop[3])67 void Operator_Ext_UPML::SetRange(const unsigned int start[3], const unsigned int stop[3])
68 {
69 for (int n=0; n<3; ++n)
70 {
71 m_StartPos[n]=start[n];
72 m_numLines[n]=stop[n]-start[n]+1;
73 }
74 }
75
Create_UPML(Operator * op,const int ui_BC[6],const unsigned int ui_size[6],string gradFunc)76 bool Operator_Ext_UPML::Create_UPML(Operator* op, const int ui_BC[6], const unsigned int ui_size[6], string gradFunc)
77 {
78 int BC[6]={ui_BC[0],ui_BC[1],ui_BC[2],ui_BC[3],ui_BC[4],ui_BC[5]};
79 unsigned int size[6]={ui_size[0],ui_size[1],ui_size[2],ui_size[3],ui_size[4],ui_size[5]};
80
81 //check if mesh is large enough to support the pml
82 for (int n=0; n<3; ++n)
83 if ( (size[2*n]*(BC[2*n]==3)+size[2*n+1]*(BC[2*n+1]==3)) >= op->GetNumberOfLines(n,true) )
84 {
85 cerr << "Operator_Ext_UPML::Create_UPML: Warning: Not enough lines in direction: " << n << ", resetting to PEC" << endl;
86 BC[2*n]=0;
87 size[2*n]=0;
88 BC[2*n+1]=0;
89 size[2*n+1]=0;
90 }
91
92 //check cylindrical coord compatiblility
93 Operator_Cylinder* op_cyl = dynamic_cast<Operator_Cylinder*>(op);
94 if (op_cyl)
95 {
96 if ((BC[0]==3) && (op_cyl->GetClosedAlpha() || op_cyl->GetR0Included()))
97 {
98 BC[0]=0;
99 size[0]=0;
100 cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in r-min direction is not possible, resetting to PEC..." << endl;
101 }
102 if ( (BC[2]==3) && (op_cyl->GetClosedAlpha()) )
103 {
104 BC[2]=0;
105 size[2]=0;
106 cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in alpha-min direction is not possible, resetting to PEC..." << endl;
107 }
108 if ( (BC[3]==3) && (op_cyl->GetClosedAlpha()) )
109 {
110 BC[3]=0;
111 size[3]=0;
112 cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in alpha-max direction is not possible, resetting to PEC..." << endl;
113 }
114 }
115
116 //check cylindrical coord compatiblility
117 if (dynamic_cast<Operator_CylinderMultiGrid*>(op))
118 {
119 if (BC[2]==3)
120 {
121 BC[2]=0;
122 size[2]=0;
123 cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in alpha direction is not possible for a cylindrical multi-grid, resetting to PEC..." << endl;
124 }
125 if (BC[3]==3)
126 {
127 BC[3]=0;
128 size[3]=0;
129 cerr << "Operator_Ext_UPML::Create_UPML: Warning: An upml in alpha direction is not possible for a cylindrical multi-grid, resetting to PEC..." << endl;
130 }
131 }
132
133
134 Operator_Ext_UPML* op_ext_upml=NULL;
135 unsigned int start[3]={0 ,0 ,0};
136 unsigned int stop[3] ={op->GetNumberOfLines(0,true)-1,op->GetNumberOfLines(1,true)-1,op->GetNumberOfLines(2,true)-1};
137
138 //create a pml in x-direction over the full width of yz-space
139 if (BC[0]==3)
140 {
141 op_ext_upml = new Operator_Ext_UPML(op);
142 op_ext_upml->SetGradingFunction(gradFunc);
143 start[0]=0;
144 stop[0] =size[0];
145 op_ext_upml->SetBoundaryCondition(BC, size);
146 op_ext_upml->SetRange(start,stop);
147 op->AddExtension(op_ext_upml);
148 }
149 if (BC[1]==3)
150 {
151 op_ext_upml = new Operator_Ext_UPML(op);
152 op_ext_upml->SetGradingFunction(gradFunc);
153 start[0]=op->GetNumberOfLines(0,true)-1-size[1];
154 stop[0] =op->GetNumberOfLines(0,true)-1;
155 op_ext_upml->SetBoundaryCondition(BC, size);
156 op_ext_upml->SetRange(start,stop);
157 op->AddExtension(op_ext_upml);
158 }
159
160 //create a pml in y-direction over the xz-space (if a pml in x-direction already exists, skip that corner regions)
161 start[0]=(size[0]+1)*(BC[0]==3);
162 stop[0] =op->GetNumberOfLines(0,true)-1-(size[0]+1)*(BC[1]==3);
163
164 if (BC[2]==3)
165 {
166 op_ext_upml = new Operator_Ext_UPML(op);
167 op_ext_upml->SetGradingFunction(gradFunc);
168 start[1]=0;
169 stop[1] =size[2];
170 op_ext_upml->SetBoundaryCondition(BC, size);
171 op_ext_upml->SetRange(start,stop);
172 op->AddExtension(op_ext_upml);
173 }
174 if (BC[3]==3)
175 {
176 op_ext_upml = new Operator_Ext_UPML(op);
177 op_ext_upml->SetGradingFunction(gradFunc);
178 start[1]=op->GetNumberOfLines(1,true)-1-size[3];
179 stop[1] =op->GetNumberOfLines(1,true)-1;
180 op_ext_upml->SetBoundaryCondition(BC, size);
181 op_ext_upml->SetRange(start,stop);
182 op->AddExtension(op_ext_upml);
183 }
184
185 //create a pml in z-direction over the xy-space (if a pml in x- and/or y-direction already exists, skip that corner/edge regions)
186 start[1]=(size[2]+1)*(BC[2]==3);
187 stop[1] =op->GetNumberOfLines(1,true)-1-(size[3]+1)*(BC[3]==3);
188
189 //exclude x-lines that does not belong to the base multi-grid operator;
190 Operator_CylinderMultiGrid* op_cyl_MG = dynamic_cast<Operator_CylinderMultiGrid*>(op);
191 if (op_cyl_MG)
192 start[0] = op_cyl_MG->GetSplitPos()-1;
193
194 if (BC[4]==3)
195 {
196 op_ext_upml = new Operator_Ext_UPML(op);
197 op_ext_upml->SetGradingFunction(gradFunc);
198 start[2]=0;
199 stop[2] =size[4];
200 op_ext_upml->SetBoundaryCondition(BC, size);
201 op_ext_upml->SetRange(start,stop);
202 op->AddExtension(op_ext_upml);
203 }
204 if (BC[5]==3)
205 {
206 op_ext_upml = new Operator_Ext_UPML(op);
207 op_ext_upml->SetGradingFunction(gradFunc);
208 start[2]=op->GetNumberOfLines(2,true)-1-size[5];
209 stop[2] =op->GetNumberOfLines(2,true)-1;
210 op_ext_upml->SetBoundaryCondition(BC, size);
211 op_ext_upml->SetRange(start,stop);
212 op->AddExtension(op_ext_upml);
213 }
214
215 BC[1]=0;
216 size[1]=0;
217 //create pml extensions (in z-direction only) for child operators in cylindrical multigrid operators
218 while (op_cyl_MG)
219 {
220 Operator_Cylinder* op_child = op_cyl_MG->GetInnerOperator();
221 op_cyl_MG = dynamic_cast<Operator_CylinderMultiGrid*>(op_child);
222 for (int n=0; n<2; ++n)
223 {
224 start[n]=0;
225 stop[n]=op_child->GetNumberOfLines(n,true)-1;
226 }
227
228 if (op_cyl_MG)
229 start[0] = op_cyl_MG->GetSplitPos()-1;
230
231 if (BC[4]==3)
232 {
233 op_ext_upml = new Operator_Ext_UPML(op_child);
234 op_ext_upml->SetGradingFunction(gradFunc);
235 start[2]=0;
236 stop[2] =size[4];
237 op_ext_upml->SetBoundaryCondition(BC, size);
238 op_ext_upml->SetRange(start,stop);
239 op_child->AddExtension(op_ext_upml);
240 }
241 if (BC[5]==3)
242 {
243 op_ext_upml = new Operator_Ext_UPML(op_child);
244 op_ext_upml->SetGradingFunction(gradFunc);
245 start[2]=op->GetNumberOfLines(2,true)-1-size[5];
246 stop[2] =op->GetNumberOfLines(2,true)-1;
247 op_ext_upml->SetBoundaryCondition(BC, size);
248 op_ext_upml->SetRange(start,stop);
249 op_child->AddExtension(op_ext_upml);
250 }
251 }
252
253 return true;
254 }
255
256
DeleteOp()257 void Operator_Ext_UPML::DeleteOp()
258 {
259 Delete_N_3DArray<FDTD_FLOAT>(vv,m_numLines);
260 vv = NULL;
261 Delete_N_3DArray<FDTD_FLOAT>(vvfo,m_numLines);
262 vvfo = NULL;
263 Delete_N_3DArray<FDTD_FLOAT>(vvfn,m_numLines);
264 vvfn = NULL;
265 Delete_N_3DArray<FDTD_FLOAT>(ii,m_numLines);
266 ii = NULL;
267 Delete_N_3DArray<FDTD_FLOAT>(iifo,m_numLines);
268 iifo = NULL;
269 Delete_N_3DArray<FDTD_FLOAT>(iifn,m_numLines);
270 iifn = NULL;
271 }
272
273
SetGradingFunction(string func)274 bool Operator_Ext_UPML::SetGradingFunction(string func)
275 {
276 if (func.empty())
277 return true;
278
279 m_GradFunc = func;
280 int res = m_GradingFunction->Parse(m_GradFunc.c_str(), "D,dl,W,Z,N");
281 if (res < 0) return true;
282
283 cerr << "Operator_Ext_UPML::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl;
284 cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n";
285 return false;
286 }
287
CalcGradingKappa(int ny,unsigned int pos[3],double Zm,double kappa_v[3],double kappa_i[3])288 void Operator_Ext_UPML::CalcGradingKappa(int ny, unsigned int pos[3], double Zm, double kappa_v[3], double kappa_i[3])
289 {
290 double depth=0;
291 double width=0;
292 for (int n=0; n<3; ++n)
293 {
294 if ((pos[n] <= m_Size[2*n]) && (m_BC[2*n]==3)) //lower pml in n-dir
295 {
296 width = (m_Op->GetDiscLine(n,m_Size[2*n]) - m_Op->GetDiscLine(n,0))*m_Op->GetGridDelta();
297 depth = width - (m_Op->GetDiscLine(n,pos[n]) - m_Op->GetDiscLine(n,0))*m_Op->GetGridDelta();
298
299 if ((m_Op_Cyl) && (n==1))
300 {
301 width *= m_Op_Cyl->GetDiscLine(0,pos[0]);
302 depth *= m_Op_Cyl->GetDiscLine(0,pos[0]);
303 }
304
305 if (n==ny)
306 depth-=m_Op->GetEdgeLength(n,pos)/2;
307 double vars[5] = {depth, width/m_Size[2*n], width, Zm, (double)m_Size[2*n]};
308 if (depth>0)
309 kappa_v[n] = m_GradingFunction->Eval(vars);
310 else
311 kappa_v[n]=0;
312 if (n==ny)
313 depth+=m_Op->GetEdgeLength(n,pos)/2;
314
315 if (n!=ny)
316 depth-=m_Op->GetEdgeLength(n,pos)/2;
317 if (depth<0)
318 depth=0;
319 vars[0]=depth;
320 if (depth>0)
321 kappa_i[n] = m_GradingFunction->Eval(vars);
322 else
323 kappa_i[n] = 0;
324 }
325 else if ((pos[n] >= m_Op->GetNumberOfLines(n,true) -1 -m_Size[2*n+1]) && (m_BC[2*n+1]==3)) //upper pml in n-dir
326 {
327 width = (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-1) - m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-m_Size[2*n+1]-1))*m_Op->GetGridDelta();
328 depth = width - (m_Op->GetDiscLine(n,m_Op->GetNumberOfLines(n,true)-1) - m_Op->GetDiscLine(n,pos[n]))*m_Op->GetGridDelta();
329
330 if ((m_Op_Cyl) && (n==1))
331 {
332 width *= m_Op_Cyl->GetDiscLine(0,pos[0]);
333 depth *= m_Op_Cyl->GetDiscLine(0,pos[0]);
334 }
335
336 if (n==ny)
337 depth+=m_Op->GetEdgeLength(n,pos)/2;
338 double vars[5] = {depth, width/(m_Size[2*n]), width, Zm, (double)m_Size[2*n]};
339 if (depth>0)
340 kappa_v[n] = m_GradingFunction->Eval(vars);
341 else
342 kappa_v[n]=0;
343 if (n==ny)
344 depth-=m_Op->GetEdgeLength(n,pos)/2;
345
346 if (n!=ny)
347 depth+=m_Op->GetEdgeLength(n,pos)/2;
348 if (depth>width)
349 depth=0;
350 vars[0]=depth;
351 if (depth>0)
352 kappa_i[n] = m_GradingFunction->Eval(vars);
353 else
354 kappa_i[n]=0;
355 }
356 else
357 {
358 kappa_v[n] = 0;
359 kappa_i[n] = 0;
360 }
361 }
362 }
363
BuildExtension()364 bool Operator_Ext_UPML::BuildExtension()
365 {
366 /*Calculate the upml coefficients as defined in:
367 Allen Taflove, computational electrodynamics - the FDTD method, third edition, chapter 7.8, pages 297-300
368 - modified by Thorsten Liebig to match the equivalent circuit (EC) FDTD method
369 - kappa is used for conductivities (instead of sigma)
370 */
371 if (m_Op==NULL)
372 return false;
373
374 DeleteOp();
375 vv = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
376 vvfo = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
377 vvfn = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
378 ii = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
379 iifo = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
380 iifn = Create_N_3DArray<FDTD_FLOAT>(m_numLines);
381
382 unsigned int pos[3];
383 unsigned int loc_pos[3];
384 int nP,nPP;
385 double kappa_v[3]={0,0,0};
386 double kappa_i[3]={0,0,0};
387 double eff_Mat[4];
388 double dT = m_Op->GetTimestep();
389
390 for (loc_pos[0]=0; loc_pos[0]<m_numLines[0]; ++loc_pos[0])
391 {
392 pos[0] = loc_pos[0] + m_StartPos[0];
393 for (loc_pos[1]=0; loc_pos[1]<m_numLines[1]; ++loc_pos[1])
394 {
395 pos[1] = loc_pos[1] + m_StartPos[1];
396 vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL);
397 for (loc_pos[2]=0; loc_pos[2]<m_numLines[2]; ++loc_pos[2])
398 {
399 pos[2] = loc_pos[2] + m_StartPos[2];
400 for (int n=0; n<3; ++n)
401 {
402 m_Op->Calc_EffMatPos(n,pos,eff_Mat,vPrims);
403 CalcGradingKappa(n, pos,__Z0__ ,kappa_v ,kappa_i);
404 nP = (n+1)%3;
405 nPP = (n+2)%3;
406 if ((kappa_v[0]+kappa_v[1]+kappa_v[2])!=0)
407 {
408 //check if pos is on PEC
409 if ( (m_Op->GetVV(n,pos[0],pos[1],pos[2]) + m_Op->GetVI(n,pos[0],pos[1],pos[2])) != 0 )
410 {
411 //modify the original operator to perform eq. (7.85) by the main engine (EC-FDTD: equation is multiplied by delta_n)
412 //the engine extension will replace the original voltages with the "voltage flux" (volt*eps0) prior to the voltage updates
413 //after the updates are done the extension will calculate the new voltages (see below) and place them back into the main field domain
414 m_Op->SetVV(n,pos[0],pos[1],pos[2], (2*__EPS0__ - kappa_v[nP]*dT) / (2*__EPS0__ + kappa_v[nP]*dT) );
415 m_Op->SetVI(n,pos[0],pos[1],pos[2], (2*__EPS0__*dT) / (2*__EPS0__ + kappa_v[nP]*dT) * m_Op->GetEdgeLength(n,pos) / m_Op->GetEdgeArea(n,pos) );
416
417
418 //operators needed by eq. (7.88) to calculate new voltages from old voltages and old and new "voltage fluxes"
419 GetVV(n,loc_pos) = (2*__EPS0__ - kappa_v[nPP]*dT) / (2*__EPS0__ + kappa_v[nPP]*dT);
420 GetVVFN(n,loc_pos) = (2*__EPS0__ + kappa_v[n]*dT) / (2*__EPS0__ + kappa_v[nPP]*dT)/eff_Mat[0];
421 GetVVFO(n,loc_pos) = (2*__EPS0__ - kappa_v[n]*dT) / (2*__EPS0__ + kappa_v[nPP]*dT)/eff_Mat[0];
422 }
423 }
424 else
425 {
426 //disable upml
427 GetVV(n,loc_pos) = m_Op->GetVV(n,pos[0],pos[1],pos[2]);
428 m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 );
429 GetVVFO(n,loc_pos) = 0;
430 GetVVFN(n,loc_pos) = 1;
431 }
432
433 if ((kappa_i[0]+kappa_i[1]+kappa_i[2])!=0)
434 {
435 //check if pos is on PMC
436 if ( (m_Op->GetII(n,pos[0],pos[1],pos[2]) + m_Op->GetIV(n,pos[0],pos[1],pos[2])) != 0 )
437 {
438 //modify the original operator to perform eq. (7.89) by the main engine (EC-FDTD: equation is multiplied by delta_n)
439 //the engine extension will replace the original currents with the "current flux" (curr*mu0) prior to the current updates
440 //after the updates are done the extension will calculate the new currents (see below) and place them back into the main field domain
441 m_Op->SetII(n,pos[0],pos[1],pos[2], (2*__EPS0__ - kappa_i[nP]*dT) / (2*__EPS0__ + kappa_i[nP]*dT) );
442 m_Op->SetIV(n,pos[0],pos[1],pos[2], (2*__EPS0__*dT) / (2*__EPS0__ + kappa_i[nP]*dT) * m_Op->GetEdgeLength(n,pos,true) / m_Op->GetEdgeArea(n,pos,true) );
443
444 //operators needed by eq. (7.90) to calculate new currents from old currents and old and new "current fluxes"
445 GetII(n,loc_pos) = (2*__EPS0__ - kappa_i[nPP]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT);
446 GetIIFN(n,loc_pos) = (2*__EPS0__ + kappa_i[n]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT)/eff_Mat[2];
447 GetIIFO(n,loc_pos) = (2*__EPS0__ - kappa_i[n]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT)/eff_Mat[2];
448 }
449 }
450 else
451 {
452 //disable upml
453 GetII(n,loc_pos) = m_Op->GetII(n,pos[0],pos[1],pos[2]);
454 m_Op->SetII(n,pos[0],pos[1],pos[2], 0 );
455 GetIIFO(n,loc_pos) = 0;
456 GetIIFN(n,loc_pos) = 1;
457 }
458 }
459 }
460 }
461 }
462 return true;
463 }
464
CreateEngineExtention()465 Engine_Extension* Operator_Ext_UPML::CreateEngineExtention()
466 {
467 Engine_Ext_UPML* eng_ext = new Engine_Ext_UPML(this);
468 return eng_ext;
469 }
470
ShowStat(ostream & ostr) const471 void Operator_Ext_UPML::ShowStat(ostream &ostr) const
472 {
473 Operator_Extension::ShowStat(ostr);
474
475 ostr << " PML range\t\t: " << "[" << m_StartPos[0]<< "," << m_StartPos[1]<< "," << m_StartPos[2]<< "] to ["
476 << m_StartPos[0]+m_numLines[0]-1 << "," << m_StartPos[1]+m_numLines[1]-1 << "," << m_StartPos[2]+m_numLines[2]-1 << "]" << endl;
477 ostr << " Grading function\t: \"" << m_GradFunc << "\"" << endl;
478 }
479