dadixy.cpp

pic 模拟程序！面向对象

							c_x2geom[i][j] = e2b/dyave/dyb;
							b_x2geom[i][j] = - ( a_x2geom[i][j] + c_x2geom[i][j] );
							return;
						}
				}

				if(i==I&&b_x2geom[i][j]!=0) {  //right hand wall, not a conductor
					dxa = X[i][j].e1()-X[i-1][j].e1();
					a_x1geom[i][j] = 2*e1a/dxa/dxa;
					c_x1geom[i][j] = 0.0;
					b_x1geom[i][j] = -( a_x1geom[i][j] + c_x1geom[i][j] );
					if(j==0) { // neumann neumann corner
						a_x2geom[i][j] = 0.0;
						c_x2geom[i][j] = e2b/dyb/dyb;
						b_x2geom[i][j] = -c_x2geom[i][j];
						return;
					} else
					if(j==J) { //neumann neumann corner
						dya = X[i][j].e2()-X[i][j-1].e2();
						c_x2geom[i][j] = 0.0;
						a_x2geom[i][j] = e1b/dya/dya;
						b_x2geom[i][j] = -a_x2geom[i][j];
						return;
					}
					else	{
						a_x2geom[i][j] = e2a/dyave/dya;
						c_x2geom[i][j] = e2b/dyave/dyb;
						b_x2geom[i][j] = - ( a_x2geom[i][j] + c_x2geom[i][j] );
						return;
					}
				}
				if(j==J&&b_x2geom[i][j]!=0) {  //we already dealt with corners, no need to here
					c_x2geom[i][j] = 0.0;
					a_x2geom[i][j] = e1b/dya/dya;
					b_x2geom[i][j] = -a_x2geom[i][j];
				}

				if(j==0&&b_x2geom[i][j]!=0) {	// already dealt with corners, no need to here
					a_x2geom[i][j] = 0.0;
					c_x2geom[i][j] = e2b/dyb/dyb;
					b_x2geom[i][j] = -c_x2geom[i][j];
				}
				break;
			}
		 case CYLINDRICAL_AXIS:
			{
				//meaningless in X-Y geometry

				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 Dadixy::get_coefficient(int j, int k) {
  return (b_x2geom[j][k]==0)? CONDUCTING_BOUNDARY:FREESPACE;
}


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

/*dadi(u_in, s, itermax, tol_test, u_x0, u_xlx, u_y0, u_yly)
int itermax;
Scalar tol_test;
Scalar **u_in, **s, *u_x0, *u_xlx, *u_y0, *u_yly;
*/
int Dadixy::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++) {
		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));

	 }

	 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= %g\n", rnorm);
#endif

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

  /*************************************************/
  if(del_t==0.0) 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 not in 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=%le\n", iter, res,del_t);
		printf("dadi: iter= %d, res = %g del_t=%e\n", iter, res,del_t);
#endif // DADI_DEBUG
    /* If the residual is less than the tolerance, SUCCESS! */
    if ((res < tol_test) && (iter))
		{
#ifdef DADI_DEBUG
		  printf("dadi: iter=%d\n", iter);
#endif// DADI_DEBUG
		  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;
#ifndef NO_STEP_ADJUST    
    /* 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];
#endif   
    
	 /* Ratio is too large. */
    if (ratio >=0.60)
    { 
      ndiscard++;
			iter--;
#ifdef DADI_DEBUG
	    //  printf("ndiscard= %d, iter=%d step=%lf\n", ndiscard, iter,del_t); 
			printf("ndiscard= %d, iter=%d step=%f\n", ndiscard, iter,del_t); 
#endif //DADI_DEBUG
      
      /* 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;
//		  if(solve(u_in,s,itermax,tol_test))
		  printf("Poisson solve FAILURE: dadi: iter= %d, ndiscard>20\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);
}


Dadixy::~Dadixy()
{
}