[d3] Spherical shell - potential boundary conditions

[csskene]

Hi
I've seen a recent commit say that test_shell_dynamo.py and test_BVP_shell.py do not seem to currently work. I believe this is in part to do with the potential boundary conditions on the inner sphere. In order to match the vector potential A to an electrically insulating BC for B I think the inner boundary should be changed (for example in test_BVP_shell.py) to

# Potential BC on A
ell_func_out = lambda ell: ell + 1
ell_func_in = lambda ell: ell
radial = lambda A: operators.RadialComponent(A)
ell_op_out = lambda A: operators.SphericalEllProduct(A, c, ell_func_out)
ell_op_in = lambda A: operators.SphericalEllProduct(A, c, ell_func_in)
A_bc_outer = radial(grad(A)(r=r_outer)) + ell_op_out(A)(r=r_outer)/r_outer
A_bc_inner = radial(grad(A)(r=r_inner)) - ell_op_in(A)(r=r_inner)/r_inner

This seems to give much less error in test_BVP_shell.py, and also prevents test_shell_dynamo.py from quickly blowing up. I think the rest of the error stems from the fact that the magnetic field B used for test_BVP_shell.py does not satisfy potential BCs for the radii given, for example the toroidal part is not zero at r=r_o and r=r_i.

[csskene]

Further to my previous message, I have made some progress with getting test_shell_dynamo.py to work. As well as the boundary condition changes I mentioned before, I also changed the initial condition from an LBVP to obtain A from B, to a direct initial condition for A

A['g'][0] = 5/8*(-3*r**2+4*r_outer*r-r_inner**4/r**2)*np.sin(theta)
A['g'][2] = 5/2*r*np.sin(np.pi*(r-r_inner))*np.cos(2*theta)

whose curl is the B given. I also made some minor changes such as changing the timestepper to SBDF2, adding dealiasing and changing the formulation to use LiftTau. Whilst this initial condition does not satisfy the Coulomb gauge condition this doesn't seem to matter and I get the magnetic energy converging to 626.303 and the kinetic energy to 30.77, which agrees with the Christensen et al. benchmark values. If it's helpful I'm happy to tidy up the script I got to work and submit a pull request?

[kburns]

Thanks! I think we're planning to move most of these examples out of this repo and into the "dedalus-spheres" repo actually, so it might make sense to wait until that's done to put in the PR. Let's follow up next week.

[csskene]

Hi, I've seen that the test_shell_dynamo.py and test_BVP_shell.py scripts have now moved from the Dedalus repo. A quick update on this issue; for the current version of Dedalus v3 I have written a script that fixes both the problems I commented on before (the initial B does satisfy the BCs so disregard what I said before). So put together, a script that solves for the initial condition for the vector potential A from B, and then solves the case1 benchmark from Christensen et. al. Would this be something that would be good to reshape into an example for Dedalus v3? I was at a dynamo workshop last week and a few people there would be interested in this. Otherwise I'm happy for the issue to be closed since these scripts are no longer in the Dedalus repo. Thanks!

[csskene]

Happy to close this issue since the scripts it concerns have been moved out of the Dedalus repo. Also just to mention, @bpbrown has working versions of similar scripts in his shell benchmark repo here. Also happy for it to be reopened in case an example similar to these scripts is wanted for the Dedalus repo.

[kburns]

Thanks, we have a separate issue #221 to remind us to add examples/docs on this.