1 / 27

How to calculate potential and electric field in 3D charge-up simulation

EECE 695. How to calculate potential and electric field in 3D charge-up simulation. S.J. Kim, H.J. Lee, and J.K. Lee. Plasma Application Modeling Lab. Department of Electronic and Electrical Engineering Pohang University of Science and Technology. Oct. 11, 2005. Ion: anisotropic

irina
Télécharger la présentation

How to calculate potential and electric field in 3D charge-up simulation

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. EECE 695 How to calculate potential and electric field in 3D charge-up simulation S.J. Kim, H.J. Lee, and J.K. Lee Plasma Application Modeling Lab. Department of Electronic and Electrical Engineering Pohang University of Science and Technology Oct. 11, 2005

  2. Ion: anisotropic Electron: isotropic Plasma Application Modeling, POSTECH What is charge-up effect? High charge-up potential Because of the electron shading effect in high aspect ratio etching, most of the ions reach the bottom of trench. High potential is generated at the bottom of trench. Trajectory of ions is changed.

  3. Plasma Application Modeling, POSTECH Simulation routine and basic assumptions Ion Electron Solve Laplace Equation Update Potential All particles arrive at boundary Move particles (by E-field) 0.2um • Due to the small size of simulation domain( ~ 1um) • Particle flight time is much shorter than the time interval of each • entering particle. • Number of space charge particle is small in the simulation • domain. • Ignore space charge effect in the E-field calculation. • Collisions in the simulation domain are neglected. • mean free path of ions or electrons(~mm) is much longer than • the simulation domain size. • The potential at the top and bottom boundary is the same. • Surface current is neglected.

  4. Plasma Application Modeling, POSTECH Flow chart

  5. Plasma Application Modeling, POSTECH Simulation domain and boundary conditions

  6. Plasma Application Modeling, POSTECH Solving electrostatic potential in 3D charge-up simulation(1) Space gradient of dielectric constant being considered, Poisson’s equation is as follows: Surface charges in the dielectric surfaces is only considered in right-hand. For solving PDE numerically, Poisson’s equation is represented as follows:

  7. Plasma Application Modeling, POSTECH Solving electrostatic potential in 3D charge-up simulation(2) Poisson’s equation is described as Matrix equation

  8. Plasma Application Modeling, POSTECH Successive OverRelaxation (SOR) • Elliptic PDE • SOR method for solving PDEs where,  = 1 : Gauss-seidel (GS) method 0<  < 1 : underrelaxation 1<  <2 : overrelaxation How to choose optimal  If Jacobi is the spectral radius of the Jacobi iteration, where,  =  : Dirichlet or Neumann boundary conditions.  = 2 : Periodic boundary conditions

  9. Plasma Application Modeling, POSTECH SOR algorithm for 2nd-order elliptic PDE From 2D finite difference model, Iterative procedure: Residual is calculated as follows: Finding optimal  in SOR with Chebyshev acceleration

  10. Plasma Application Modeling, POSTECH Solving electric field in 3D charge-up simulation Electric field is calculated by using Gauss’s law. i) Cases without changes in  ii) Cases with changes in  x-directional electric field in the plane of yz1 : where,

  11. Plasma Application Modeling, POSTECH Main.c /***********************************************************/ /* The main physics loop */ void XGMainLoop() { int isp, phi_flag; it++; phi_flag=0; /* to determine if it is necessary to solve Laplace's Equaton */ for(isp=0; isp<nsp; isp++) { (*move_ptr)(isp); boundary(isp); phi_flag+=sp[isp].np; } if (theRunWithXFlag) history2(); if (phi_flag==0) { /* if no particle is left */ t+=DT; fields(); load(); if (theRunWithXFlag) history1(); } } /***************************************************************/ void display_title() { puts("\nEtch-3d Version 1.0"); puts("Made by Hae June Lee"); puts("Plasma Application Modeling Group"); puts("Pohang University of Science and Technology\n"); } /***************************************************************/

  12. Main.c #define SIDE_ONN_C 14 #define SIDE_POP_C 15 #define SIDE_PON_C 16 #define SIDE_NOP_C 17 #define SIDE_NON_C 18 #define CORN_PPP_C 19 #define CORN_PNP_C 20 #define CORN_NPP_C 21 #define CORN_NNP_C 22 #define CORN_PPN_C 23 #define CORN_PNN_C 24 #define CORN_NPN_C 25 #define CORN_NNN_C 26 #define SIDE_PPO_V 27 #define SIDE_PNO_V 28 #define SIDE_NPO_V 29 #define SIDE_NNO_V 30 #define SIDE_OPP_V 31 #define SIDE_OPN_V 32 #define SIDE_ONP_V 33 #define SIDE_ONN_V 34 #define SIDE_POP_V 35 #define SIDE_PON_V 36 #define SIDE_NOP_V 37 #define SIDE_NON_V 38 #define CORN_PPP_V 39 #define CORN_PNP_V 40 #define CORN_NPP_V 41 #define CORN_NNP_V 42 #define CORN_PPN_V 43 #define CORN_PNN_V 44 #define CORN_NPN_V 45 #define CORN_NNN_V 46 // F means flat in xz surface #define CORN_PPP_F 47 #define CORN_NPP_F 48 #define CORN_PPN_F 49 #define CORN_NPN_F 50 // U means there is a side in up-direction #define CORN_PPP_U 51 #define CORN_NPP_U 52 #define CORN_PPN_U 53 #define CORN_NPN_U 54 /**************************************/ /* defining internal structures */ #define ACCUMULATION 0 #define INJECTION 1 #define OPEN 3 /******************************************/ /* def.h */ #define EPS0 8.8542e-12 /* (F/m) */ #define NperTORR 8.3221e20 #define HISTMAX 2048 #define HISTMAX2 4096 #define NSMAX 2 #define QUIET_START 0 #define RANDOM 1 #define GROUNDED 0 #define VOLTAGE_D 1 #define CURRENT_D 2 #define DIRICHLET 0 #define NEUMANN 1 #define PERIODIC 2 #define SYMMETRIC 1 #define LEFT 0 #define RIGHT 1 #define UP 2 #define DOWN 3 #define SIDES 2 /* ---- PNO means positive, negative, and ordinary in each directions, and C and V means concave and convex. --*/ #define NO_FACE 0 #define FACE_X_UP 1 #define FACE_X_DOWN 2 #define FACE_Y_UP 3 #define FACE_Y_DOWN 4 #define FACE_Z_UP 5 #define FACE_Z_DOWN 6 #define SIDE_PPO_C 7 #define SIDE_PNO_C 8 #define SIDE_NPO_C 9 #define SIDE_NNO_C 10 #define SIDE_OPP_C 11 #define SIDE_OPN_C 12 #define SIDE_ONP_C 13

  13. Plasma Application Modeling, POSTECH field.c (1) #include "et3d.h" #define my_div(x, y) ((fabs(y) <= 1e-25) ? 0 : x/y) void set_voltage(); float hdx, hdy, hdz; /***************************************************************/ void fields() { register int isp, i, j, k, n, ip1, im1, jp1, jm1, kp1, km1; float ex_mag, ey_mag, ez_mag, e_mag, temp, epsi1, epsi2; static int init_field_flag=1, stat_counter=0; static float ***epsx, ***epsy, ***epsz; if(init_field_flag) { hdx = 0.5*idx; hdy = 0.5*idy; hdz = 0.5*idz; epsx = tensor3(ngx, ngy, ngz); epsy = tensor3(ngx, ngy, ngz); epsz = tensor3(ngx, ngy, ngz); for(i=0; i<ngx; i++) for(j=0; j<ngy; j++) for (k=0; k<ngz; k++) eps_array[i][j][k] *= EPS0; for(i=0; i<ngx; i++) for(j=0; j<ngy; j++) for (k=0; k<ngz; k++) { if (bc_x==PERIODIC) im1= (i)?i-1:ncx-1; else im1= (i)?i-1:0; jm1= (j)?j-1:0; km1= (k)?k-1:ncz-1; epsx[i][j][k]=0.25*(eps_array[i][j][k] + eps_array[i][jm1][k] + eps_array[i][j][km1] + eps_array[i][jm1][km1]); epsy[i][j][k]=0.25*(eps_array[i][j][k] + eps_array[im1][j][k] + eps_array[i][j][km1] + eps_array[im1][j][km1]); epsz[i][j][k]=0.25*(eps_array[i][j][k] + eps_array[im1][j][k] + eps_array[i][jm1][k] + eps_array[im1][jm1][k]); } set_poisson_coefficient(); init_field_flag=0; } Setting initial values Symmetric Symmetric (y-direction) Periodic (z-direction) Call function which sets coefficients for solving poisson’s eq.

  14. Plasma Application Modeling, POSTECH field.c (2) /**************************************************/ /*** Get the charge densities.. ****/ for (i=0; i<= ncx; i++) for (j=0; j<= ncy; j++) for (k=0; k<ngz; k++) { sigma[i][j][k] = 0.0; source[i][j][k]= 0; } for (isp=0; isp< nsp; isp++) { for (i=0; i< ngx; i++) for (j=0; j< ngy; j++) for (k=0; k<ngz; k++) { source[i][j][k]-= sp[isp].q_per_cell*sp[isp].sigma[i][j][k]; sigma[i][j][k]+= sp[isp].q*sp[isp].sigma[i][j][k]; /* 아직 area로 나눠지지 않았음. */ } } /*******************************************/ /*** Solve Poisson's eqn with SOR method */ set_voltage(); sor_num=sor_cheb(aa, bb, cc, dd, ee, ff, gg, source, phi, 1.0); Set =0 at the xzs plane.

  15. Plasma Application Modeling, POSTECH field.c (3) /**********************************************/ /* Calculation of surface charge of conductor */ /* Divergence theorem 을 이용해서 (i,j,k) grid를 둘러싼 Cell에서의 Q = epsilon* \int E ds - \int rho dV 로 여분의 charge를 구할 수 있다. 이것을 area로 나누면 균일한 surface charge를 갖게 된다. */ for (i=0; i<ngx; i++) for (j=1; j<ngy; j++) for (k=0; k<ngz; k++) { if(grid_mask[i][j][k]==3 && face[i][j][k]) { if (bc_x==PERIODIC) { ip1= (i==ncx)?1:i+1; im1= (i)?i-1:ncx-1; } else { ip1= (i==ncx)?ncx:i+1; im1= (i)?i-1:0; } jp1= (j==ncy)?ncy:j+1; kp1= (k==ncz)?1:k+1; jm1= (j)?j-1:0; km1= (k)?k-1:ncz-1; sigma[i][j][k] = dz*dy*(phi[i][j][k]-phi[ip1][j][k])*idx* 0.25*(eps_array[i][j][k]+eps_array[i][jm1][k]+ eps_array[i][j][km1]+eps_array[i][jm1][km1]); sigma[i][j][k]+= dz*dy*(phi[i][j][k]-phi[im1][j][k])*idx* 0.25*(eps_array[im1][j][k]+eps_array[im1][jm1][k]+ eps_array[im1][j][km1]+eps_array[im1][jm1][km1]); sigma[i][j][k]+= dz*dx*(phi[i][j][k]-phi[i][jp1][k])*idy* 0.25*(eps_array[i][j][k]+eps_array[im1][j][k]+ eps_array[i][j][km1]+eps_array[im1][j][km1]); sigma[i][j][k]+= dz*dx*(phi[i][j][k]-phi[i][jm1][k])*idy* 0.25*(eps_array[i][jm1][k]+eps_array[im1][jm1][k]+ eps_array[i][jm1][km1]+eps_array[im1][jm1][km1]); sigma[i][j][k]+= dx*dy*(phi[i][j][k]-phi[i][j][kp1])*idz* 0.25*(eps_array[i][j][k]+eps_array[i][jm1][k]+ eps_array[im1][j][k]+eps_array[im1][jm1][k]); sigma[i][j][k]+= dx*dy*(phi[i][j][k]-phi[i][j][km1])*idz* 0.25*(eps_array[i][j][km1]+eps_array[i][jm1][km1]+ eps_array[im1][j][km1]+eps_array[im1][jm1][km1]); } if (sigma[i][j][k]) sigma[i][j][k] = my_div(sigma[i][j][k], area[i][j][k]); } Note Note ngx=ncx+1 ngy=ncy+1 ngz=ncz+1 Note ip1(i+1)=1, 2, 3, …, ncx, 1 Grid mask = 0 : empty 1 : inside dielectric 2 : dielectric surface 3 : conductor 4 : simulation boundary im1(i-1)=ncx-1, 0, 1, …, ncx-1 (symmetic boundary) ip1(i+1)=1, 2, 3, …, ncx, ncx im1(i-1)=0, 0, 1, …, ncx-1

  16. Plasma Application Modeling, POSTECH field.c (4) /*******************************************/ /*** calculate the E field... ***/ for(i=0; i<=ncx; i++) for(j=0; j<=ncy; j++) for(k=0; k<=ncz; k++) { if (bc_x==PERIODIC) { ip1=(i==ncx) ? 1:i+1; im1=(i) ? i-1:ncx-1; } else { ip1=(i==ncx) ? ncx:i+1; im1=(i) ? i-1:0; } jp1= (j==ncy)?ncy:j+1; kp1= (k==ncz)?1:k+1; jm1= (j)?j-1:0; km1= (k)?k-1:ncz-1; if (j==0) { ex[i][j][k]=ez[i][j][k]=0; ey[i][j][k] = (phi[i][0][k]-phi[i][1][k])*idy ; } else if (j==ncy) { ex[i][j][k]=ez[i][j][k]=0; ey[i][j][k] = (phi[i][jm1][k]-phi[i][j][k])*idy; } else if (grid_mask[i][j][k]==3) { /* biased conductor의 경우, surface charge가 외부 회로에 의해 영향을 받기 때문에 이 방법을 쓸 수 없다. */ ex[i][j][k] = ey[i][j][k] = ez[i][j][k] = 0; if (face[i][j][k] ==FACE_Z_UP) ez[i][j][k] = sigma[i][j][k] /epsz[i][j][k]; else if (face[i][j][k]==FACE_Z_DOWN) ez[i][j][k] = -sigma[i][j][k] /epsz[i][j][km1]; Note

  17. Plasma Application Modeling, POSTECH field.c (5) else if (face[i][j][k]==FACE_Y_UP) ey[i][j][k] = sigma[i][j][k] /epsy[i][j][k]; else if (face[i][j][k]==FACE_Y_DOWN) ey[i][j][k] = -sigma[i][j][k] /epsy[i][jm1][k]; else if (face[i][j][k]==FACE_X_UP) ex[i][j][k] = sigma[i][j][k] /epsx[i][j][k]; else if (face[i][j][k]==FACE_X_DOWN) ex[i][j][k] = -sigma[i][j][k] /epsx[im1][j][k]; else if (face[i][j][k]==CORN_PPP_F || face[i][j][k]==CORN_NPP_F || face[i][j][k]==CORN_PPN_F || face[i][j][k]==CORN_NPN_F) ey[i][j][k] = sigma[i][j][k] /epsy[i][j][k]; else if (face[i][j][k]==SIDE_PPO_V) { epsi2=(eps_array[im1][j][km1]+eps_array[im1][j][k]+ eps_array[i][jm1][km1]+eps_array[i][jm1][k]+ eps_array[i][j][km1]+eps_array[i][j][k])/6.0; ex_mag = (phi[i][j][k] -phi[ip1][j][k])*idx; ey_mag = (phi[i][j][k] -phi[i][jp1][k])*idy; e_mag = sqrt(ex_mag*ex_mag + ey_mag*ey_mag); ex[i][j][k] = fabs(sigma[i][j][k])*my_div(ex_mag, e_mag)/epsi2; ey[i][j][k] = fabs(sigma[i][j][k])*my_div(ey_mag, e_mag)/epsi2; } else if (face[i][j][k]==SIDE_NPO_V) { epsi2=(eps_array[im1][jm1][km1]+eps_array[im1][jm1][k]+ eps_array[im1][j][km1]+eps_array[im1][j][k]+ eps_array[i][j][km1]+eps_array[i][j][k])/6.0; ex_mag = (phi[im1][j][k] -phi[i][j][k])*idx; ey_mag = (phi[i][j][k] -phi[i][jp1][k])*idy; e_mag = sqrt(ex_mag*ex_mag + ey_mag*ey_mag); ex[i][j][k] = fabs(sigma[i][j][k])*my_div(ex_mag, e_mag)/epsi2; ey[i][j][k] = fabs(sigma[i][j][k])*my_div(ey_mag, e_mag)/epsi2; } else if (face[i][j][k]==SIDE_OPP_V) { epsi2=(eps_array[i][j][km1]+eps_array[im1][j][km1]+ eps_array[i][j][k]+eps_array[im1][j][k]+ eps_array[i][jm1][k]+eps_array[im1][jm1][k])/6.0; ez_mag = (phi[i][j][k] -phi[i][j][kp1])*idz; ey_mag = (phi[i][j][k] -phi[i][jp1][k])*idy; e_mag = sqrt(ez_mag*ez_mag + ey_mag*ey_mag); ez[i][j][k] = fabs(sigma[i][j][k])*my_div(ez_mag, e_mag)/epsi2; ey[i][j][k] = fabs(sigma[i][j][k])*my_div(ey_mag, e_mag)/epsi2; }

  18. Plasma Application Modeling, POSTECH field.c (6) else if (face[i][j][k]==SIDE_OPN_V) { epsi2=(eps_array[i][j][km1]+eps_array[im1][j][km1]+ eps_array[i][j][k]+eps_array[im1][j][k]+ eps_array[i][jm1][km1]+eps_array[im1][jm1][km1])/6.0; ez_mag = (phi[i][j][km1] -phi[i][j][k])*idz; ey_mag = (phi[i][j][k] -phi[i][jp1][k])*idy; e_mag = sqrt(ez_mag*ez_mag + ey_mag*ey_mag); ez[i][j][k] = fabs(sigma[i][j][k])*my_div(ez_mag, e_mag)/epsi2; ey[i][j][k] = fabs(sigma[i][j][k])*my_div(ey_mag, e_mag)/epsi2; } } /* end of conductor surface */ else if (grid_mask[i][j][k]==2) { ex[i][j][k] = hdx*(phi[im1][j][k] - phi[ip1][j][k]); ey[i][j][k] = hdy*(phi[i][jm1][k] - phi[i][jp1][k]); ez[i][j][k] = hdz*(phi[i][j][km1] - phi[i][j][kp1]); if (face[i][j][k]==FACE_X_UP) { /* For the case that epsilon doesn't change in y,z directions */ epsi1=epsx[im1][j][k]; epsi2=epsx[i][j][k]; ex[i][j][k]=sigma[i][j][k] + epsi1*( (phi[im1][j][k]-phi[i][j][k])*idx - 0.5*dx* ( idy*(2*phi[i][j][k]-phi[i][jp1][k]-phi[i][jm1][k])*idy +idz*(2*phi[i][j][k]-phi[i][j][kp1]-phi[i][j][km1])*idz) ); ex[i][j][k]/=epsi2; } else if (face[i][j][k]==FACE_X_DOWN) { epsi2=epsx[im1][j][k]; epsi1=epsx[i][j][k]; ex[i][j][k]=sigma[i][j][k] + epsi1*( (phi[ip1][j][k]-phi[i][j][k])*idx - 0.5*dx* ( idy*(2*phi[i][j][k]-phi[i][jp1][k]-phi[i][jm1][k])*idy +idz*(2*phi[i][j][k]-phi[i][j][kp1]-phi[i][j][km1])*idz) ); ex[i][j][k]/=-epsi2; } else if (face[i][j][k]==FACE_Y_UP) { epsi1=epsy[i][jm1][k]; epsi2=epsy[i][j][k]; ey[i][j][k]=sigma[i][j][k] + epsi1*( (phi[i][jm1][k]-phi[i][j][k])*idy - 0.5*dy* ( idx*(2*phi[i][j][k]-phi[ip1][j][k]-phi[im1][j][k])*idx +idz*(2*phi[i][j][k]-phi[i][j][kp1]-phi[i][j][km1])*idz) ); ey[i][j][k]/=epsi2; } else if (face[i][j][k]==FACE_Y_DOWN) { epsi2=epsy[i][jm1][k]; Note

  19. Plasma Application Modeling, POSTECH field.c (7) epsi1=epsy[i][j][k]; ey[i][j][k]=sigma[i][j][k] + epsi1*( (phi[i][jp1][k]-phi[i][j][k])*idy - 0.5*dy* ( idx*(2*phi[i][j][k]-phi[ip1][j][k]-phi[im1][j][k])*idx +idz*(2*phi[i][j][k]-phi[i][j][kp1]-phi[i][j][km1])*idz) ); ey[i][j][k]/=-epsi2; } else if (face[i][j][k]==FACE_Z_UP) { epsi1=epsz[i][j][km1]; epsi2=epsz[i][j][k]; ez[i][j][k]=sigma[i][j][k] + epsi1*( (phi[i][j][km1]-phi[i][j][k])*idz - 0.5*dz* ( idx*(2*phi[i][j][k]-phi[ip1][j][k]-phi[im1][j][k])*idx +idy*(2*phi[i][j][k]-phi[i][jp1][k]-phi[i][jm1][k])*idy) ); ez[i][j][k]/=epsi2; } else if (face[i][j][k]==FACE_Z_DOWN) { epsi2=epsz[i][j][km1]; epsi1=epsz[i][j][k]; ez[i][j][k]=sigma[i][j][k] + epsi1*( (phi[i][j][kp1]-phi[i][j][k])*idz - 0.5*dz* ( idx*(2*phi[i][j][k]-phi[ip1][j][k]-phi[im1][j][k])*idx +idy*(2*phi[i][j][k]-phi[i][jp1][k]-phi[i][jm1][k])*idy) ); ez[i][j][k]/=-epsi2; } } else { ex[i][j][k] = hdx*(phi[im1][j][k] - phi[ip1][j][k]); ey[i][j][k] = hdy*(phi[i][jm1][k] - phi[i][jp1][k]); ez[i][j][k] = hdz*(phi[i][j][km1] - phi[i][j][kp1]); } } } /***************************************************************/ void set_voltage() { register int i, k; for (i=0;i<ngx;i++) for (k=0;k<ngz;k++) { phi[i][0][k]= 0; phi[i][ncy][k]=0; } }

  20. Plasma Application Modeling, POSTECH SOR3d.c(1) #include "et3d.h" /* Successive Overrelaxation (SOR) with Chebyshev acceleration */ int sor_cheb(float ***a, float ***b, float ***c, float ***d, float ***e, float ***f, float ***g, float ***s, float ***u, double omega) { register int i,j,k; int im1,ip1,jm1,jp1,km1,kp1; int iterations; float residue, change, max_change; int i_order,j_order,k_order; int i_max2, j_max2, k_max2; float anorm, old_anorm=0.0; /* rjac is input as the spectral radius of the Jacobi iteration, or an estimate of it */ if (bc_x == SYMMETRIC) i_max2 = ncx; else if (bc_x == PERIODIC) i_max2 = ncx-1; if (bc_y == SYMMETRIC) j_max2 = ncy; else if (bc_y == PERIODIC) j_max2 = ncy-1; if (bc_z == SYMMETRIC) k_max2 = ncz; else if (bc_z == PERIODIC) k_max2 = ncz-1; for(iterations=1;iterations<=SOR_MAX;iterations++) { anorm = 0.0; max_change = 0; for (i_order = 0; i_order < 2; i_order++) { k_order = i_order; for (i=0;i<=i_max2;i++) { j_order = k_order; for(k=0;k<=k_max2;k++) { /* if i_order = 0, i is even => j_order = 0, even j is swept i is odd => j_order = 1, odd j is swept if i_order = 1, i is even => j_order = 1, odd j is swept i is odd => j_order = 0, even j is swept */

  21. Plasma Application Modeling, POSTECH SOR3d.c(2) for(j=j_order;j<=j_max2;j+=2) { if(grid_mask[i][j][k] !=3 && j>0 && j<ncy) { im1 = i - 1; ip1 = i + 1; jm1 = j - 1; jp1 = j + 1; km1 = k - 1; kp1 = k + 1; if (bc_x == SYMMETRIC) { if(im1 == -1) im1 = 1; if(ip1 == ncx+1) ip1 = ncx-1; } else if (bc_x == PERIODIC) { if(im1 == -1) im1 = ncx-1; if(ip1 == ncx) ip1 = 0; } if (bc_y == SYMMETRIC) { if(jm1 == -1) jm1 = 1; if(jp1 == ncy+1) jp1 =ncy-1; } else if (bc_y == PERIODIC) { if(jm1 == -1) jm1 = ncy-1; if(jp1 == ncy) jp1 = 0; } if (bc_z == SYMMETRIC) { if(km1 == -1) km1 = 1; if(kp1 == ncz+1) kp1 =ncz-1; } else if (bc_z == PERIODIC) { if(km1 == -1) km1 = ncz- 1; if(kp1 == ncz) kp1 = 0; } /* residual */ residue = a[i][j][k]*u[im1][j][k] + b[i][j][k]*u[ip1][j][k] + c[i][j][k]*u[i][jm1][k] + d[i][j][k]*u[i][jp1][k] + e[i][j][k]*u[i][j][km1] + f[i][j][k]*u[i][j][kp1] + g[i][j][k]*u[i][j][k] - s[i][j][k]; change = omega*residue/g[i][j][k]; if (i_order == 0) anorm += change*change; u[i][j][k] -= change; Note Residual (2D): Note In 2D case,

  22. Plasma Application Modeling, POSTECH SOR3d.c(3) if (bc_x == PERIODIC && i == 0) u[ncx][j][k] = u[0][j][k]; if (bc_y == PERIODIC && j == 0) u[i][ncy][k] = u[i][0][k]; if (bc_z == PERIODIC && k == 0) u[i][j][ncz] = u[i][j][0]; max_change = max(max_change, fabs(change)); } } /* 0 , 1 swap */ j_order = 1 - j_order; } /* 0 , 1 swap */ k_order = 1 - k_order; } if (iterations==1 && k_order==0) omega=1.0/(1.0-0.5*rjac_square); else omega=1.0/(1.0-0.25*rjac_square*omega); } if (bc_x == PERIODIC && bc_z == PERIODIC) for (j=0;j<ngy;j++) u[ncx][j][ncz]=u[0][j][0]; if (max_change < tol_pois) return iterations; #if 0 if (!cheb_flag) { if (iterations > 1) { opt_omega=obtain_opt_omega(anorm, old_anorm, omega); } old_anorm = anorm; } #endif } /* End of Iteration Loop */ fprintf(stderr, " Max. # of Iteration exceeded in SOR\n"); return -1; } Note Finding optimal  in SOR with Chebyshev acceleration

  23. Plasma Application Modeling, POSTECH SOR3d.c(4) /* rjac is input as the spectral radius of the Jacobi iteration, or an * estimate of it */ void set_rjac_square() { double rjac,factor; factor = dx*dx/dy/dy; if (bc_x == SYMMETRIC && bc_y == SYMMETRIC) rjac = (cos(M_PI/ncx)+factor*cos(M_PI/ncy)); else if (bc_x == PERIODIC && bc_y == PERIODIC) rjac = (cos(2.*M_PI/ncx)+factor*cos(2.*M_PI/ncy)); else if (bc_x == PERIODIC && bc_y == SYMMETRIC) rjac = (cos(2.*M_PI/ncx)+factor*cos(M_PI/ncy)); else if (bc_x == SYMMETRIC && bc_y == PERIODIC) rjac = (cos(M_PI/ncx)+factor*cos(2.*M_PI/ncy)); rjac /= (1.+factor); rjac_square = rjac*rjac; } double obtain_opt_omega(double anorm, double old_anorm, double omega) { double mu_W, mu_J, opt_omega; mu_W = sqrt(anorm/old_anorm); mu_J = (mu_W+omega-1.0)/(omega*sqrt(mu_W)); opt_omega = 2.0/(1.0+sqrt(1.0-mu_J*mu_J)); return opt_omega; } void set_poisson_coefficient() { register int i, j, k; int ii, im1, jm1, jj, km1, kk, isp; for (i=0;i<ngx;i++) for (j=0;j<ngy;j++) for (k=0;k<ngz;k++) { if (grid_mask[i][j][k] !=3) { ii = adjust_x(i); jj = adjust_y(j); kk = adjust_z(k); im1 = adjust_x(i-1); jm1 = adjust_y(j-1); km1 = adjust_z(k-1); Note where,  =  : Dirichlet or Neumann boundary conditions.  = 2 : Periodic boundary conditions Note

  24. Plasma Application Modeling, POSTECH SOR3d.c(5) aa[i][j][k]= eps_array[im1][jj][kk] + eps_array[im1][jm1][kk] + eps_array[im1][jj][km1]+ eps_array[im1][jm1][km1]; bb[i][j][k]= eps_array[ii][jj][kk] + eps_array[ii][jm1][kk] + eps_array[ii][jj][km1]+ eps_array[ii][jm1][km1]; cc[i][j][k]= eps_array[ii][jm1][kk] + eps_array[im1][jm1][kk] + eps_array[ii][jm1][km1]+ eps_array[im1][jm1][km1]; dd[i][j][k]= eps_array[ii][jj][kk] + eps_array[im1][jj][kk] + eps_array[ii][jj][km1]+ eps_array[im1][jj][km1]; ee[i][j][k]= eps_array[ii][jj][km1] + eps_array[im1][jj][km1] + eps_array[ii][jm1][km1]+ eps_array[im1][jm1][km1]; ff[i][j][k]= eps_array[ii][jj][kk] + eps_array[im1][jj][kk] + eps_array[ii][jm1][kk]+ eps_array[im1][jm1][kk]; } aa[i][j][k]*=idx2; bb[i][j][k]*=idx2; cc[i][j][k]*=idy2; dd[i][j][k]*=idy2; ee[i][j][k]*=idz2; ff[i][j][k]*=idz2; gg[i][j][k]= -(aa[i][j][k] + bb[i][j][k] + cc[i][j][k] + dd[i][j][k] + ee[i][j][k] + ff[i][j][k]); } } int adjust_x(int i) { if (bc_x == SYMMETRIC) { /* not valid for all integer range */ if (i < 0) i = (-i-1) % ncx; else if (i >= ncx) i = ncx - 1 - i % ncx; } else if (bc_x == PERIODIC) { if (i < 0) i = ncx - 1 - (-i-1) % ncx; else if (i >= ncx) i = i % ncx; } return i; } Note idx2=0.25*idx*idx

  25. Plasma Application Modeling, POSTECH SOR3d.c(6) int adjust_y(int j) { if (bc_y == SYMMETRIC) { if (j < 0) j = (-j-1) % ncy; else if (j >= ncy) j = ncy - 1 - j % ncy; } else if (bc_y == PERIODIC) { if (j < 0) j = ncy - 1 - (-j-1) % ncy; else if (j >= ncy) j = j % ncy; } return j; } int adjust_z(int k) { if (bc_z == SYMMETRIC) { if (k < 0) k = (-k-1) % ncz; else if (k >= ncz) k = ncz - 1 - k % ncz; } else if (bc_z == PERIODIC) { if (k < 0) k = ncz - 1 - (-k-1) % ncz; else if (k >= ncz) k = k % ncz; } return k; }

  26. Plasma Application Modeling, POSTECH Motions of electrons and ions Electron Saturated Initial Ion Initial Saturated

  27. B C A C B A Plasma Application Modeling, POSTECH Potential profiles Aspect ratio = 7 AR=3, long z-direction AR=3

More Related