Commit fd94fb85 authored by Siamak Pourkeivannour's avatar Siamak Pourkeivannour
Browse files

DC is working

AC converges to a wrong result
parents
*.txt
*.mph
*.png
*.lock
Ccore GMSH/
From Andrea/
2D_MH_lin.py
import nutils, numpy
from matplotlib import pyplot as plt
topo, geom = nutils.mesh.rectilinear([numpy.linspace(0, 2, 10)])
ns = nutils.function.Namespace()
ns.x = geom
ns.basis = topo.basis('spline', degree=1)
ns.u = 'basis_n ?lhs_n'
sqr = topo.boundary['left'].integral('(u - 0)^2 d:x' @ ns, degree=1)
cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)
res = topo.integral(' (basis_n,i u_,i ) d:x' @ ns, degree=0)
res -= topo.boundary['right'].integral('basis_n d:x' @ ns, degree=1)
lhs = nutils.solver.solve_linear('lhs', res, cons)
bezier = topo.sample('bezier', 32)
nanjoin = lambda array, tri: numpy.insert(array.take(tri.flat, 0).astype(float), slice(tri.shape[1], tri.size, tri.shape[1]), numpy.nan, axis=0)
sampled_x = nanjoin(bezier.eval(ns.x[0]), bezier.tri)
def plot_line(func, **arguments):
plt.plot(sampled_x, nanjoin(bezier.eval(func, **arguments), bezier.tri))
plt.xlabel('$x_0$')
plt.xticks(numpy.linspace(0, 2, 11))
plot_line(ns.basis)
plt.show()
plot_line(ns.u, lhs=lhs)
plt.show()
import nutils
from matplotlib import pyplot as plt
import numpy as np
L_copper = 0.06
L_air = 0.03
nelems = 1, 3, 1, 3, 1
lengths = L_air, L_copper, L_air, L_copper, L_air
subnames = 'Air1', 'Copper1', 'Air2', 'Copper2', 'Air3'
subverts = np.cumsum([0, *lengths])
verts = np.concatenate([np.linspace(l, r, n+1)[:-1] for l, r, n in zip(subverts[:-1], subverts[1:], nelems)]+[subverts[-1:]])
domain, geom = nutils.mesh.rectilinear([verts])
hverts = np.cumsum([0,*nelems])
domain = domain.withgroups(vgroups={k: domain[l:r] for k, l, r in zip(subnames, hverts[:-1], hverts[1:])})
ns = nutils.function.Namespace()
ns.x = geom
bezier = domain.sample('bezier', 2)
nanjoin = lambda array, tri: np.insert(array.take(tri.flat, 0).astype(float), slice(tri.shape[1], tri.size, tri.shape[1]), np.nan, axis=0)
sampled_x = nanjoin(bezier.eval(ns.x[0]), bezier.tri)
x = bezier.eval('x_0' @ ns)
print(x.take(bezier.tri, 0))
domain.boundary['left']
def plot_line(func, **arguments):
plt.plot(sampled_x, nanjoin(bezier.eval(func, **arguments), bezier.tri))
plt.xlabel('x_0')
plt.xticks(np.linspace(0, 0.222, 10))
ns.basis = domain.basis('spline', degree=1)
plot_line(ns.basis)
plt.show()
import nutils, numpy
import matplotlib.pyplot as plt
import matplotlib.collections
from nutils import cli, mesh, export
import time
__matrix__ = nutils.matrix.MKL() or nutils.matrix.Scipy() or nutils.matrix.Numpy()
__matrix__.__enter__()
def main():
pi = numpy.pi
btype = 'std'
degree = 2
m = 2 # number of winding layers
L = 50e-3
nelems = 10
nelems_x = nelems,nelems,nelems,nelems,nelems
nelems_y = nelems,nelems,nelems
lengths_x = L, L, L, L, L
lengths_y = L, L, L
subnames_x = 'Air', 'CondX1', 'Air', 'CondX2', 'Air'
subnames_y = 'Air', 'CondY', 'Air'
subverts_x = numpy.cumsum([0, *lengths_x])
verts_x = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_x[:-1], subverts_x[1:], nelems_x)] + [
subverts_x[-1:]])
subverts_y = numpy.cumsum([0, *lengths_y])
verts_y = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_y[:-1], subverts_y[1:], nelems_y)] + [
subverts_y[-1:]])
domain, geom = nutils.mesh.rectilinear([verts_x, verts_y])
yverts = numpy.cumsum([0, *nelems_y])
domain = domain.withgroups(vgroups={k: domain[:, l:r] for k, l, r in zip(subnames_y, yverts[:-1], yverts[1:])})
xverts = numpy.cumsum([0, *nelems_x])
domain = domain.withgroups(
vgroups={k: domain['CondY'][l:r] for k, l, r in zip(subnames_x, xverts[:-1], xverts[1:])})
ns = nutils.function.Namespace()
ns.x = geom
# with export.mplfigure('geometry.png') as fig:
# # Add geometry
# ax = fig.add_subplot(1, 1, 1)
# #ax.set_aspect('equal')
# smpl_gauss = domain.sample('gauss',degree*2)
# smpl_bezier = domain.sample('bezier',15)
# gauss = smpl_gauss.eval(geom)
# bezier = smpl_bezier.eval(geom)
# ax.add_collection(matplotlib.collections.LineCollection(bezier[smpl_bezier.hull], colors='k', linewidths=0.25, alpha=.25))
# # Add Gauss points
# #for i, q in enumerate(gauss):
# # ax.text(q[0] - 1e-3, q[1] + 0.25e-3, i, fontsize=4)
# # ax.plot(q[0], q[1], 'k.', markersize=1)
#
# ax.set_xlim([0, 175e-3])
# ax.set_ylim([0, 150e-3])
# Define helper function for creating values that are constant per subdomain.
subdomconst = lambda *, default=0, **values: sum(v * domain.indicator(k) for k, v in values.items()) + default
start_time = time.time()
freq = 10e3
µ_0 = 4 * pi * 1e-7
σ_C = 5.998e7 # Copper conductivity
ns.SigmaCond = σ_C
omega = 2 * pi * freq
I = 1 #coil current
Cond_Area = L * L
J_dc = I/Cond_Area
ns.µ0 = µ_0
ns.ϑ0 = 1 / µ_0
ns.σC = σ_C
ns.omega = omega
ns.J = J_dc
ns.mu = subdomconst(default=1 * ns.µ0)
ns.nu = subdomconst(default=1 / ns.µ0)
ns.sigma = subdomconst(CondX1=ns.σC, CondX2=ns.σC, default=0)
ns.Jdc = subdomconst( CondX1=ns.J, CondX2=-ns.J, default=0)
spline_km_y = numpy.ones([sum(nelems_y) + 1], int)
spline_km_y[numpy.cumsum([0, *nelems_y])] = degree
ns.ucbasis, ns.usbasis = nutils.function.chain([
domain.basis('spline', degree=degree, knotvalues=[verts_x, verts_y], knotmultiplicities=[None, spline_km_y]),
domain.basis('spline', degree=degree, knotvalues=[verts_x, verts_y], knotmultiplicities=[None, spline_km_y])])
# Setup magnetics problem
ns.eps = numpy.array([[0, 1], [-1, 0]]) # 2D
ns.curlucbasis_ni = 'eps_ij ucbasis_n,j' # <<< ?
ns.curlusbasis_ni = 'eps_ij usbasis_n,j'
ns.uc = 'ucbasis_n ?lhs_n'
ns.us = 'usbasis_n ?lhs_n'
ns.curluc_i = 'eps_ij uc_,j'
ns.curlus_i = 'eps_ij us_,j'
res = domain.integral('( curlucbasis_ni curluc_i + curlusbasis_ni curlus_i ) d:x' @ ns, degree=degree * 2)
res -= domain.integral('( mu sigma ucbasis_n us - mu sigma usbasis_n uc) d:x' @ ns, degree=degree * 2)
res += domain.integral('(mu ucbasis_n Jdc ) d:x' @ ns, degree=degree * 2) # careful with parenthesis
# Solve the problem
lhs = nutils.solver.solve_linear('lhs', res)
ns.A = ' sqrt( uc uc + us us )'
bezierC = domain.sample('bezier', 5)
xC, J_dc, sigma, Az_c, Az_s, A = bezierC.eval(['x_i', 'Jdc', 'sigma', 'uc', 'us', 'A'] @ ns, lhs=lhs)
print("--- %s seconds ---" % (time.time() - start_time))
# Plot
plt.figure(1)
plt.tripcolor(xC[:, 0] * 1000, xC[:, 1] * 1000, bezierC.tri, A, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Magnetic Vector Potential [Wb/m]')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
matplotlib.rcParams["figure.dpi"] = 0.1
plt.show()
# # Define a 2d list or numpy matrix with coordinates
# x_end = 5*L
# y_end = 3*L
# x_coordinates = numpy.linspace(0,x_end,250)
# y_coordinates = numpy.linspace(0,y_end,150)
# coordinates = numpy.array(numpy.meshgrid(x_coordinates, y_coordinates)).T.reshape(-1, 2)
# # Define a sampling object
# located_points = domain.locate(geom, coordinates, tol=1e-12)
# # Evaluate the sampling object for a desired quantity
# Az = numpy.reshape(located_points.eval('u' @ ns, lhs=lhs), (-1,150)).T
# print(len(Az))
# Az_comsol = numpy.loadtxt('Az_COMSOL.txt')
# print(len(Az_comsol))
#
# error = 100*numpy.true_divide((Az - Az_comsol),Az_comsol)
# error_max = numpy.max(error)
# print(error_max)
# error_min = numpy.min(error)
# print(error_min)
#
# fig, ax = plt.subplots()
# c = ax.pcolormesh(x_coordinates, y_coordinates, error, cmap='RdBu', vmin=error_min, vmax=error_max)
# ax.set_title('Error (%)')
# # set the limits of the plot to the limits of the data
# ax.axis([x_coordinates.min(), x_coordinates.max(), y_coordinates.min(), y_coordinates.max()])
# fig.colorbar(c, ax=ax)
# plt.show()
if __name__ == '__main__':
nutils.cli.run(main)
import nutils, numpy
import matplotlib.pyplot as plt
import matplotlib.collections
from nutils import cli, mesh, export
import time
__matrix__ = nutils.matrix.MKL() or nutils.matrix.Scipy() or nutils.matrix.Numpy()
__matrix__.__enter__()
def main():
pi = numpy.pi
btype = 'std'
degree = 2
m = 2 # number of winding layers
L = 50e-3
nelems = tuple([100])
nelems_x = tuple([100])
nelems_y = tuple([100])
lengths_x = tuple([50e-3])
lengths_y = tuple([50e-3])
subnames_x = 'CondX'
subnames_y = 'CondY'
subverts_x = numpy.cumsum([0, *lengths_x])
verts_x = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_x[:-1], subverts_x[1:], nelems_x)] + [
subverts_x[-1:]])
subverts_y = numpy.cumsum([0, *lengths_y])
verts_y = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_y[:-1], subverts_y[1:], nelems_y)] + [
subverts_y[-1:]])
domain, geom = nutils.mesh.rectilinear([verts_x, verts_y])
#
# yverts = numpy.cumsum([0, *nelems_y])
# domain = domain.withgroups(vgroups={k: domain[:, l:r] for k, l, r in zip(subnames_y, yverts[:-1], yverts[1:])})
# xverts = numpy.cumsum([0, *nelems_x])
# domain = domain.withgroups(
# vgroups={k: domain['CondY'][l:r] for k, l, r in zip(subnames_x, xverts[:-1], xverts[1:])})
ns = nutils.function.Namespace()
ns.x = geom
# with export.mplfigure('geometry.png') as fig:
# # Add geometry
# ax = fig.add_subplot(1, 1, 1)
# #ax.set_aspect('equal')
# smpl_gauss = domain.sample('gauss',degree*2)
# smpl_bezier = domain.sample('bezier',15)
# gauss = smpl_gauss.eval(geom)
# bezier = smpl_bezier.eval(geom)
# ax.add_collection(matplotlib.collections.LineCollection(bezier[smpl_bezier.hull], colors='k', linewidths=0.25, alpha=.25))
# # Add Gauss points
# #for i, q in enumerate(gauss):
# # ax.text(q[0] - 1e-3, q[1] + 0.25e-3, i, fontsize=4)
# # ax.plot(q[0], q[1], 'k.', markersize=1)
#
# ax.set_xlim([0, 50e-3])
# ax.set_ylim([0, 50e-3])
# Define helper function for creating values that are constant per subdomain.
# subdomconst = lambda *, default=0, **values: sum(v * domain.indicator(k) for k, v in values.items()) + default
start_time = time.time()
freq = 10e3
µ_0 = 4 * pi * 1e-7
σ_C = 5.87 * 1e7 # Copper conductivity
omega = 2 * pi * freq
I = 1 #coil current
Cond_Area = L * L
J_dc = I/Cond_Area
ns.µ0 = µ_0
ns.ϑ0 = 1 / µ_0
ns.σC = σ_C
ns.omega = omega
ns.J = J_dc
ns.mu = µ_0
ns.nu = 1 / ns.µ0
ns.sigma = ns.σC
ns.Jdc = ns.J
ns.ucbasis, ns.usbasis = nutils.function.chain([
domain.basis(btype, degree=degree),
domain.basis(btype, degree=degree)])
# Setup magnetics problem
ns.eps = numpy.array([[0, 1], [-1, 0]]) # 2D
ns.curlucbasis_ni = 'eps_ij ucbasis_n,j' # <<< ?
ns.curlusbasis_ni = 'eps_ij usbasis_n,j'
ns.uc = 'ucbasis_n ?lhs_n'
ns.us = 'usbasis_n ?lhs_n'
ns.curluc_i = 'eps_ij uc_,j'
ns.curlus_i = 'eps_ij us_,j'
res = domain.integral('( curlucbasis_ni curluc_i + curlusbasis_ni curlus_i ) d:x' @ ns, degree=degree * 2)
res += domain.integral('( mu sigma omega ucbasis_n us - mu sigma omega usbasis_n uc) d:x' @ ns, degree=degree * 2)
res -= domain.integral('(mu ucbasis_n Jdc ) d:x' @ ns, degree=degree * 2) # careful with parenthesis
sqr = domain.boundary['top'].integral('(uc^2 + us^2) d:x' @ ns, degree=degree*2)
sqr += domain.boundary['right'].integral('(uc^2 + us^2) d:x' @ ns, degree=degree*2)
sqr += domain.boundary['left'].integral('(uc^2 + us^2) d:x' @ ns, degree=degree*2)
cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)
# Solve the problem
lhs = nutils.solver.solve_linear('lhs', res,constrain=cons)
print("--- %s seconds ---" % (time.time() - start_time))
bezierC = domain.sample('bezier', 8)
xC, J_dc, Az_c, Az_s = bezierC.eval(['x_i', 'Jdc', 'uc', 'us'] @ ns, lhs=lhs)
# Plot
plt.figure()
plt.tripcolor(xC[:, 0] * 1000, xC[:, 1] * 1000, bezierC.tri, Az_s, shading='gouraud', cmap='viridis')
plt.colorbar(orientation="vertical")
plt.title('Magnetic Vector Potential [Wb/m]')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
matplotlib.rcParams["figure.dpi"] = 0.1
plt.show()
if __name__ == '__main__':
nutils.cli.run(main)
import nutils, numpy
import matplotlib.pyplot as plt
import matplotlib.collections
from nutils import cli, mesh, export
__matrix__ = nutils.matrix.MKL() or nutils.matrix.Scipy() or nutils.matrix.Numpy()
__matrix__.__enter__()
def main():
pi = numpy.pi
btype = 'std'
degree = 2
m = 2 # number of winding layers
clear_core_LV = 20e-3
clear_core_HV = 20e-3
clear_LV_HV = 20e-3
clear_h = 20e-3
copper_w = 0.5e-3
copper_h = 50e-3
ins_w = 0.1e-3
nelems_Clear_core_LV = 50
nelems_Clear_core_HV = 50
nelems_Clear_LV_HV = 50
nelems_clear_h = 50
nelems_copper_w = 5
nelems_copper_h = 100
nelems_ins_w = 2
# Build tuple of number of elements
def build_geom():
nelems_x = []
nelems_x.append(nelems_Clear_core_LV)
ielems = 0
while ielems < m:
nelems_x.append(nelems_copper_w)
if ielems != m - 1:
nelems_x.append(nelems_ins_w)
ielems += 1
nelems_x.append(nelems_Clear_LV_HV)
ielems = 0
while ielems < m:
nelems_x.append(nelems_copper_w)
if ielems != m - 1:
nelems_x.append(nelems_ins_w)
ielems += 1
nelems_x.append(nelems_Clear_core_HV)
nelems_x = tuple(nelems_x)
nelems_y = nelems_clear_h, nelems_copper_h, nelems_clear_h
# Build tuple of length of elements
lengths_x = []
lengths_x.append(clear_core_LV)
ielems = 0
while ielems < m:
lengths_x.append(copper_w)
if ielems != m - 1:
lengths_x.append(ins_w)
ielems += 1
lengths_x.append(clear_LV_HV)
ielems = 0
while ielems < m:
lengths_x.append(copper_w)
if ielems != m - 1:
lengths_x.append(ins_w)
ielems += 1
lengths_x.append(clear_core_HV)
lengths_x = tuple(lengths_x)
lengths_x = clear_core_LV, copper_w, ins_w, copper_w, clear_LV_HV, copper_w, ins_w, copper_w, clear_core_HV
lengths_y = clear_h, copper_h, clear_h
# Build tuple of name of elements
subnames_x = []
subnames_x.append('Air')
ielems = 0
while ielems < m:
subnames_x.append('CondP' + str(ielems + 1))
if ielems != m - 1:
subnames_x.append('Air')
ielems += 1
subnames_x.append('Air')
ielems = 0
while ielems < m:
subnames_x.append('CondN' + str(ielems + 1))
if ielems != m - 1:
subnames_x.append('Air')
ielems += 1
subnames_x.append('Air')
subnames_x = tuple(subnames_x)
subnames_y = 'Air', 'CondY', 'Air'
return nelems_x, nelems_y, lengths_x, lengths_y, subnames_x, subnames_y
nelems_x, nelems_y, lengths_x, lengths_y, subnames_x, subnames_y = build_geom()
subverts_x = numpy.cumsum([0, *lengths_x])
verts_x = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_x[:-1], subverts_x[1:], nelems_x)] + [
subverts_x[-1:]])
subverts_y = numpy.cumsum([0, *lengths_y])
verts_y = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_y[:-1], subverts_y[1:], nelems_y)] + [
subverts_y[-1:]])
domain, geom = nutils.mesh.rectilinear([verts_x, verts_y])
yverts = numpy.cumsum([0, *nelems_y])
domain = domain.withgroups(vgroups={k: domain[:, l:r] for k, l, r in zip(subnames_y, yverts[:-1], yverts[1:])})
xverts = numpy.cumsum([0, *nelems_x])
domain = domain.withgroups(
vgroups={k: domain['CondY'][l:r] for k, l, r in zip(subnames_x, xverts[:-1], xverts[1:])})
ns = nutils.function.Namespace()
ns.x = geom
# with export.mplfigure('geometry.png') as fig:
# # Add geometry
# ax = fig.add_subplot(1, 1, 1)
# #ax.set_aspect('equal')
# smpl_gauss = domain.sample('gauss',degree*2)
# smpl_bezier = domain.sample('bezier',15)
# gauss = smpl_gauss.eval(geom)
# bezier = smpl_bezier.eval(geom)
# ax.add_collection(matplotlib.collections.LineCollection(bezier[smpl_bezier.hull], colors='k', linewidths=0.25, alpha=.25))
# # Add Gauss points
# #for i, q in enumerate(gauss):
# # ax.text(q[0] - 1e-3, q[1] + 0.25e-3, i, fontsize=4)
# # ax.plot(q[0], q[1], 'k.', markersize=1)
#
# ax.set_xlim([0, 175e-3])
# ax.set_ylim([0, 150e-3])
# Define helper function for creating values that are constant per subdomain.
subdomconst = lambda *, default=0, **values: sum(v * domain.indicator(k) for k, v in values.items()) + default
freq = 10e3
µ_0 = 4 * pi * 1e-7
σ_C = 5.998 * 1e7 # Copper conductivity
omega = 2 * pi * freq
I = 1 #coil current
Cond_Area = copper_w * copper_h
J_dc = I/Cond_Area
ns.µ0 = µ_0
ns.ϑ0 = 1 / µ_0
ns.σC = σ_C
ns.omega = omega
ns.J = J_dc
ns.mu = subdomconst(default=1 * ns.µ0)
ns.nu = subdomconst(default=1 / ns.µ0)
ns.sigma = subdomconst(CondP1=ns.σC, CondP2=ns.σC, CondN1=ns.σC, CondN2=ns.σC, default=0)
ns.Jdc = subdomconst(CondP1=ns.J, CondP2=ns.J, CondN1=-ns.J, CondN2=-ns.J, default=0)