Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Problem with solving a Poisson equation

More
1 year 3 weeks 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:
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")
From the code above I get an incorrect solution. The exact solution for this equation is (x2y2)/2. What am I doing wrong?
More
1 year 3 weeks ago - 1 year 3 weeks 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

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 3 weeks ago by hvwahl. Reason: error in code print statement
The following user(s) said Thank You: Ikaruga
Time to create page: 0.134 seconds

Login