spbound.cpp

pic 模拟程序！面向对象

#ifdef MPI_VERSION
  cout << "Received stripe:" << endl;
  cout << "E1 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lE[kk].e1();
  cout << endl;
  cout << "E2 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lE[kk].e2();
  cout << endl;
  cout << "E3 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lE[kk].e3();
  cout << endl;
  cout << "B1 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lB[kk].e1();
  cout << endl;
  cout << "B2 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lB[kk].e2();
  cout << endl;
  cout << "B3 =";
  for(int kk=0; kk<=k2-k1; kk++) cout << " " << lB[kk].e3();
  cout << endl;
#endif
*/

  switch(dir){

// Fields coming from right (right moving window)
    case 1:{
      assert(j1 == j2);
// On KC's authority, one always has k1 <= k2
      for(int k=k1; k<=k2; k++) {
        intEdl[j1-1][k].set_e1(lE[k-k1].e1());
        intEdl[j1][k].set_e2(lE[k-k1].e2());
        intEdl[j1][k].set_e3(lE[k-k1].e3());
        intBdS[j1][k].set_e1(lB[k-k1].e1());
        intBdS[j1-1][k].set_e2(lB[k-k1].e2());
        intBdS[j1-1][k].set_e3(lB[k-k1].e3());
        I[j1-1][k].set_e1(lJ[k-k1].e1());
        I[j1][k].set_e2(lJ[k-k1].e2());
        I[j1][k].set_e3(lJ[k-k1].e3());
      }
    }
    break;

// Fields coming from left (left moving window)
    case 2:{
      assert(j1 == j2);
// On KC's authority, one always has k1 <= k2
      for(int k=k1; k<=k2; k++) {
        intEdl[j1][k] = lE[k-k1];
        intBdS[j1][k] = lB[k-k1];
        I[j1][k] = lJ[k-k1];
      }
    }
    break;

// Fields coming from bottom (from greater k, downward moving window)
    case 3:{
      assert(k1 == k2);
// On KC's authority, one always has k1 <= k2
      for(int j=j1; j<=j2; j++) {
        intEdl[j][k1].set_e1(lE[j-j1].e1());
        intEdl[j][k1-1].set_e2(lE[j-j1].e2());
        intEdl[j][k1].set_e3(lE[j-j1].e3());
        intBdS[j][k1-1].set_e1(lB[j-j1].e1());
        intBdS[j][k1].set_e2(lB[j-j1].e2());
        intBdS[j][k1-1].set_e3(lB[j-j1].e3());
        I[j][k1].set_e1(lJ[j-j1].e1());
        I[j][k1-1].set_e2(lJ[j-j1].e2());
        I[j][k1].set_e3(lJ[j-j1].e3());
      }
    }
    break;

// Fields coming from top (from smaller k, upward moving window)
    case 4:{
      assert(k1 == k2);
// On KC's authority, one always has k1 <= k2
      for(int j=j1; j<=j2; j++) {
        intEdl[j][k1] = lE[j-j1];
        intBdS[j][k1] = lB[j-j1];
        I[j][k1] = lJ[j-j1];
      }
    }
    break;

  }

}

void SpatialRegionBoundary::SRBadvanceB(Scalar t, Scalar dt) {
   Scalar dt2 = dt/2;
   Vector3 **intBdS = fields->getIntBdS();
   Vector3 **intEdl = fields->getIntEdl();
   Vector3 **I = fields->getI();
   Vector3 **iC = fields->getiC();
   Vector3 **iL = fields->getiL();

   int kl = MIN(k1,k2);
   int kh = MAX(k1,k2);
   int jl = MIN(j1,j2);
   int jh = MAX(j1,j2);
   int K = fields->getK();
   int J = fields->getJ();

   //  Advance B 1/2 step
   if(k1==k2) {  //horizontal boundary 
      for(int j=jl;j<jh;j++) {
	 //The regular field solve handles B2[j][k]
	 if(normal == 1) { //  the virtual region is below us
	    lB[j-jl] -= dt2 * Vector3 ( intEdl[j][k1].e3() - lE[j-jl].e3(),
				       0,  //handled by regular field solve
				       lE[j-jl].e1() - intEdl[j][k1].e1() +
				       lE[j-jl+1].e2() - lE[j-jl].e2());
	 }
	 else {    // the virtual region is above us
	    lB[j-jl] -= dt2 * Vector3 ( -intEdl[j][k1].e3() + lE[j-jl].e3(),
				       0,
				       -lE[j-jl].e1() + intEdl[j][k1].e1() +
				       lE[j-jl+1].e2() - lE[j-jl].e2());
	 }
      }
   }
   else {  //vertical boundary
      for(int k=kl;k<kh;k++) {
	 int kk = k-kl;
	 // write the local E fields to lEa.  The other side is not "authoritative"
	 // and if we don't, we get an instability!!
	 lEa[kk].set_e2(intEdl[j1][k].e2());
	 lEa[kk].set_e3(intEdl[j1][k].e3());
	 //The regular field solve handles B1[j][k]
	 if(normal == 1) {  // lB handles the virtual fields on the left
	    /* new way of doing this:  use the LOCAL E's where possible. */
	    lBa[kk] = intBdS[j1][k];  // These are the authoratitive fields

	    lBb[kk] -= dt2 * Vector3( NEW_B1( lEb[kk], lEb[kk+1]),
				     NEW_B2( lEb[kk], intEdl[j1][k]),
				     NEW_B3( lEb[kk], intEdl[j1][k], lEb[kk+1]));
	    
	    lBc[kk] -= dt2 * Vector3( NEW_B1( lEc[kk], lEc[kk+1]),
				     NEW_B2( lEc[kk], lEb[kk]),
				     NEW_B3( lEc[kk], lEb[kk],lEc[kk+1]));
	 }
	 else {  //lB handles tp *intEdl[4]@4he virtual fields on the right
	    /* new way of doing this:  we use the OTHER side's version of E3, E2 to make sure we're consistent */
	    lBa[kk] -= dt2 * Vector3 ( NEW_B1( lEa[kk], lEa[kk+1]), //These have correct  E3, E2, not E1   
				      NEW_B2( lEa[kk], lEb[kk]),  
				      NEW_B3( lEa[kk], lEb[kk], lEa[kk+1]));

	    lBb[kk] -= dt2 * Vector3 ( NEW_B1( lEb[kk],lEb[kk+1]),
				      NEW_B2( lEb[kk], lEc[kk]),
				      NEW_B3( lEb[kk],lEc[kk],lEb[kk+1]));

	    lBc[kk] -= dt2 * Vector3 ( NEW_B1( lEc[kk], lEc[kk+1]),
				      0, // cannot calculate:  insufficient info
				      0); // cannot calculate:  insufficient info
	    // no update of intBds necessary, lB's are all we need, lB1 done properly by reg.
	    // field solve.
	 }
      }
   }
}


void SpatialRegionBoundary::SRBadvanceE(Scalar t, Scalar dt) {
   Scalar dt2 = dt/2;
   Vector3 **intBdS = fields->getIntBdS();
   Vector3 **intEdl = fields->getIntEdl();
   Vector3 **I = fields->getI();
   Vector3 **iC = fields->getiC();
   Vector3 **iL = fields->getiL();

   int kl = MIN(k1,k2);
   int kh = MAX(k1,k2);
   int jl = MIN(j1,j2);
   int jh = MAX(j1,j2);
   int K = fields->getK();
   int J = fields->getJ();

   if(k1==k2) {  //horizontal boundary 
      for(int j=jl;j<jh;j++) {
	 
	 if(normal==1) {  //real region is below
	    //first update the current
	    I[j][k1]+=Vector3(lJ[j-jl].e1(),0,lJ[j-jl].e3());
	    

	    //update the E fields
	    intEdl[j][k1]+=dt * iC[j][k1-1].jvMult(Vector3(
							   ((j<J)?lB[j-jl].e3() * liL[j-jl].e3() - intBdS[j][k1-1].e3() * iL[j][k1-1].e3():0),
							   0,  //This field isn't my problem here, fields takes care of it
							   +intBdS[j][k1-1].e1() * iL[j][k1-1].e1() - lB[j-jl].e1() * liL[j-jl].e1()
							   +((j<J)?intBdS[j][k1].e2() * iL[j][k1].e2():0) - ((j>0)?intBdS[j-1][k1].e2()*iL[j-1][k1].e2():0))
						   -Vector3(I[j][k1].e1(),0,I[j][k1].e3()));

	    // need to update this component for to-nodes.
	    //		  lE[k-kl]+= dt * iC[j1+1][k].jvMult(Vector3(lB[k-kl].e3()*iL[j1+1][k].e3()
	    //						-((k>0)?lB[k-kl-1].e3()*iL[j1+1][k-1].e3():0) - lJ[k-kl].e1(),0,0));
	    
	 }
	 else {  // real region is above.
	    //first update the current
	    I[j][k1]+=Vector3(lJ[j-jl].e1(),0,lJ[j-jl].e3());

	    //update the E fields
	    intEdl[j][k1]+=dt * iC[j][k1-1].jvMult(Vector3(
							   -((j<J)?lB[j-jl].e3() * liL[j-jl].e3() + intBdS[j][k1].e3() * iL[j][k1].e3():0),
							   0,  //This field isn't my problem here, fields takes care of it
							   -intBdS[j][k1-1].e1() * iL[j][k1-1].e1() + lB[j-jl].e1() * liL[j-jl].e1()
							   +((j<J)?intBdS[j][k1].e2() * iL[j][k1].e2():0) - ((j>0)?intBdS[j-1][k1].e2()*iL[j-1][k1].e2():0))
						   -Vector3(I[j][k1].e1(),0,I[j][k1].e3()));

	    
	    
	 }
      }
      
   }
   else {  //vertical boundary
      for(int k=kl;k<=kh;k++) {
	 int kk = k - kl;
	 // First, update the current
	 I[j1][k]+=Vector3(0,lJa[kk].e2(),lJa[kk].e3());
	 lJa[kk] = Vector3(lJa[kk].e1(),I[j1][k].e2(),I[j1][k].e3());  // make sure we got the right value in lJa too

	 if(normal==1) { //real region is to the right.
	    // ASSUME:  UNIFORM MESH (assume iL)
	    // ASSUME:  UNIFORM DIELECTRIC CONSTANT (assume iC)
	    // ASSUME:  NO BOUNDARY CONDITION ON OPPOSITE SIDE  (assume iC)
	    // To remove these assumptions, iC and iL must be passed and recalculated.

	    lEa[kk] += dt * Vector3(0,  // Reg. field solve calculates this.
				    NEW_E2( iC[j1+1][k], iL[j1+1][k], iL[j1+1][k], intBdS[j1][k], lBb[kk], I[j1][k]),
				    NEW_E3( iC[j1+1][k], iL[j1+1][k], iL[j1+1][MAX(k-1,0)], iL[j1+1][k], intBdS[j1][k], ((kk>0)? intBdS[j1][k-1]: Vector3(0)),
					   lBb[kk], I[j1][k]));

	    lEb[kk] += dt * Vector3( NEW_E1( iC[j1+1][k], iL[j1+1][k], iL[j1+1][MAX(k-1,0)], lBb[kk], ((kk>0)? lBb[kk-1]: Vector3(0)),
					    lJb[kk]),
				    NEW_E2( iC[j1+1][k], iL[j1+1][k], iL[j1+1][k], lBb[kk], lBc[kk], lJb[kk]),
				    NEW_E3( iC[j1+1][k], iL[j1+1][k], iL[j1+1][MAX(k-1,0)], iL[j1+1][k], lBb[kk], ((kk>0)? lBb[kk-1]: Vector3(0)),
					   lBc[kk], lJb[kk]));

	    lEc[kk] += dt * Vector3( NEW_E1( iC[j1+1][k], iL[j1+1][k], iL[j1+1][MAX(k-1,0)], lBc[kk], ((kk>0)? lBc[kk-1]: Vector3(0)),
					    lJc[kk]),
				    0,  //  Cannot calculateNEW_E2( iC[j1+1][k], iL[j1+1][k], iL[j1+1][k], lBb[kk], lBc[kk], lJb[kk]),
				    0);// cannot calculate NEW_E3( iC[j1+1][k], iL[j1+1][k], iL[j1+1][k], lBb[kk], ((kk>0)? lBb[kk-1]: Vector3(0)),
	    //	 lBc[kk], lJb[kk]));
	 }
	 else {  // real region to the left.
	    // ASSUME:  UNIFORM MESH (assume iL)
	    // ASSUME:  UNIFORM DIELECTRIC CONSTANT (assume iC)
	    // ASSUME:  NO BOUNDARY CONDITION ON OPPOSITE SIDE  (assume iC)
	    // To remove these assumptions, iC and iL must be passed and recalculated.
	    lEa[kk] += dt * Vector3( NEW_E1( iC[j1-1][k], iL[j1-1][k], iL[j1-1][MAX(k-1,0)], lBa[kk], ((kk>0)? lBa[kk-1]: Vector3(0)),
					    lJa[kk]),
				    NEW_E2( iC[j1-1][k], iL[j1-1][k], iL[j1-1][k], lBa[kk], intBdS[j1-1][k], I[j1][k]),
				    NEW_E3( iC[j1-1][k], iL[j1-1][k],iL[j1-1][MAX(k-1,0)], iL[j1-1][k], lBa[kk], ((kk>0)? lBa[kk-1]: Vector3(0)),
					   intBdS[j1-1][k], I[j1][k]));

	    lEb[kk] += dt * Vector3( NEW_E1( iC[j1-1][k], iL[j1-1][k], iL[j1-1][MAX(k-1,0)], lBb[kk], ((kk>0)? lBb[kk-1]: Vector3(0)),
					    lJb[kk]),
				    NEW_E2( iC[j1-1][k], iL[j1-1][k], iL[j1-1][k], lBb[kk], lBa[kk], lJb[kk]),
				    NEW_E3( iC[j1-1][k], iL[j1-1][k],iL[j1-1][MAX(k-1,0)], iL[j1-1][k], lBb[kk], ((kk>0)? lBb[kk-1]: Vector3(0)),
					   lBa[kk], lJb[kk]));

	    lEc[kk] += dt * Vector3( NEW_E1( iC[j1-1][k], iL[j1-1][k], iL[j1-1][MAX(k-1,0)], lBc[kk], ((kk>0)? lBc[kk-1]: Vector3(0)),
					    lJc[kk]),
				    0,  //  Cannot calculateNEW_E2( iC[j1+1][k], iL[j1+1][k], iL[j1+1][k], lBb[kk], lBc[kk], lJb[kk]),
				    0);// cannot calculate NEW_E3( iC[j1+1][k], iL[j1+1][k], iL[j1+1][k], lBb[kk], ((kk>0)? lBb[kk-1]: Vector3(0)),
	    //	 lBc[kk], lJb[kk]));


	    
	 }
	 // write back the updated E fields
	 intEdl[j1][k] = Vector3(intEdl[j1][k].e1(),lEa[kk].e2(),lEa[kk].e3());
      }
   }
}

/**
 * Ask for a stripe of NGDs for this boundary
 */
void SpatialRegionBoundary::askNGDStripe() {
  assert(theLink);
  theLink->askShiftedNGDs();
}
/**
 * Send a stripe of the data structures related to 
 * the Neutral Gas Density  for this boundary
 *
 * Again, only the right moving window is implemented
 * for now: case 1 in the switch statement. 
 *
 * @param dir the direction integer
 */
void SpatialRegionBoundary::sendNGDStripe(int dir) 
  throw(Oops){
 
#ifdef  MPI_VERSION
  extern int MPI_RANK;
#endif

  switch(dir){

    // NGDs sending to left (right moving window)
    case 1:{
      assert(j1 == j2);
      // On KC's authority, one always has k1 <= k2
      // note that for the NGDs we loop to k<k2; dad 05/23/01.
      int cntr = 0;
      NGD* ptrNGD;
      oopicListIter<NGD> NGDIter(*(ptrMapNGDs->getPtrNGDList()));
      for ( NGDIter.restart(); !NGDIter.Done(); NGDIter++ ) {
        ptrNGD = NGDIter.current(); 
        for(int k=k1; k<k2; k++) {
          try{
            lNGDbufSend[cntr++] = ptrNGD->getNGD(j1, k);
            if ( ptrNGD->getisTIOn() ) {
              lNGDbufSend[cntr++] = ptrNGD->getNGD_excessNumIons(j1, k);
              lNGDbufSend[cntr++] = ptrNGD->getNGD_TI_np2c(j1, k);
            }
          }
          catch(Oops& oops){
            oops.prepend("SpatialRegionBoundary::sendNGDStripe: Error: \n"); //OK
            throw oops;
          }
        }
      }
      /*
        #ifdef  MPI_VERSION
        cout << "MPI_RANK = " << MPI_RANK << 
        ", SpatialRegionBoundary::sendStripe sending NGDs." << endl;
        #endif
      */
      try{
        theLink->sendShiftedNGDs( lNGDbufSend, cntr );
      }
      catch(Oops& oops){
        oops.prepend("SpatialRegionBoundary::sendNGDStripe: Error:\n");  //OK
        throw oops;
      }


      /*
        #ifdef  MPI_VERSION
        cout << "MPI_RANK = " << MPI_RANK << 
	", SpatialRegionBoundary::sendStripe sent NGDs." << endl;
        #endif
      */
    }
    break;
  default:{
    stringstream ss (stringstream::in | stringstream::out);
    ss << "SpatialRegionBoundary::sendNGDStripe: Error: \n"<<
        "Only right moving window (dir = 1) " << endl
           << "is implemented for now! You have requested dir = " << dir << endl;
         
      std::string msg;
      ss >> msg;
      Oops oops(msg);
      throw oops;    // exit()MapNGDs::sendNGDsAndShift

    }
    break;
  }
}
/**
 * Receive a stripe of NGDs from this boundary.
 */
void SpatialRegionBoundary::recvNGDStripe(int dir) 
  throw(Oops){

  assert(theLink);
#ifdef  MPI_VERSION
  extern int MPI_RANK;
#endif
  
  theLink->waitShiftedNGDs();  //  do a synchronize

  switch(dir){

    //
    // NGDs coming from right (right moving window)
    // Only this case is impleneted for now. 
    // dad 05/24/01.
    //
    case 1:{
      assert(j1 == j2);
// On KC's authority, one always has k1 <= k2
      int cntr = 0;
      int j = j1 - 1;
      NGD* ptrNGD;
      oopicListIter<NGD> NGDIter(*(ptrMapNGDs->getPtrNGDList()));
      for ( NGDIter.restart(); !NGDIter.Done(); NGDIter++ ) {
        ptrNGD = NGDIter.current(); 
        for(int k = k1; k < k2; k++) {
          ptrNGD->set(j, k, lNGDbuf[cntr++]);
          if ( ptrNGD->getisTIOn() ) {
            try{
              ptrNGD->setNGD_excessNumIons(j, k, lNGDbuf[cntr++]);
              ptrNGD->setNGD_TI_np2c(j, k, lNGDbuf[cntr++]);
            }
            catch(Oops& oops){
              oops.prepend("SpatialRegionBoundary::recvNGDStripe: Error: \n"); //OK
              throw oops;
            }
          }
        }
      }
    }
    break;
    //
    // the default case is to print an error message and exit
    //
  default:{
    stringstream ss (stringstream::in | stringstream::out);
    ss <<"SpatialRegionBoundary::recvNGDStripe: Error: \n"<<
        "Only right moving window (dir = 1)\n" <<
        "is implemented for now! You have requested dir = " << dir << endl;
      std::string msg;
      ss >> msg;
      Oops oops(msg);
      throw oops;    // exit() MapNGDs::sendNGDsAndShift: 

    }
    break;
  }
}

/**
 * A helper function to allocate memory for the 
 * data structures needed for the send/recv of
 * the NGD stuff. 
 */
void SpatialRegionBoundary::allocMemNGDbuffers() {
  if ( ptrMapNGDs ) {
    int numNGDs = ptrMapNGDs->getNumNGDs();
    if ( numNGDs ) {
#ifdef MPI_VERSION
      extern int MPI_RANK;
      cout << "MPI rank " << MPI_RANK 
           << " in void SpatialRegionBoundary::allocMemNGDbuffers()" 
           << endl;
#endif
      cout << "Allocating memory for NGD buffers..." << endl;
      int length = MAX(abs(k1-k2),abs(j1-j2));
      //
      // currently the number of elements in the NGD buffers
      // is determined by the length times the number of NGDs
      // instantiated times at most 3 (there are currently three data 
      // sets of size length that have to be communicated; 
      // these are the NGDs, the np2c's, and the excess number
      // of ions, see the NGD class). dad 01/18/02. 
      // 
      numElementsNGDbuf = 0; 
      NGD* ptrNGD;
      oopicListIter<NGD> NGDIter(*(ptrMapNGDs->getPtrNGDList()));
      for ( NGDIter.restart(); !NGDIter.Done(); NGDIter++ ) {
        ptrNGD = NGDIter.current(); 
        if ( ptrNGD->getisTIOn() )
          numElementsNGDbuf += 3*length;
        else
          numElementsNGDbuf += length;
      }
      lNGDbuf = new Scalar[numElementsNGDbuf];
      lNGDbufSend = new Scalar[numElementsNGDbuf];    
      //
      // zero the memory
      // 
      memset(lNGDbuf,     0, numElementsNGDbuf*sizeof(Scalar));
      memset(lNGDbufSend, 0, numElementsNGDbuf*sizeof(Scalar));
    }
  }
}