Is there a bug in resistive MHD?

[abransgrove]

Hi everyone,

I am playing around with a simple Harris sheet in resistive MHD, and I noticed some strange behavior of the current sheet near the boundaries. This motivated me to understand how the non-ideal MHD works at the edge of the grid. When I read diffusivity.cpp the ohmic diffusivity is set as:

  if (pfdif->eta_ohm > 0.0) { // Ohmic resistivity is turned on
    for (int k=ks; k<=ke; k++) {
      for (int j=js; j<=je; j++) {
#pragma omp simd
        for (int i=is; i<=ie; i++)
          pfdif->etaB(FieldDiffusion::DiffProcess::ohmic, k,j,i) =
              pfdif->eta_ohm;
      }
    }
  }

As far as I can tell, the diffusivity is cell centered, while the electric fields are defined at the cell edges. When the resistive electric field (E=eta*J) is calculated, eta is interpolated to the cell edges. However, as you can see from the loop bounds above, eta is not defined in the ghost cells. Does this mean that the interpolation of eta is wrong at the edge of the grid?

[felker]

At first this appeared to be a bug to me, but the function is called in field_diffusion.cpp in a way that extends the loop limits to all ghost cells:

void FieldDiffusion::SetDiffusivity(const AthenaArray<Real> &w,
                                    const AthenaArray<Real> &bc) {
  MeshBlock *pmb = pmy_block;
  int il = pmb->is - NGHOST; int jl = pmb->js; int kl = pmb->ks;
  int iu = pmb->ie + NGHOST; int ju = pmb->je; int ku = pmb->ke;
  if (pmb->pmy_mesh->f2) {
    jl -= NGHOST; ju += NGHOST;
  }
  if (pmb->pmy_mesh->f3) {
    kl -= NGHOST; ku += NGHOST;
  }

  for (int k=kl; k<=ku; ++k) {
    for (int j=jl; j<=ju; ++j) {
#pragma omp simd
      for (int i=il; i<=iu; ++i) {
        Real Bsq = SQR(bc(IB1,k,j,i)) + SQR(bc(IB2,k,j,i)) + SQR(bc(IB3,k,j,i));
        bmag_(k,j,i) = std::sqrt(Bsq);
      }
    }
  }
  // set diffusivities
  CalcMagDiffCoeff_(this, pmb, w, bmag_, il, iu, jl, ju, kl, ku);

  return;
}

Did you print out the interpolated eta in the EMF calculations to confirm that the constant diffusivity is messed up? Can you share your test results?

[abransgrove]

Thanks Kyle, I just realized this as well. I have not tried printing the interpolated eta, but this is a good suggestion and I will try it (although I am now convinced eta is correct).