Initial condition is not correct

[hedlofr]

I modified the example for solving the heat equation with Neumann BC's found here such that the initial condition everywhere in the domain is $u(0, x) = 300$, $u_x(t, 0) = -5\times 10^4$, and $u_x(t, 1) = 0.0$. When the solutions are plotted at the end, I get this:

Why is the value of $u$ at $t = 0$ and $x = 0$ not equal to 300 (c.f., blue curve)? It appears to be over $500$ which is inconsistent with the initial condition provided.

I'm using Julia 1.8.2 on RHEL8 with MethodOfLines v0.9.6, OrdinaryDiffEq v6.44.1, ModelingToolkit v8.69.1, DomainSets v0.6.7

[hedlofr]

FWIW, I dumped the solution after a couple of timesteps and it appears as though this incorrect value for the initial condition is being used for the next timestep; it's not just an artifact. Also, the value at the first node appears to be related to the value used for the Neumann BC at $x=0$; as the magnitude of the Neumann BC goes up, so does the value of $u(0,0)$

[xtalax]

Can you post the whole system definition?

[hedlofr]

@xtalax, thanks for getting to this so quickly! Here are the contents of the script

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
# Method of Manufactured Solutions: exact solution
#u_exact = (x,t) -> exp.(-t) * cos.(x)

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq  = Dt(u(t, x)) ~ Dxx(u(t, x))
bcs = [u(0, x) ~ 300.0,
        Dx(u(t, 0)) ~ -5e4,
        Dx(u(t, 1)) ~ 0.0e4]

# Space and time domains
domains = [t ∈ Interval(0.0, 2.0e-3),
        x ∈ Interval(0.0, 1.0)]

# PDE system
@named pdesys = PDESystem(eq, bcs, domains,[t, x],[u(t, x)])

# Method of lines discretization
# Need a small dx here for accuracy
dx = 0.01
order = 2
discretization = MOLFiniteDifference([x => dx],t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization)

# Solve ODE problem
using OrdinaryDiffEq
#sol = solve(prob, Tsit5(), saveat=0.2, dt=1e-4)
sol = solve(prob, Euler(), dt=1e-5, saveat=0.2)

# Plot results and compare with exact solution
discrete_x = sol[x]
discrete_t = sol[t]

solu = sol[u(t, x)]

using Plots
plt = plot()

for i in 1:length(discrete_t)
    plot!(discrete_x, solu[i, :], label="Numerical, t=$(discrete_t[i])")
#    scatter!(discrete_x, u_exact(discrete_x, discrete_t[i]), label="Exact, t=$(discrete_t[i])")
end
display(plt)

[xtalax]

I'm wondering if the Neumann condition on the left boundary is forcing this somehow

[hedlofr]

@xtalax, yes I believe that's the case. If you print out solu[0,:], it looks like it's only the single node on the boundary that is affected; everything else is at 300.0 like it should be. Seems like somewhere in the backend the code is updating this boundary value prior to actually doing the solve. Happens with both Euler() and Tsit5(), and it occurs if you flip the domain and apply the Neumann BC to the other side also.

[xtalax]

The boundary value has to be locked high for both the initial condition and the boundary condition to be true. This isn't a bug, its just the way the system has been specified

[hedlofr]

@xtalax Sorry to say it, but that makes no physical sense. Say you're solving the heat equation in a hypothetical solid that has $\rho = 1~\mathrm{kg/m^3}$, $C_p = 1.0~\mathrm{J/(kg\cdot K)}$, and $\kappa = 1~\mathrm{W/(m\cdot K)}$ and your solid starts at $T = 300~\mathrm{K}$. Further, imagine that there is no thermal transfer into or out of the solid for the first, say, 1 nanosecond. Then you apply a constant energy flux $\Gamma^\epsilon = 5 \times 10^4~\mathrm{W/m^2}`$ to the surface (by say a laser, plasma, space heater, etc.). The code above models this physical system, except that there is no delay in imparting the energy in the code.

After the laser starts to hit the solid, it doesn't magically jump to $T > 500~\mathrm{K}$, right? It will take incremental steps to the value of $T = 500~\mathrm{K}$. This is definitely a bug.

From a numerical perspective, I don't think satisfying your initial condition and a Neumann BC on either/both sides of the domain are related.

[hedlofr]

@xtalax, FWIW, there's an analytical solution for this system if you model the solid as being semi-infinite:

T(t,x) = T_0 + \frac{2 \Gamma^\epsilon \sqrt{\frac{\alpha t}{\pi}}}{\kappa} \exp\bigg(-\frac{x^2}{4\alpha t}\bigg) - \frac{\Gamma^\epsilon x}{\kappa} \mathrm{erfc}\bigg(\frac{x}{2\sqrt{\alpha t}}\bigg),

where $\alpha = \kappa / (\rho C_p)$ (Frank P. Incropera and David P. DeWitt. Fundamentals of Heat and Mass Transfer. John Wiley & Sons, New York, 2001.)

If you make the domain big enough that the heat doesn't diffuse to the opposite side of the domain from the boundary, you can use this equation as a verification problem/regression test for your code so that future updates don't break it again...

[xtalax]

This is how finite differencing works I'm afraid, the boundary points are effectively ghost points and may take unphysical values because it is required for the discretization you requested to be true. If this bothers you, ignore the ghost cells, the interior will be accurate.

Thanks for showing me the analytic solution though, I will add it to the tests.

[hedlofr]

@xtalax, you can take any code out there (COMSOL, ANSYS, Sierra) except yours or just solve it via finite differences and you won't get the same behavior. The end points aren't ghost cells, they're part of the domain. The ghost cells are what you use to define the boundary values but retain the order of the discretization. For the left side (x = 0, index j = 1), the boundary condition is

\frac{\partial u}{\partial x} = \frac{\Gamma^\epsilon}{\kappa}

Numerically, at the surface of the boundary at index $j=1$ this means

\frac{u_2 - u_0}{2 \Delta x} = \frac{\Gamma^\epsilon}{\kappa}

From there, you solve that for $u_0$ and plug it back into your discretized heat equation on the boundary. If you're using explicit-Euler you get this

u^{n+1}_1 = u^n_1 + \frac{\kappa \Delta t}{\rho C_p \Delta x^2} (u_2^n - u_1^n) - 2 \Delta x \Gamma^\epsilon

[xtalax]

Yes that's the weak form discretization, it becomes much more difficult to apply in general which is what we are trying to do here, but I may try to detect and move to it where possible