Specifying inconsisten initial and boundary conditions

Main Forum Old-Forum - do not post here
1

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
2

Hi,

Set zeros-out the GridFunction where nothing is set, i.e., when only setting on a boundary. The way to get around this ist to use two GridFunctions, Set one to the boundary condition and one to the initial condition. You can then loop over the interior/free degrees of freedoms and copy the initial values into the GridFunction with the correct boundary condition:

gfu, gfu_int = GridFunction(fes), GridFunction(fes)
gfu.Set(values_list, definedon=mesh.Boundaries(dirichlet_boundaries))
gfu_int.Set(1-x**2)

free_dofs = fes.FreeDofs()

for dof in range(fes.ndof):
    if free_dofs[dof]:
        gfu.vec.data[dof] = gfu_int.vec[dof]

Best wishes,
Henry

3

This works perfectly, Henry - Thank you very much!