Directly specifying electric field boundary condition

[kyriakosTapinou]

Has anyone implemented an electric field boundary condition through directly specifying electric field values e.g. in corners within the CT algorithm? There seems to be direct access to the electric field register i.e. pmb->pfileld->e.x1e etc. but I am wondering if this is overridden by the CT algorithm? I couldn't see a flow diagram including CT in the programmer guide. I suppose it is done as part of the flux solve routine and so it will override whatever is put in the boundary function?

[c-white]

Short answer: The values will be overwritten.

The MHD diagrams here are applicable: https://github.com/PrincetonUniversity/athena/wiki/Code-Structure

The face-centered electric fields are calculated via the Riemann solver in CalculateFluxes (now named CalculateHydroFlux). These are next moved to edges in CalculateEMF. Then they are used to update the face-centered magnetic fields in FieldIntegrate (now named IntegrateField). The boundary function call, PhysicalBoundary, never falls in between any of these steps.

The diagrams are a bit out of date (though it doesn't affect this issue). You can always check the actual dependencies of the tasks in src/task_list/time_integrator.cpp. There you'll find wrapper functions named CalculateHydroFlux, CalculateEMF, IntegrateField, and PhysicalBoundary, which call more involved member functions of the appropriate classes. These have corresponding shorthand task names CALC_HYDFLX, CALC_FLDFLX, INT_FLD, and PHY_BVAL. If you search the file for these, you can see their immediate dependencies in the second argument to AddTask. For example, AddTask(INT_FLD,RECV_FLDFLX); means magnetic field integration can only be executed once electric fields have been received from and synced with any neighboring MeshBlocks.

(Crawling through the dependencies can be a bit tedious. This is perhaps why no one has had the time to write a robust way of diagramming the task flow automatically, and I haven't manually remade the diagrams in quite a long time.)

Modifying electric fields, like modifying hydrodynamical fluxes, is not particularly supported by the code. That said, if you really want, this should be doable, and perhaps has been done by others, but it will require writing some code outside the usual safe places for user modifications.

[msbc]

What is your pressure BC doing?

[kyriakosTapinou]

@msbc I've tried a couple different standard options from const, lin extrap, zero grad etc. I think the problem may be fundamental issue with the scenario simulated but I was wondering if there was something in the chain the of magnetic field update (face -> cell -> e_edge -> e_corner->face) and total energy conservation that I may have overlooked.

[c-white]

It sounds like this is the same issue that one would have if one were to overwrite the boundary fluxes -- no matter what you put there, conservation laws are satisfied, but you might force the system to evolve to an unphysical state (by, e.g., giving a cell lots of momentum but no energy).

When you overwrite E, the code will guarantee that div B doesn't change and in fact that the evolution of B satisfies the induction equation. Normally, however, E is related to B-fluxes across faces, which are inseparably joined to energy and momentum fluxes via the Riemann solver. By overwriting E, you might force B to evolve inconsistent with total energy evolution, even though both are evolving realistically in isolation.

Put another way, externally applying an EMF also means externally applying an energy flux (and maybe a momentum flux? certainly in relativity it does), and those should be accounted for. Perhaps this has a simple form.

A more tedious and untested approach that occurs to me is this: Consider the set of ideal MHD fluxes and try to find a set of them that is self-consistent (i.e. can be derived from a valid primitive state) while obeying the electric fields you impose in the vicinity.

$$ \begin{align} F^x(d) & = \rho v^x \\ F^x(e) & = (\rho \vec{v}^2 / 2 + \Gamma p / (\Gamma - 1) + \vec{B}^2) v^x - (\vec{B} \cdot \vec{v}) B^x \\ F^x(m^i) & = \rho v^i v^x + (p + \vec{B}^2 / 2) \delta^{ix} - B^i B^x \\ -E^z = F^x(B^y) & = B^y v^x - v^y B^x \\ +E^y = F^x(B^z) & = B^z v^x - v^z B^x \end{align} $$

One could imagine fixing $F^x(d)$, $E^y$, and $E^z$ (the latter set by averaging the adjacent desired EMFs to the face) and solving for $F^x(e)$ (and possibly $F^x(m^i)$ ) by fixing sufficiently many primitives on the face (almost certainly $B^x$, maybe $\rho$, $v^y$, and $v^z$ as needed?).

[kyriakosTapinou]

@c-white Yes I agree.

For context, the method (which works for isothermal, pointing me to energy flux issues with my present adiabatic implementation) works as follows.

Edit in the first case when testing isothermal I only set magnetic field and velocity values in the ghost cells, I wasn't messing with fluxes.

1. Read external HDF5 files at initialization - these hold the B and E information
2. Apply in the ix3 boundary conditions the relevant properties - for example the face centered magnetic field and the cell centered velocity that "drives" the simulation, as well a generic values of rho and p.
3. In hydro/calculate_fluxes.cpp

• overwrite values of (a) face centered magnetic field values for the z-boundary ks to the magnetic field desired i.e. what the soon to be applied E values will generate,
• (b) overwrite the cell-centered magnetic field values in the firs ghost cell (ks-1) and the first interior cell (ks) to 0.5*(b.x3f(k) + b.x3f(k-1) )
• for ks-1, use DonorCellX3 reconstruction so the values from ghost cell are applied at wr and wl reconstructions,
• for ks, use DonorCellX3 reconstruction, BUT then overwrite the wr_ register with the ks-1 cell values for rho, v, p, and overwrite the wr_ IBY and IBZ values to IB1 and IB2 from our ks-1 cell-centered values.

4. Within calculate_corner_e.cpp I overwrite the ex and ey values on the z-boundary i.e. ks=0 values of edge centered e1 and e2

I think this covers all the bases in terms of which values the ghost cell driving would affect it were not using E overwrite (i.e. I'm trying to make the rest of the code like the energy flux operate with what should be the consistent magnetic field values).

@c-white Are you suggesting calculating these flux values (that would correspond to the electric field being forced) and overwrite the riemann solve that otherwise assigns them?

[c-white]

Yes, I suppose I am suggesting overwriting the fluxes. The Riemann solver will only give the correct energy fluxes if given the correct inputs, and it's not being given the correct exterior state (i.e. the one responsible for the EMFs you are imposing). This is essentially figuring out what Riemann problem would give you the EMFs you want, and then asking what energy flux such a problem would have yielded.

[kyriakosTapinou]

@msbc do you recall if you had any issues with AMR when you were overwriting the electric field to zero?

When I have been using AMR with (especially with small block sizes) I am seeing what look like block boundary artifacts. Plotting the elecrtric field values from the register (e.x1e etc.) and comparing to my "desired" values, there are deviations even though the former is specified directly from the latter with the exact same data extents.

@c-white I am thinking this could be something to do with block synchronisation, or refinement and derefinement. Is there any danger that edge-centered electric field values specified according to Matt's suggestion (below) could be altered after the ComputeCornerE? There doesn't seem to be anything odd going on in the TimeIntegratorTaskList::IntegrateField or the Field::CT functions, but I'm unfamiliar with the stage weights and how they affect the field update on shared boundaries with different refinements.

`if (pmb->loc.lx1==0) {

for (int k=ks; k<=ke; ++k) {

  for (int j=js; j<=je+1; ++j) {

    e3(k,j,is) = 0.0;

  }

}

for (int k=ks; k<=ke+1; ++k) {

  for (int j=js; j<=je; ++j) {

    e2(k,j,is) = 0.0;

  }

}

}`