Parameterized Airfoil Flow Network Won't Converge

I am making a PINN with Modulus Sym to predict the flow around an airfoil. The airfoil is represented as a parametric curve. I have been training my model but the loss hovers around 10^7, sometimes spiking but never falling below 10^7. Here is the code for my PINN:

Making the network

airflow_net = FullyConnectedArch(
        input_keys = [
            Key("B"),
            Key("T"),
            Key("P"),
            Key("C"),
            Key("E"),
            Key("R"),
            Key("alpha"),
            Key("u_flight"),
            Key("Re"),
            Key("x"),
            Key("y")
        ],
        output_keys = [
            Key("u"),
            Key("v"),
            Key("p")
        ],
        nr_layers = 5,
        layer_size = 512
    )

    ns = NavierStokes()

    nodes = ns.make_nodes() + [airflow_net.make_node(name="airflow_net")]

    #Set up parameters
    input_variables = ["B", "T", "P", "C", "E", "R", "alpha", "u_flight", "Re"]
    x, y = Symbol("x"), Symbol("y")

Geometry & Constraints

#Freestream conditions
    airflow_domain = Domain()

    farfield = 5 # Farfield is defined as 5c above and below
    inlet = -5 # Inlet is 5c ahead of airfoil
    outlet = 10 # Outlet is 10c behind airfoil

    # Pressure at 30k feet altitude = 1395Pa
    atm = 1395

    computational_domain = Rectangle((inlet, -farfield), (outlet, farfield))

    freestream = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=computational_domain,
        outvar={"u": symbols["u_flight"] ,"v": 0, "p": 1395.0},
        batch_size=cfg.batch_size.Freestream,
        parameterization=parameterization
    )

    airflow_domain.add_constraint(freestream, "freestream")

    # Airfoil noslip condition

    #Parametric equations defining geometry
    def X(theta):
        return 0.5 + 0.5 * (
                sym.Abs(
                    sym.cos(theta)
                ) ** symbols["B"]
                / sym.cos(theta))

    #Generate Y value from theta
    def Y(theta):
        x = X(theta)

        y = symbols["T"] / 2
        y *= sym.Abs(sym.sin(theta)) ** symbols["B"] / sym.sin(theta)
        y *= 1 - x ** symbols["P"]
        y += symbols["C"] * sym.sin(sym.pi * x ** symbols["E"])
        y += symbols["R"] * sym.sin(2 * sym.pi * x)

        return y

    theta = Symbol("theta")
    theta_range = (1e-3, 2 * math.pi)

    airfoil_x = X(theta)
    airfoil_y = Y(theta)

    functions = {
            "x": airfoil_x,
            "y": airfoil_y,
            "normal_x": 1e-3, 
            "normal_y": 1
    }

    print("Making airfoil curve")
    airfoil_curve = SympyCurve(
        functions=functions,
        area=symbols["T"] + 2.05, #Approx area
        parameterization=theta_parameterization
    )

My SDF function

def sdf(*args, **kwargs):
        x_arr = args[0]["x"].flatten()
        y_arr = args[0]["y"].flatten()

        B_arr = args[1]["B"].flatten()
        T_arr = args[1]["T"].flatten()
        P_arr = args[1]["P"].flatten()
        C_arr = args[1]["C"].flatten()
        E_arr = args[1]["E"].flatten()
        R_arr = args[1]["R"].flatten()

        @np.vectorize
        def _sdf(x, y, B, T, P, C, E, R):
            cm = C * np.sin(np.pi * x ** E) + R * np.sin(2 * np.pi * x)

            if x < 0:
                distance = np.sqrt(x ** 2 + y ** 2)
            elif x > 1:
                distance = np.sqrt((x - 1) ** 2 + y ** 2)
            else:
                # Get theta
                if 0 < x < 0.5:
                    neg_radical = -(-(x-0.5)/0.5) ** (1 / (B - 1))
                    if y > cm:
                        theta = np.arccos(neg_radical)
                    elif y < cm:
                        theta = 2 * np.pi - np.arccos(neg_radical)
                elif 0.5 < x < 1:
                    pos_radical = ((x-0.5)/0.5) ** (1 / (B - 1))
                    if y > cm:
                        theta = np.arccos(pos_radical)
                    elif y < cm:
                        theta = 2 * np.pi - np.arccos(pos_radical)

                boundary_x = 0.5 + 0.5 * (
                    np.abs(
                        np.cos(theta)
                    ) ** B
                    / np.cos(theta))


                boundary_y = T / 2
                boundary_y *= np.abs(np.sin(theta)) ** B / np.sin(theta)
                boundary_y *= 1 - boundary_x ** P
                boundary_y += C * np.sin(np.pi * boundary_x ** E)
                boundary_y += R * np.sin(2 * np.pi * boundary_x)

                distance = np.sign(boundary_y - cm) * (y - boundary_y)

            return distance

        return {"sdf": _sdf(x_arr, y_arr, B_arr, T_arr, P_arr, C_arr, E_arr, R_arr).reshape(args[0]["x"].shape)}