- Thank you received: 0
Volume integration
4 years 9 months ago #2373
by andrew22
Volume integration was created by andrew22
I've written a script (based on the examples) to get the electric field between conductors: a cylindrical rod and a surrounding larger cylindrical shield. I would like to integrate the square of the field over the volume mesh to give the stored energy, from which I can calculate capacitance. I'm interested to see how accurate the FEM value will be.
I integrate with:
func= E*E
integ = Integrate(func, mesh, VOL)
print (integ)
but this gives NaN.
I'm probably making an elementary mistake and would appreciate any help.
Andrew.
I integrate with:
func= E*E
integ = Integrate(func, mesh, VOL)
print (integ)
but this gives NaN.
I'm probably making an elementary mistake and would appreciate any help.
Andrew.
Attachments:
- mneunteufel
-
- Offline
- Premium Member
-
Less
More
- Thank you received: 59
4 years 9 months ago #2375
by mneunteufel
Replied by mneunteufel on topic Volume integration
Hi Andrew,
your definition of epsl
yields a singular matrix a as epsl is only defined on the first material domain. The sparsecholesky solver does not say if the system is singular, but integrating over the resulting u gives already NaN ( the UMFPACK solver complains that it is singular).
Defining
solves the problem. If you want to have different values of epsl on the domains you have to write
Best
Michael
your definition of epsl
Code:
epsl = CoefficientFunction( [1.0006] )
yields a singular matrix a as epsl is only defined on the first material domain. The sparsecholesky solver does not say if the system is singular, but integrating over the resulting u gives already NaN ( the UMFPACK solver complains that it is singular).
Defining
Code:
epsl = 1.0006
solves the problem. If you want to have different values of epsl on the domains you have to write
Code:
epsl = CoefficientFunction( [value first domain, value second domain, ...] )
Best
Michael
4 years 9 months ago #2378
by andrew22
Replied by andrew22 on topic Volume integration + adaptive meshing
Hi Michael,
Thank you for the help. I'm very impressed by the speed and accuracy (the calculated capacitance compares well to results I've calculated with three other programs).
If I can trouble you further, I would like to make the adaptive meshing work. It is not actually needed for this problem, but it would help me to gain understanding. I've followed the recipe set out in the tutorial - but it does not quite work. See attachment (same .vol file needed).
Regards,
Andrew.
Thank you for the help. I'm very impressed by the speed and accuracy (the calculated capacitance compares well to results I've calculated with three other programs).
If I can trouble you further, I would like to make the adaptive meshing work. It is not actually needed for this problem, but it would help me to gain understanding. I've followed the recipe set out in the tutorial - but it does not quite work. See attachment (same .vol file needed).
Regards,
Andrew.
Attachments:
- mneunteufel
-
- Offline
- Premium Member
-
Less
More
- Thank you received: 59
4 years 9 months ago #2379
by mneunteufel
Replied by mneunteufel on topic Volume integration + adaptive meshing
Hi Andrew,
you need the autoupdate=True flag also for die GridFunction u1
Then, you also need to set the boundary data in every solve stage
Best
Michael
you need the autoupdate=True flag also for die GridFunction u1
Code:
u1 = GridFunction(fes, name="extension", autoupdate=True)
Then, you also need to set the boundary data in every solve stage
Code:
def SolveBVP():
a.Assemble()
f = u.vec.CreateVector()
u1.Set(extendone, BND)
f.data = a.mat * u1.vec
u.vec.data -= a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky") * f
Redraw (blocking=True)
Best
Michael
Attachments:
4 years 9 months ago #2385
by andrew22
Replied by andrew22 on topic Volume integration
Hi Michael,
I've edited the script as you suggested, but I'm still having problems with the refinement. I would appreciate your help again. The expected capacitance is about 2.326 pF for the dimensions given.
I get:
maxerr = 0.026982948684526564
maxerr = 0.9045386846912055
maxerr = 0.7627369747838466
maxerr = 0.9376416794932843
Capacitance = 24.875 pF
The revised script imports a module (meshcoaxcap.py) which generates the .vol file (line 14). If I comment out the line (#MakeCoaxMesh(..) ) and use last .vol file I get different output.
maxerr = 0.0269829486845266
maxerr = nan
maxerr = nan
etc.
That's rather strange, because line 14 is at the top of the script and is called only once.
I note that the outer boundary shown in Netgen changes from red (1V) to blue (0V) after the first refinement. This suggests that there might be a problem with refinement on the conductor surfaces. Or perhaps some of the metadata associated with the mesh (surfnr etc.) is not correct?
I have 6.2.2001-0~ubuntu18.04.1 installed.
Regards,
Andrew.
I've edited the script as you suggested, but I'm still having problems with the refinement. I would appreciate your help again. The expected capacitance is about 2.326 pF for the dimensions given.
I get:
maxerr = 0.026982948684526564
maxerr = 0.9045386846912055
maxerr = 0.7627369747838466
maxerr = 0.9376416794932843
Capacitance = 24.875 pF
The revised script imports a module (meshcoaxcap.py) which generates the .vol file (line 14). If I comment out the line (#MakeCoaxMesh(..) ) and use last .vol file I get different output.
maxerr = 0.0269829486845266
maxerr = nan
maxerr = nan
etc.
That's rather strange, because line 14 is at the top of the script and is called only once.
I note that the outer boundary shown in Netgen changes from red (1V) to blue (0V) after the first refinement. This suggests that there might be a problem with refinement on the conductor surfaces. Or perhaps some of the metadata associated with the mesh (surfnr etc.) is not correct?
I have 6.2.2001-0~ubuntu18.04.1 installed.
Regards,
Andrew.
Attachments:
- mneunteufel
-
- Offline
- Premium Member
-
Less
More
- Thank you received: 59
4 years 9 months ago #2390
by mneunteufel
Replied by mneunteufel on topic Volume integration
Hi andrew,
the problem with the different boundary condition after the first refinement was that you first set u to u1 and then set the right boundary condition for u1. You have to switch these lines to
I also observed different behavior when commented out the mesh generation function. I am not an expert of merging meshing, so I don't really have a clue where this behavior comes from.
However, for this geometry, you can simplifying the mesh generation by only using CSG geometries
Is there a reason why you mesh also the inner part? The solution seems to be always zero there.
With this, you can now use strings to identify your boundaries and you can also Curve your mesh.
With the attached code I get the following output:
ndof = 10423
maxerr = 0.13384372761117735 , toterr = 2.267302450023096
ndof = 15180
maxerr = 0.027621924443630222 , toterr = 1.809503014134304
ndof = 26375
maxerr = 0.007188620406700055 , toterr = 1.3764282981704608
ndof = 52185
maxerr = 0.007488492893794685 , toterr = 0.8460445676985052
ndof = 53249
maxerr = 0.001808892090986848 , toterr = 0.8097602634799458
ndof = 83743
maxerr = 0.0005241508456320649 , toterr = 0.64365894648496
ndof = 141903
maxerr = 0.00018470303727374993 , toterr = 0.5030248548701931
ndof = 245948
maxerr = 9.058668270776927e-05 , toterr = 0.37027924480707936
ndof = 326343
maxerr = 4.618238944138729e-05 , toterr = 0.30930491094782586
ndof = 418213
maxerr = 2.0765179180659515e-05 , toterr = 0.2665229929125358
ndof = 594250
maxerr = 1.1146830833658861e-05 , toterr = 0.22296841251806954
Capacitance = 2.3274 pF
which is quite near to your expected capacitance of 2.326 pF
Best,
Michael
the problem with the different boundary condition after the first refinement was that you first set u to u1 and then set the right boundary condition for u1. You have to switch these lines to
Code:
def SolveBVP():
a.Assemble()
#first set u1
u1.Set(voltage, definedon=mesh.Boundaries("outer"))# Boundary sources
u.vec.data = u1.vec
I also observed different behavior when commented out the mesh generation function. I am not an expert of merging meshing, so I don't really have a clue where this behavior comes from.
However, for this geometry, you can simplifying the mesh generation by only using CSG geometries
Code:
cube = OrthoBrick( Pnt(-outer_rad,-outer_rad,-leng/2), Pnt(outer_rad,outer_rad,2.0*quarterleng) )
cylb = Cylinder ( Pnt(0.0, 0.0, -leng), Pnt(0.0, 0.0, 0.0), outer_rad)
cube2 = OrthoBrick( Pnt(-inner_rad,-inner_rad,-leng/4), Pnt(inner_rad,inner_rad,quarterleng) )
cyli = Cylinder ( Pnt(0.0, 0.0, -leng), Pnt(0.0, 0.0, 0.0), inner_rad)
inner = (cyli*cube2).bc('inner')
outer = (cylb*cube-inner).bc('outer')
geo = CSGeometry()
geo.Add(outer)
#geo.Add(inner)
mesh = Mesh(geo.GenerateMesh (maxh=h))
Is there a reason why you mesh also the inner part? The solution seems to be always zero there.
With this, you can now use strings to identify your boundaries and you can also Curve your mesh.
With the attached code I get the following output:
ndof = 10423
maxerr = 0.13384372761117735 , toterr = 2.267302450023096
ndof = 15180
maxerr = 0.027621924443630222 , toterr = 1.809503014134304
ndof = 26375
maxerr = 0.007188620406700055 , toterr = 1.3764282981704608
ndof = 52185
maxerr = 0.007488492893794685 , toterr = 0.8460445676985052
ndof = 53249
maxerr = 0.001808892090986848 , toterr = 0.8097602634799458
ndof = 83743
maxerr = 0.0005241508456320649 , toterr = 0.64365894648496
ndof = 141903
maxerr = 0.00018470303727374993 , toterr = 0.5030248548701931
ndof = 245948
maxerr = 9.058668270776927e-05 , toterr = 0.37027924480707936
ndof = 326343
maxerr = 4.618238944138729e-05 , toterr = 0.30930491094782586
ndof = 418213
maxerr = 2.0765179180659515e-05 , toterr = 0.2665229929125358
ndof = 594250
maxerr = 1.1146830833658861e-05 , toterr = 0.22296841251806954
Capacitance = 2.3274 pF
which is quite near to your expected capacitance of 2.326 pF
Best,
Michael
Attachments:
The following user(s) said Thank You: andrew22
Time to create page: 0.112 seconds