dadirz.cpp

pic 模拟程序！面向对象

	 case DIELECTRIC_BOUNDARY:
		{
			if(i==0&&b_x2geom[i][j]!=0) {  //Don't override a dirichlet with a Neumann
					Scalar xa = X[i][j].e1();
					Scalar xave = 0.5*(2*xa + dxib);

					a_x1geom[i][j] = 0.0;
					c_x1geom[i][j] = e1b/(xave-xa)/dxib;
					b_x1geom[i][j] = -( a_x1geom[i][j] + c_x1geom[i][j] );
					if(j==0) {
						Scalar dria = X[i][1].e2() - X[i][0].e2();
						a_x2geom[i][0] = 0.0;
						b_x2geom[i][0] = -e2b/dria/dria;
						c_x2geom[i][0] = -b_x2geom[i][0];
					}
					if(j==J) { //neumann neumann corner
						c_x2geom[i][j] = 0.0;
						a_x2geom[i][j] = e2a*2*ra2/(dr2*dra);
						b_x2geom[i][j] = -a_x2geom[i][j];
					}

			 }
			if(i==I&&b_x2geom[i][j]!=0) {
					Scalar xb = X[i][j].e1();
					Scalar xave = 0.5*(X[i][j].e1() + X[i-1][j].e1());
					
					a_x1geom[i][j] = e1a/(xb-xave)/dxia;
					c_x1geom[i][j] = 0.0;
					b_x1geom[i][j] = -( a_x1geom[i][j] + c_x1geom[i][j] );
				if(j==0&&b_x2geom[i][j]!=0) { //neumann neumann corner
				  Scalar dria = X[i][j+1].e2() - X[i][j].e2();
				  a_x2geom[i][j] = 0.0;
				  c_x2geom[i][j] = e2b/dria/dria;
				  b_x2geom[i][j] = -c_x2geom[i][j];
				}
				if(j==J&&b_x2geom[i][j]!=0) { //neumann neumann corner
				  c_x2geom[i][j] = 0.0;
				  a_x2geom[i][j] = e2a*2*ra2/(dr2*dra);
				  b_x2geom[i][j] = -a_x2geom[i][j];
				}
			}
			
			if(j==J && i!=0&&b_x2geom[i][j]!=0 ) { 
				Scalar dria = X[i][j].e2()-X[i][j-1].e2();

				a_x2geom[i][nc2] = 2 * e2a * ra2/ ( dria * ( r*r - ra2*ra2));
				b_x2geom[i][nc2] = -a_x2geom[i][nc2];
				c_x2geom[i][nc2] = 0.0;
			}
			break;
		}
	 case CYLINDRICAL_AXIS:
		{ if(i==0||i==I) break;  // don't want cyl_axis in corners
		  Scalar dria = X[i][1].e2() - X[i][0].e2();
		  a_x2geom[i][0] = 0.0;
		  b_x2geom[i][0] = -e2b/dria/dria;
		  c_x2geom[i][0] = -b_x2geom[i][0];

		  a_x1geom[i][j] = e1a/(dxia*dxave);
		  c_x1geom[i][j] = e1b/(dxib*dxave);
		  b_x1geom[i][j] = -( a_x1geom[i][j] + c_x1geom[i][j] );

		  break;
		}
	 default: 
		{
		  cout << "DADI doesn't know about this boundary condition type!\n;" << endl;
		  cout << "(was either PERIODIC_BOUNDARY or SPATIAL_REGION_BOUNDARY)\n" << endl;
		  //e_xit(1);
		}
	 }
 }

BCTypes Dadirz::get_coefficient(int j, int k) {
  return (b_x2geom[j][k]==0)? CONDUCTING_BOUNDARY:FREESPACE;
}

/**********************************************************************/

int Dadirz::solve(Scalar **u_in, Scalar **s, int itermax,Scalar tol_test)
{
  register int i, j;
  int iter, ndiscard;
  static Scalar del_t=0.0;
  Scalar del_td=0, tptop=0, tpbot=0, ratio=0;
  Scalar rnorm=0, rsum=0, res=0, errchk=0, dxdxutrm=0, dydyutrm=0;

  rnorm=rsum=0.0;
  for(i=0; i< ng1; i++)
    for(j=0; j< ng2; j++) {

		/* Residual normalization.  */
		// dirichlet conditions don't add to rnorm
		rnorm +=( (b_x2geom[i][j]==0)? 0:s[i][j]*s[i][j]) ;
      
  /*  copy u_in to u for working purposes.  */  
		u[i][j]=(Scalar)u_in[i][j];

		//calculate an initial estimate of the residual
	/* Use the residual as the absolute error and if it is bigger 
	   than the stored maximum absolute error, record new maximum 
	   absolute error and location.  */

		if(i>0&&j>0&&i<nc1&&j<nc2) {
		  dxdxutrm=a_x1geom[i][j]*u_in[i-1][j] +b_x1geom[i][j]*u_in[i][j] +
			 c_x1geom[i][j]*u_in[i+1][j];
		  dydyutrm=a_x2geom[i][j]*u_in[i][j-1] +b_x2geom[i][j]*u_in[i][j] +
			 c_x2geom[i][j]*u_in[i][j+1];
		}
		  /* only include points which are not in structures. */
		  errchk = ((b_x2geom[i][j]==0)? 0:dxdxutrm +dydyutrm -s[i][j]);
	
	/* Residual sums. */
		  rsum += errchk*errchk;
	 }

  // If rnorm is zero, we must deal with it...
  if (rnorm == 0.0) {

	 // check dirichlet conditions
	 for(i=0;i<ng1;i++) for(j=0;j<ng2;j++) {
		// check left
		rnorm += sqr(((i>0&&b_x2geom[i-1][j]!=0)?c_x1geom[i-1][j]*u[i][j]:0));
		// check right
		rnorm += sqr(((i<nc1&&b_x2geom[i+1][j]!=0)?a_x1geom[i+1][j]*u[i][j]:0));
		// check up
		rnorm += sqr(((j>0&&b_x2geom[i][j-1]!=0)?c_x2geom[i][j-1]*u[i][j]:0));
		// check right
		rnorm += sqr(((j<nc2&&b_x2geom[i][j+1]!=0)?c_x2geom[i][j+1]*u[i][j]:0));

	 }
	 //rnorm = sqrt(rnorm);
	 if(rnorm == 0) { //still zero, we don't need to iterate
		for(i=0;i<ng1;i++) for(j=0;j<ng2;j++) u_in[i][j]=0;
		return 0;
	 }
  }
  rnorm=sqrt(rnorm);
  res = sqrt(rsum)/rnorm;
#ifdef DADI_DEBUG
  printf("dadi: res = %g\n",res);
  printf("dadi: rnorm= %lg\n", rnorm);
#endif

  if(res<tol_test) return 0;  // no need to iterate

  /*************************************************/
  if(del_t<1.0e-10*del_t0) del_t=del_t0; else del_t/=4;
  del_td= 2.0*del_t;
  ndiscard=0;

  /********************/
  /* Begin iteration. */
  
  for(iter=0; iter<itermax; iter++)
  {
    /*************************************************/
    /* Copy u into the work array and storage array. */
    
    for(i=0; i< ng1; i++)
      for(j=0; j< ng2; j++) uwork[i][j] = ustor[i][j] = u[i][j];
    
    /************************************/
    /* Two advances of u via ADI at del_t. */
    
    adi(u, s, del_t);
    adi(u, s, del_t);
    
    /*****************************************/
    /* One advance of uwork via ADI at 2*del_t. */

    adi(uwork, s, del_td);
    
    /*******************************************************/
    /* Calculate test parameter and normalized error.
       For Dirichlet BCs, no need to worry about boundary
       points since u,uwork, and ustor should be the same. */
    
    tptop= tpbot= rsum = 0.0;

    for(i=1; i< nc1; i++)
      for(j=1; j< nc2; j++)
      {
		  /* Test paramter sums. */
		  tptop += (u[i][j]-uwork[i][j])*(u[i][j]-uwork[i][j]);
		  tpbot += (u[i][j]-ustor[i][j])*(u[i][j]-ustor[i][j]);
	
	/* Residual terms. */

		  dxdxutrm=a_x1geom[i][j]*u[i-1][j] +b_x1geom[i][j]*u[i][j] +c_x1geom[i][j]*u[i+1][j];
		  dydyutrm=a_x2geom[i][j]*u[i][j-1] +b_x2geom[i][j]*u[i][j] +c_x2geom[i][j]*u[i][j+1];

	/* Use the residual as the absolute error and if it is bigger 
	   than the stored maximum absolute error, record new maximum 
	   absolute error and location.  */
	/*  only include points which are outside of structures. */
		  errchk = ((b_x2geom[i][j]==0)? 0:dxdxutrm +dydyutrm -s[i][j]);
	
	/* Residual sums. */
		  rsum += errchk*errchk;
      }

    /* Calculate normalized residual. */
    res = sqrt(rsum)/rnorm;
#ifdef DADI_DEBUG
    printf("dadi: iter= %d, res = %lg del_t=%lf\n", iter, res,del_t);
#endif
    /* If the residual is less than the tolerance, SUCCESS! */
    if ((res < tol_test) && (iter))
		{
#ifdef DADI_DEBUG
		  printf("dadi: iter=%d\n", iter);
#endif
		  for(i=0;i<ng1;i++)
			 for(j=0;j<ng2;j++)
				u_in[i][j]=(Scalar)u[i][j];

		  return(0);
		}

    /* Determine ratio used to find the time step change.  If tpbot 
       is zero but tptop is finite, consider this a case of a large 
       ratio and act accordingly.  DWH does about the same thing 
       except he does NOT discard the solution if tpbot=0. */
    
    if (tpbot > 0.0) ratio=tptop/tpbot;
    if (tpbot== 0.0) ratio=1.0;
    
    /* Get next time step. */
    if (ratio < 0.02) del_t *= 8.0;
    if (ratio >=0.02 && ratio < 0.05) del_t *= 4.0;
    if (ratio >=0.05 && ratio < 0.10) del_t *= 2.0;
    if (ratio >=0.10 && ratio < 0.30) del_t *= 0.80;
    if (ratio >=0.30 && ratio < 0.40) del_t *= 0.50;
    if (ratio >=0.40 && ratio < 0.60) del_t *= 0.25;
	 for(i=0;i<ng1;i++)
		for(j=0;j<ng2;j++)
		  u_in[i][j]=(Scalar)u[i][j];
   
    
    /* Ratio is too large. */
    if (ratio >=0.60)
    { 
      ndiscard++;
#ifdef DADI_DEBUG
      printf("ndiscard= %d, iter=%d step=%lf\n", ndiscard, iter,del_t); 
#endif      
      /* Check if too many discards. */
      if (ndiscard > 20)
      {
		  for(i=0;i<ng1;i++)
			 for(j=0;j<ng2;j++)
				u_in[i][j]=(Scalar)u[i][j];
		  del_t = del_t0;
		  printf("Poisson solve FAILURE: dadi: iter= %d, ndiscard>10\n", iter);
		  return 1;
      }	
      /* Discard by replacing u with what we started with. */
      for(i=0; i< ng1; i++)
		  for(j=0; j< ng2; j++) u[i][j] = ustor[i][j];
      
      /* Reduce del_t. */
      del_t /= 8.0;
		//del_t = del_t0;
    }
    del_td=2*del_t;
  }
  /* Fail if used up the maximum iterations. */
  printf("Poisson solve FAILURE:dadi:  iter>= %d\n",itermax);
  for(i=0;i<ng1;i++)
	 for(j=0;j<ng2;j++)
		u_in[i][j]=(Scalar)u[i][j];
 
  return(2);
}



Dadirz::~Dadirz()
{

}