The forum is in read only mode.

#
Specifying inconsisten initial and boundary conditions

*Specifying inconsisten initial and boundary conditions*was created by

*mcmcclur*

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

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

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:

Then,

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

Code:

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

Code:

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:

Code:

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,

Code:

%%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

Last edit: 3 years 5 months ago by mcmcclur.

Replied by

*hvwahl*on topic*Specifying inconsisten initial and boundary conditions*
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:

Best wishes,

Henry

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:

Code:

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

Time to create page: 0.151 seconds