Project 291 - John Benjamin

Two dimensional simulation of electric field magnitude for dielectric device

In this project the propegation behaviour of electric fields were investigated in a system which has a dielectric and absorbing boundary. The investigated system was an infinite rectangular rod, using such a construct, the dimensionality of the governing equations (Maxwell's equations) can be reduced. The infinite rod is shown in the figure below:

1. The procedure and approach to implementation of Maxwell's equation

The approach in the finite difference time domain (FDTD) is determined by the style of materials used. Since perfectly matched layers (PML) is included in the model to absorb and remove energy from the system, the Maxwell equations has to be considered in the frequency domain. PML is strictly a frequency domain property, and has to be moved to the time domain using inverse Laplace transforms. Using inverse Laplace transforms the frequency domain is then moved to the time domain in the following manner:

$$j \omega [\,\,\,] = \frac{\partial}{\partial t}[\,\,\,]$$

and,

$$\frac{1}{j \omega} [\,\,\,] = \int_{0}^{t} [\,\,\,]dt'$$

The PML functions as an absorbing boundary placed at the edge of the problem space to simulate waves leaving the system. This is done by matching electromagnetic impedance with the problem space while gradually increasing the conductivity as the wave propegate deeper into the PML. The interal structure is shown below:

The dielectric has the following properties: Diagonally anisotropic permittivity $$[ \epsilon_{r} ]$$, anisotropic permeability $$[ \mu_{r} ]$$, and anisotropic conductivity $$[ \sigma ]$$. The brackets represent that the quantities are forming a matrix when interacting with the field vectors of Maxwells equation.

1.1 Maxwell's equation from 3D to 2D

The Maxwell's curl equations (and constitutive relations) in frequency space with PML matrix $$[ S ]$$ is given by:

$$\nabla \times E(\omega) = -j \omega \mu_{0} [\mu_{r}][S] H(\omega)$$ $$\nabla \times H(\omega) = [\sigma]E(\omega) + j \omega [S] D(\omega)$$ $$D(\omega) = \epsilon_{0} [\epsilon_{r}] E(\omega)$$

where the Electric field is given by E, the magnetic field by H, the displacement field by D, the vacuum permittivity $$\epsilon_{0}$$, and the vacuum permeability $$\mu_{0}$$. Before reducing the dimensions, the electric field is normalized: $$\tilde{E}=\sqrt{ \frac{\epsilon_{0}}{\mu_{0}} } E=\frac{1}{\eta_{0}} E,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \tilde{D}=\sqrt{ \frac{1}{ \epsilon_{0} \mu_{0}} }=c_{0}D$$ Using a Cartesian coordinate system we note that the rod is extended infinitely far in the z-direction. This implies that all the functions are constant when displaced in the z-direction, i.e. $$\partial / \partial z = 0$$. Expanding the Maxwell's equation in terms of its cartesian coordinate system using the vectors x,y,z:

$$C_{x}^{E}=\frac{\partial \tilde{E}_{z}}{\partial y}=-j\omega \big( \frac{\mu_{rx}}{c_{0}} S_{y} H_{x} \big)$$ $$C_{y}^{E}=-\frac{\partial \tilde{E}_{z}}{\partial x}=-j\omega \big( \frac{\mu_{ry}}{c_{0}} S_{y}^{-1} H_{y} \big)$$ $$C_{z}^{H}=\frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial x}=(\eta_{0} \sigma_{z})\tilde{E}_{z} + j\omega \big( \frac{1}{c_{0}} S_{y} \tilde{D}_{z} \big)$$

and,

$$C_{x}^{H}=\frac{\partial H_{z}}{\partial y}=(\eta_{0}\sigma_{x}) \tilde{E}_{x} + j\omega \big( \frac{1}{c_{0}} S_{y} \tilde{D}_{x} \big)$$ $$C_{y}^{H}=-\frac{\partial H_{z}}{\partial x}=(\eta_{0}\sigma_{y}) \tilde{E}_{y} + j\omega \big( \frac{1}{c_{0}} S_{y}^{-1} \tilde{D}_{y} \big)$$ $$C_{z}^{E}=\frac{\partial \tilde{E}_{y}}{\partial x}-\frac{\partial \tilde{E}_{x}}{\partial x}=-j\omega \big( \frac{\mu_{rz}}{c_{0}} S_{y} H_{z} \big)$$

for the PML oriented such that it absorbes radiation in the y-direction:

$$[ S ] =\begin{bmatrix} S_{y} & 0 & 0\\0 & S_{y}^{-1} & 0\\0 &0 &S_{y}\end{bmatrix},\,\,\,\,\,\,\,\, S_{y}=1 + \frac{1}{j \omega} \frac{\sigma'(y)}{\epsilon_{0}},\,\,\,\,\,\,\,\, \sigma'(y)= \frac{\epsilon_{0}}{2 \Delta t} \big( \frac{y'}{L_{y}} \big)^{3}$$

where $$L_{y}$$ is the length of the PML, and $$y'$$ is the current depth within the PML. Outside of the PML, $$[S] =1$$. The choice of the max value of $$[S]$$ is usually the one shown above, but will depend on the exact parameters (source intensity, geometry of the chamber, material distribution, etc...), and is usually tweaked untill it produces a stable absorption. Selecting the wrong value of the max value can make the system numerically unstable. From the expanded curl equations it is noticed system becomes seperated into two independent sub-system of differential equations. To investigate the wave properties one of these modes can then be selected (given that the source is polarized). We select the mode corresponding to $$\{E_{x} (x,y),E_{y}(x,y),H_{z}(x,y)\}$$. Moving the above equations to the time domain using the inverse Laplace transform:

$$C_{z}^{E}|_{t}=-\big( \frac{\mu_{ry}}{c_{0}} \big)\frac{\partial H_{z}}{\partial t}|_{t} - \big( \frac{\mu_{ry}\sigma'}{c_{0}\epsilon_{0}} \big) H_{z}|_{t}$$ $$C_{x}^{H}|_{t}=(\eta_{0}\sigma_{x})\tilde{E}_{x}|_{t} + \big( \frac{1}{c_{0}} \big)\frac{\partial \tilde{D}_{x}}{\partial t}|_{t} + \big( \frac{\sigma'}{c_{0}\epsilon_{0}} \big)\tilde{D}_{x}|_{t}$$ $$C_{y}^{H}|_{t} + \big( \frac{\sigma'}{\epsilon_{0}} \big) \int_{0}^{t} C_{y}^{H} dt = (\eta_{0}\sigma_{y})\tilde{E}_{y}|_{t} + \big( \frac{\eta_{0}\sigma_{y}\sigma'}{\epsilon_0} \big) \int_{0}^{t}\tilde{E}_{y}|_{t} dt + \big( \frac{1}{c_{0}} \big) \frac{\partial \tilde{D}_{y}}{\partial t}|_{t}$$

1.2a From a continuous model to a discrete model: Temporal discretization

The electric field and the magnetic field is temporally displaced relative to eachother. Using the electric field at $$t_{n}$$, the magnetic field defined at the next half step in time can then be calculated $$t_{n}+\frac{\Delta t}{2}$$. Using this magnetic field the electric field can be at the time step $$t_{n}+\Delta t =t_{n+1}$$. This is how the temporal displacements work and how the fields evolve with time when in discrete form. Using finite differences the differentials and the integrals can be substituted with differences and sums, and from the differentials the next time-steps can be evaluated. The update equations are then given by:

$$H_{z}|_{t+\frac{\Delta t}{2}}=-h_{z1}C_{z}^{E}|_{t} + h_{z2}H_{z}|_{t-\frac{\Delta t}{2}}$$ $$\tilde{D}_{x}|_{t+\Delta t}= d_{x1}C_{x}^{H}|_{t+\frac{\Delta t}{2}} + d_{x1}C_{x}^{H}|_{t-\frac{\Delta t}{2}} - d_{x2}\tilde{E}_{x}|_{t} + \tilde{D}_{x}|_{t-\Delta t} - d_{x3}\tilde{D}_{x}|_{t}$$ $$\tilde{D}_{y}|_{t+\Delta t}= d_{y1}C_{y}^{H}|_{t+\frac{\Delta t}{2}} + d_{y1}C_{y}^{H}|_{t-\frac{\Delta t}{2}} + d_{y2} \sum_{t_{i}=\frac{\Delta t}{2}}^{t-\frac{\Delta t}{2}}C_{y}^{H}|_{t_{i}} - d_{x3}\tilde{E}_{y}|_{t} -d_{y4} \sum_{t_{i}=0}^{t}\tilde{E}_{y}|_{t_{i}} + \tilde{D}_{y}|_{t-\Delta t}$$

Note that linear interpolation has to be used on the time steps that are not defined. The full derivation of these update equations are given here. The constant parameters shown in the update equation above are given by:

$$h_{z1}=\frac{2c_{0}\epsilon_{0}\Delta t}{2\mu_{ry}\epsilon_{0}+\mu_{ry}\sigma'\Delta t}, \,\,\,\, h_{z2}=\frac{2\mu_{ry}\epsilon_{0}-\mu_{ry}\sigma'\Delta t}{2\mu_{ry}\epsilon_{0}+\mu_{ry}\sigma'\Delta t}$$ $$d_{x1}=c_{0}\Delta t, \,\,\,\,\,\,\,\,\,\,\, d_{x2}=2c_{0}\eta_{0}\sigma_{x}\Delta t, \,\,\,\,\,\,\,\,\,\,\, d_{3}=\frac{2\sigma'\Delta t}{\epsilon_{0}}$$ $$d_{y1}=2c_{0}\Delta t\big(\frac{1}{2} + \frac{\sigma' \Delta t}{4 \epsilon_{0}} \big), \,\,\,\,\,\,\,\,\,\,\,\,\,\, d_{y2}=2\sigma' \eta_{0}(c_{0}\Delta t)^{2}$$ $$d_{y3}=2\sigma_{y} \eta_{0}c_{0}\Delta t, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, d_{y4}=2\sigma_{y}\sigma' (\eta_{0}c_{0}\Delta t)^{2}$$

1.2b From a continuous model to a discrete model: Spatially discretization

To discreteze space a Yee grid was used. A Yee grid is a spacial discretization where the electric fields and the magnetic fields are not defined in the same position, and thus spatially displaced from eachother. The motivation for using a Yee grid is that (1) the divergence equations are intrinsically fulfilled, (2) the material boundary conditions are intrinsically fulfilled. The Yee grid for the 2D system is show below:

The curl terms, $$\{C_{z}^{E},C_{x}^{H},C_{y}^{H}\}$$ can then be written as:

$$C_{z}^{E}(i_{x},i_{y})=\frac{\tilde{E}_{y}(i_{x}+1,i_{y})-\tilde{E}_{y}(i_{x},i_{y})}{\Delta x}-\frac{\tilde{E}_{x}(i_{x},i_{y}+1)-\tilde{E}_{x}(i_{x},i_{y})}{\Delta y}$$ $$=c_{1}\tilde{E}_{y}(i_{x}+1,i_{y})-c_{1}\tilde{E}_{y}(i_{x},i_{y})-c_{2}\tilde{E}_{x}(i_{x},i_{y}+1)+c_{2}\tilde{E}_{x}(i_{x},i_{y})$$ $$C_{x}^{H}(i_{x},i_{y})=\frac{H_{z}(i_{x},i_{y})-H_{z}(i_{x},i_{y}-1)}{\Delta y}=c_{2}H_{z}(i_{x},i_{y})-c_{2}H_{z}(i_{x},i_{y}-1)$$ $$C_{y}^{H}(i_{x},i_{y})=-\frac{H_{z}(i_{x},i_{y})-H_{z}(i_{x}-1,i_{y})}{\Delta x}=-c_{1}H_{z}(i_{x},i_{y})+c_{1}H_{z}(i_{x}-1,i_{y})$$

where, $$c_{1}=\frac{1}{\Delta x}$$ and $$c_{2}=\frac{1}{\Delta y}$$.

1.3 The spatial and the temporal resolution

The general spatial resolution $$ds$$ is given by the following:

$$ds:=\frac{\lambda_{0}}{10}=dx=dy,$$ where $$\lambda_{0}$$ is the wavelengh produced by the source while in the vacuum. The selection of the value for the divisor is the main component for a high resolution result, however, as seen later also influces some visual errors. The temporal resolution $$dt$$ is then given by the time light uses to propegate (since the propagation is restricted to 2D the fastest propegation is diagonally, which means it has to be scaled by $$(\sqrt{2})^{-1}$$): $$dt=\frac{ds}{\sqrt{2}c_{0}}.$$

But from certain books is also written as,

$$dt=\frac{1}{c_{0}\sqrt{\big(\frac{1}{\Delta x}\big)^{2}+\big(\frac{1}{\Delta y}\big)^{2}}}$$

1.4 Summary of the iteration procedure

The iteration process is shown in the figure below, and uses the formulas shown above.

For the full derivation of the theory see the following PDF form: [Theory.pdf]

2. The code constructed by C++ and ROOT

The simulation code was constructed using C++ and Cern ROOT. This code is shown below:

  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 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 #include #include #include using namespace std; double** allocMatrix(int dimX,int dimY); void initializeArrayZero(double **arr,double dimX,double dimY); void initializeArrayOne(double **arr,double dimX,double dimY); double findMaxfromArray(double **arr,int dimX,int dimY); void Simulation(){ // Input parameters: double x_max = 1.0; // Size of the chamber, x-direction double y_max = 1.0; // Size of the chamber, y-direction double E0 = 1.0; // (Norm.) Electric field amplitude double frequency = 2.4 * pow(10.0,9.0); // Frequency of the source double diameter_dielectric = 0.4; // Diameter dielectric, ratio of chamber double source_diameter = 0.5; // Diameter EM source, ratio of chamber double source_xpos = 0.5; // Center position EM source, ratio of chamber double epsilonrx = 2; // Relative permittivity of dielectric, x-direction double epsilonry = 2; // Relative permittivity of dielectric, y-direction double mury = 1.0; // Relative permeability of dielectric, y-direction double sigmax = 0.00; // Conductivity of dielectric, x-direction double sigmay = 0.00; // Conductivity of dielectric, y-direction double ntime = 500; // Desired time steps double PML_length_frac = 0.4; // Thickness of absorbing boundary relative to the heigh double PML_strength = 1.0; // Strength of absorbing boundary. 0 = No PML // The larger PML strength is, the more absorption occurs, however it is less numerically stable const TString kGifName = "PML.gif"; // Used for GIF. The name of the gif that shows the result. Is placed in the same directory as the script const TString kGifSleep = "3"; // Used for GIF. Time between each frame in the GIF. Given in centiseconds // (1.0) Constructing the constants: double epsilon0 = 8.8541878128 * pow(10.0,-12.0); // Free space permittivity double mu0 = 1.25663706212 * pow(10.0,-6.0); // Free space permeability double c0 = sqrt(1/(epsilon0*mu0)); // Free space propegation velocity (used in the resolution and normlization) double eta0 = sqrt(mu0/epsilon0); // Free space impedence (used in the normalization) double cMat = 0; if (epsilonrx > epsilonry){ cMat = c0 * sqrt(1/(epsilonrx*mury)); } else{ cMat = c0 * sqrt(1/(epsilonry*mury)); } double lambdaMat = cMat/frequency; // Material wavelength (used in the resolution) double H0 = (eta0)/(c0)*E0; // For the given (Norm.) E-Field, the associated propagating magnetic field has magnitude H0 // (1.1) Resolution of the model: dt,dx,dy double ds = lambdaMat/10; // Spatial resolution of the model // Note: the following is the round function. My version does not support round, so this is a workaround given the input is positive int nx = int(std::floor(x_max/ds + 0.5)); // Getting the most snug integer fit with the partitions given by the length ds. int ny = int(std::floor(y_max/ds + 0.5)); // ------------------\\-------------------- x_max = nx*ds; // Redefining the size of the chamber such that they form a snug integer fit for the grid y_max = ny*ds; // Redefining the size of the chamber such that they form a snug integer fit for the grid double dx = x_max/nx; double dy = y_max/ny; double dt = 1/(c0*sqrt(pow(dx,-2)+pow(dy,-2))); // Courant's stability criterion // (1.2) The key parameters related to the (integer) indices of the PML, source, and the dielectric: // For the PML: int nPML = int(std::floor(PML_length_frac*ny + 0.5)); // For the Dielectric: int ny_dielectric = int(std::floor(0.5*ny + 0.5)); // Setting the dielectric in the middle of the chamber, y-direction int nx_dielectric = int(std::floor(0.5*nx + 0.5)); // Setting the dielectric in the middle of the chamber, x-direction double n_dielectric_radius = 0; double r_dielectric_radius = 0; if (dx > dy){ n_dielectric_radius = 0.5*diameter_dielectric*nx; r_dielectric_radius = dx*n_dielectric_radius; // Used for the figure } else{ n_dielectric_radius = 0.5*diameter_dielectric*ny; r_dielectric_radius = dy*n_dielectric_radius; // Used for the figure } double n_dielectric_radius_square = pow(n_dielectric_radius,2); // For the Source: int iblocks_source_radius; double ix_source_pos; int ix_source_start; int ix_source_end; if (nx % 2 ==1){ // The if and the else statement center the source properly, such that if source_xpos is 0.5 it is symmetric regardless of if nx is odd or even iblocks_source_radius = int(std::floor(0.5*source_diameter*nx + 0.5)); // Approximating the radius by number of index steps ix_source_pos = int(std::floor(source_xpos*nx + 0.5)); // Finding the index where the center of the source is located ix_source_start = ix_source_pos-iblocks_source_radius; // Start position of the source ix_source_end = ix_source_pos+iblocks_source_radius; // End position of the source } else{ iblocks_source_radius = int(std::floor(0.5*source_diameter*nx + 0.5)); // Approximating the radius by number of index steps ix_source_pos = source_xpos*(nx+1); // Finding the index where the center of the source is located ix_source_start = int(std::floor(ix_source_pos-iblocks_source_radius+0.5 +0.5))+1; // Start position of the source ix_source_end = int(std::floor(ix_source_pos+iblocks_source_radius-0.5 +0.5))+1; // End position of the source } // (1.3) The iteration constants that are not needed in matrix form. double c1 = 1/dx; double c2 = 1/dy; double dx1 = c0*dt; // (2.0) Allocating the iteration constants (for those needed in matrix form). int dimX = nx+2; // Size of the problem space in the x-dimension +2 Due to the boundary on both sides int dimY = ny+nPML+2; // Size of the problem space in the y-dimension +2 Due to the boundary on both sides double** erxx = allocMatrix(dimX,dimY); // Relative permittivity in the x-direction for all spatial coordinates double** eryy = allocMatrix(dimX,dimY); // Relative permittivity in the y-direction for all spatial coordinates double** muryy = allocMatrix(dimX,dimY); // Relative permeability in the y-direction for all spatial coordinates double** sigmaxx = allocMatrix(dimX,dimY); // Conductivity of the medium in every spatial coordinate (anisotropy:x-direction) double** sigmayy = allocMatrix(dimX,dimY); // Conductivity of the medium in every spatial coordinate (anisotropy:y-direction) double** sigmaYYprime = allocMatrix(dimX,dimY); // Fictitious conductivity causing loss (absorption) in the PML region in the top of the chamber // (2.1) Setting the standard values based on free space initializeArrayOne(erxx, dimX, dimY); // Is the relative permittivity (x-dirr.) is 1 in vacuum initializeArrayOne(eryy, dimX, dimY); // Is the relative permittivity (y-dirr.) is 1 in vacuum initializeArrayOne(muryy, dimX, dimY); // Is the relative permeability is 1 in vacuum initializeArrayZero(sigmaxx, dimX, dimY); // The conductivity of vacuum is equal to 0 initializeArrayZero(sigmayy, dimX, dimY); // The conductivity of vacuum is equal to 0 initializeArrayZero(sigmaYYprime, dimX, dimY); // The ficticious conductivity of the PML is zero except in the PML region // (2.2) Filling in the values when in the dielectric: // There are alternative and quicker ways to create a circular dielectric element, but we will do it the intuitive way. // Using the circle formula. for (int iy=1; iyUnlink(kGifName.Data()); // Used for GIF. Delete existing GIF TPaveText *t = new TPaveText(0.02,0.92,0.24,0.99,"brNDC"); // Used for GIF. Textbox TH2D *hist = new TH2D("hist","Norm. E-field magnitude",dimX,0,x_max+2*dx,dimY,0,y_max+dy*(2+nPML)); // Used for plot. The main graphical object used for plotting the E-field magnitude hist->GetXaxis()->SetTitle("x-direction, [m]"); hist->GetXaxis()->SetTitleOffset(1.1); // Sets the label of the x-axis with a pleasant spacing from the axis hist->GetYaxis()->SetTitle("y-direction, [m]"); hist->GetYaxis()->SetTitleOffset(1.1); // Sets the label of the y-axis with a pleasant spacing from the axis // Drawing the PLM double dx2cube = dx*(nx/2); double dy2cube = dy*(nPML/2); // The half width of the PML rectangle, in both directions (xy) double forigin[2]; forigin[0] = dx*(nx+2)/2; forigin[1] = dy*(1+ny)+dy*(nPML/2); // The center point of the PML rectangle in the problemspace auto pml = new TBox(forigin[0]-dx2cube,forigin[1]-dy2cube,forigin[0]+dx2cube,forigin[1]+dy2cube); // Creating the rectangle object representing the position of the PML pml->SetFillStyle(3001); // Pseudo transparrent due to the fill style //Drawing the dielectric double rmax = r_dielectric_radius; // The radius of a circle describing the dielectric double forigin2[2]; forigin2[0] = dx*(nx+2)/2; forigin2[1] = dy*((ny/2)+1); // Setting the center point of the dielectric auto dielectric = new TEllipse(forigin2[0],forigin2[1],rmax); // Creating the circular object representing the position of the dielectric dielectric->SetFillStyle(3001); // Pseudo transparrent due to the fill style // Drawing the source double ysource_start = dy*(ny-11-0.5); double ysource_end = dy*(ny-11+0.5); // The start and the end position of the source auto source = new TBox(dx*(ix_source_start),ysource_start,dx*(ix_source_end),ysource_end); // Constructing the source object source->SetFillStyle(3001); // Pseudo transparrent due to the fill style // (4) The Iteration loop for (int it = 0; it < ntime+1; it++){ string itString = Form("i_{t} = %d",it); // Used for plot. Will be inserted into the textbox created by TPaveText double timer = it * dt; // Current simulation time string timeString = Form("t_{i} = %e s",timer); // Used for plot. Will be inserted into the textbox created by TPaveText // (4.1) Implementing the source for (int ix = ix_source_start; ix <= ix_source_end; ix++){ Ex[ix][ny-11]=E0*sin(2*(TMath::Pi())*frequency*(dt*it)); Hz[ix][ny-11]=H0*sin(2*(TMath::Pi())*frequency*(dt*it)); } // (4.2) Iterating the fields. First time zero is plotted // Calculating the curl of the electric field at every point for (int iy=1; iyReset("ICESM"); // clearing the histogram hist->SetMaximum(1.0*E0); hist->SetMinimum(-0.01*E0); for (int iy=0; iyFill(xc,yc,magEc); } } if(it==0){ hist->Draw("colz"); // Used for plot. 2D Color plot hist->SetStats(0); // Remove the histogram statbox } else{ c->Update(); // Used for plot. Updates the canvas with the new values t->Clear(); // Used for plot. Clear the textbox from the previous iteration } t->AddText(itString.c_str()); // Used for plot. Adding the current iteration number to the textbox object t->AddText(timeString.c_str()); // Used for plot. Adding the current time to the textbox object t->AddText(maxmagEString.c_str()); // Used for plot. Adding the max E-field to the textbox object t->Draw(); // Used for plot. Drawing the textbox with the values above pml->Draw("SAME"); // Drawing the rectangular PML dielectric->Draw("SAME"); // Drawing the cicular dielectric source->Draw("SAME"); // Drawing the rectangular source c->Print(kGifName+"+"+kGifSleep); // Used for GIF. Adds another frame with a hold time } c->Print(kGifName+"++"+kGifSleep); //Used for GIF. Makes the last frame of the iteration loop to the beginning } // Allocates space in memory for a 2D matrix of sizeX, sizeY double** allocMatrix(int sizeX, int sizeY){ double** matrix = new double*[sizeX]; for (int ix = 0; ix < sizeX; ix++){ matrix[ix] = new double[sizeY]; } return matrix; } // Fills the memory allocated 2D matrix with zeros void initializeArrayZero(double **arr,double sizeX,double sizeY){ for (int i = 0; i < sizeX; i++){ for (int j = 0; j < sizeY; j++){ arr[i][j] = 0; } } } // Fills the memory allocated 2D matrix with ones void initializeArrayOne(double **arr,double sizeX,double sizeY){ for (int i = 0; i < sizeX; i++){ for (int j = 0; j < sizeY; j++){ arr[i][j] = 1; } } } // Returns the maximum value of a 2D matrix double findMaxfromArray(double **arr,int sizeX,int sizeY){ double max = arr[0][0]; for (int i = 0; i < sizeX; i++){ for (int j = 0; j < sizeY; j++){ if (arr[i][j]>max){ max=arr[i][j]; } } } return max; } 

2.1 The results

Here are some examples with different input parameters showing some of the output that works well, as well as some graphical bugs and other odd behaviors. The roots of these odd behaviors were not found. The simplest case is an isotropic system with no dielectric and no PML with $$\frac{\lambda}{10}$$ defining the grid lengths dx and dy:

As seen above, there is one huge graphical bug present which is causing the weird patterns. This bug seem to be associated with the resolution of the system as the following sorted out the graphical bug: by adjusting the resolution governed by $$ds=\frac{\lambda}{10}$$ into the following, $$ds=\frac{\lambda}{19}$$, causes the pattern to almost completely go away:

To resolve this bug, different definitions of the grid resolution has been invoked and tested as seen above, however all of them produced some graphical issues (to a large or small extent), but most importantly they depend on the other parameters as well such as the dielectric constant, as seen later. However, we can notice further issues: (1) Larger than expected field amplitudes, (2) no divergence around the source. The latter actually arrising from the fact that the fields does not spread propperly from one cell to the next in the x-direction, and will cause issues in the continuity of the wave equation to be broken when we add the dielectric. However, wavebehavior is attained as predicted: (1) reflections from the boundary condition, (2) production of standing waves due these reflections. Since the Yee grid should set the governing divergence equations, the mistake had to be with the Yee Grid. Taking a second look at the derivation and other books on FDTD electrodynamics no mistake was found by the Author. It is still interesting to see how the system performs when adding dielectric elements as this might give some further clues. For the next test we add a dielectric element: $$\epsilon_{r}=2, \mu_{r}=1, \sigma=0, ds=\frac{\lambda}{19}$$:

This as in the first picture, we attain a lot of junk. Just as previously, we change the grid resolution, $$ds=\frac{\lambda}{25}$$ as a temporary solution:

As mentioned earlier, here we see that there is some discontinuity that we should not be attaining from the Maxwell's equation on the boundary of the dielectric, which seems to stem from the inability to propegate along the x-direction. Otherwise we also have no standing waves beneath the dielectric. There thus seems to be two mains bugs, one bug related to the divergence, and one bug related to the divergence and wave behavior in the x-direction that makes waves exhibit odd behavior. Since the PML was also included in the derivation, it would be a shame to not introduce it. Including the PML:

The PML does not seem to absorb the waves as it should be doing, even when manually increasing the strength of the PML to 20 times the original strength it does not abosorb the waves properly. However, the code that simulates free space has issues the field amplitudes, and these large field values might mask the effect of the PML. What is left is two definite bugs that has been characterized, but unsucessfully resolved: The error with the divergence, and the resolution error.

Regarding the author

Created by John Benjamin in association with the course PHYS291 at the University of Bergen.
Email: jlo016@uib.no