Implementing the periodic boundary condition manually

[rahulgaur104]

Hi everyone,

I am simulating a rapidly rotating plasma in a cylinder in presence of a background magnetic field. To understand the setup, please look at the snapshot below:

The magnetic field lines are roughly along the contours of constant v2 velocity. The viscosity and resistivity are turned on with values of 0.01 and 0.0002, respectively.

I want to implement reflecting boundary conditions on the v2 velocity in the marked(black curly bracket) region and periodic boundary conditions on the rest of the variables everywhere on the z-boundaries. Our argument is that the field lines intersecting the top wall should be frozen-in as the wall is a perfect conductor. Hence, the plasma cannot move perpendicular to the field lines.

This is what my custom boundary condition function looks like.

`

void OuterZ(MeshBlock *pmb, Coordinates *pco,
                         AthenaArray<Real> &prim,
                         FaceField &b, Real time, Real dt,
            		 int il, int iu, int jl, int ju, int kl, int ku, int ngh) {


  for (int k=1; k<=ngh; k++) {
    for (int j=jl; j<=ju; j++) {
      for (int i=il; i<=iu+1; i++) {
	Real r  = pco->x1f(i);
	Real th = pco->x2f(j);
	Real z2 = pco->x3f(pmb->ks+NGHOST+k);
	Real z1 = pco->x3f(pmb->ks+NGHOST+k-1);
        // setting periodic boundary condition on the plasma generated B_r 
	// Attempt1
        //b.x1f(ku+k,j,i) =  -(A2(r, th, z2) - A2(r, th, z1))/pco->dx3f(kl+k) + b.x1f(kl+k, j, i);
	// Attempt2
        b.x1f(ku+k,j,i) =  -(A2(r, th, z2) - A2(r, th, z1))/pco->dx3f(ku+k) + b.x1f(pmb->ks+NGHOST+k, j, i);
      }
    }
  }
  for (int k=1; k<=ngh; k++) {
    for (int j=jl; j<=ju+1; j++){
      for (int i=il; i<=iu; i++) {
        b.x2f(ku+k,j,i) = b.x2f(pmb->ks+NGHOST+k,j,i);
      }
    }
  }
  for (int k=1; k<=ngh; k++){
    for (int j=jl; j<=ju; j++) {
      for (int i=il; i<=iu; i++) {
	Real r2 = pco->x1f(i);
	Real r1 = pco->x1f(i-1);
	Real r4 = pco->x1f(i+1);
	Real r3 = pco->x1f(i);
	Real th = pco->x2f(j);
	Real z  = pco->x2f(ku+k+1);

        // Periodic BC on the plasma-generated field
        b.x3f(ku+k+1,j,i) =  b.x3f(kl+k,j,i) +  (((A2(r4, th, z) - A2(r3, th, z))/pco->dx1f(i) + A2(r3, th, z)/r3) );
      }
    }
  }

  for (int k=1; k<=ngh; k++) {
    for (int j=jl; j<=ju; j++) {
      for (int i=il; i<=iu; i++) {

	// Attempt 1 Trying periodic BC only
        //prim(IDN,ku+k,j,i) = prim(IDN,kl+k-1,j,i);
        //prim(IM1,ku+k,j,i) = prim(IM1,kl+k-1,j,i);
        //prim(IM3,ku+k,j,i) = prim(IM3,kl+k-1,j,i);
        //prim(IM2,ku+k,j,i) = prim(IM2,kl+k-1,j,i);
	
	// Attempt 2 Trying to apply no-slip in the marked region and periodic BC elsewhere
	////if (pco->x1f(i) >= 0.46) {
        ////	prim(IM2,ku+k,j,i) = -prim(IM2,ku-k,j,i);
	////}
	////else {
        ////	prim(IM2,ku+k,j,i) = prim(IM2,kl+k,j,i);
	////}
	
	// Attempt 3 trying periodic BC --- doesn't work possibly because of meshblocks
        prim(IDN,ku+k,j,i) = prim(IDN,pmb->ks+NGHOST+k,j,i);

        prim(IM1,ku+k,j,i) = prim(IM1,pmb->ks+NGHOST+k,j,i);
        prim(IM3,ku+k,j,i) = prim(IM3,pmb->ks+NGHOST+k,j,i);
        prim(IM2,ku+k,j,i) = prim(IM2,pmb->ks+NGHOST+k,j,i);

      }
    }
  }
  return;
}

namespace {

//----------------------------------------------------------------------------------------
//! \fn Real A2(const Real x1,const Real x2,const Real x3)
//  \brief A2: phi-component of the vector potential generated by two coils at z = +-z0
//---------------------------------------------------------------------------------------
Real A2(const Real x1, const Real x2, const Real x3) {
  Real d1 =  std::pow((x1 + a0),2) + std::pow((x3 - z0),2) ;
  Real d2 =  std::pow((x1 + a0),2) + std::pow((x3 + z0),2) ;
  Real d3 =  4*a0*x1;
  Real sd31 = std::sqrt(d3/d1);
  Real sd32 = std::sqrt(d3/d2);

  Real Aphi = B0*x1 + I0*std::sqrt(d1)*((2-d3/d1)*std::tr1::comp_ellint_1(sd31) - 2*std::tr1::comp_ellint_2(sd31))/(d3) +I0*std::sqrt(d2)*((2-d3/d2)*std::tr1::comp_ellint_1(sd32) - 2*std::tr1::comp_ellint_2(sd32))/(d3);

  return Aphi;
  }
}

`
but it doesn't work. I suspect that the indices of the boundaries are different(not kl, ku, ks or ke) when we use mesh blocks. Also, I couldn't find any example pgen files where one imposes the periodic boundary condition manually. I would appreciate any help.

Best regards,
Rahul

[felker]

Definitely don't use pmb->ks or pmb->ke or even NGHOST; use the passed-in function parameters kl, ku, ngh. Look to polarwedge_cc.cpp and other files in src/bvals/ for guidance.

[rahulgaur104]

Thanks, Kyle! I'll try and understand some more of the pgen files and get back to you.

[tomidakn]

Your implementation copies the boundaries within one MeshBlock, so it does not work when you have more than one MeshBlocks in the periodic direction even if you fix Kyle's point. When you have more than one MeshBlocks, Periodic boundary conditions are implemented as inter-MeshBlock communications and cannot be implemented as a user-defined boundary function.

If you really want to implement what you are thinking of, you have to hack into the code. Perhaps the easiest way is to apply v2 reflecting condition after normal periodic boundary communications are completed. It is not impossible, but not trivial. At least the code does not provide such an interface. I would rather suggest you to reconsider if you really need this boundary condition - for example, can't you do that just with reflecting or outflowing boundary conditions, or whatever you can set locally without communication?

[felker]

Good point @tomidakn, I didnt think about the inability of user-defined BCs to handle any periodicity.

[rahulgaur104]

Thanks for the reply, Tomida-san and Kyle! I think it should be possible to get around this exotic boundary condition. In fact, when I set resistivity to 0, the velocity in the bracketed region does go to zero as shown in the GIF below.
(https://imgur.com/a/mvW0LMM)
Therefore, I don't need this boundary condition anymore. This strange boundary condition is one of many ideas that I am using to debug my original problem. The simulation that I am running doesn't reach a steady state. As you can see in the GIF, the rate of change of the flow seems to decrease before increasing again.