- Thank you received: 0
boundary integral
5 years 7 months ago #1517
by noname
boundary integral was created by noname
Hello everyone,
I have a very simple question. I am doing flux calculations and want to calculate line integral on the "outlet" boundary (the right boundary) of a square mesh. I know my solution i.e u_h is correct because it converges but it seems like I am doing something wrong with flux calculation since it is not giving me the results I expect.
The analytical flux I expect is grad(u_analy)*n and I choose u_analytical to be (x^2+y^2,0) in two dimensions. Since I am in right boundary of a square mesh my n is (1,0). So I expect flux analytical to be (2x , 0). The square has dimensions of [-1,1]^2 so the analytical result is 0. But mu numerical result is far from that.. I am doing the boundary integration below...
Am I doing it wrong ?
I have a very simple question. I am doing flux calculations and want to calculate line integral on the "outlet" boundary (the right boundary) of a square mesh. I know my solution i.e u_h is correct because it converges but it seems like I am doing something wrong with flux calculation since it is not giving me the results I expect.
The analytical flux I expect is grad(u_analy)*n and I choose u_analytical to be (x^2+y^2,0) in two dimensions. Since I am in right boundary of a square mesh my n is (1,0). So I expect flux analytical to be (2x , 0). The square has dimensions of [-1,1]^2 so the analytical result is 0. But mu numerical result is far from that.. I am doing the boundary integration below...
Code:
VH1= H1(mesh, order=2, dirichlet="wall|inlet")
gradu00 = GridFunction(VH1)
gradu01 = GridFunction(VH1)
gradu10 = GridFunction(VH1)
gradu11 = GridFunction(VH1)
gradu00.Set(grad(u_h.components[0])[0])
gradu01.Set(grad(u_h.components[0])[1])
gradu10.Set(grad(u_h.components[1])[0])
gradu11.Set(grad(u_h.components[1])[1])
et = GridFunction(VH1)
et.Set(CoefficientFunction(1.0), definedon=mesh.Boundaries("outlet"))
flux_x = Integrate(et*( gradu00*n[0] + gradu01*n[1] ), mesh, BND)
flux_y = Integrate(et*( gradu10*n[0] + gradu11*n[1] ), mesh, BND)
print ("flux_x:", flux_x)
print ("flux_y:", flux_y)
Am I doing it wrong ?
- mneunteufel
-
- Offline
- Premium Member
-
Less
More
- Thank you received: 59
5 years 7 months ago #1518
by mneunteufel
Replied by mneunteufel on topic boundary integral
Hello noname,
on the right boundary I expect x to be constant 1 and therefore the analytical flux to be 4, not 0.
What numerical value do you obtain?
Best,
Michael
on the right boundary I expect x to be constant 1 and therefore the analytical flux to be 4, not 0.
What numerical value do you obtain?
Best,
Michael
5 years 7 months ago #1519
by noname
Replied by noname on topic boundary integral
Hi Michael,
I get almost 4 numerically... but I really do not follow why you expect 4 analytically ...
I get almost 4 numerically... but I really do not follow why you expect 4 analytically ...
5 years 7 months ago - 5 years 7 months ago #1520
by noname
Replied by noname on topic boundary integral
With the same logic, on the circular boundary defined by level set function x^2+y^2-0.5^2=0 (using the same analytical solution) i calculate n to be <2x/sqrt(4x^2+4y^2) , 2y/sqrt(4x^2+4y^2)> and I calculate grad(u_analy)=<2x , 2y>
that makes analytical flux around the circle which is grad(u_analy).n to be
<4x^2/sqrt(4x^2+4y^2) + 4y^2/sqrt(4x^2+4y^2)> = 1. Am I also wrong here ?
that makes analytical flux around the circle which is grad(u_analy).n to be
<4x^2/sqrt(4x^2+4y^2) + 4y^2/sqrt(4x^2+4y^2)> = 1. Am I also wrong here ?
Last edit: 5 years 7 months ago by noname.
- mneunteufel
-
- Offline
- Premium Member
-
Less
More
- Thank you received: 59
5 years 7 months ago #1522
by mneunteufel
Replied by mneunteufel on topic boundary integral
1) Let the domain [-1,1]^2 and u = (x^2+y^2,0). Then grad(u)= (2x, 2y; 0, 0). The right boundary is given by the set { (1,y) : -1<=y<=1} and the according normal vector is n=(1,0). Therefore, grad(u)*n = (2x, 0) on the right boundary. Further x=1 on the right boundary => grad(u)*n = (2, 0) on the right boundary. Integrating this flux over the right boundary (it has length 2) yields int(2, y=-1..1) = 4.
This is what I also observe numerically.
2) Let the domain given by the level-set function x^2+y^2-0.5^2=0. The normal vector is then given by n= (2x,2y) (simplifying your expression by using the level-set equation). Thus, grad(u)*n = (4x^2 + 4y^2, 0) = (1,0). Integrating over the boundary of the circle yields int(1) = 2*r*pi= pi.
The numerical result is also near pi.
Attached you find the according python-code for both examples.
Best,
Michael
This is what I also observe numerically.
2) Let the domain given by the level-set function x^2+y^2-0.5^2=0. The normal vector is then given by n= (2x,2y) (simplifying your expression by using the level-set equation). Thus, grad(u)*n = (4x^2 + 4y^2, 0) = (1,0). Integrating over the boundary of the circle yields int(1) = 2*r*pi= pi.
The numerical result is also near pi.
Attached you find the according python-code for both examples.
Best,
Michael
Attachments:
The following user(s) said Thank You: noname
Time to create page: 0.114 seconds