Problem with solving a Poisson equation
- Ikaruga
-
Topic Author
- New Member
-
Less
More
1 year 7 months ago #4723
by Ikaruga
Problem with solving a Poisson equation was created by Ikaruga
I'm trying to solve the Poisson equation given as follows (sorry, I don't how to use Latex here):
delta u = x2 + y2
with the boundary conditions:
u(x,0) = 0
u(0,y) = 0
ux(0,y) = 0
uy(x,0) = 0
Here is my code:
From the code above I get an incorrect solution. The exact solution for this equation is (x2y2)/2. What am I doing wrong?
delta u = x2 + y2
with the boundary conditions:
u(x,0) = 0
u(0,y) = 0
ux(0,y) = 0
uy(x,0) = 0
Here is my code:
Code:
from ngsolve import *
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
source = x**2 + y**2
fes = H1(mesh, order=2, dirichlet='bottom|left')
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes, symmetric=True)
a += -grad(u) * grad(v) * dx
a.Assemble()
f = LinearForm(fes)
f += source * v * dx
f.Assemble()
gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec
Draw(gfu, mesh, "numerical_solution")
1 year 7 months ago - 1 year 7 months ago #4724
by hvwahl
Replied by hvwahl on topic Problem with solving a Poisson equation
Hi Ikaruga,
one small Issue: You should consider -Delta u = f. The main issue is, that your solution does not have homogeneous Neumann values (-Grad(u) * n ) on the right and top boundaries, which you imply in your code. You therefore either need to consider inhomogeneous Dirichlet conditions on these boundaries, or implement the inhomogeneous Neumann BC in the forcing term f
Best wishes,
Henry
one small Issue: You should consider -Delta u = f. The main issue is, that your solution does not have homogeneous Neumann values (-Grad(u) * n ) on the right and top boundaries, which you imply in your code. You therefore either need to consider inhomogeneous Dirichlet conditions on these boundaries, or implement the inhomogeneous Neumann BC in the forcing term f
Best wishes,
Henry
Code:
from ngsolve import *
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
source = x**2 + y**2
u_ex = - x**2 * y**2 / 2
fes = H1(mesh, order=2, dirichlet='bottom|left')
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes, symmetric=True)
a += grad(u) * grad(v) * dx
a.Assemble()
f = LinearForm(fes)
f += source * v * dx
f += - y**2 * v * ds(definedon=mesh.Boundaries('right'))
f += - x**2 * v * ds(definedon=mesh.Boundaries('top'))
f.Assemble()
gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec
Draw(gfu, mesh, "numerical_solution")
Draw(u_ex, mesh, "exact")
Draw(gfu - u_ex, mesh, "error")
print('L2 Error: ', sqrt(Integrate(InnerProduct(gfu - u_ex, gfu - u_ex) * dx, mesh)))
Last edit: 1 year 7 months ago by hvwahl. Reason: error in code print statement
The following user(s) said Thank You: Ikaruga
Time to create page: 0.093 seconds