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;}