spbound.cpp

pic 模拟程序！面向对象

		 ENode[j][k].set_e1(0.5* (intEdl[j][k].e1() + lEb[kk].e1())/grid->dl1(j+1,k));

		 if(k>0 && k<K) {
		    ENode[j][k].set_e2(w2p*intEdl[j][k].e2()/grid->dl2(j, k) 
				       + w2m*intEdl[j][k-1].e2()/grid->dl2(j, k-1));
		    BNode[j][k].set_e1(intBdS[j][k].e1()*w2p/grid->dS1(j, k)
				       + intBdS[j][k-1].e1()*w2m/grid->dS1(j, k-1));

		    BNode[j][k].set_e3(intBdS[j][k].e3()*w1p*w2p/grid->dS3(j+1, k) // uniform mesh
				       + lBb[kk].e3()*w1m*w2p/grid->dS3(j+1, k)  //uniform mesh
				       + intBdS[j][k-1].e3()*w1p*w2m/grid->dS3(j+1, k-1) // uniform mesh
				       + lBb[kk-1].e3()*w1m*w2m/grid->dS3(j+1, k-1)); //uniform mesh
		 }

		 if(k > 0 || !grid->axis()) {
		    BNode[j][k].set_e2(intBdS[j][k].e2()*(w1p/grid->dS2(j+1, k))  // uniform mesh 
				       + lBb[kk].e2()*w1m/grid->dS2(j+1, k));  //uniform mesh
		    ENode[j][k].set_e3(intEdl[j][k].e3()/grid->dl3(j+1,k)); // uniform mesh
		 }		       
		 if(k==K) {
		    ENode[j][k].set_e2(intEdl[j][k-1].e2()/grid->dl2(j+1, k-1));
		    BNode[j][k].set_e1(2*(intBdS[j][k-1].e1()/grid->dS1(j+1,k-1))
				       -BNode[j][k-1].e1());
		    BNode[j][k].set_e3(2* ( intBdS[j][k-1].e3()*w1p/grid->dS3(j+1,k-1) +
					   lBb[kk-1].e3()*w1m/grid->dS3(j+1,k-1))-BNode[j][k-1].e3()); //uniform mesh
		 }


	      }
	      else {  //normal == -1 , real region on left
		 // ASSUME:  Uniform mesh

		 Scalar w2m = grid->dl2(j,k) / (grid->dl2(j,k) + grid->dl2(j,MAX(0,k-1)));
		 Scalar w2p = 1-w2m;
		 Scalar w1m = 0.5;  // uniform mesh assumption
		 Scalar w1p = 0.5;  // uniform mesh assumption

		 ENode[j][k].set_e1(0.5* (intEdl[j-1][k].e1() + lEa[kk].e1())/grid->dl1(j-1,k));

		 if(k>0 && k<K) {
		    ENode[j][k].set_e2(w2p*intEdl[j][k].e2()/grid->dl2(j, k) 
				       + w2m*intEdl[j][k-1].e2()/grid->dl2(j, k-1));
		    BNode[j][k].set_e1(intBdS[j][k].e1()*w2p/grid->dS1(j, k)
				       + intBdS[j][k-1].e1()*w2m/grid->dS1(j, k-1));

		    BNode[j][k].set_e3(lBa[kk].e3()*w1p*w2p/grid->dS3(j-1, k) // uniform mesh
				       + intBdS[j-1][k].e3()*w1m*w2p/grid->dS3(j-1, k)  //uniform mesh
				       + lBa[kk-1].e3()*w1p*w2m/grid->dS3(j-1, k-1) // uniform mesh
				       + intBdS[j-1][k-1].e3()*w1m*w2m/grid->dS3(j-1, k-1)); //uniform mesh
		 }


		 if(k > 0 || !grid->axis()) {
		    BNode[j][k].set_e2(intBdS[j-1][k].e2()*(w1p/grid->dS2(j-1, k))  // uniform mesh 
				       + lBa[kk].e2()*w1m/grid->dS2(j-1, k));  //uniform mesh
		    ENode[j][k].set_e3(intEdl[j][k].e3()/grid->dl3(j-1,k)); // uniform mesh
		 }		       

		 if(k==K) {
		    ENode[j][k].set_e2(intEdl[j][k-1].e2()/grid->dl2(j-1, k-1));
		    BNode[j][k].set_e1(2*(intBdS[j][k-1].e1()/grid->dS1(j-1,k-1))
				       -BNode[j][k-1].e1());
		    BNode[j][k].set_e3(2* ( intBdS[j-1][k-1].e3()*w1p/grid->dS3(j-1,k-1) +
					   lBa[kk-1].e3()*w1m/grid->dS3(j-1,k-1))-BNode[j][k-1].e3()); //uniform mesh
		 }
		 

	      }
	   }

	   // This needs to be AFTER the for loop because it depends on BNode higher up than the for
	   // loop gets to!!!
	   if(normal==1) {
	      k = k1;
	      int kk = k - k1;
	      if(k==0) {
		 BNode[j][k].set_e1(2*(intBdS[j][k].e1()/grid->dS1(j+1,k))
				 -BNode[j][k+1].e1());
		 if(grid->axis()) {
		    ENode[j][k].set_e2(0);  // none of these fields may be nonzero on axis.
		    BNode[j][k].set_e2(0);
		    BNode[j][k].set_e3(0);

		 } else {
		    ENode[j][k].set_e2(intEdl[j][k].e2()/grid->dl2(j+1,k));  // uniform mesh
		    // Refer to boundary.cpp for this one.
		    BNode[j][k].set_e3(2*(intBdS[j][k].e3()*w1p/grid->dS3(j+1,k)
					  +lBb[kk].e3()*w1m/grid->dS3(j+1,k))-BNode[j][1].e3()); // uniform mesh
		 }
	      }
	   }
	   else {
	      k = k1;
	      int kk = k - k1;
	      if(k==0) {
		 BNode[j][k].set_e1(2*(intBdS[j][k].e1()/grid->dS1(j-1,k))
				    -BNode[j][k+1].e1());
		 if(grid->axis()) {
		    ENode[j][k].set_e2(0);  // none of these fields may be nonzero on axis.
		    BNode[j][k].set_e2(0);
		    BNode[j][k].set_e3(0);
		 } else {
		    ENode[j][k].set_e2(intEdl[j][k].e2()/grid->dl2(j-1, k));
		    BNode[j][k].set_e3(2*(intBdS[j-1][k].e3()*w1p/grid->dS3(j-1,k)
					  +lBa[kk].e3()*w1m/grid->dS3(j-1,k))-BNode[j][1].e3()); // uniform mesh
		 }
	      }
	   }
	}
     else if (k1==k2)                
	{
	   // NIY
	}
  }

}

  

SpatialRegionBoundary::~SpatialRegionBoundary() {
  delete [] lJ;
  delete [] lE;
  delete [] lB;
  lJ=lE=lB=0;
  //
  // delete the memory for the NGD data buffers only if
  // it was allocated
  // 
  if ( lNGDbuf )
    delete [] lNGDbuf; 
  if ( lNGDbufSend )
    delete [] lNGDbufSend;
}

//
//  This is the function which packages and sends particles over 
// to the paired SRB.
//
	 
void SpatialRegionBoundary::passParticles() {
  theLink->sendParticles(number_Particles_passing,ptclPassArray);
}

void SpatialRegionBoundary::collect(Particle &p, Vector3& dxMKS) {

// check and see if we've exceeded the passing array memory
  while(number_Particles_passing+1 > ptclPassArraySize) {
    particlePassDat *array = new particlePassDat[2*ptclPassArraySize];
    memcpy(array, ptclPassArray,
	number_Particles_passing * sizeof(particlePassDat));
    delete[] ptclPassArray;
    ptclPassArray = array;
    ptclPassArraySize*=2;
  }

// Add this particle to the array for passing
  ptclPassArray[number_Particles_passing].u = p.get_u();
  if ( p.isVariableWeight() )
       ptclPassArray[number_Particles_passing].vary_np2c = 1.;
  else ptclPassArray[number_Particles_passing].vary_np2c = 0.;
  ptclPassArray[number_Particles_passing].np2c = p.get_np2c();
  ptclPassArray[number_Particles_passing].speciesID = (Scalar) p.get_speciesID();
  ptclPassArray[number_Particles_passing].dxMKS = dxMKS;
  ptclPassArray[number_Particles_passing].x =
     (j1==j2) ? p.get_x().e2() : p.get_x().e1();

// Diagnostics
/*
#ifdef MPI_VERSION
  if ( number_Particles_passing == 0) {
    extern int MPI_RANK;
    cout << "MPI rank " << MPI_RANK << " collected a particle with\n";
    cout << "\n";
  }
#endif
*/

  delete &p;
  number_Particles_passing++;

}

// Emit the particles coming across the boundary.  In this case
// the particles come in and must be weighted in the incoming region.
// This method assumes that calling routine/class will take
// care of emptying the returned ParticleList.  Thus, multiple
// calls to recvParticles accumulate particles.

ParticleList &SpatialRegionBoundary::emit(Scalar t, Scalar dt, 
	Species *species) {

  int i,n_in;
  particlePassDat *inParticles = theLink->waitParticles(&n_in);
  Boundary *bPtr;

  number_Particles_passing=0;  // our send has completed here.
#ifdef DEBUG
  printf("\nSRB:emit:  emitting %d particles",n_in);
#endif

  for(i=0;i<n_in;i++) {

// Get the particle's coordinate
    Vector2 x( (j1==j2)? j1* (1 +normal * 1e-6) +
	normal * 1e-20:inParticles[i].x, 
	(j1==j2) ? inParticles[i].x : k1*(1+normal*1e-6)+normal * 1e-20);

	 // figure out if it is a variable-weight particle
	 BOOL _vary_np2c = FALSE;
	 if(inParticles[i].vary_np2c > 0.) _vary_np2c = TRUE;


    Particle * p = new Particle(x,inParticles[i].u, 
	speciesArray[(int)(inParticles[i].speciesID+0.1)],
	inParticles[i].np2c, _vary_np2c);

// finish the push
    x = p->get_x();
    if((bPtr = fields->translateAccumulate(x, inParticles[i].dxMKS, 
		p->get_q()/fields->get_dt()))) {
      p->set_x(x);
// Send this particle to boundary
      bPtr->collect(*p,inParticles[i].dxMKS);

    }
    else {
      p->set_x(x);
      particleList.add(p);
    }

  }

  return particleList;
}

/**
 * Collect the particle into the shiftedParticles array.
 * This is a handoff, SpatialRegionBoundary is subsequently
 * responsible for deletion of the particle.
 *
 * @param p the particle to be collected
 */
void SpatialRegionBoundary::collectShiftedParticle(Particle* p){
  shiftedParticles.push_back(p);
}

/**
 * Pass the shifted particles over the boundary
 *
 * Because we are responsible for deleting the particles,
 * we must delete after sending.
 */
void SpatialRegionBoundary::passShiftedParticles(){

// Send the particles to the link
/*
#ifdef MPI_VERSION
  extern int MPI_RANK;
  cerr << "MPI_RANK of " << MPI_RANK << ": sending " << 
	shiftedParticles.size() << " particles.\n";
#endif
*/

  if(theLink) theLink->sendShiftedParticles(shiftedParticles);

// Clear out the particles
  for(int i=0; (unsigned)i<shiftedParticles.size(); i++) delete shiftedParticles[i];
  shiftedParticles.clear();

}

/**
 * Receive the shifted particles from the boundary
 */
ParticleList& SpatialRegionBoundary::recvShiftedParticles(){

#ifdef MPI_VERSION
  extern int MPI_RANK;
#endif

// Must have a link
  assert(theLink);

// Return the particles from the link
  // cout << "MPI_RANK = " << MPI_RANK << 
	// ": receiving particles in recvShiftedParticles." << endl;
  ParticleList& pl = theLink->recvShiftedParticles(speciesArray);
  // cout << "MPI_RANK = " << MPI_RANK << 
	// ": particles received in recvShiftedParticles." << endl;

/*
#ifdef MPI_VERSION
  extern int MPI_RANK;
  cerr << "MPI_RANK of " << MPI_RANK << ": received " << 
	pl.nItems() << " particles.\n";
#endif
*/

// Shift the particles over to the right
  oopicListIter<Particle> piter(pl);
  Particle* pp;
  Vector2 x;
  bool prtd = false;
  for(piter.restart(); !piter.Done(); piter++){
    pp = piter.current();
    // if(!prtd){
      // cout << "Before shift, first particle received has x1 = " << 
	// pp->get_x().e1() << endl;
    // }
    x = pp->get_x();
    x.set_e1(x.e1() + j2);
    pp->set_x(x);

    // if(!prtd){
      // cout << "After shift, first particle received has x1 = " << 
	// pp->get_x().e1() << endl;
    // }
    prtd = true;
  }

  return pl;

}

// Ask the link to get the stripe of fields
//
void SpatialRegionBoundary::askStripe(){
  assert(theLink);
  theLink->askShiftedFields();
}


// This is like putCharge_and_Current, except that it sends over
// a pure stripe as is needed for shifting fields
//
void SpatialRegionBoundary::sendStripe(int dir){

#ifdef  MPI_VERSION
  extern int MPI_RANK;
#endif

  Vector3 **intBdS = fields->getIntBdS();
  Vector3 **intEdl = fields->getIntEdl();
  Vector3 **I = fields->getI();

  switch(dir){

// Fields sending to left (right moving window)
    case 1:{
      assert(j1 == j2);
// On KC's authority, one always has k1 <= k2
      for(int k=k1; k<=k2; k++) {
	lEsend[k-k1] = Vector3( intEdl[j1][k].e1(), intEdl[j1+1][k].e2(),
		intEdl[j1+1][k].e3() );
	lBsend[k-k1] = Vector3( intBdS[j1+1][k].e1(), intBdS[j1][k].e2(),
		intBdS[j1][k].e3() );
	lJsend[k-k1] = Vector3( I[j1][k].e1(), I[j1+1][k].e2(),
		I[j1+1][k].e3() );
      }
    }
    break;

// Fields sending to right (left moving window)
    case 2:{
      assert(j1 == j2);
// On KC's authority, one always has k1 <= k2
      for(int k=k1; k<=k2; k++) {
        lEsend[k-k1] = intEdl[j1-1][k];
        lBsend[k-k1] = intBdS[j1-1][k];
        lJsend[k-k1] = I[j1-1][k];
      }
    }
    break;

// Fields sending up (to smaller 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++) {
	lEsend[j-j1] = Vector3( intEdl[j][k1+1].e1(), intEdl[j][k1].e2(),
		intEdl[j][k1+1].e3() );
	lBsend[j-j1] = Vector3( intBdS[j][k1].e1(), intBdS[j][k1+1].e2(),
		intBdS[j][k1].e3() );
	lJsend[j-j1] = Vector3( I[j][k1+1].e1(), I[j][k1].e2(),
		I[j][k1+1].e3() );
      }
    }
    break;

// Fields sending down (to greater 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++) {
        lEsend[j-j1] = intEdl[j][k1-1];
        lBsend[j-j1] = intBdS[j][k1-1];
        lJsend[j-j1] = I[j][k1-1];
      }
    }
    break;

  }

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

// Send the fields through the link
/*
#ifdef  MPI_VERSION
  cout << "MPI_RANK = " << MPI_RANK << 
	", SpatialRegionBoundary::sendStripe sending fields. j1,j2,k1,k2 = " <<
	j1 << "," << j2 << "," << k1 << "," << k2 << endl;
#endif
*/

  theLink->sendShiftedFields(lEsend,lBsend,lJsend);

/*
#ifdef  MPI_VERSION
  cout << "MPI_RANK = " << MPI_RANK << 
	", SpatialRegionBoundary::sendStripe sent fields." << endl;
#endif
*/

}

void SpatialRegionBoundary::recvStripe(int dir){

  Vector3 **intBdS = fields->getIntBdS();
  Vector3 **intEdl = fields->getIntEdl();
  Vector3 **I = fields->getI();

  assert(theLink);
#ifdef  MPI_VERSION
  extern int MPI_RANK;
#endif

/*
#ifdef  MPI_VERSION
  cout << "MPI_RANK = " << MPI_RANK << 
	", SpatialRegionBoundary::recvStripe awaiting fields. j1,j2,k1,k2 = " <<
	j1 << "," << j2 << "," << k1 << "," << k2 << endl;
#endif
*/

  theLink->waitShiftedFields();  //  do a synchronize

/*
#ifdef  MPI_VERSION
  cout << "MPI_RANK = " << MPI_RANK << 
	", SpatialRegionBoundary::recvStripe received fields." << endl;
#endif
*/

/*