fields.cpp

pic 模拟程序！面向对象

      for(int j1=0; j1<=J; j1++)  cout << " " << I[j1][kout].e1();
      cout << endl;

      #ifdef MPI_VERSION
      }
      #endif
      */

   int	j, k;
   local_dt = dt;  // so that objects needing to manipulate
   // fields directly can get_dt() easily
   Scalar local_t = t;

   //  get any charge on any dielectrics into rho, put any currents generated by loops
   // this current weighting does not subcycle properly

   //ApplyToList(putCharge_and_Current(t),*boundaryList,Boundary);

   if(ElectrostaticFlag) 
      ElectrostaticSolve(t,dt);

   else 
      {
         if(DivergenceCleanFlag) {
            updateDivDerror();
            DivergenceClean();
            updateDivDerror();
         }
         
         if(CurrentWeightingFlag)  InodeToI();  //  the bilinear weighting needs this
         if (grid->axis()) radialCurrentCorrection(); // correct current at 0 and K

         for (int i = 0; i < FieldSubFlag; i++, local_t += dt)
            {
               if (EMdamping>0)
                  intEdl = intEdlPrime;

               if (i==0) advanceB(0.5*dt);		   //	advance B^(n-1) -> B^(n-1/2)
               else advanceB(dt);		   			//	advance B^(n-3/2) -> B^(n-1/2)

               if (EMdamping>0)
                  {
                     intEdl = intEdlBasePtr;
                     for (j=0; j<=J; j++)
                        for (k=0; k<=K; k++) {
                           intEdlBar[j][k] *= EMdamping;
                           intEdlBar[j][k] += (1-EMdamping)*intEdl[j][k];
                           intEdlPrime[j][k] = .5*((1-EMdamping)*intEdlBar[j][k] -intEdl[j][k]);
                        }
                  }
               
               advanceE(dt);                   //	advance E^(n-1) -> E^n

               //	Apply Boundary Conditions at t
               try{
                 ApplyToList(applyFields(local_t,dt),*boundaryList,Boundary);
               }
               catch(Oops& oops3){
                 oops3.prepend("Fields::advance: Error: \n");
                 throw oops3;
               }
               if (EMdamping>0)
                  for (j=0; j<=J; j++)
                     for (k=0; k<=K; k++) 
                        intEdlPrime[j][k] += (1+.5*EMdamping)*intEdl[j][k];
            }

         if (EMdamping>0)
            intEdl = intEdlPrime;

         advanceB(0.5*dt);
         
         if (EMdamping>0)
            intEdl = intEdlBasePtr;
      }
   
   if(FieldSubFlag==1)   //ask peter about this
      if(MarderIter>0) MarderCorrection(dt);

   /*
      #ifdef MPI_VERSION
      extern int MPI_RANK;
      if(MPI_RANK == 0){
      cout << "MPI_RANK = " << MPI_RANK << ", leaving Fields::advance.  E1 =";
      #else
      cout << "Leaving Fields::advance.  E1 =";
      #endif
      for(int i1=0; i1<=J; i1++)  cerr << " " << intEdl[i1][kout].e1();
      cout << endl;
      cout << "E1 =";
      for(int j1=0; j1<=J; j1++)  cout << " " << intEdl[j1][kout].e2();
      cout << endl;
      cout << "I1 =";
      for(int j1=0; j1<=J; j1++)  cout << " " << I[j1][kout].e1();
      cout << endl;

      #ifdef MPI_VERSION
      }
      #endif
      */
   
}

void Fields::ElectrostaticSolve(Scalar t, Scalar dt)
  throw(Oops){
   if(BoltzmannFlag)	{
      oopicListIter<Boundary> nextb(*boundaryList);
      for (nextb.restart(); !nextb.Done(); nextb++){
        try{
          nextb.current()->applyFields(t,dt);
        }
        catch(Oops& oops3){
          oops3.prepend("Fields::ElectrostaticSolve: Error: \n");
          throw oops3;
        }
      }
      theBoltzmann->updatePhi(rho,Phi,t,dt);
   }
   else updatePhi(rho,Phi,t,dt);
   setEfromPhi();
}

//--------------------------------------------------------------------
//	Advances B by the time step specified.  This function is called twice
//	per time step by advance() with a half timestep each time.

void Fields::advanceB(Scalar dt2)
{
   int	j, k;					//	indices
   for (j=0; j<J; j++)
      {
	 for (k=0; k<K; k++)
	    {
	       intBdS[j][k] -= dt2*Vector3(intEdl[j][k+1].e3() - intEdl[j][k].e3(),
					   -intEdl[j+1][k].e3() + intEdl[j][k].e3(),
					   intEdl[j][k].e1() - intEdl[j][k+1].e1()
					   + intEdl[j+1][k].e2() - intEdl[j][k].e2());
	    }
	 //	get intBdS[j][K].e2 only
	 intBdS[j][K] -= dt2*Vector3(0, intEdl[j][K].e3() - intEdl[j+1][K].e3(), 0);
      }
   for (k=0; k<K; k++)				//	intBdS[J][k].e1
      intBdS[J][k] -= dt2*Vector3(intEdl[J][k+1].e3()	- intEdl[J][k].e3(), 0, 0);
   
}

void Fields::advanceE(Scalar dt) {
   int j,k;

   Vector3*	iL0 = &iL[0][0];
   intEdl[0][0] += dt*iC[0][0].jvMult(Vector3(
					      intBdS[0][0].e3()*iL0->e3(),
					      -intBdS[0][0].e3()*iL0->e3(),
					      -intBdS[0][0].e1()*iL0->e1() + intBdS[0][0].e2()*iL0->e2()) 
				      - I[0][0]);
   iL0 = &iL[0][K];
   Vector3*	iLk1 = &iL[0][K-1];
   intEdl[0][K] += dt*iC[0][K].jvMult(Vector3(
					      -intBdS[0][K-1].e3()*iLk1->e3(),
					      0,											 
					      intBdS[0][K-1].e1()*iLk1->e1()	+ intBdS[0][K].e2()*iL0->e2())	
				      - I[0][K]);
   iL0 = &iL[J][0];
   Vector3*	iLj1 = &iL[J-1][0];
   intEdl[J][0] += dt*iC[J][0].jvMult(Vector3(
					      0,
					      intBdS[J-1][0].e3()*iLj1->e3(),
					      -intBdS[J][0].e1()*iL0->e1() - intBdS[J-1][0].e2()*iLj1->e2()) 
				      - I[J][0]);
   iL0 = &iL[J][K];
   iLk1 = &iL[J][K-1];
   iLj1 = &iL[J-1][K];
   intEdl[J][K] += dt*iC[J][K].jvMult(Vector3(
					      0,
					      0,
					      intBdS[J][K-1].e1()*iLk1->e1() - intBdS[J-1][K].e2()*iLj1->e2()) 
				      - I[J][K]);
   for (k=1; k<K; k++)
      {
	 iL0 = &iL[0][k];
	 iLk1 = &iL[0][k-1];
	 intEdl[0][k] += dt*iC[0][k].jvMult(Vector3(
						    intBdS[0][k].e3()*iL0->e3() - intBdS[0][k-1].e3()*iLk1->e3(),
						    -intBdS[0][k].e3()*iL0->e3(),
						    intBdS[0][k-1].e1()*iLk1->e1() - intBdS[0][k].e1()*iL0->e1()
						    + intBdS[0][k].e2()*iL0->e2()) 
					    - I[0][k]);
	 iL0 = &iL[J][k];
	 iLk1 = &iL[J][k-1];
	 iLj1 = &iL[J-1][k];
	 intEdl[J][k] += dt*iC[J][k].jvMult(Vector3(
						    //								 intBdS[J][k].e3()*iL0->e3() - intBdS[J][k-1].e3()*iLk1->e3(),
						    0,   // component 1 outside system
						    intBdS[J-1][k].e3()*iLj1->e3(),
						    intBdS[J][k-1].e1()*iLk1->e1() - intBdS[J][k].e1()*iL0->e1()
						    - intBdS[J-1][k].e2()*iLj1->e2()) 
					    - I[J][k]);
      }
   Vector3 *Bk1, *Bj1, *B0;
   for (j=1; j<J; j++)
      {
	 iL0 = &iL[j][0];
	 iLj1 = &iL[j-1][0];
	 B0 = &intBdS[j][0];
	 Bj1 = &intBdS[j-1][0];
	 intEdl[j][0] += dt*iC[j][0].jvMult(Vector3(
						    intBdS[j][0].e3()*iL0->e3(),
						    -intBdS[j][0].e3()*iL0->e3() + intBdS[j-1][0].e3()*iLj1->e3(),
						    -intBdS[j][0].e1()*iL0->e1() + intBdS[j][0].e2()*iL0->e2()
						    - intBdS[j-1][0].e2()*iLj1->e2()) 
					    - I[j][0]);
	 for (k=1; k<K; k++)
	    {
	       iLk1 = iL0;
	       iLj1 = &iL[j-1][k];
	       iL0 = &iL[j][k];
	       Bk1 = B0;
	       Bj1 = &intBdS[j-1][k];
	       B0 = &intBdS[j][k];
	       intEdl[j][k] += dt*iC[j][k].jvMult(Vector3(
							  B0->e3()*iL0->e3() - Bk1->e3()*iLk1->e3(),
							  -B0->e3()*iL0->e3() + Bj1->e3()*iLj1->e3(),
							  Bk1->e1()*iLk1->e1() - B0->e1()*iL0->e1()
							  + B0->e2()*iL0->e2() - Bj1->e2()*iLj1->e2()) 
						  - I[j][k]);
	    }
	 iL0 = &iL[j][K];
	 iLk1 = &iL[j][K-1];
	 intEdl[j][K] += dt*iC[j][K].jvMult(Vector3(
						    -intBdS[j][K-1].e3()*iLk1->e3(),
						    0,
						    intBdS[j][K-1].e1()*iLk1->e1()
						    + intBdS[j][K].e2()*iL0->e2() - intBdS[j-1][K].e2()*iLj1->e2())
					    - I[j][K]);
      }

}

//--------------------------------------------------------------------
//	Translates to new locations by incrementally locating cell crossings
//	and accumulates current along the trajectory.  Initial position passed
//	in x is updated to final position on return.  Return value is a
//	pointer to a boundary if encountered, NULL otherwise.
//	qOverDt = q/dt*particleWeight for this particle.

Boundary* Fields::translateAccumulate(Vector2& x, Vector3& dxMKS, Scalar
                                      qOverDt)
{
   Boundary*	boundary = NULL;
   int	which_dx = (dxMKS.e1() != 0.0) ? 1 : 2;
   Scalar	prev_dxMKS = (which_dx == 1) ? dxMKS.e1() : dxMKS.e2();
   Scalar	frac_v3dt;

   while (fabs(dxMKS.e1()) + fabs(dxMKS.e2()) > 1E-25){
      Vector2 	xOld = x;
      boundary = grid->translate(x, dxMKS);

      //	If the particle translates in x1 or x2, apply the fraction of
      //	rotation for this cell only; otherwise it spins only and all
      //	current is collected at this location.
      if (prev_dxMKS != 0) frac_v3dt = dxMKS.e3()*(prev_dxMKS -
                                                   ((which_dx==1) ? dxMKS.e1() : dxMKS.e2()))/prev_dxMKS;
      else frac_v3dt = dxMKS.e3();
      accumulate(xOld, x, frac_v3dt, qOverDt);
      prev_dxMKS = (which_dx == 1) ? dxMKS.e1() : dxMKS.e2();
      if (boundary) break;
   }
   return boundary;
}

//--------------------------------------------------------------------
//	Weight intEdl and intBdS to nodes to obtain Enode and Bnode in
//	physical units.  This routine has been modified to set only the 
// interior points; edges and internal boundaries are set by calling
// Boundary::toNodes().

void Fields::toNodes(Scalar t) {

   int j, k;
   Scalar w1p, w1m, w2p, w2m;

   // all interior points in 1-direction
   for (k=j=1; j<J; j++) {

      // Calculate weight factors for 1-directions
      w1m = grid->dl1(j, k)/(grid->dl1(j-1, k) + grid->dl1(j, k));
      w1p = 1 - w1m;

      // Not sure what this is doine
      BNode[j][0]=BNode[j][K]=0;  // unbounded edges.

      // All interior points in the 2-direction

      // JRC: Could this be accelerated through use of smaller loops?
      for (k=1; k<K; k++) {

         // Get weights in 2-direction
         w2m = grid->dl2(j, k)/(grid->dl2(j, k-1) + grid->dl2(j, k));
         w2p = 1 - w2m;

         // Weight each of the fields

         // E1 defined on 2-node values, but weight in 1-direction
         ENode[j][k].set_e1(w1p*intEdl[j][k].e1()/grid->dl1(j, k) 
                            + w1m*intEdl[j-1][k].e1()/grid->dl1(j-1, k));

         // E2 defined on 1-node values, but weight in 2-direction
         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));

         // E3 defined at the nodes
         ENode[j][k].set_e3(intEdl[j][k].e3()/grid->dl3(j,k));

         // B1 defined on 1-node values, but weight in 2-direction
         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));

         // B2 defined on 2-node values, but weight in 1-direction
         BNode[j][k].set_e2(intBdS[j][k].e2()*(w1p/grid->dS2(j, k)) 
                            + intBdS[j-1][k].e2()*w1m/grid->dS2(j-1, k));

         // B3 weighted in both directions (a cell center quantity)
         BNode[j][k].set_e3(intBdS[j][k].e3()*w1p*w2p/grid->dS3(j, k)
                            + intBdS[j-1][k].e3()*w1m*w2p/grid->dS3(j-1, k)
                            + intBdS[j][k-1].e3()*w1p*w2m/grid->dS3(j, k-1)
                            + intBdS[j-1][k-1].e3()*w1m*w2m/grid->dS3(j-1, k-1));
      }
   }

   // Not sure what this is doing.  Shouldn't this be set by boundaries?
   for (k=0; k<=K; k++)	BNode[0][k] = BNode[J][k] = 0;  // unbounded edges.

   // boundaries interpolate coincident fields (assumes interior previously set):
   ApplyToList(toNodes(),*boundaryList,Boundary);

   // Add in the static magnetic field
   for (j=0; j<=J; j++)
      for (k=0; k<=K; k++) {
         BNodeDynamic[j][k] = BNode[j][k];
         BNode[j][k] += BNodeStatic[j][k];
      }

}

//--------------------------------------------------------------------
//	Zero the current accumulator array, I

void Fields::clearI(int i)
{
   Vector3** Iptr;
   if (ElectrostaticFlag) return;
   if (i == -1) Iptr = I;
   else Iptr = Ispecies[i];
   if (!Iptr) return;
   for (int j=0; j<=J; j++){
      if (CurrentWeightingFlag!=0) memset(Inode[j], 0, (K+1)*sizeof(Inode[j][0]));
      memset(Iptr[j], 0, (K+1)*sizeof(Iptr[j][0]));
   }
}

// Accumulate current for species I:

void Fields::setAccumulator(int i)
{
   if (Ispecies[i]) accumulator = Ispecies[i];
   else if (CurrentWeightingFlag) accumulator = Inode;
   else accumulator = I;
}

void Fields::addtoI(int i)
{
   if (!Ispecies[i]) return;
   for (int j=0; j<=J; j++)
      for (int k=0; k<=K; k++)
         I[j][k] += Ispecies[i][j][k]; 
}

void Fields::Excite() //ntg8-16 excites a field component at t = 0
{
   this->setIntEdl(J/2, K/2, 1, 1.0);
   this->setIntEdl(J/2, K/2, 2, -1.0);
   this->setIntEdl(J/2, K/2+1, 1, -1.0);
   this->setIntEdl(J/2+1, K/2, 2, 1.0);
   /*
      Scalar freq = 20E9;
      Scalar z = getMKS(J/2, 0).e1();
      Scalar outerradius = getMKS(J/2, K).e2();
      for(int k = 0; k <= K; k++)
      {Scalar r = getMKS(J/2, k).e2();
      Scalar Ezval = AnalyticEz(z, r, freq, 0.0, outerradius); 
      this->setIntEdl(J/2, k, 1, Ezval);}
      */
}

//--------------------------------------------------------------------
//	Seed the field emitter with a large negative field to get it
//	started emitting.

void Fields::SeedFieldEmitter()
{
   for (int j=0; j<J; j++)
      for (int k=0; k<=K; k++)
         setIntEdl(j, k, 1, -1000);
}

//--------------------------------------------------------------------
//	Set initial uniform magnetic field in x1. 
//	Here, Bz0 is in Tesla.

BOOL Fields::setBNodeStatic(Vector3 B0,Scalar betwig, Scalar zoff, char* BFname,const ostring &B01analytic,const ostring &B02analytic, const ostring &B03analytic)
{
   FILE *openfile;
   int status;
   ostring B01a = B01analytic;
   ostring B02a = B02analytic;
   ostring B03a = B03analytic;
   Vector3 Btemp;
   if (strcmp(BFname, "NULL")) {
      if ((openfile = fopen (BFname, "r"))  != NULL) {
         // read in
         Scalar dum1, dum2, B1, B2, B3; // dummies, remove later
         cout << "starting to read B field file" << endl;
         for (int j=0; j<=J; j++) {
            for (int k=0; k<=K; k++) {
#ifdef SCALAR_DOUBLE