Solving for a FRC

[jtpaquet]

I'm trying to solve for the poloidal flux for a FRC. For this, I have three coils with constant current that create a vertical magnetic field. I was trying to define some isoflux and x-points according to this drawing hoping that it would find an O-point at the specified position.

I tried following the tutorial with my own dimensions but because the psi boundary is set to 0, the equilibrium is found in the center of the domain.

I also get some warnings when I try to evalutae psi at r=0. I can get rid of some of them by specifying my radial domain from 0.001 to 1. However, I think it would be important to include 0 because the separatrix would be along the axis r=0.

I can also add a plasma current using a coil in the opposite direction inside the chamber to get psi profiles that somewhat look like an FRC.
coils.append( ("Ip", Coil(R/2, 0, current=-I/2, turns=N_turns)) ) # Plasma Current
Here the lein R=0 is nan.

How could I change my script to get adequate equilibrium and model an FRC?

Here is what I tried:

from freegs.machine import Coil, ShapedCoil, Wall, Machine, Circuit
import freegs

import numpy as np
import matplotlib.pyplot as plt

### Machine and coil dimensions
R = 0.30 # m
L = 0.60 # m
I = 10 # A
c_w = 0.02 # Coil width, m
N_turns = 5

coils = []
for i in [-1,0,1]:
    x0, z0 = np.sqrt(4/7)*R, np.sqrt(3/7)*R * i
    if i == 0:
        x0 = R
    z_label = {-1:'D', 0:'C', 1:'U'} # Down, Center, Up
    coil = ("VC"+z_label[i], Coil( x0, z0, turns=N_turns, control=True, current=I) )
    coils.append(coil)

### Show magnetic field
# Make a 2D grid of R, Z values
# Note: Number of cells 65 = 2^n + 1 is useful later
R, Z = np.meshgrid(np.linspace(0.0, 1.0, 65), np.linspace(-0.5, 0.5, 65), indexing='ij')

psi = np.zeros(R.shape)
B_r = np.zeros(R.shape)
B_z = np.zeros(R.shape)

for coil in coils:
    psi += coil[1].psi(R,Z)
    B_r += coil[1].Br(R,Z)
    B_z += coil[1].Bz(R,Z)

# Added nan_to_num because the half opposite to the coil is nan

B_p = np.sqrt(B_r**2 + B_z**2)

plt.contour(R, Z, psi, 40)
plt.contourf(R, Z, B_p, 40)
plt.xlabel("Major radius R [m]")
plt.ylabel("Height Z [m]")
plt.gca().set_aspect('equal')
plt.show()

#########################################
# Create the machine, which specifies coil locations
# and equilibrium, specifying the domain to solve over

frc_chamber = Machine(coils)

eq = freegs.Equilibrium(tokamak=frc_chamber,
                 Rmin=0.0, Rmax=1.0, # Radial domain
                 Zmin=-0.5, Zmax=0.5, # Height range
                 nx=65, ny=65, # Number of grid points
                 boundary= freegs.boundary.freeBoundaryHagenow)

#########################################
# Plasma profiles

# profiles = freegs.jtor.ConstrainPaxisIp(eq, paxis=1e4, Ip=1e6, fvac=0.0)
profiles = freegs.jtor.ConstrainBetapIp(eq, betap=0.5, Ip=1e6, fvac=0.0)

#########################################
# Coil current constraints
#
# Specify locations of the X-points
# to use to constrain coil currents

xpoints = [(0.0, -R/2), (0.0, R/2)] # (R,Z) locations of X-points
isoflux = [(0.0,0, R/2,0)] # (R1,Z1, R2,Z2) pairs
constrain = freegs.control.constrain(xpoints=xpoints, isoflux=isoflux)

#########################################
# Solve

freegs.solve(eq, # The equilibrium to adjust
             profiles, # The toroidal current profile function
             constrain, # Constraint function to set coil currents
             psi_bndry=0.0, # Because no X-points, specify the separatrix psi
             show=True)