multigrid.cpp

pic 模拟程序！面向对象

		error = errorest();
		for (i = 0; (i < iter2) && (error > error0*tol); i++){
			mgrelax();
			for (int foo =0; foo<3; foo++)
	 			GSRB(soln[level],rhs[level],length[level],GSRBCoeff[level], PeriodicFlagX1, PeriodicFlagX2); 
			error = errorest();
#ifdef MGRID_DEBUG
			printf("%d "     "  %g\n",i,error/error0);
#endif
		}
#ifdef MGRID_DEBUG
//		printf("%d "     "  %g\n",i,error/error0);
#endif
		if (error > error0*tol)
			cout << "Multigrid did not reach tolerance.  Residual = "<< error/error0 << endl;
	}
	for (i = 0; i <= length[0].e1(); i++)   
		for (j = 0; j <= length[0].e2(); j++)
			phi[i][j] = oldphi[i][j] = soln[0][i][j];
	
	return 0;
}

/*void Multigrid::init_solve(Grid *grid);
	{

	}
	*/

void Multigrid::set_coefficient(int j, int k, BCTypes type, Grid *grid)
{
// not used
//	level = 0;	
//	MGSetCoeff(j,k,type,level);

}

void Multigrid::MGSetCoeff(int j, int k, BCTypes type, int level)
{
	Scalar iCoeff;
	
	Grid* grid = MultiGrid[level];
	int J,K;
	J=grid->getJ();
	K=grid->getK();

	Scalar* resCjk, *GSRBCjk;
	resCjk = resCoeff[level][j][k];
	GSRBCjk = GSRBCoeff[level][j][k];

	if((type==DIELECTRIC_BOUNDARY) && (((0<k)&&(k<K))&&((0<j)&&(j<J)))) type = FREESPACE;


	switch(type)
		{
		case FREESPACE:
			{
				resCjk[NORTH] = (epsi_local[level][j][k]+epsi_local[level][j-1][k])/2*
					grid->dS2Prime(j,k)/grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
				resCjk[EAST] = (epsi_local[level][j][k]+epsi_local[level][j][k-1])/2*
					grid->dS1Prime(j,k)/grid->dl1(j,k)/grid->get_halfCellVolumes()[j][k];
				resCjk[SOUTH] = (epsi_local[level][j][k-1]+epsi_local[level][j-1][k-1])/2*
					grid->dS2Prime(j,k-1)/grid->dl2(j,k-1)/grid->get_halfCellVolumes()[j][k];
				resCjk[WEST] = (epsi_local[level][j-1][k]+epsi_local[level][j-1][k-1])/2*
					grid->dS1Prime(j-1,k)/grid->dl1(j-1,k)/grid->get_halfCellVolumes()[j][k];
				resCjk[SOURCE] = -resCjk[NORTH]- resCjk[EAST]
					-resCjk[SOUTH] - resCjk[WEST];
				
				iCoeff = -1/resCjk[SOURCE];
				GSRBCjk[NORTH] = resCjk[NORTH]*iCoeff;
				GSRBCjk[EAST] = resCjk[EAST]*iCoeff;
				GSRBCjk[SOUTH] = resCjk[SOUTH]*iCoeff;
				GSRBCjk[WEST] = resCjk[WEST]*iCoeff;
				GSRBCjk[SOURCE] = -iCoeff;
			}
			break;
		case CONDUCTING_BOUNDARY:
			{
				resCjk[NORTH] = 0.0;
				resCjk[EAST] = 0.0;
				resCjk[SOUTH] = 0.0;
				resCjk[WEST] = 0.0;
				resCjk[SOURCE] = 0.0;

				GSRBCjk[NORTH] = 0.0;
				GSRBCjk[EAST] = 0.0;
				GSRBCjk[SOUTH] = 0.0;
				GSRBCjk[WEST] = 0.0;
				GSRBCjk[SOURCE] = 0.0;
				break;
			}
		case PERIODIC_BOUNDARY:
			{
				int jm, km, jp;
				jm = j-1;
				jp = j;
				km = k-1;
				if (PeriodicFlagX1&&PeriodicFlagX2){
					if (j==0)
						jm = J-1;
					if (k==0) 
						km = K-1;
					if (j==J)
						j=0;
					if (k==K)
						k=0;
					
					resCjk[NORTH] = (epsi_local[level][j][k]+epsi_local[level][jm][k])/2*
						grid->dS2Prime(j,k)/grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
					resCjk[EAST] = (epsi_local[level][j][k]+epsi_local[level][j][km])/2*
						grid->dS1Prime(j,k)/grid->dl1(j,k)/grid->get_halfCellVolumes()[j][k];
					resCjk[SOUTH] = (epsi_local[level][j][km]+epsi_local[level][jm][km])/2*
						grid->dS2Prime(j,km)/grid->dl2(j,km)/grid->get_halfCellVolumes()[j][k];
					resCjk[WEST] = (epsi_local[level][jm][k]+epsi_local[level][jm][km])/2*
						grid->dS1Prime(jm,k)/grid->dl1(jm,k)/grid->get_halfCellVolumes()[j][k];
				}
				else if (PeriodicFlagX1){
					if (j==0){
						jm = J-1;
						jp = 0;
					}
					if (j==J){
						jm = J-1;
						jp = 0;
					}
					if (k==0){
						resCjk[NORTH] = (epsi_local[level][jp][k]+epsi_local[level][jm][k])/2*
							grid->dS2Prime(jp,k)/grid->dl2(jp,k)/grid->get_halfCellVolumes()[jp][k];
						resCjk[SOUTH] = 0.0;
					}
					else if (k==K){
						resCjk[NORTH] = 0.0;
						resCjk[SOUTH] = (epsi_local[level][jp][km]+epsi_local[level][jm][km])/2*
							grid->dS2Prime(jp,km)/grid->dl2(jp,km)/grid->get_halfCellVolumes()[jp][k];
					}
					else{
						resCjk[NORTH] = (epsi_local[level][jp][k]+epsi_local[level][jm][k])/2*
							grid->dS2Prime(jp,k)/grid->dl2(jp,k)/grid->get_halfCellVolumes()[jp][k]; 
						resCjk[SOUTH] = (epsi_local[level][jp][km]+epsi_local[level][jm][km])/2*
							grid->dS2Prime(jp,km)/grid->dl2(jp,km)/grid->get_halfCellVolumes()[jp][k];
					}

					if (k==K){
						resCjk[WEST] = epsi_local[level][jm][km]*
							grid->dS1Prime(jm,k)/grid->dl1(jm,k)/grid->get_halfCellVolumes()[jp][k];
						resCjk[EAST] = epsi_local[level][jp][km]*
							grid->dS1Prime(jp,k)/grid->dl1(jp,k)/grid->get_halfCellVolumes()[jp][k]; 
					}
					else{ 
						resCjk[WEST] = (epsi_local[level][jm][k]+epsi_local[level][jm][km])/2*
							grid->dS1Prime(jm,k)/grid->dl1(jm,k)/grid->get_halfCellVolumes()[jp][k];
						resCjk[EAST] = (epsi_local[level][jp][k]+epsi_local[level][jp][km])/2*
							grid->dS1Prime(jp,k)/grid->dl1(jp,k)/grid->get_halfCellVolumes()[jp][k];
					}
				}
				else if (PeriodicFlagX2){
					if (k==0) 
						km = K-1;
					if (j==0){
						resCjk[EAST] = (epsi_local[level][j][k]+epsi_local[level][j][km])/2*
							grid->dS1Prime(j,k)/grid->dl1(j,k)/grid->get_halfCellVolumes()[j][k];
						resCjk[WEST] = 0.0;
					}
					else if (j==J){
						resCjk[EAST] = 0.0;
						resCjk[WEST] = (epsi_local[level][jm][k]+epsi_local[level][jm][km])/2*
							grid->dS1Prime(jm,k)/grid->dl1(jm,k)/grid->get_halfCellVolumes()[j][k];
					}
					else {
						resCjk[EAST] = (epsi_local[level][j][k]+epsi_local[level][j][km])/2*
							grid->dS1Prime(j,k)/grid->dl1(j,k)/grid->get_halfCellVolumes()[j][k];
						resCjk[WEST] = (epsi_local[level][jm][k]+epsi_local[level][jm][km])/2*
							grid->dS1Prime(jm,k)/grid->dl1(jm,k)/grid->get_halfCellVolumes()[j][k];
					}

					if (j==0){
						resCjk[NORTH] = epsi_local[level][j][k]*
							grid->dS2Prime(j,k)/grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
						resCjk[SOUTH] = epsi_local[level][j][km]*
							grid->dS2Prime(j,km)/grid->dl2(j,km)/grid->get_halfCellVolumes()[j][k];
					}
					else if (j==J){
						resCjk[NORTH] = epsi_local[level][jm][k]*
							grid->dS2Prime(j,k)/grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
						resCjk[SOUTH] = epsi_local[level][jm][km]*
							grid->dS2Prime(j,km)/grid->dl2(j,km)/grid->get_halfCellVolumes()[j][k];
					}
					else {
						resCjk[NORTH] = (epsi_local[level][j][k]+epsi_local[level][jm][k])/2*
							grid->dS2Prime(j,k)/grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
						resCjk[SOUTH] = (epsi_local[level][j][km]+epsi_local[level][jm][km])/2*
							grid->dS2Prime(j,km)/grid->dl2(j,km)/grid->get_halfCellVolumes()[j][k];
					}
				}

				resCjk[SOURCE] = -resCjk[NORTH]- resCjk[EAST]
					-resCjk[SOUTH] - resCjk[WEST];	
				iCoeff = -1/resCjk[SOURCE];
				GSRBCjk[NORTH] = resCjk[NORTH]*iCoeff;
				GSRBCjk[EAST] = resCjk[EAST]*iCoeff;
				GSRBCjk[SOUTH] = resCjk[SOUTH]*iCoeff;
				GSRBCjk[WEST] = resCjk[WEST]*iCoeff;
				GSRBCjk[SOURCE] = -iCoeff;

				break;
			}
		case DIELECTRIC_BOUNDARY:
		case CYLINDRICAL_AXIS:	
			{
				if ((j==J)&&(k==K)){
					resCjk[NORTH] = 0.0;
					resCjk[EAST] = 0.0;
					resCjk[SOUTH] = epsi_local[level][j-1][k-1]*grid->dS2Prime(j,k-1)/
						grid->dl2(j,k-1)/grid->get_halfCellVolumes()[j][k];
					resCjk[WEST] = epsi_local[level][j-1][k-1]*grid->dS1Prime(j-1,k)/
						grid->dl1(j-1,k)/grid->get_halfCellVolumes()[j][k];
				}				
				else if ((j==0)&&(k==0)){
					resCjk[NORTH] = epsi_local[level][j][k]*grid->dS2Prime(j,k)/
						grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
					resCjk[EAST] = epsi_local[level][j][k]*grid->dS1Prime(j,k)/
						grid->dl1(j,k)/grid->get_halfCellVolumes()[j][k];
					resCjk[SOUTH] = 0.0;
					resCjk[WEST] = 0.0;
				}				
				else if ((j==0)&&(k==K)){
					resCjk[NORTH] = 0.0;
					resCjk[EAST] = epsi_local[level][j][k-1]*grid->dS1Prime(j,k)/
						grid->dl1(j,k)/grid->get_halfCellVolumes()[j][k];
					resCjk[SOUTH] = epsi_local[level][j][k-1]*grid->dS2Prime(j,k-1)/
						grid->dl2(j,k-1)/grid->get_halfCellVolumes()[j][k];
					resCjk[WEST] = 0.0;
				}
				else if ((j==J)&&(k==0)){
					resCjk[NORTH] = epsi_local[level][j-1][k]*grid->dS2Prime(j,k)/
						grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
					resCjk[EAST] = 0.0;
					resCjk[SOUTH] = 0.0;
					resCjk[WEST] = epsi_local[level][j-1][k]*grid->dS1Prime(j-1,k)/
						grid->dl1(j-1,k)/grid->get_halfCellVolumes()[j][k];
				}
				else if (k==K){
					resCjk[NORTH] = 0.0;
					resCjk[EAST] = epsi_local[level][j][k-1]*grid->dS1Prime(j,k)/
						grid->dl1(j,k)/grid->get_halfCellVolumes()[j][k];
					resCjk[SOUTH] = (epsi_local[level][j-1][k-1]+epsi_local[level][j][k-1])/2*
						grid->dS2Prime(j,k-1)/
							grid->dl2(j,k-1)/grid->get_halfCellVolumes()[j][k];
					resCjk[WEST] = epsi_local[level][j-1][k-1]*grid->dS1Prime(j-1,k)/
						grid->dl1(j-1,k)/grid->get_halfCellVolumes()[j][k];
				}
				else if (k==0){
					resCjk[NORTH] = (epsi_local[level][j][k]+epsi_local[level][j-1][k])/2*
						grid->dS2Prime(j,k)/
							grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
					resCjk[EAST] = epsi_local[level][j][k]*grid->dS1Prime(j,k)/
						grid->dl1(j,k)/grid->get_halfCellVolumes()[j][k];
					resCjk[SOUTH] = 0.0;
					resCjk[WEST] = epsi_local[level][j-1][k]*grid->dS1Prime(j-1,k)/
						grid->dl1(j-1,k)/grid->get_halfCellVolumes()[j][k];
				}				
				else if (j==0){
					resCjk[NORTH] = epsi_local[level][j][k]*
						grid->dS2Prime(j,k)/
							grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
					resCjk[EAST] = (epsi_local[level][j][k]+epsi_local[level][j][k-1])/2*
						grid->dS1Prime(j,k)/
							grid->dl1(j,k)/grid->get_halfCellVolumes()[j][k];
					resCjk[SOUTH] = epsi_local[level][j][k-1]*grid->dS2Prime(j,k-1)/
						grid->dl2(j,k-1)/grid->get_halfCellVolumes()[j][k];
					resCjk[WEST] = 0.0;
				}
				else if (j==J){
					resCjk[NORTH] = epsi_local[level][j-1][k]*
						grid->dS2Prime(j,k)/
							grid->dl2(j,k)/grid->get_halfCellVolumes()[j][k]; 
					resCjk[EAST] = 0.0;
					resCjk[SOUTH] = epsi_local[level][j-1][k-1]*grid->dS2Prime(j,k-1)/
						grid->dl2(j,k-1)/grid->get_halfCellVolumes()[j][k];
					resCjk[WEST] = (epsi_local[level][j-1][k]+epsi_local[level][j-1][k])/2*
						grid->dS1Prime(j-1,k)/
							grid->dl1(j-1,k)/grid->get_halfCellVolumes()[j][k];
				}
				else 
					cout << "MultiGrid edge error" << endl;
				

				resCjk[SOURCE] = -resCjk[NORTH]- resCjk[EAST]
					-resCjk[SOUTH] - resCjk[WEST];
				
				iCoeff = -1/resCjk[SOURCE];
				GSRBCjk[NORTH] = resCjk[NORTH]*iCoeff;
				GSRBCjk[EAST] = resCjk[EAST]*iCoeff;
				GSRBCjk[SOUTH] = resCjk[SOUTH]*iCoeff;
				GSRBCjk[WEST] = resCjk[WEST]*iCoeff;
				GSRBCjk[SOURCE] = -iCoeff;
		
				break;
			}
		default:
			{
//			  cout << "Multigrid doesn't know about this boundary condition type!\n;" << endl;
//			  cout << "(was it  SPATIAL_REGION_BOUNDARY?)\n" << endl;
//			  e_xit(1);
			}
			
		}
}

Vector2** Multigrid::GridCoarsen(Grid* grid, int ratio1, int ratio2)
{	
	Vector2** X;
	int j;
	X = new Vector2*[length[level].e1()+1];
	for (j=0; j<=length[level].e1(); j++){
		X[j] = new Vector2[length[level].e2()+1];
		for (int k=0; k<=length[level].e2(); k++)
			X[j][k] = grid->getMKS(ratio1*j,ratio2*k);
	}
	
	return X;
}

Scalar Multigrid::Resid(Scalar**rho, Scalar**phi)
{
	Scalar norm;
	for (int i = 0; i <= length[0].e1(); i++) 
		for (int j = 0; j <= length[0].e2(); j++) {
			if (GSRBCoeff[0][i][j][SOURCE]==0.0)
				soln[0][i][j] = phi[i][j];
			else
				soln[0][i][j] = 0.0;
			rhs[0][i][j] = -rho[i][j];
		}
	Scalar error0 = errorest();
	Residual(phi, rhs[0], res[0],length[0],resCoeff[0], PeriodicFlagX1, PeriodicFlagX2);
	norm = Norm(res[0],length[0],MultiGrid[0]);
	norm /= error0;
	return norm;
}

void Multigrid::PSolveBoltzCoeff(const Scalar& n0, const Scalar& qbyT, Scalar** brho, const Scalar& MinPot)
{

	Scalar temp;
	Scalar iCoeff;
	int J,K;
	Scalar* resCjk, *GSRBCjk;

	for (int i = 0; i <= length[0].e1(); i++) 
		for (int j = 0; j <= length[0].e2(); j++) 
			rhs[0][i][j] = brho[i][j];
		

	for (level=0; level < (tree_size-1); level++)
		Average(rhs[level], rhs[level+1], length[level], length[level+1], 
						PeriodicFlagX1, PeriodicFlagX2);

	for (level=0; level < tree_size; level++){
		J = length[level].e1();
		K = length[level].e2();
		for(int j=0;j<=J; j++)
			for(int k=0;k<=K; k++){
				GSRBCjk = GSRBCoeff[level][j][k];
				if (GSRBCjk[SOURCE])
					if(MultiGrid[level]->PlasmaRegion(j,k)){
						resCjk = resCoeff[level][j][k];
						temp = -rhs[level][j][k];
						if (n0)
							resCjk[SOURCE] = -resCjk[NORTH]- resCjk[EAST]
								-resCjk[SOUTH] - resCjk[WEST] + temp*qbyT;
						else
							resCjk[SOURCE] = -resCjk[NORTH]- resCjk[EAST]
								-resCjk[SOUTH] - resCjk[WEST];
						
						//this code chould be sped up a lot resCoeff_local is 
						//stored in the boltzmann object.

						iCoeff = -1/resCjk[SOURCE];
						GSRBCjk[NORTH] = resCjk[NORTH]*iCoeff;
						GSRBCjk[EAST] = resCjk[EAST]*iCoeff;
						GSRBCjk[SOUTH] = resCjk[SOUTH]*iCoeff;
						GSRBCjk[WEST] = resCjk[WEST]*iCoeff;
						GSRBCjk[SOURCE] = -iCoeff;
					}
			}
	}
	return;
}