fields.cpp

pic 模拟程序！面向对象

               // status = fscanf (openfile, "%lg %lg %lg %lg %lg", &dum1, &dum2, &B2, &B3, &B1);               
               status = fscanf (openfile, "%lg %lg %lg %lg %lg", &dum1, &dum2, &B1, &B2, &B3);
#else
               // status = fscanf (openfile, "%g %g %g %g %g", &dum1, &dum2, &B2, &B3, &B1);
               status = fscanf (openfile, "%g %g %g %g %g", &dum1, &dum2, &B1, &B2, &B3);
#endif
               Btemp.set_e1(B1*1e-4);
               Btemp.set_e2(B2*1e-4);
               Btemp.set_e3(B3*1e-4);
               BNodeStatic[j][k] = Btemp;
            }
            //				status = fscanf(openfile, "%e", &dum1);
         }
         cout << "finished reading B field file" << endl;
         fclose (openfile);
      }
      else {
         cout << "open failed on magnetic field file Bf in Control group";
         return FALSE;
      }
   }
   else if (B0.e1()!=0 || B0.e2()!=0 || B0.e3()!=0) {
      for (int j=0; j<=J; j++)
         for (int k=0; k<=K; k++)
            {Bz0 = (grid->query_geometry() == ZXGEOM) ? B0.e1() : B0.e1()
                + B0.e2()*(cos(betwig*(getMKS(j,k).e1()-zoff)));
             Bx0 = (grid->query_geometry() == ZXGEOM) ? B0.e2() :
                + .5*B0.e2()*betwig*getMKS(j,k).e2()*(sin(betwig*(getMKS(j,k).e1()-zoff)));
             Btemp.set_e1(Bz0);
             Btemp.set_e2(Bx0);
             Btemp.set_e3(B0.e3());
             BNodeStatic[j][k] = Btemp;
          }
   }
   else {
      Scalar x1,x2;  // These can now be scalar -- RT, 2003/12/09
      adv_eval->add_indirect_variable("x1",&x1);
      adv_eval->add_indirect_variable("x2",&x2);
      B01a+='\n';
      B02a+='\n';
      B03a+='\n';
      for(int j=0;j<=J;j++)
         for(int k=0;k<=K;k++) {
            x1 = grid->getMKS(j,k).e1();
            x2 = grid->getMKS(j,k).e2();
            Btemp.set_e1(adv_eval->Evaluate(B01a.c_str()));
            Btemp.set_e2(adv_eval->Evaluate(B02a.c_str()));
            Btemp.set_e3(adv_eval->Evaluate(B03a.c_str()));
            BNodeStatic[j][k] = Btemp;
         }
   }
   return TRUE;
}

void Fields::setBeamDump(int _j1, int _j2)
{
   int j1 = MIN(_j1, _j2);
   int j2 = MAX(_j1, _j2);
   if (j1 < 0 || j2 > J) return;			// default values or error
   Scalar x1 = getMKS(j1).e1();
   Scalar x2 = getMKS(j2).e1();
   Scalar beta = M_PI/(2*(x2 - x1));
   Scalar geomMult = (grid->query_geometry() == ZXGEOM) ? 1 : 0.5;
   int j,k;
   for (j=j1; j<=j2; j++)
      for (k=0; k<=K; k++)
         BNodeStatic[j][k] = Bz0*Vector3(cos(beta*(getMKS(j,k).e1()-x1)),
                                         geomMult*beta*getMKS(j,k).e2()*sin(beta*(getMKS(j,k).e1()-x1)),
                                         BNodeStatic[j][k].e3());
   for (j=j2+1; j<=J; j++)
      for (k=0; k<=K; k++)
         BNodeStatic[j][k] = Vector3(0, 0, 0);
}

//--------------------------------------------------------------------
//	Update the DivD array of divergences.  This is not intended to
//      be kept updated but rather to be updated by request.

void Fields::updateDivDerror(void)
{ int j,k;
  Vector3 D1,D2;
  Scalar Div,ic;
  //This section of the code computes the divergence of D everywhere.

  for(j=1;j<=J-1;j++)   // bounds are temporary, need to check for metals eventually
     for(k=1 ;k<=K-1;k++) {

        if( j>0&&(ic=iC[j-1][k].e1())!=0.0)
           D1.set_e1( intEdl[j-1][k].e1()/ic );
        else
           D1.set_e1(0);

        if( k>0&&(ic=iC[j][k-1].e2())!=0.0)
           D1.set_e2( intEdl[j][k-1].e2()/ic );
        else 
           D1.set_e2(0);

        if( j<J&& (ic=iC[j][k].e1())!=0.0)
           D2.set_e1( intEdl[j][k].e1()/ic );
        else
           D2.set_e1(0);
        if( k<K&&(ic=iC[j][k].e2()) !=0.0)
           D2.set_e2( intEdl[j][k].e2()/ic );
        else
           D2.set_e2(0);

        //		D1.set_e2(D2.e2()); //from the previous
        //D1.set_e1( intEdl[j][k].e1()*iC[j][k].e1() );
        //D2.set_e1( intEdl[j+1][k].e1()*iC[j+1][k].e1() );
        //D2.set_e2( intEdl[j][k+1].e2()*iC[j][k+1].e2() );

        Div = D2.e1() - D1.e1() + D2.e2() - D1.e2();
        Div /= grid->get_halfCellVolumes()[j][k];
        DivDerror[j][k]=-Div + rho[j][k];
     }
}


//--------------------------------------------------------------------
//        compute the potential everywhere using Fields::rho as a 
//	  source term


void Fields::updatePhi(Scalar **source,Scalar **dest,Scalar t,Scalar dt)
  throw(Oops){  
   int j,k;

   //	presiduetol_test=1e-4;  //tolerance for the poisson solver
   static int itermax=100;  //maximum # of iterations for psolve
   oopicListIter<Boundary> nextb(*boundaryList);
   
   // set up the boundary arrays if necessary
   //  they will not have been set up if u_z0 = 0
   if(minusrho==0) { 
      minusrho = new Scalar* [J+1];
      for(j=0;j<=J;j++) {
         minusrho[j]=new Scalar[K+1];
         memset(minusrho[j],0,(K+1)*sizeof(Scalar));
      }
      itermax = 100;
   }

   //  Scale the source term by -1
   for(j=0;j<=J;j++) for(k=0;k<=K;k++) 
      minusrho[j][k]=-(Scalar)source[j][k];

   //actually DO the solve
   psolve->solve(PhiP,minusrho,itermax,presidue);

   // To add all the LaPlace solutions, we first need to start
   //  with the Poisson solution.
   for(j=0;j<=J;j++) for(k=0;k<=K;k++) 
      dest[j][k]=PhiP[j][k];

   // apply any laplace boundary conditions to phi
   for (nextb.restart(); !nextb.Done(); nextb++)
   {
     try{
         nextb.current()->applyFields(t,dt);
     }
     catch(Oops& oops3){
       oops3.prepend("Fields::updatePhi: Error: \n");
       throw oops3;
     }
   }
   //  assign the solution to dest
   for(j=0;j<=J;j++) for(k=0;k<=K;k++) 
      dest[j][k]=Phi[j][k];


}

// -----------------------------
//  complete those tasks necessary to initialize the poisson
//  solver--symmetry boundary at zero, dirichlet otherwise

void Fields::initPhi() 
throw (Oops){

#ifndef PETSC
   if(grid->getPeriodicFlagX1()&&ElectrostaticFlag == DADI_SOLVE) 
      ElectrostaticFlag = PERIODIC_X1_DADI; 
   switch(ElectrostaticFlag) {
    default:
       fprintf(stderr,"WARNING, poisson solve defaulting to Multigrid\n");

     case ELECTROMAGNETIC_SOLVE:  
        psolve=new Multigrid(epsi,grid);
       break;
     case DADI_SOLVE:

        switch(grid->query_geometry()) {
         case ZRGEOM:
            psolve=new Dadirz(epsi,grid);
            break;
          case ZXGEOM:
             psolve=new Dadixy(epsi,grid);
            break;
         }
       break;

       // morgan
       // begin Krylov inverters
     case CG_SOLVE:
       try{
         set_up_inverter(grid);
       }
       catch(Oops& oops3){
         oops3.prepend("Fields::initPhi: Error: \n");
         throw oops3;
       }

       //Conjugate_Gradient *cg = new Conjugate_Gradient(dom,poisson);
       psolve = new Conjugate_Gradient(dom,poisson);
       break;
     case GMRES_SOLVE:
       try{
         set_up_inverter(grid);
       }
       catch(Oops& oops3){
         oops3.prepend("Fields::initPhi: Error: \n");  //OK
         throw oops3;
       }

       psolve = new GMRES(dom,poisson);
       break;
       // end Krylov inverters
     case MGRID_SOLVE:
        psolve=new Multigrid(epsi,grid);
       break;
     case PERIODIC_X1_DADI:
       try{
         psolve=new DadixyPer(epsi,grid);
       }
       catch(Oops& oops3){
         oops3.prepend("Fields::initPhi: Error: \n");//OK
         throw oops3;
       }
     break;
       
       
    }
#else // ndef-PETSC
  try{
    psolve = new ParallelPoisson(epsi,grid);
  }
  catch(Oops& oops3){
    oops3.prepend("Fields::initPhi: Error: \n");//OK
    throw oops3;
  }

#endif // ndef-PETSC
}



//--------------------------------------------
//  this collects the charge from the particle group
//  into rho[][]  If the group == 0 then it clears rho[][]

void Fields::collect_charge_to_grid(ParticleGroup *group,Scalar **_rho)
{
   int i,j,k;
   Scalar q,jx,kx;
   Vector2 *x;

   if(group==0) {  // reset flag - zero this rho
      for(j=0;j<=J;j++) memset(_rho[j],0,(K+1)*sizeof(Scalar));
      return;
   }

   q = group->get_q();  //charge * np2c
   x = group->get_x();
   for(i=0; i<group->get_n(); i++) {
      jx = x[i].e1();
      kx = x[i].e2();
      j = (int)jx; k=(int)kx;  //get the integer part
      jx -= j; kx -= k;  //get the fractional part only
      if (group->vary_np2c()) q = group->get_q(i);

      _rho[j][k] += (1-jx)*(1-kx)*q;
      _rho[j][k+1] += (1-jx)*(kx)*q;
      _rho[j+1][k] += (jx)*(1-kx)*q;
      _rho[j+1][k+1] += (jx)*(kx)*q;
   }
}

//----------------------------------------------------------------------
//function iterates through the particle group list given it
//to sum up all the charge using collect_charge_to_grid()
//

void Fields::compute_rho(ParticleGroupList **list,int nSpecies) 
{
   int j;
   
   if (rho_species==0){ // not initialized, we need to allocate and initialize
      int J=grid->getJ();
      int K=grid->getK();
      rho_species = new Scalar **[nSpecies];
      for(int i=0;i<nSpecies;i++) {
         rho_species[i] = new Scalar *[J+1];
         for(j=0;j<J+1;j++) {
            rho_species[i][j]=new Scalar[K+1];
            //zero the memory
            memset(rho_species[i][j],0,(K+1)*sizeof(Scalar));
         }
      }
   }

#ifndef PHI_TEST
   collect_charge_to_grid(0,rho);  //reset the global rho
   for (int i=0; i<nSpecies; i++){
      collect_charge_to_grid(0,rho_species[i]);  //reset the species rho
      oopicListIter<ParticleGroup> pgscan(*list[i]);
      for(pgscan.restart();!pgscan.Done();pgscan++)
         collect_charge_to_grid(pgscan(),rho_species[i]);
   }
   charge_to_charge_density(nSpecies);
   grid->binomialFilter(rho, nSmoothing);
   addBGRho();

#else
   
   for(j=0;j<J;j++)
      for(int k=0;k<K;k++)
         rho[j][k] = analytic_rho(j,k);
#endif
}

void Fields::ZeroCharge(void)
{
   for(int j=0;j<=J;j++)
      for(int k=0;k<=K;k++)
         Charge[j][k] = 0.0;
}
void Fields::charge_to_charge_density(int nSpecies)
{ 
   int i,j,k;
   Scalar **CellVolumes=grid->get_halfCellVolumes();
   
   for(i=0;i<nSpecies;i++)
      for(j=0;j<=J;j++) 
         for(k=0;k<=K;k++) {
            rho_species[i][j][k]/=CellVolumes[j][k];
            rho[j][k]+=rho_species[i][j][k];
         }
   if (grid->getPeriodicFlagX1()){
      for (k=0; k<=K; k++) {
         rho[0][k] += rho[J][k]; // note: CellVol along these edges is doubled
         rho[J][k] = rho[0][k];
      }
      for (i=0; i<nSpecies; i++)
         for (k=0; k<=K; k++) {
            rho_species[i][0][k] += rho_species[i][J][k];
            rho_species[i][J][k] = rho_species[i][0][k];
         }
   }
   if (grid->getPeriodicFlagX2()){
      for (j=0; j<=J; j++) {
         rho[j][0] += rho[j][K];
         rho[j][K] = rho[j][0];
      }
      for (i=0; i<nSpecies; i++)
         for (j=0; j<=J; j++) {
            rho_species[i][j][0] += rho_species[i][j][K];
            rho_species[i][j][K] = rho_species[i][j][0];
         }
   }
}

void Fields::addBGRho()
{
   if(BGRhoFlag)
      for(int j=0;j<=J;j++) 
         for(int k=0;k<=K;k++) 
            rho[j][k] += backGroundRho[j][k];
}
//----------------------------------------------------------------------
//  Given a z-r solution of poisson\'s equation from updatePhi in
//  Scalar **Phi, make the electric fields consistent with it.
//  This is an initialization of the E-fields.  It will destroy any
//  previous E fields.  This function cannot touch e3.

void Fields::setEfromPhi()
{ 
   int j,k;
   
   //  phi solution is for phi at nodes.  This makes it convenient to
   //  get E at half-cells.  Fortunately, this is what intEdl needs.

   for(j=0;j<J;j++)
      for(k=0;k<K;k++) {
         setIntEdl(j,k,1,  (Phi[j][k]-Phi[j+1][k])   );
         setIntEdl(j,k,2,  (Phi[j][k]-Phi[j][k+1])   );
      }

   for(j=0;j<J;j++) setIntEdl(j,K,1, (Phi[j][K]-Phi[j+1][K])  );
   for(k=0;k<K;k++) setIntEdl(J,k,2, (Phi[J][k]-Phi[J][k+1])  );

}


//----------------------------------------------------------------------
//  Use the current rho and a poisson solver to correct the errors in the
//  irrotational E fields.  compute_rho and updateDivDerror should be called first.
//  Also initPhi should be called.
//  

void Fields::DivergenceClean(void)
{ int j,k;

  if(!delPhi) {  // get space if this is the first call.
     delPhi= new Scalar *[J +1];
     for(j=0;j<=J;j++) 
        delPhi[j] = new Scalar [K+1];
     for(j=0;j<=J;j++) for(k=0;k<=K;k++) delPhi[j][k]=0;
  }

  // get the phi we need
  psolve->solve(delPhi,DivDerror,5,presidue); // use the Poisson solver directly,
  //delPhi is actually returned as -delPhi

  // Now we need to perform the actual correction to the intEdl everywhere.