I’d like to solve heat problems with inhomogeneous Dirichlet conditions and a non-zero initial condition that need not be consistent with the boundary conditions.
As a simple example, we might have a square whose initial temperature distribution is given by u(x,y,0)=1-x^2. At time zero, the bottom is set to temp 0 and the top is set to temp 1 while the sides are insulated. How can I model this problem?
Following the Parabolic Model tutorial, I see that I can use syntax such as
gfu = GridFunction(fes)
gfu.Set(1-x**2)
to specify the initial condition and then iteratively redefine gfu via an implicit Euler method that respects any stated Dirichlet conditions. I’m also aware that I can specify specify inhomogendous Dirichlet conditions via code like
gfu.Set(values_list, definedon=mesh.Boundaries(dirichlet_boundaries))
where dirichlet_boundaries is a regular expression specifying the boundaries and values_lists contains the values.
I can’t seem to find a way to accomplish both of these however. Perhaps there’s a way to set gfu to values_list on the boundary but 1-x**2 in the interior? Or, perhaps there’s another way to accomplish this?
My code to this point (which is mostly based on the parabolic tutorial) is below:
from netgen import gui
from math import pi
from ngsolve import *
from netgen.geom2d import SplineGeometry
time = 0.0
dt = 0.001
tstep = 4
geo = SplineGeometry()
geo.AddRectangle( (-1, -1), (1, 1),
bcs = ("bottom", "right", "top", "left"))
mesh = Mesh( geo.GenerateMesh(maxh=0.1))
fes = H1(mesh, order=3, dirichlet="bottom|top")
gfu = GridFunction(fes)
## This simple block for a single initial condition
#gfu.Set(1-x**2)
## This block for inhomogenous boundary conditions
boundary_conditions = {'bottom':'0', 'top':'1'}
dirichlet_boundaries = '|'.join(boundary_conditions.keys())
for k in boundary_conditions.keys():
boundary_conditions[k] = eval(boundary_conditions[k])
values_list = [
boundary_conditions[bc] if bc in boundary_conditions.keys() else 0
for bc in mesh.GetBoundaries()
]
gfu.Set(values_list, definedon=mesh.Boundaries(dirichlet_boundaries))
res = gfu.vec.CreateVector()
Draw(gfu,mesh,"init")
u,v = fes.TnT()
a = BilinearForm(fes, symmetric=False)
a += 0.5*grad(u)*grad(v)*dx
a.Assemble()
m = BilinearForm(fes, symmetric=False)
m += u*v*dx
m.Assemble()
mstar = m.mat.CreateMatrix()
mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()
invmstar = mstar.Inverse(freedofs=fes.FreeDofs())
Then,
%%time
t_intermediate=0 # time counter within one block-run
while t_intermediate < tstep - 0.5 * dt:
res.data = - dt * a.mat * gfu.vec
gfu.vec.data += invmstar * res
t_intermediate += dt
print("\r",time+t_intermediate,end="")
Redraw(blocking=True)
print("")
time+=t_intermediate