1D diffusion equation: wildly different results when dividing by diffusion constant

[OthmanElHammouchi]

Hi, I’m trying to diagnose an issue with my own code using Fenics, but I’m running into the same problem with it. I need to solve the very elementary 1D diffusion equation

-k \frac{d^2T}{dx^2} = q

on the interval [5, 108] where k and q are piecewise constant functions whose values I have in a vector. The equation is subject to a Dirichlet boundary condition at x = 5 and the natural boundary condition on the left endpoint (i.e. k \frac{dT}{dx} = 0, which is the default if I understand correctly). Depending on whether I divide both sides by k, I get wildly different answers:

from fenics import *

mesh = IntervalMesh(n_elems, vec[0], vec[-1])
V = FunctionSpace(mesh, "P", 3)


def boundary_L(x, on_boundary):
    r = on_boundary and np.abs(x - vec[0]) < 1e-4
    return r


left_val = Constant(Ta)
bc = DirichletBC(V, left_val, boundary_L)

T = TrialFunction(V)
v = TestFunction(V)

class q_fen(UserExpression):
    def eval(self, value, x):
        value[0] = q(x[0])

class k_fen(UserExpression):
    def eval(self, value, x):
        value[0] = k(x[0])

class f(UserExpression):
    def eval(self, value, x):
      value[0] = q(x[0]) / k(x[0])

q_fen = q_fen(element = V.ufl_element())
k_fen = k_fen(element = V.ufl_element())
f = f(element = V.ufl_element())

a = k_fen*dot(grad(T), grad(v)) * dx
b = dot(grad(T), grad(v)) * dx

L1 = q_fen * v * dx
L2 = f * v * dx

T1 = Function(V)
T2 = Function(V)
solve(a == L1, T1, bc)
solve(b == L2, T2, bc)

I’ve been breaking my head trying to understand why. Any help would be welcome.

[Ekrem_Ekici]

Please share the full code as we do not know “vec”

[OthmanElHammouchi]

My apologies, vec is just the mesh

array([  5.,  12.,  18.,  24.,  30.,  37.,  39.,  44.,  50.,  53.,  57.,
        63.,  70.,  76.,  83.,  90.,  96., 102., 108.])

determining the values of k and q, i.e. they take different values on [5, 12], [12, 18], etc.

[dparsons]

In one system you’re setting L using q, while in the other you’re setting L using q/k. It’s not clear to me why you’d expect the two solutions to be the same.

[OthmanElHammouchi]

Because \frac{d^2T}{dx^2} = -\frac{q}{k} should be equivalent to -k\frac{d^2T}{d^2x} = q.

[OthmanElHammouchi]

Oh I see what you mean, my mistake, I used the a twice. The post is now corrected, but the same problem remains: they give different results.

[dparsons]

You haven’t given the value of k. Is it constant?

[OthmanElHammouchi]

No, it’s piecewise constant, so between every two points of vec it has a different value (and the same goes for q, see graphs below). Both are implemented as function which give the correct value at any x.

(orange is q)

[dparsons]

I suggest checking the derivation of the weak formulations from first principles. Relocating k in the strong formulation (the differential equation) will affect the application of integration by parts used to reduce the Laplacian operator down to the gradient operator. Keep in mind the weak formulation is an integral equation, so you can’t simply move k directly from one integrand to the other if it’s not constant.

[dokken]

First of all, you should make sure that your code can reproduce the result.

If you use a constant k, as @dparsons mentions, you can differentiate by it on both sides of the equation.
i.e.
an MWE would be:

    from fenics import *
    import numpy as np
    n_elems = 10
    vec = np.array([5.,  12.,  18.,  24.,  30.,  37.,  39.,  44.,  50.,  53.,  57.,
                    63.,  70.,  76.,  83.,  90.,  96., 102., 108.])

    mesh = IntervalMesh(n_elems, vec[0], vec[-1])
    V = FunctionSpace(mesh, "P", 3)


    def boundary_L(x, on_boundary):
        r = on_boundary and np.abs(x - vec[0]) < 1e-4
        return r


    Ta = 0
    left_val = Constant(Ta)
    bc = DirichletBC(V, left_val, boundary_L)

    T = TrialFunction(V)
    v = TestFunction(V)

    x = SpatialCoordinate(mesh)
    def k(x): return 0.1


    def q(x):
        return 0.3


    class q_fen(UserExpression):
        def eval(self, value, x):
            value[0] = q(x[0])


    class k_fen(UserExpression):
        def eval(self, value, x):
            value[0] = k(x[0])


    class f(UserExpression):
        def eval(self, value, x):
            value[0] = q(x[0]) / k(x[0])


    q_fen = q_fen(element=V.ufl_element())
    k_fen = k_fen(element=V.ufl_element())
    f = f(element=V.ufl_element())

    a = k_fen*dot(grad(T), grad(v)) * dx
    b = dot(grad(T), grad(v)) * dx

    L1 = q_fen * v * dx
    L2 = f * v * dx

    T1 = Function(V)
    T2 = Function(V)
    solve(a == L1, T1, bc)
    solve(b == L2, T2, bc)
    print(T1.vector().get_local())
    print(T2.vector().get_local())
    print(T1.vector().get_local()-T2.vector().get_local())

Usually, if k is not constant, the diffusion equation you solve is -\frac{\partial d}{\partial x} \left(k\frac{\partial }{\partial x} T \right) = q
which after integration by parts yields
\int k\frac{\partial}{\partial x} T \frac{\partial}{\partial x} v ~\mathrm{d} x = \int q v ~\mathrm{d} x
which means that you cannot divide by k.

[OthmanElHammouchi]

Thank you for the feedback, I am new here so I will try to keep it in mind. My reasoning was that it would be ok to divide by k if it is piecewise constant, because you’re just dividing up the domain and dividing on each piece:

-k_i\frac{dT}{dx} \bigg \rvert_{[x_{i - 1}, x_i]} = q_i \Rightarrow \frac{dT}{dx} \bigg \rvert_{[x_{i - 1}, x_i]} = -q_i / k_i \,.

Could you perhaps explain why this thinking is incorrect?

[dparsons]

That relationship is correct as far as it goes. The problem is that it applies only at a specific point x, while the weak formulation is an integral equation, as dokken showed, defined over the whole space,
\int k\frac{\partial}{\partial x} T \frac{\partial}{\partial x} v ~\mathrm{d} x = -\int q v ~\mathrm{d} x
(adding a minus sign here, I haven’t formally verified the sign).

So your relationship transforms that into
-\int q v ~\mathrm{d} x = \int k\frac{\partial}{\partial x} T \frac{\partial}{\partial x} v ~\mathrm{d} x = \int k \left( -\frac{q}{k}\right) v ~\mathrm{d} x = -\int q v ~\mathrm{d} x
which is the tautology

0 = 0

i.e. it doesn’t directly help construct the weak formulation

[Edit] Wait, no, the transformation would be

-\int q v ~\mathrm{d} x = \int k\frac{\partial}{\partial x} T \frac{\partial}{\partial x} v ~\mathrm{d} x = \int k \left( -\frac{q}{k}\right) \frac{\partial}{\partial x} v ~\mathrm{d} x = -\int q v ~\mathrm{d} x
which in any case is not the correction formulation.

[dparsons]

Your premise does not seem to be correct, though. The differential equation (strong form) is
- \frac{d}{dx} \left( k(x) \frac{dT}{dx} \right) = q
That does not mean that
-k\frac{dT}{dx} = q

[OthmanElHammouchi]

I think you are right. Thank you very much for the help!