Magnetostatic Equation Solver for Circular Domain with Copper Coil

I am currently working on solving the magnetostatic equation within a circular domain filled with air, containing a copper coil with a given current density, J_coil. The challenge arises from the differing permeabilities of air and copper, requiring me to adapt my PDE implementation for accurate results.

My current implementation includes two approaches, both of which have encountered specific challenges.
The first approach involves defining a piecewise or Heaviside function inside the PDE that evaluates to 1 within the copper coil and 0 outside. In this approach, I multiply the function by J_coil * permeability of copper inside the coil region, and outside, only the Laplacian term remains as J_noncoil = 0. However, I have observed that the model evaluates this condition only once during the creation of the class object. Consequently, the results I obtain are equivalent to having the entire region filled with copper (which is accurate when compared to FEMM Simulation data for the entire domain made of copper)

The second approach I explored involves defining two separate domains for the copper coil and air. Although this approach entails creating two networks, I have encountered stability issues, and the model fails to converge. On the other hand, the single-domain model converges rapidly and provides accurate results.
Below is a sample code snippet that demonstrates how I am currently implementing the magnetostatic equation within a single domain:

class MagnetostaticsEquation2D(PDE):
    name = "MagnetostaticsEquation2D"
    def __init__(self, nu_air, nu_copper, J_copper):
        # coordinates
        x = Symbol("x")
        y = Symbol("y")

        # make input variables
        input_variables = {"x": x, "y": y}

        # Magnetic vector potential (A) output
        A = Function("A")(*input_variables)
        # Set the equation based on the permeability and current density
        self.equations = {}
        self.equations["magnetostatics_equation"] = A.diff(x, 2) + A.diff(y, 2) + nu_copper * J_copper * heaviside_square(x, y)

I have also introduced non-dimensionalization and scaling, although I have not included those details in this post.

I am reaching out to this community for guidance and advice on resolving these challenges.