Hello together,
the new timesteppers are great. But wouldn’t it be nice if the u.dt term evaluates to zero in case of classical use of bilinear form and solver.
Consider e.g. following problem:
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
shape = Rectangle(1,1).Face()
shape.edges.Max(X).name="right"
shape.edges.Min(X).name="left"
shape.edges.Max(Y).name="top"
shape.edges.Min(Y).name="bot"
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.05))
fes = H1(mesh, order=4, dirichlet="left|right")
u,v = fes.TnT()
eq = ( grad(u) * grad(v) + u.dt * v) * dx
gfu = GridFunction(fes)
def init_gfu(gfu):
t_left = 10
t_right = 20
t_ref = 15
gfu.Interpolate(15)
gfu.Interpolate(mesh.BoundaryCF({"left": t_left, "right": t_right}, default=t_ref),definedon=mesh.Boundaries("left|right"))
works as expected in case of timestepper usage:
init_gfu(gfu)
ts = timestepping.ImplicitEuler(equation=eq, dt=1e-1)
scene = Draw(gfu)
def callback(t, gfu):
scene.Redraw()
ts.Integrate(gfu, start_time=0, end_time=2, callback=callback)
But if we plug eq directly into a bfa and solve the same problem, the u.dt evaluates to u and we get a different result instead of the quasi-static solution:
init_gfu(gfu)
a = BilinearForm(fes, symmetric=True)
a += eq.Compile()
solvers.Newton(a=a, u=gfu, inverse="sparsecholesky", printing=False)
scene = Draw (gfu)
import matplotlib.pyplot as plt
import numpy as np
X = np.linspace(0, 1, num=100)
Y = np.ones_like(X) * 0
plt.plot(X, gfu(mesh(X, Y)))
plt.show()
