[WIP] Add Vloop BC

[theo-brown]

Adding Vloop boundary condition on the psi equation, so that the current can be evolved self-consistently given a rate of change of edge flux.

I have:

• Added an optional config argument, config["profile_conditions"]["Vloop_bound_right"] for the potentially time-varying Vloop at LCFS.
  Made config["profile_conditions"]["Ip"] and config["profile_conditions"]["Vloop_bound_right"] mutually incompatible (with unit test)
• Added a right_face_constraint to the Psi element of compute_boundary_conditions, set by Vloop_bound_right

This necessitated the following changes:

• Pass dt into compute_boundary_conditions
• Minor renamings for clarity (eg t - dt, t, and t + dt now all get used at various points)

[theo-brown]

I have not yet written a test case for this, as I wanted to run the pattern past one of the team first. Let me know your initial thoughts and then I can continue!

[jcitrin]

Some general thoughts.

1. What's missing is a consideration of how to initialize psi to begin with, for all cases. If Ip is None, then how do we set the right boundary condition at the first timestep, for cases where we initialize based on plasma current? Recall that there are several options currently (which can be further extended later):

i) profile_conditions['initial_psi_from_j'] = False, i.e. Initialize psi based on the input geometry file. Two suboptions:
a) geometry['Ip_from_parameters'] = False.
Here's it's clear what to do, psi_right = psi_geometry[-1]
b) geometry['Ip_from_parameters'] = True. For this we need Ip to rescale psi, so if Ip=None we have a problem

ii) profile_conditions['initial_psi_from_j'] = True, i.e. Initialize psi based on the "nu" current profile model. This needs normalization to a total current, i.e. Ip. So if Ip=None we have a problem. Note that at the initial condition, the Vloop_bound_right needs to be such that Psi.face_grad()[-1] is consistent with Ip, according to the same formula as the standard boundary condition. See fvm.cell_variable.face_grad, that face_grad()[-1] is just a standard derivative between the dx/2 distance between the last cell point and the rightmost face.

To solve this, I recommend changing the name of Ip to Ip_input, and for cases where the Vloop boundary condition is used, then it's used (if necessary) to initialize psi. That means that Vloop_bound_right and Ip can coexist. The pattern should be a bit different, e.g. an extra boolean "use_Vloop_boundary_condition" or something like that. Note that Ip_input can still be a TimeInterpolatedInput, but for the Vloop case only the initial condition at t=t_initial will be used.

2. We need to adjust to the Vloop boundary case where profile_conditions['Ip_input'] is not actually total current, but rather the initial condition. For example, the output Ip is taken from there, and also in plotruns_lib.PlotData.i_total. For the code to work with both boundary condition cases, we should take Ip (total current) from elsewhere. A candidate could be to extend physics.calc_jtot_from_psi to output I_tot as well (can change the name to Ip_profile), and then save this quantity in CoreProfiles.Currents . Total current is then Ip_profile[-1]. This can be a separate PR. This way, we also get information on the total current evolution in the Vloop boundary condition case.

[theo-brown]

Great points, and this is exactly why I wanted to run it through you! I think because for the type of simulation I've learnt how to do we always prescribe psi from a geometry I hadn't really thought about the consequences for any situation other than 1ia).

Noting down for myself - I need to consider the cases:

	Psi from...	Ip
1ai	Geo	Self-consistent
1aii	Geo	Prescribed
2	Initial j profile	Prescribed

and their interactions with the Vloop BCs

[jcitrin]

#492 was done with this PR in mind