How to write piecewise function in the source of a Poisson equation?

my equation is :

-(u_xx + u_yy + u_zz) = f(r)

where r = sqrt(x**2 + y**2 + z**2)

and

f(r) = 1 for r < 0.5
f(r) = 0 for r > 0.5

I have tried:

class Poisson(PDE):
    name = "Poisson"

    def __init__(self):
        # coordinates
        x, y, z = Symbol("x"), Symbol("y"), Symbol("z")
        rr = Symbol("w")
        # make input variables
        input_variables = {"x": x, "y": y, "z": z}

        # make u function
        u = Function("u")(*input_variables)

        # source term is $f = \left\{ \begin{array}{ll}
        #   1 & \text{if } r < 0.5\\
        #   0 & \text{if } r > 0.5
        # \end{array} \right. \text{. }$
        r = (x ** 2 + y ** 2 + z ** 2) ** 0.5
        f = Piecewise((Number(1), rr < Number(0.5)), (Number(0), rr > Number(0.5)))

        # set equations
        self.equations = {}
        self.equations["poisson_u"] = - (u.diff(x, 2) + u.diff(y, 2) + u.diff(z, 2)) - f.subs(rr,r)

but failed:

could not unroll graph!
This is probably because you are asking to compute a value that is not an output of any node
####################################
invar: [x, y, z, sdf, area]
requested var: [poisson_u]
computable var: [x, y, z, sdf, area, u]
####################################
Nodes in graph: 
node: Sympy Node: poisson_u
evaluate: SympyToTorch
inputs: [Piecewise]
derivatives: [u__x__x, u__y__y, u__z__z]
outputs: [poisson_u]
optimize: False
node: Arch Node: poisson_network
evaluate: FullyConnectedArch
inputs: [x, y, z]
derivatives: []
outputs: [u]
optimize: True
####################################

I have searched the keyword piecewise on the examples but found nothing helpful.

I am a very beginner so I dont know if there is some method to write such a source term.

Thank you very much.

Now I have tried to use Heaviside to solve this problem.
But I still dont know how to write the code if this piecewise function is not easy to be described by Heaviside function.

Hi @Zhao-ZC

Interesting that the piecewise function does not work. But you should be able to convert it into a Heaviside representation (e.g. this example here).

Another option is to have two constraints, one for each part of the domain in your piecewise function. I.e have one constraint with a geometry that only contains the portion with the source term 1 and another constraint with a geometry that only contains the portion with source term 2. Not the cleanest, but would allow you to address the PDE source term in a isolated forms.

Multi-constraints is a good solution!
Thank you very much!

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.