Add vertical field vector potential and remove unify vector potential parameters - [merged] #537

Closed jmstone closed 3 months ago

In GitLab by @gnwong on Jul 28, 2023, 17:57

Merges feature/verticalfield -> master

This PR provides an implementation of vertical (and vertical-adjoint) magnetic field vector potentials in the gr_torus problem generator.

It also removes the is_sane and is_mad parameters in favor of a unified always-available set of "vector potential tuning" parameters. Standardized values for these parameters are available in the comments of the pgen/gr_torus.cpp file near the parameter-parsing code.

Physics

The code first defines the cylindrical radius cyl_radius = r * sinvartheta_ks

For cylindrical radii within r_edge (i.e., within the inner edge of the torus), the vector potential is identically zero Aphi = 0.

Outside of r_edge

Aphi = pow(cyl_radius / r_edge, potential_r_pow) * exp(-cyl_radius / potential_falloff) + offset

where offset is chosen to continuously match Aphi around r_edge.

When the potential_rho_pow parameter is not zero, the vector potential everywhere is modified as Aphi *= pow(rho/rho_max, potential_rho_pow);

Example parameters

The following set of parameters seems to launch a reasonable flow from the Chakrabarti initial condition that eventually becomes MAD. In my experience, it takes about 10,000 GM/c^3 for the normalized magnetic flux to saturate at the MAD value given the usual disk size parameters.

# chmad_vertical_2_80.athinput

<problem>
user_hist = true      # enroll user-defined history function
rho_min   = 1.0e-5    # background on rho given by rho_min ...
rho_pow   = -1.5      # ... * r^rho_pow
pgas_min  = 0.333e-7  # background on p_gas given by pgas_min ...
pgas_pow  = -2.5      # ... * r^pgas_pow
rho_max    = 1.0     # if > 0, rescale rho to have this peak; rescale pres by same factor
k_adi      = 1.0     # adiabat K, p_gas = K * rho^Gamma
r_edge     = 15.0    # radius of inner edge of disk
r_peak     = 58.0    # radius of pressure maximum; use l instead if negative
l          = 0.0     # constant ang. mom. per unit mass u^t u_phi; only used if r_peak < 0
tilt_angle = 0.0     # angle (deg) to incl disk spin axis relative to BH spin in dir of x
pert_amp   = 2.0e-2  # perturbation amplitude
mad = true                  # vector potential for magnetic field (SANE or MAD config.) currently unused for vertical fields.
vertical_field = true       # turn on vertical field functionality
potential_rho_pow  = 0.0    # modifies behavior when non-zero (default=0; see the rest of this PR)
chakrabarti_torus = true    # use chakrabarti initial condition for fluid
potential_beta_min = 100.0  # ratio of gas to magnetic pressure at maxima (diff locations)
potential_falloff = 80.0    # falloff for exponential in vector potential
potential_r_pow    = 2.0    # vector potential proportional to this power of radius

Known issues to address before merging

[x] The tilted torus machinery should be checked. The implementation has been cribbed from the original vector potential implementation, which defines ks and bl angles, applies appropriate coordinate transforms so that sin_vartheta_bl and sin_vartheta_ks are their "effective" versions, i.e., for a coordinate system that is aligned with the (potentially) tilted torus. The LogHAux(...)function operates on sin_vartheta_bl as in the original, but I derive a cylindrical radius from sin_vartheta_ks. At the end of the day this is just a choice, but it may be worth unifying choices.
[x] New output has been added during problem generation. The code that finds the "outer edge" of the torus in the pgen file writes "Found outer edge: ###" to the console. This output may need to be removed to be consistent with code and output style choices.
[x] New parameters I: The code introduces several new parameters for the vector potential including is_vertical_field and potential_rho_pow. The former switches the vertical magnetic field on/off . As of this writing, when is_vertical_field is true, the code ignores is_mad or is_sane, but since there could be a future where the vertical field could be configured to MAD/SANE flows, I have not made the flags all mutually exclusive. This is ultimately a design choice.
[x] New parameters II: The code also introduces the new parameter potential_rho_pow. When this parameter is nonzero (default = 0), the vector potential is zeroed outside of the torus (based on LogHAux(...)) and Aphi is rescaled by pow(normalized_rho, potential_rho_pow) everywhere within the disk. The motivation for this parameter is to describe the vector potential for an easy "enclosed field" geometry. Whether or not this feature "works as desired" for obvious parameter values is not yet clear.

In GitLab by @gnwong on Jul 28, 2023, 17:58

added 1 commit

d4e0289d - Clean up vertical field code in gr_torus.cpp a bit more.

Compare with previous version

In GitLab by @pdmullen on Aug 10, 2023, 09:06

Commented on src/pgen/gr_torus.cpp line 755

Is the vertical field outside the torus imperative to have the thing go MAD?

In GitLab by @pdmullen on Aug 10, 2023, 09:09

Commented on src/pgen/gr_torus.cpp line 256

potentially add a max_iter and iter that is incremented in your while loops so that a user isn't sitting around forever trying to construct a torus that has no realizable solution.

In GitLab by @pdmullen on Aug 10, 2023, 09:09

Commented on src/pgen/gr_torus.cpp line 438

??? Is it possible that the maximum gas (+ radiation) pressure is ever not inside the torus?

In GitLab by @pdmullen on Aug 10, 2023, 09:10

Please make your changes compliant with our C++ linter.

In GitLab by @pdmullen on Aug 10, 2023, 09:11

Can you drop a png of a tilted system with a vertical field (magnetic field drawn with streamlines)? We should be able to eyeball if it looks correct, but I will also take a closer look at the source code.

In GitLab by @pdmullen on Aug 10, 2023, 09:15

As discussed over Slack:

I think we should take the opportunity in this MR to remove the is_mad and is_sane parameters. These parameters were originally designed so that a user could just run a SANE and MAD via is_sane=true and is_mad=true, respectively. As you discovered, just because you specify is_mad doesn't mean that any torus (e.g., Chakrabarti) will actually go MAD. I think we should just revert to specifying the type of field configuration with tunable knobs with accompanying comments that can recommend values to achieve MAD or SANE for a variety of standard FM or standard Chakrabarti initial conditions.

In GitLab by @gnwong on Aug 17, 2023, 14:08

added 3 commits

d4e0289d...5f8e1b52 - 2 commits from branch master
d848c0f7 - Merge remote-tracking branch 'origin/master' into feature/verticalfield

Compare with previous version

In GitLab by @gnwong on Aug 20, 2023, 03:47

Commented on src/pgen/gr_torus.cpp line 755

I am not 100% certain that the field must extend beyond the torus, but that's the only class of configuration that I was able to easily and consistently go MAD. There are hints that other configurations (i.e., those with enclosed fields) may go MAD if you wait long enough.

In GitLab by @gnwong on Aug 20, 2023, 03:52

Commented on src/pgen/gr_torus.cpp line 438

This set of code (which mirrors the max_bsq_intorus code) is primarily for parallel code structure. It seemed prudent to separate out both of the "_intorus" variables since the normalization is written as the ratio of the two pressures. The alternative "ptotmax/0.5/bsqmax_intorus" felt unnecessarily heterogeneous.

In GitLab by @gnwong on Aug 21, 2023, 02:14

So I've thought about this more since we last spoke, and although I agree about removing the is_sane and is_mad parameters, I'm not convinced that this is the right place to do it. That modification could/would/will be much more disruptive to users and is, in my mind, ideologically distinct vs. vertical fields. I'd be happy to start another MR for removing those parameters (and work on getting together a solution that agrees with the plan), but I'd prefer it not be a part of this branch. Thoughts?

In GitLab by @pdmullen on Aug 21, 2023, 09:20

Commented on src/pgen/gr_torus.cpp line 755

I would be concerned that huge magnetic tension forces at the interface of torus and background would cause numerical instability, but after all, we are invoking relativistic dynamics....

In GitLab by @pdmullen on Aug 21, 2023, 09:25

Commented on src/pgen/gr_torus.cpp line 256

I'd prefer different variable names over extra braces if they were introduced due to a scope issue.

In GitLab by @pdmullen on Aug 21, 2023, 09:27

Commented on src/pgen/gr_torus.cpp line 438

I guess...

For the case that the maximum pressure is inside the torus, then we aren't getting too much more extra information from your extra reduction. How about we reduce to find the maximum pressure inside the torus and maximum pressure outside the torus?

Then outside the kernel and after the MPI reduction you can do something like

const Real pmax_tot = (is_vertical_field) ? std::max(pmax_interior, pmax_exterior) : pmax_interior

In GitLab by @pdmullen on Aug 21, 2023, 09:37

MRs take time... I think we should at least include some fix right now.

I propose the simplest fix:

Keep is_sane and is_mad parameters for now, but when is_chakrabarti && is_mad, can you point the code to a configuration known to go MAD (i.e., is_vertical_field = true). Similarly, make some choice about where the code defaults for is_sane (albeit there are choices are here...)

Also, can you cleanup the is_mad and is_sane associations with FM tori? Currently, you must specify is_mad or is_sane BUT ALSO potential_r_pow and potential_falloff. This is bad design (and my fault!).

I trust these are simple fixes though... (and the inputs/grmhd/ mad and sane inputs define the appropriate designations for the standard FM tori).

In GitLab by @gnwong on Aug 21, 2023, 16:18

Commented on src/pgen/gr_torus.cpp line 256

changed this line in version 4 of the diff

In GitLab by @gnwong on Aug 21, 2023, 16:18

added 1 commit

21514c4d - change while for for in code that finds the outer edge of the torus to avoid...

Compare with previous version

In GitLab by @gnwong on Aug 21, 2023, 16:20

Commented on src/pgen/gr_torus.cpp line 256

The first part of this (iter and max_iter) are addressed by 21514c4d6629673c85aca414f9986c3e6ce6dac4. The braces are not for a scope issue but rather conceptual separation, although I agree it might be more in line with code style to remove them. Will do in the next commit.

In GitLab by @gnwong on Aug 21, 2023, 16:23

Commented on src/pgen/gr_torus.cpp line 755

I see your point, but this (having external field) is at least not novel. It appears in other initial conditions for vertical fields, including, e.g., Scepi+ 2023 (2302.10226).

In GitLab by @gnwong on Aug 21, 2023, 16:26

Commented on src/pgen/gr_torus.cpp line 438

I would be happy with something that gives pmax_interior and pmax_exterior, but I'm hesitant to then set pmax_tot based on something like is_vertical_field since the differences between interior/exterior might be significantly influenced by other creative fluid initial conditions. Thoughts?

In GitLab by @pdmullen on Aug 21, 2023, 18:51

Commented on src/pgen/gr_torus.cpp line 438

I don't have a strong opinion. You can change it to something like !(in_torus_beta_init_mask) but that is starting to feel super clumsy.

In GitLab by @gnwong on Aug 22, 2023, 11:09

added 1 commit

b75599ca - remove braces around gr_torus edge-finding routine and lint code

Compare with previous version

In GitLab by @gnwong on Aug 22, 2023, 11:10

Commented on src/pgen/gr_torus.cpp line 256

Braces removed with commit b75599ca0a791340ec443aa037a185c02c05306e. Should be ready to go (from the perspective of this thread).

In GitLab by @gnwong on Aug 23, 2023, 15:07

added 1 commit

54a26c26 - remove is_mad and is_sane flags from gr_torus parameter file.

Compare with previous version

In GitLab by @gnwong on Aug 23, 2023, 15:13

Adding in defaults for is_sane and is_mad in the vertical field configuration would just make more work for the future so I removed them entirely instead as requested at the top of this thread. I also added adjusted the "regular" field configuration to accept the potential_rho_pow parameter. I'm also modifying the original text of the MR and its title.

Standard defaults for the SANE/MAD vector potential (for FM, which is established) are included near the parameter parsing code. The only nuance is that the standard SANE vector potential does not include an exponential cutoff, so the potential_falloff parameter has to be set to "very large" (or alternatively another parameter must be introduced, which would be even uglier in my opinion).

Thoughts?

Addressed in commit 54a26c26d29beb5a67129be45952b176caf22860

In GitLab by @gnwong on Aug 23, 2023, 15:28

Commented on src/pgen/gr_torus.cpp line 438

changed this line in version 7 of the diff

In GitLab by @gnwong on Aug 23, 2023, 15:28

added 1 commit

c6dc4591 - remove ptotmax_intorus in favor of single ptotmax

Compare with previous version

In GitLab by @gnwong on Aug 23, 2023, 15:31

Commented on src/pgen/gr_torus.cpp line 438

I've just removed ptotmax_intorus as originally requested. I'm not particularly enthusiastic about this change, but I agree it is the minimal (i.e., least redundant code) solution given the current implementation of the torus initial conditions.

It will probably have to be revisited in the future, but until those models actually exist it probably won't be possible to cleanly design for them.

Addressed in commit c6dc4591f393b2898e832d9cc0d61cf7ea58fe05 and ready for your review again, Patrick.

In GitLab by @pdmullen on Aug 23, 2023, 20:56

Commented on src/pgen/gr_torus.cpp line 1202

this is highly technical, and it has been a while since I looked at this... what is the concern here? i.e., is there a suspicion that the original tilt infrastructure was wrong?

In GitLab by @pdmullen on Aug 23, 2023, 21:01

Commented on src/pgen/gr_torus.cpp line 1249

orus.potential_falloff at 1.e30 seems weird for configurations where we don't want to invoke a potential falloff. Not sure my suggestion is any better?

torus.potential_falloff = pin->GetOrAddReal("problem", "potential_falloff", 0.0);

scaling_param *= (torus.potential_falloff > 0.0) ? exp(-r/pgen.potential_falloff)) : 1.0;

and similar elsewhere...

In GitLab by @pdmullen on Aug 23, 2023, 21:04

I am happy with this change, modulo using large values for the potential falloff... I left some comments elsewhere about this.

In GitLab by @pdmullen on Aug 23, 2023, 21:04

Commented on src/pgen/gr_torus.cpp line 1245

Jim disallows capital variable names.

In GitLab by @gnwong on Aug 25, 2023, 07:27

Commented on src/pgen/gr_torus.cpp line 1202

I think I am confused by the conventions for tilt. As I understand it, tilted tori are handled by defining a sin_vartheta value, which represents "sin(theta)" with theta defined relative to the intended polar axis of the tilted system.

In the fluid initialization code (e.g., where LogHAux is called), CKS x,y,z -> BL r,theta,phi, and the tilt is then defined in terms of BL. Then the torus is tilted about the line defined by theta=pi/2 with phi_bl=pi/2 or 3pi/2. Notice that this axis is different compared to if phi_bl -> phi_ks.

When computing the vector potential and when we rely on things like rho, we certainly want to be computing them at a point (x,y,z) in the same way that they're computed above for the fluid initialization. So then certainly sin_vartheta_bl should be passed into (e.g.) LogHAux (as it is ... see https://gitlab.com/theias/hpc/jmstone/athena-parthenon/athenak/-/blob/5f8e1b52dd32300450c57df11438a850aac647a3/src/pgen/gr_torus.cpp#L1159). But then shouldn't we also compute the other parts of the potential in terms of vartheta_bl (cf https://gitlab.com/theias/hpc/jmstone/athena-parthenon/athenak/-/blob/5f8e1b52dd32300450c57df11438a850aac647a3/src/pgen/gr_torus.cpp#L1170)?

In the code linked above, it seems like the fluid is rotated about phi_bl = [fixed] in the midplane whereas the vector potential component that does not depend on the fluid is rotated about phi_ks = [fixed]. So you're rotating the field around an axis that differs from the fluid rotation axis with some radial dependence.

I guess my default here would've been to rotate about phi_ks = [fixed] since phi_bl = [fixed] doesn't makes sense as you approach the horizon anyway. And then you would rotate the field around phi_ks = [fixed] to be consistent.

But am I missing something?

In GitLab by @gnwong on Aug 25, 2023, 07:28

I'm holding off on doing this until the other issue below is resolved (w.r.t. rotating about phi_bl or phi_ks = fixed in the midplane).

In GitLab by @gnwong on Aug 25, 2023, 07:37

Commented on src/pgen/gr_torus.cpp line 1249

changed this line in version 8 of the diff

In GitLab by @gnwong on Aug 25, 2023, 07:37

Commented on src/pgen/gr_torus.cpp line 1245

changed this line in version 8 of the diff

In GitLab by @gnwong on Aug 25, 2023, 07:37

added 1 commit

73f10761 - update variable style, use more sensible defaults for potential parameters,...

Compare with previous version

In GitLab by @gnwong on Aug 25, 2023, 07:38

Commented on src/pgen/gr_torus.cpp line 1249

Addressed in 73f107619dd7d24e093d6f42ae271e01d22e121a. Using default value of potential_falloff==0.0 now turns off the exponential cutoff. Ready for review.

In GitLab by @gnwong on Aug 25, 2023, 07:38

Commented on src/pgen/gr_torus.cpp line 1245

Addressed in commit 73f107619dd7d24e093d6f42ae271e01d22e121a. Ready for review.

In GitLab by @gnwong on Aug 25, 2023, 07:38

Addressed in commit 73f107619dd7d24e093d6f42ae271e01d22e121a and commented on in the other thread. Ready for review.

In GitLab by @pdmullen on Aug 28, 2023, 10:01

Commented on src/pgen/gr_torus.cpp line 1202

@c-white, can you take a look?

In GitLab by @pdmullen on Aug 28, 2023, 10:02

this is still failing @gnwong

In GitLab by @gnwong on Aug 28, 2023, 10:06

Yes, this is intentional :) The only (persistent) failure is a too-long comment about the field rotation part of the code that I've left in to make sure the CI fails until we resolve everything.

In GitLab by @gnwong on Aug 28, 2023, 18:42

Commented on src/pgen/gr_torus.cpp line 1202

Hey @c-white, happy to provide any more missing context. But I think it boils down to:

1) Is there a compelling reason to rotate around fixed BL phi? 2) Is it intended that the vector potential be rotated around fixed KS phi while the fluid is rotated around fixed BL phi (or maybe I am missing something)?

In GitLab by @c-white on Aug 29, 2023, 14:08

Commented on src/pgen/gr_torus.cpp line 1202

There was no particular reason for any of the choices, including explicitly using any BL coordinates (other than FM being written in BL).

For A tracing rho and/or pgas, the fluid and the field really should be tilted in the same way, whatever that is. For A specified independent of the torus, they should probably be the same.

I never worried too much for two reasons.

(1) Tilted tori are already out of equilibrium, and will undergo rapid adjustment once the simulation starts.

(2) All this is approximate and ambiguous. Does tilting look like the Euclidean definition in SKS, SBL, CKS, or something else? Does tilt preserve the constancy of u_ph u_BL^t from an untilted FM, or the constancy of -u_phi / u_t from AJS? Does a constant tilt angle of i result in constant i = arccos(L^z/L) (and for which L?), or a constant latitude i of fluid elements' antinodes, or a constant angle i = arctan(u^phi / u^theta) along the line of nodes? All these are reasonable expectations, but the options in each set are mutually exclusive.

So with all that, I'm okay with whatever changes can be made to make the conventions more internally consistent.

In GitLab by @gnwong on Aug 30, 2023, 06:30

Commented on src/pgen/gr_torus.cpp line 1202

changed this line in version 9 of the diff

In GitLab by @gnwong on Aug 30, 2023, 06:30

added 1 commit

ed97ac18 - unify vector potential coordinates with fluid coordinates for tilted tori (BL...

Compare with previous version

In GitLab by @gnwong on Aug 30, 2023, 06:32

After fighting with VisIt for several hours, I've produced the following two screenshots of field lines in tilt=0 and tilt=30 Chakrabarti tori with the "vertical field" configuration as specified above and below. Looks good to me, so I'd count this as resolved if you agree @pdmullen

Screen_Shot_2023-08-30_at_6.29.13_AM

Screen_Shot_2023-08-30_at_6.29.00_AM

In GitLab by @gnwong on Aug 30, 2023, 06:34

Commented on src/pgen/gr_torus.cpp line 1202

Thanks both. I agree with Chris about the initial configuration (which only differs close to the hole) not mattering so much. The minimal change here was to use (tilted) BL coordinates for both the geometric and fluid components of the vector potential, so that's what I've done.

Addressed in commit ed97ac1818b7a91cf043e0873559fd9425a79345 and ready for review.

In GitLab by @gnwong on Aug 30, 2023, 06:35

marked this merge request as ready

IAS-Astrophysics / athenak