Specifying inconsisten initial and boundary conditions

  • mcmcclur
  • Offline
  • New Member
  • New Member
  • Professor of Mathematics at UNC - Asheville
More
3 years 9 months ago - 3 years 9 months ago #3505 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
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.
More
3 years 9 months ago #3507 by hvwahl
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:
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
  • mcmcclur
  • Offline
  • New Member
  • New Member
  • Professor of Mathematics at UNC - Asheville
More
3 years 9 months ago #3508 by mcmcclur
This works perfectly, Henry - Thank you very much!
Time to create page: 0.102 seconds