BH curve
- lahoussaine
- Topic Author
- New Member
Less
More
1 year 8 months ago - 1 year 8 months ago #4629
by lahoussaine
Replied by lahoussaine on topic BH curve
Here is the code :
u, v = fes.TnT()
mu_r = {"air": 1, "magnet": 1.1, "fer1": 1000, "fer2": 1000}
# Definition des permeabilites
def nuf(u):
return IfPos(B-2.28,(1/(-0.0039954722*B**3+0.00595660*B**2+0.00646251*B + 0.001925625)), ((1902.211579*B+-4266.42343)/0.016694669))
B = Norm(grad(u))
nu_test = [nuf( if mat=="fer1" else 1.0/(mu0*mu_r[mat]) for mat in mesh.GetMaterials()]
nu = CoefficientFunction(nu_test)
a = BilinearForm(fes, symmetric=False)
a += nu *(grad(u)*grad(v))*dx
...
u, v = fes.TnT()
mu_r = {"air": 1, "magnet": 1.1, "fer1": 1000, "fer2": 1000}
# Definition des permeabilites
def nuf(u):
return IfPos(B-2.28,(1/(-0.0039954722*B**3+0.00595660*B**2+0.00646251*B + 0.001925625)), ((1902.211579*B+-4266.42343)/0.016694669))
B = Norm(grad(u))
nu_test = [nuf( if mat=="fer1" else 1.0/(mu0*mu_r[mat]) for mat in mesh.GetMaterials()]
nu = CoefficientFunction(nu_test)
a = BilinearForm(fes, symmetric=False)
a += nu *(grad(u)*grad(v))*dx
...
Last edit: 1 year 8 months ago by lahoussaine.
- mrambausek
- Offline
- New Member
Less
More
- Thank you received: 3
1 year 8 months ago #4630
by mrambausek
Replied by mrambausek on topic BH curve
You still use Norm(grad(u)) in the bilinear form instead of nu.
Is this an oversight?
Did you call a.Assemble or a.AssembleLinearization further doen behind the dots?
Is this an oversight?
Did you call a.Assemble or a.AssembleLinearization further doen behind the dots?
- lahoussaine
- Topic Author
- New Member
Less
More
1 year 8 months ago - 1 year 8 months ago #4631
by lahoussaine
Replied by lahoussaine on topic BH curve
Here is the reste of code
a += nu*(grad(u)*grad(v))*dx
# f= LinearForm(fes)
j = CoefficientFunction(10*1e3)
a+=- (j*v) * dx("magnet")
def SimpleNewtonSolve(gfu, a, tol=1e-8, maxits=30): res = gfu.vec.CreateVector()
du = gfu.vec.CreateVector()
fes = gfu.space
for it in range(maxits):
print("Iteration {} ".format(it), end="")
a.Apply(gfu.vec, res)
a.AssembleLinearization(gfu.vec)
du.data = a.mat.Inverse(fes.FreeDofs()) * res
gfu.vec.data -= du
# stopping criteria
stopcritval = sqrt(abs(InnerProduct(du, res)))
print("<A u", it, ", A u", it, ">_{-1}^0.5 = ", stopcritval)
if stopcritval < tol:
break
# Solve the linear system
gfa = GridFunction(fes)
SimpleNewtonSolve(gfa, a)
a += nu*(grad(u)*grad(v))*dx
# f= LinearForm(fes)
j = CoefficientFunction(10*1e3)
a+=- (j*v) * dx("magnet")
def SimpleNewtonSolve(gfu, a, tol=1e-8, maxits=30): res = gfu.vec.CreateVector()
du = gfu.vec.CreateVector()
fes = gfu.space
for it in range(maxits):
print("Iteration {} ".format(it), end="")
a.Apply(gfu.vec, res)
a.AssembleLinearization(gfu.vec)
du.data = a.mat.Inverse(fes.FreeDofs()) * res
gfu.vec.data -= du
# stopping criteria
stopcritval = sqrt(abs(InnerProduct(du, res)))
print("<A u", it, ", A u", it, ">_{-1}^0.5 = ", stopcritval)
if stopcritval < tol:
break
# Solve the linear system
gfa = GridFunction(fes)
SimpleNewtonSolve(gfa, a)
Last edit: 1 year 8 months ago by lahoussaine.
- mrambausek
- Offline
- New Member
Less
More
- Thank you received: 3
1 year 8 months ago #4632
by mrambausek
Replied by mrambausek on topic BH curve
Okay, there a couple of issues:
First of all, using grad(u) and grad(v) suggests that you are using the (auxiliary) scalar magnetic potential (H=-grad(u)). However, everywhere else you work like you'd use the vector potential (B=curl(A)), i.e. using nu=1/mu and j*v as loading term. So, to me, the whole formulation seems to be inconsistent. Could you recheck, please?
Second, when using a.AssembleLinearization(...), ngsolve will linearize the expressions it finds in the bilinear form in some way. However, the nu*grad(u)*grad(v) is a linearized form already.
When trying to linearize this term, you get NaN's because of a division by zero (or a very small number) resulting from the norm (at B=0).
When you work with CFs representing already "linearized" forms,
they must not depend on trial functions but grid functions only.
So nu must not depend on "u" but "gfu" in this case.
Please note that I did not try to fix the example, so there may be some things that I have overlooked (also in my own explanations). I hope they are still useful though.
Best,
Matthias
First of all, using grad(u) and grad(v) suggests that you are using the (auxiliary) scalar magnetic potential (H=-grad(u)). However, everywhere else you work like you'd use the vector potential (B=curl(A)), i.e. using nu=1/mu and j*v as loading term. So, to me, the whole formulation seems to be inconsistent. Could you recheck, please?
Second, when using a.AssembleLinearization(...), ngsolve will linearize the expressions it finds in the bilinear form in some way. However, the nu*grad(u)*grad(v) is a linearized form already.
When trying to linearize this term, you get NaN's because of a division by zero (or a very small number) resulting from the norm (at B=0).
When you work with CFs representing already "linearized" forms,
they must not depend on trial functions but grid functions only.
So nu must not depend on "u" but "gfu" in this case.
Please note that I did not try to fix the example, so there may be some things that I have overlooked (also in my own explanations). I hope they are still useful though.
Best,
Matthias
The following user(s) said Thank You: lahoussaine
- lahoussaine
- Topic Author
- New Member
Less
More
1 year 8 months ago - 1 year 8 months ago #4633
by lahoussaine
Replied by lahoussaine on topic BH curve
Thank you for your time, I will just try to clarify my problem one last time.
You will find below the code of the Linear solution.
Concerning the linear solution, the weak formulation has been simplified, so the equation might be not clear, sorry for that. However, this simulation has been validated both with a commercial solver and another FEA api.
My problem is not about the weak formulation but the way to carry on a NonLinear simulation, the NonLinearity coming from the nu expression which depend on B and thus on u.
The Nonlinear simulation is thus the same as the linear one but with nu(
The question is : How would you perform a Non linear algorithm iterating on nu ?
You'll find below the linear code that worked for me ...
stator = Circle(center=(0, 0), radius=10, bc="bc_rect_ext" )
rotor = Circle(center=(0, 0), radius=6)
circle_airgap = Circle(center=(0, 0), radius=7.5)
magnet1 = Circle(center=(0, 0), radius=7)
magnet_f = magnet1-rotor
airgap = circle_airgap-(magnet1)
stator_f = stator-circle_airgap
# change domain name and maxh
magnet_f.Mat("magnet").Maxh(0.5)
rotor.Mat("fer1").Maxh(0.5)
stator_f.Mat("fer2").Maxh(0.5)
airgap.Mat("air").Maxh(0.5)
###
geo.Add(magnet_f)
geo.Add(rotor)
geo.Add(stator_f)
geo.Add(airgap)
mesh = Mesh(geo.GenerateMesh(maxh=0.1))
Draw(mesh)
# FEM
fes = H1(mesh, order=3, dirichlet='bc_rect_ext')
u, v = fes.TnT()
mu_r = {"air": 1, "magnet": 1.1, "fer1": 1000, "fer2": 1000}
# Definition des permeabilites
mu0 = 4*pi*1e-7
nu_coef = [1/(mu0*mu_r[mat]) for mat in mesh.GetMaterials()]
nu = CoefficientFunction(nu_coef)
a = BilinearForm(fes)
a += nu*(grad(u)*grad(v))*dx
M=80000
f = LinearForm(fes)
j = CoefficientFunction(10*1e3)
Mx = CoefficientFunction(M*cos(atan2(y,x))*1e0)
My = CoefficientFunction(0)
f += (Mx*grad(v)[1]-My*grad(v)[0]) * dx("magnet")
with TaskManager():
a.Assemble()
f.Assemble()
# Solve the linear system
gfa = GridFunction(fes)
gfa.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
You will find below the code of the Linear solution.
Concerning the linear solution, the weak formulation has been simplified, so the equation might be not clear, sorry for that. However, this simulation has been validated both with a commercial solver and another FEA api.
My problem is not about the weak formulation but the way to carry on a NonLinear simulation, the NonLinearity coming from the nu expression which depend on B and thus on u.
The Nonlinear simulation is thus the same as the linear one but with nu(
The question is : How would you perform a Non linear algorithm iterating on nu ?
You'll find below the linear code that worked for me ...
stator = Circle(center=(0, 0), radius=10, bc="bc_rect_ext" )
rotor = Circle(center=(0, 0), radius=6)
circle_airgap = Circle(center=(0, 0), radius=7.5)
magnet1 = Circle(center=(0, 0), radius=7)
magnet_f = magnet1-rotor
airgap = circle_airgap-(magnet1)
stator_f = stator-circle_airgap
# change domain name and maxh
magnet_f.Mat("magnet").Maxh(0.5)
rotor.Mat("fer1").Maxh(0.5)
stator_f.Mat("fer2").Maxh(0.5)
airgap.Mat("air").Maxh(0.5)
###
geo.Add(magnet_f)
geo.Add(rotor)
geo.Add(stator_f)
geo.Add(airgap)
mesh = Mesh(geo.GenerateMesh(maxh=0.1))
Draw(mesh)
# FEM
fes = H1(mesh, order=3, dirichlet='bc_rect_ext')
u, v = fes.TnT()
mu_r = {"air": 1, "magnet": 1.1, "fer1": 1000, "fer2": 1000}
# Definition des permeabilites
mu0 = 4*pi*1e-7
nu_coef = [1/(mu0*mu_r[mat]) for mat in mesh.GetMaterials()]
nu = CoefficientFunction(nu_coef)
a = BilinearForm(fes)
a += nu*(grad(u)*grad(v))*dx
M=80000
f = LinearForm(fes)
j = CoefficientFunction(10*1e3)
Mx = CoefficientFunction(M*cos(atan2(y,x))*1e0)
My = CoefficientFunction(0)
f += (Mx*grad(v)[1]-My*grad(v)[0]) * dx("magnet")
with TaskManager():
a.Assemble()
f.Assemble()
# Solve the linear system
gfa = GridFunction(fes)
gfa.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Last edit: 1 year 8 months ago by lahoussaine.
- mrambausek
- Offline
- New Member
Less
More
- Thank you received: 3
1 year 8 months ago #4636
by mrambausek
Replied by mrambausek on topic BH curve
In my previous post, I might have had a confusion of the meaning of nu. Is the your nu "H(" or "dH/dB" evaluated for some B? I am used to work with "incremental" coefficients, ie "dH/dB". So *if you do not so*, then my second "issue" in my previous response is not valid.
Assuming that your nu if such that H = nu B, then a.AssembleLinearization(gfa.vec) will linearized it for you.
However, linearizing Norm( when B is zero, results in a division by zero. You can work around that as follows (see comments in the code below, in particular around "def nuf"):
fes = H1(mesh, order=3, dirichlet='bc_rect_ext')
u, v = fes.TnT()
mu_r = {"air": 1, "magnet": 1.1, "fer1": 1000, "fer2": 1000}
# Definition des permeabilites
mu0 = 4*pi*1e-7
def nuf(u):
# You need to avoid division by zero in the linearization! Either by perturbing the norm of by providing
# the linearized norm for the norm being close to zero.
# Perturbation
#B = Norm(grad(u)+CF((1e-6, 1e-6)))
# manual linearization
B = Norm(grad(u))
return IfPos(B-2.28,
(1/(-0.0039954722*B**3+0.00595660*B**2+0.00646251*B + 0.001925625)),
IfPos(B - 1e-8,
(1902.211579*B+-4266.42343)/0.016694669,
1902.211579/0.016694669) # linearized nu (dh/dB) for small B
)
nu_coef = [nuf(u) if mat=="fer1" else 1.0/(mu0*mu_r[mat]) for mat in mesh.GetMaterials()]
nu = CoefficientFunction(nu_coef)
a = BilinearForm(fes)
a += nu*(grad(u)*grad(v))*dx
M=80000
f = LinearForm(fes)
j = CoefficientFunction(10*1e3)
Mx = CoefficientFunction(M*cos(atan2(y,x))*1e0)
My = CoefficientFunction(0)
f += (Mx*grad(v)[1]-My*grad(v)[0]) * dx("magnet")
### "SIMULATING" the first Newton step only
# define gfa upfront
gfa = GridFunction(fes)
with TaskManager():
# check that AssembleLinearization() does not give NaNs
a.AssembleLinearization(gfa.vec)
f.Assemble()
# Solve the linear system
gfa.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Hope that helps!
Assuming that your nu if such that H = nu B, then a.AssembleLinearization(gfa.vec) will linearized it for you.
However, linearizing Norm( when B is zero, results in a division by zero. You can work around that as follows (see comments in the code below, in particular around "def nuf"):
fes = H1(mesh, order=3, dirichlet='bc_rect_ext')
u, v = fes.TnT()
mu_r = {"air": 1, "magnet": 1.1, "fer1": 1000, "fer2": 1000}
# Definition des permeabilites
mu0 = 4*pi*1e-7
def nuf(u):
# You need to avoid division by zero in the linearization! Either by perturbing the norm of by providing
# the linearized norm for the norm being close to zero.
# Perturbation
#B = Norm(grad(u)+CF((1e-6, 1e-6)))
# manual linearization
B = Norm(grad(u))
return IfPos(B-2.28,
(1/(-0.0039954722*B**3+0.00595660*B**2+0.00646251*B + 0.001925625)),
IfPos(B - 1e-8,
(1902.211579*B+-4266.42343)/0.016694669,
1902.211579/0.016694669) # linearized nu (dh/dB) for small B
)
nu_coef = [nuf(u) if mat=="fer1" else 1.0/(mu0*mu_r[mat]) for mat in mesh.GetMaterials()]
nu = CoefficientFunction(nu_coef)
a = BilinearForm(fes)
a += nu*(grad(u)*grad(v))*dx
M=80000
f = LinearForm(fes)
j = CoefficientFunction(10*1e3)
Mx = CoefficientFunction(M*cos(atan2(y,x))*1e0)
My = CoefficientFunction(0)
f += (Mx*grad(v)[1]-My*grad(v)[0]) * dx("magnet")
### "SIMULATING" the first Newton step only
# define gfa upfront
gfa = GridFunction(fes)
with TaskManager():
# check that AssembleLinearization() does not give NaNs
a.AssembleLinearization(gfa.vec)
f.Assemble()
# Solve the linear system
gfa.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Hope that helps!
The following user(s) said Thank You: lahoussaine
Time to create page: 0.113 seconds