parpoi.cpp

pic 模拟程序！面向对象

  if(Idown >=0) { /* valid index */
	  PetscMatSetValue(mpiCoef,Icenter,Idown,downc,INSERT_VALUES);
  }
  if(Ileft >=0) { /* valid index */
	  PetscMatSetValue(mpiCoef,Icenter,Ileft,leftc,INSERT_VALUES);
  }
  if(Iright >=0) { /* valid index */
	  PetscMatSetValue(mpiCoef,Icenter,Iright,rightc,INSERT_VALUES);
  }	
}

// Set the matrix coefficients to reflect the local epsilon and
// boundary conditions, for XY geometry.

void ParallelPoisson::set_coefficient_rz(int j,int k, BCTypes type, Grid *grid) {
	Vector2 **X = grid->getX();
	int J = grid->getJ();
	int K = grid->getK();

	int Iup,Idown,Ileft,Iright,Icenter;  //  global indices...
	Iup = grid->getGlobalIndex(j,k,UP);
	Idown = grid->getGlobalIndex(j,k,DOWN);
	Ileft = grid->getGlobalIndex(j,k,LEFT);
	Iright = grid->getGlobalIndex(j,k,RIGHT);
	Icenter = grid->getGlobalIndex(j,k,HERE);

	if(type == PERIODIC_BOUNDARY)  { //need to diddle the indices a bit
		Iup = grid->getGlobalIndex(j,k,UPWRAP);
		Idown = grid->getGlobalIndex(j,k,DOWNWRAP);
		Ileft = grid->getGlobalIndex(j,k,LEFTWRAP);
		Iright = grid->getGlobalIndex(j,k,RIGHTWRAP);
		//Icenter handled already.
	}
  
	// These calculations in the z direction are the same
	// as for the x-y.
	Scalar dxa =  X[j][k].e1()  -X[MAX(0,j-1)][k].e1();
	Scalar dxb =  X[MIN(J,j+1)][k].e1()-X[j][k].e1();
	Scalar dxave= 0.5*(dxa + dxb);

	// Calculations peculiar to r-z geometry.  radial, in other words.
	
	Scalar ra = X[j][MAX(0,k-1)].e2();
	Scalar r  = X[j][k].e2();
	Scalar rb = X[j][MIN(K,k+1)].e2();
	

	// Need to do some special things at spatial region boundaries!
	// we just correct some things which weren't right for the 
	// normal FREESPACE or DIELECTRIC
	if(j==0 && Ileft != -1) { // we're on the left hand side
		dxave=dxa=dxb;  // uniform mesh assumption.
	}
	if(j==J && Iright != -1) { // we're on the righthand side
		dxave=dxb=dxa;  // uniform mesh assumption.
	}
	if(k==0 && Iup != -1) { // we're on the top side
		// uniform mesh assumption.
		rb = r + (r - ra);
	}
	if(k==K && Idown != -1) { // we're on the top side
		// uniform mesh assumption.
		ra = r - (rb - r);
	}

	Scalar ra2= 0.5*(ra + r);
	Scalar rb2= 0.5*(r +rb);
	Scalar dr2 = rb2*rb2-ra2*ra2;
	Scalar drb = rb -  r;
	Scalar dra = r  - ra;


	Scalar da1a,da1b,da2aa,da2ab,da2ba,da2bb;  // Areas of surfaces, for epsilons
	
	Scalar e1a,e1b,e2a,e2b;  /* epsilon (i-1/2, j), (i+1/2,j),(i,j-1/2),(i,j+1/2) */
	
	// the area weightings to use for the epsilons
	da1a = sqr(r) - sqr(ra2);
	da1b = sqr(rb2)-sqr(r);
	da2aa= dxa;
	da2ab= dxb;
	da2ba= dxa;
	da2bb= dxb;

	int jm = MAX(j-1,0);
	int km = MAX(k-1,0);
	e1a = (epsi[jm][km]*da1a   + epsi[jm][MIN(k,K-1)]*da1b ) / ( da1a  + da1b);
	e1b = (epsi[MIN(j,J-1)][km]*da1a + epsi[MIN(j,J-1)][MIN(k,K-1)]*da1b  ) / ( da1a  + da1b);
	e2a = (epsi[jm][km]*da2aa  + epsi[MIN(j,J-1)][km] * da2ab  ) / ( da2aa + da2ab);
	e2b = (epsi[jm][MIN(k,K-1)]*da2ba  + epsi[MIN(j,J-1)][MIN(k,K-1)] * da2bb  ) / ( da2aa + da2ab);


	Scalar upc,downc,rightc,leftc,centerc;  //  the coefficients for the stencil.

	// recast as freespace any dielectric found in the interior
	if(type==DIELECTRIC_BOUNDARY && (j>0 && j<J && k>0 && k<K))
		type = FREESPACE;
	// similarly recast as freespace any dielectric found on an edge where
	// there is a SRB
	if(type==DIELECTRIC_BOUNDARY && (
			(j==0 && Ileft >=0) ||
			(j==J && Iright>=0) ||
         (k==0 && Iup >=0)   ||
         (k==K && Idown >=0))) 
		type=FREESPACE;
	
	switch(type) {
	 case CONDUCTING_BOUNDARY:// easy case first
		 {
			 upc=downc=leftc=rightc=0;
			 centerc = 1;
			 break;
		 }
	  case SPATIAL_REGION_BOUNDARY:
	  case PERIODIC_BOUNDARY:
	  case FREESPACE: //  no boundaries, center of system.
		  {
			  centerc = upc = downc = leftc = rightc = 0;
			  // cartesian part of the coefficients.
			  if(Ileft >=0) leftc= e1a/(dxa * dxave);
			  if(Iright >=0) rightc = e1b/(dxb * dxave);
			  centerc -= leftc + rightc;

			  if(Iup >=0) upc = e2a*2*ra2/(dr2 * dra);
			  if(Idown >=0) downc = e2b*2*rb2/(dr2 * drb);
			  centerc -= upc + downc;
			  break;
		  }
	  case DIELECTRIC_BOUNDARY:
		  {
			  upc=downc=leftc=rightc=0;
			  if(j==0 && Ileft==-1) { // left wall Neumann
				  // left right coefficients
				  leftc = 0;
				  rightc = 2 * e1b / (dxb * dxb);
				  centerc -= rightc;
				  // r coefficients
				  if(k==0 && Iup == -1) { // neuman/neumann corner bottom left
					  Scalar dria = X[j][1].e2() -X[j][0].e2();
					  upc = 0;
					  downc =  e2b / (dria * dria);
					  centerc -= upc + downc;
				  } else
				  if(k==K && Idown == -1) {  //neumann/naumann corner top left
					  downc = 0;
					  upc = e2a* 2*ra2/(dr2*dra);
					  centerc -= upc + downc;
				  } else { // on left face
					  upc =  e2a*2*ra2/(dr2 * dra);
					  downc = e2b*2*rb2/(dr2 * drb);
					  centerc -= upc + downc;
				  }
			  }
			  if(j==J&&Iright == -1) { // right wall Neumann
				  rightc = 0;
				  leftc = 2 * e1a / (dxa * dxa);
				  centerc -= rightc + leftc;
				  if(k==0&&Iup == -1) {  //neumann/neumann corner
					  Scalar dria =  X[j][k+1].e2() - X[j][k].e2();
					  upc = 0;
					  downc = e2b / (dria * dria);
					  centerc -= upc + downc;
				  } else
				  if(k==K & Idown == -1) { //neumann/neumann corner
					  downc = 0;
					  upc = e2a*2*ra2/(dr2*dra);
					  centerc -= upc + downc;
				  } else { //on right face
					  upc =  e2a*2*ra2/(dr2 * dra);
					  downc = e2b*2*rb2/(dr2 * drb);
					  centerc -= upc + downc;
				  }

			  }

			  if(k==K && Idown == -1) { // top face, not a corner.
				  Scalar dria = X[j][k].e2()-X[j][k-1].e2();
				  upc = 2 * e2a * ra2/(dria * (r*r-ra2 * ra2));
				  downc = 0;
				  leftc = e1a / (dxa * dxave);
				  rightc = e1b / (dxb * dxave);
				  centerc -= upc +downc +leftc +rightc;
			  }
			  break;
		  }
	  case CYLINDRICAL_AXIS:
		  {
			  centerc = upc = downc = rightc = leftc = 0;
			  if(j==0&&Ileft == -1 || j==J && Iright == -1) break;  // case handled above.
			  Scalar dria = X[j][1].e2() - X[j][0].e2();
			  upc = 0;
			  downc = e2b/(dria*dria);
			  leftc = e1a/(dxa*dxave);
			  rightc = e1b/(dxb*dxave);
			  centerc -= upc + downc + leftc + rightc;
			  
			  break;
		  }
	  default: {
		  cout << "Error, unhandled boundary type in ParallelPoisson-xy\n" << endl;
		  break;
	  }
		 
	 }
	// All cases have to set up their coefficients.
	PetscMatSetValue(mpiCoef,Icenter,Icenter,centerc,INSERT_VALUES);
	if(Iup >=0) { /* valid index */
		PetscMatSetValue(mpiCoef,Icenter,Iup,upc,INSERT_VALUES);
	}
	if(Idown >=0) { /* valid index */
		PetscMatSetValue(mpiCoef,Icenter,Idown,downc,INSERT_VALUES);
	}
	if(Ileft >=0) { /* valid index */
		PetscMatSetValue(mpiCoef,Icenter,Ileft,leftc,INSERT_VALUES);
	}
	if(Iright >=0) { /* valid index */
		PetscMatSetValue(mpiCoef,Icenter,Iright,rightc,INSERT_VALUES);
	}	

}


// because this solver uses a different matrix solver
int ParallelPoisson::laplace_solve(Scalar **u_in, Scalar **s, int itermax, Scalar tol_test) {
	int j,k;

	int J = grid->getJ();
	int K = grid->getK();
	for(j=0;j<=J;j++) {
		for(k=0;k<=K;k++) {
			s[j][k] = u_in[j][k];
		}
	}
	return solve(u_in,s,itermax,tol_test);
}

// perform the actual parallel solve.
	
int ParallelPoisson::solve(Scalar **u_in, Scalar **s, int itermax, Scalar tol_test) {

  double *s_petsc;
  double *u_petsc;
  int ierr;
  int J = grid->getJ();
  int K = grid->getK();
  int j,k;
  int iternum;
  Viewer RhoViewer;

  ierr= ViewerASCIIOpen(PETSC_COMM_WORLD,"RhoViewer",&RhoViewer);

  CHKERRA(ierr);
  ierr = ViewerPushFormat(RhoViewer,VIEWER_FORMAT_ASCII_INDEX,"Rho");CHKERRA(ierr);
  // get the array from PetSC-- it's almost certainly quicker
  //to set the values than have PetSC do it:  it likes to send messages.

  ierr=(VecGetArray(mpiRho,&s_petsc));
  CHKERRA(ierr);
  // write the s_data into the local PETSC array
  for(j=0;j<=J;j++) {
	 for(k=0;k<=K;k++) {
		 s_petsc[getVecLocalIndex(j,k)]=s[j][k];
//		 s_petsc[getVecLocalIndex(j,k)]=getVecLocalIndex(j,k);
	 }
  }
  ierr=(VecRestoreArray(mpiRho,&s_petsc));
  VecView(mpiRho,RhoViewer);

  CHKERRA(ierr);
  // Actually call the solve
  ierr=(SLESSolve(petsc_solver,mpiRho,mpiPhi,&iternum));
  CHKERRA(ierr);

  // Copy the data back into our own datastructures.
  ierr=(VecGetArray(mpiPhi,&u_petsc));
  CHKERRA(ierr);
  for(j=0;j<=J;j++) {
	 for(k=0;k<=K;k++) {
		u_in[j][k]=u_petsc[getVecLocalIndex(j,k)];
	 }
  }
  ierr=(VecRestoreArray(mpiPhi,&u_petsc));
  CHKERRA(ierr); 
  return 0;
}

ParallelPoisson::~ParallelPoisson() {
  
  CHKERRA(SLESDestroy(petsc_solver));

}

void ParallelPoisson::init_solve(Grid *grid,Scalar **epsi)
{
	// All the initialization work is done in parpoi.cpp

}

#endif  /* PETSC */