Specifying inconsisten initial and boundary conditions
3 years 9 months ago - 3 years 9 months ago #3505
by mcmcclur
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 9 months ago by mcmcclur.
3 years 9 months ago #3507
by hvwahl
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
The following user(s) said Thank You: mcmcclur
Time to create page: 0.102 seconds