[Getdp] Magnetodynamic 3D - Electromagnetic circuit element
Florin CIUPRINA
florin at elmat.pub.ro
Tue Apr 22 15:51:34 CEST 2003
Hello Christophe,
>You should try to increase Nb_Fill: this will make the incomplete LU
>decomposition (used as a preconditioner) be closer and closer to the
>complete LU. This should improve convergence a lot.
Thanks. Indeed, this has improved the convergence a lot.
I tried several resolutions for different frequencies between 10 and 1.e5 Hz.
Looking at the results I noticed that, when the frequency increases, the J
vectors
"flow" closer to one L-shaped surface, so looking along X axis there is an
asymmetry,
the J field lines preferring the L-shape surface at z=0. But from the
symmetry point
of view, J should have the same behaviour near z=zmax.
With other words I have a symmetrical problem with respect to the geometry
and boundary conditions (a=0 on all the boundary of the conductor,
and v=1 on one terminal and v=0 on the other terminal) and the distribution
of J vectors is
asymmetrical when the frequency increases. I refined the mesh for having
element size lower
than the skin depth, but the results are still asymmetrical.
I have no idea where this assymmetry comes from, therefore I attach again
the problem,
hoping that you will find a second to take a look at it. Your opinion would
help me a lot.
Best regards,
Florin
Florin Ciuprina, Ph.D.
Associate Professor
----------------------------------------------------------------------------
POLITEHNICA University of Bucharest
Electrical Engineering Faculty
Electrical Materials Laboratory
Spl.Independentei 313, 77206, Bucharest, Romania
Tel: +40.21.402 92 91
Fax: +40.21.410 43 55
e-mail: florin at elmat.pub.ro
----------------------------------------------------------------------------
-------------- next part --------------
/*
Matrix_Format (Integer):
- 1 Sparse
- 2 Full
- default : 1
Matrix_Printing (Integer): Disk write ('fort.*')
- 1 matrix (csr)
- 2 preconditioner (msr)
- 3 both
- default : 0
Matrix_Storage (Integer): Disk Write or Read in internal format
- 0 none
- 1 write matrix (sparse)
- 2 read matrix (sparse)
- default : 0
Scaling (Integer): Scale system
- 0 no
- 1 on basis of diagonal elements (no loss of possible symmetry)
- 2 on basis of inf. norm of first rows and then columns (asymmetric)
- 3 on basis of norm 1 of first rows and then columns (asymmetric)
- 4 on basis of norm 2 of first rows and then columns (asymmetric)
- default : 0
Renumbering_Technique (Integer):
- 0 No renumbering
- 1 Reverse Cuthill-Mc Kee
- default : 1
Preconditioner (Integer):
- 0 NONE No Factorization
- 1 ILUT Incomplete LU factorization with dual truncation strategy
- 2 ILUTP ILUT with column pivoting
- 3 ILUD ILU with single dropping + diagonal compensation (~MILUT)
- 4 ILUDP ILUD with column pivoting
- 5 ILUK level-k ILU
- 6 ILU0 simple ILU(0) preconditioning
- 7 MILU0 MILU(0) preconditioning
- 8 DIAGONAL
- default : 2
Preconditioner_Position (Integer):
- 0 No Preconditioner
- 1 Left Preconditioner
- 2 Right Preconditioner
- 3 Both Left and Right Preconditioner
- default : 2
Nb_Fill (Integer):
- ILUT/ILUTP : maximum number of elements per line
of L and U (except diagonal element)
- ILUK : each element whose fill-in level is greater than NB_FILL
is dropped.
- default : 20
Permutation_Tolerance (Real): Tolerance for column permutation in ILUTP/ILUDP.
At stage i, columns i and j are permuted if
abs(a(i,j))*PERMUTATION_TOLERANCE > abs(a(i,i)).
- 0 no permutations
- 0.001 -> 0.1 classical
- default : 0.05
Dropping_Tolerance (Real):
- ILUT/ILUTP/ILUK: a(i,j) is dropped if
abs(a(i,j)) < DROPPING_TOLERANCE * abs(diagonal element in U).
- ILUD/ILUDP : a(i,j) is dropped if
abs(a(i,j)) < DROPPING_TOLERANCE * [weighted norm of line i].
Weighted norm = 1-norm / number of nonzero elements on the line.
- default : 0
Diagonal_Compensation (Real): ILUD/ILUDP: the term 'DIAGONAL_COMPENSATION * (sum
of all dropped elements of the line)' is added to the diagonal element in U
- 0 ~ ILU with threshold
1 ~ MILU with threshold.
- default : 0
Re_Use_ILU (Integer): Reuse ILU decomposition (and renumbering if any)
- 0 no
- 1 yes
- default : 0
Algorithm (Integer):
- 1 CG Conjugate Gradient
- 2 CGNR CG (Normal Residual equation)
- 3 BCG Bi-Conjugate Gradient
- 4 DBCG BCG with partial pivoting
- 5 BCGSTAB BCG stabilized
- 6 TFQMR Transpose-Free Quasi-Minimum Residual
- 7 FOM Full Orthogonalization Method
- 8 GMRES Generalized Minimum RESidual
- 9 FGMRES Flexible version of GMRES
- 10 DQGMRES Direct versions of GMRES
- 11 LU LU Factorization
- 12 PGMRES Alternative version of GMRES
- default : 8
Krylov_Size (Integer): Krylov subspace size
- default : 40
IC_Acceleration (Real): IC accelerator
- default : 1
Re_Use_LU (Integer): Reuse LU decomposition
- 0 no
- 1 yes
- default : 0
Iterative_Improvement (Integer): Iterative improvement of the solution obtained by a LU
- default : 0
Nb_Iter_Max (Integer): Maximum number of iterations
- default : 1000
Stopping_Test (Real): Target relative residual
- default : 1e-10
*/
Matrix_Format 1
Matrix_Printing 0
Matrix_Storage 0
Scaling 4
Renumbering_Technique 1
Preconditioner 2
Preconditioner_Position 2
Nb_Fill 60
Permutation_Tolerance 0.05
Dropping_Tolerance 0
Diagonal_Compensation 0
Re_Use_ILU 0
Algorithm 10
Krylov_Size 40
IC_Acceleration 1
Re_Use_LU 0
Iterative_Improvement 0
Nb_Iter_Max 20000
Stopping_Test 1e-10
-------------- next part --------------
/* -------------------------------------------------------------------
File "L3Dind.pro"
This file defines the problem dependent data structures for
an Magnetodynamic 3D problem.
------------------------------------------------------------------- */
Group {
/* Let's start by defining the interface (i.e. Elementary groups)
between the mesh file and GetDP (no mesh object is defined, so
the default mesh will be assumed to be in GMSH format and located
in "L3D.msh") */
Conductor = Region[100];
Ground = Region[111]; Hot = Region[110];
Surf_Lat = Region[112];
/* We can then define a global group (used in "MagDyn_av.pro",
the file containing the function spaces and formulations) */
DomainC_Mag = Region[{Conductor}];
Surf = Region[{Surf_Lat, Ground, Hot}];
Electrode_Mag = Region[{Hot}] ; // only used with global quantities
}
Function {
/* The relative permeability (needed in the formulation) is piecewise
defined in Elementary groups */
sigma[Conductor] = 5.e7;
mur[Conductor] = 1.;
mu0 = 4.e-7 * Pi;
nu[] = 1. / (mu0 * mur[]);
Freq1 = 1.e1;
Freq2 = 1.e2;
Freq3 = 1.e3;
Freq4 = 1.e4;
Freq5 = 1.e5;
Tc = 1.;
Mag_Time0 = 0.; Mag_TimeMax = 1.*Tc; Mag_DTime[] = Tc/20.;
Mag_Theta[] = 1.;
}
Constraint {
/*
A 1-form function space containing vector potentials can be associated with a gauge condition,
which can be defined as a constraint, e.g. a zero value is fixed for all circulations along edges
of a tree (EdgesOfTreeIn) built in the mesh (Domain), having to be complete on certain
boundaries (StartingOn Surf)
*/
{ Name GaugeCondition_a_Mag_3D; Type Assign;
Case {
{ Region DomainC_Mag ; SubRegion Surf; Value 0.; }
}
}
/* UNUSED
{ Name MagneticVectorPotential_3D;
Case {
{ Region Surf; Value 0.; }
}
}
*/
/* Now, some Dirichlet conditions are defined. The name
'ElectricScalarPotential' refers to the constraint name given in
the function space */
{ Name ElectricScalarPotential; Type Assign;
Case {
{ Region Ground; Value 0.; }
}
}
{ Name GlobalElectricPotential ;
Case {
{ Region Hot ; Value 1.; }
}
}
}
/* The formulation used and its tools, considered as being
in a black box, can now be included */
Include "Jacobian_Lib.pro"
Include "Integration_Lib.pro"
Include "MagDyn_av.pro"
/* Finally, we can define some operations to output results */
PostOperation Map_vf_v UsingPost MagDyn_avf {
Print[ v, OnElementsOf DomainC_Mag, File "Map_v.pos" ] ;
}
PostOperation Map_vf_j UsingPost MagDyn_avf {
Print[ j, OnElementsOf DomainC_Mag, File "Map_j.pos" ] ;
}
PostOperation Map_vf_h UsingPost MagDyn_avf {
Print[ h, OnElementsOf DomainC_Mag, File "Map_h.pos" ] ;
}
PostOperation Map_vf_a UsingPost MagDyn_avf {
Print[ a, OnElementsOf DomainC_Mag, File "Map_a.pos" ] ;
}
PostOperation Map_vf_jg UsingPost MagDyn_avf {
Print [ j, OnPlane {{0,1.5,0}{1,1.5,0}{0,1.5,1}} {10,10}, File "Cut_j" ];
}
PostOperation Val_V UsingPost MagDyn_avf {
Print[ V, OnRegion Electrode_Mag, Format Table, File "VIZ.txt" ] ;
Print[ I, OnRegion Electrode_Mag, Format Table, File >> "VIZ.txt" ] ;
Print[ Z, OnRegion Electrode_Mag, Format Table, File >> "VIZ.txt" ] ;
}
// For test: post-calculation of a global quantity
// - can be used if the GlobalTerm with I is not included in the Formulation).
// - a function of type vf (in a SubSpace) must be defined.
// - OnRegion is for BasisFunction entity number.
// - the whole domain (DomainC_Mag) is the domain of integration.
PostOperation Ipos UsingPost MagDyn_avf {
Print[ Ipos[DomainC_Mag], OnRegion Hot, Format Table, File "Ipos.txt" ] ;
}
PostOperation Power UsingPost MagDyn_avf {
Print [ roj2[Conductor], OnGlobal, Format Table, File "Power.txt" ];
}
-------------- next part --------------
/* -------------------------------------------------------------------
File "MagDyn_v.pro"
Magnetodynamic - Magnetic vector potential and electric scalar
potential a-v formulation
with floating electric potential
-------------------------------------------------------------------
I N P U T
---------
Global Groups : (Extension '_Mag' is for Magnetic problem)
-------------
Domain_Mag Whole magnetic domain
DomainC_Mag Conducting regions
Electrode_Mag Electrodes (surfaces)
Surf Lateral surfaces (Surf_Lat) + Electrodes
Function :
--------
mur[] Relative permeability
sigma[] Conductivity
Constraint :
----------
MagneticVectorPotential Fixed magnetic vector potential
(boundary condition)
ElectricScalarPotential Fixed electric scalar potential
(classical boundary condition)
GlobalElectricPotential Fixed global electric scalar potential
(for perfect electrodes)
GlobalElectricCurrent Fixed global electric current
(for perfect electrodes)
/* O U T P U T
-----------
PostQuantities :
--------------
v : Electric scalar potential
a : Magnetic vector potential
h : Magnetic field
b : Magnetic flux density
j : Electric current density
V : Voltage
I : Current
Z : Impedance
*/
/* --------------------------------------------------------------------------*/
Group {
DefineGroup[ Domain_Mag, DomainC_Mag, Electrode_Mag, Surf ] ;
}
Function {
DefineFunction[ mur ];
DefineFunction[ sigma ];
DefineFunction[ nu ];
DefineVariable[ Freq1, Freq2, Freq3, Freq4, Freq5 ];
DefineVariable[Mag_Time0, Mag_TimeMax, Mag_DTime, Mag_Theta];
}
/* --------------------------------------------------------------------------*/
/* Edge finite element space with gauge condition */
FunctionSpace {
{ Name Hcurl_a_Gauge; Type Form1;
BasisFunction {
// a = a s
// e e
{ Name se; NameOfCoef ae; Function BF_Edge;
Support DomainC_Mag ; Entity EdgesOf[All]; }
}
Constraint {
//gauge condition ae = 0 on tree edges and on some boundaries
{ NameOfCoef ae;
EntityType EdgesOfTreeIn; EntitySubType StartingOn;
NameOfConstraint GaugeCondition_a_Mag_3D; }
}
}
/* Edge finite element space */
/* UNUSED
FunctionSpace {
{ Name Hcurl_a_Gauge; Type Form1;
BasisFunction {
// a = a s
// e e
{ Name se; NameOfCoef ae; Function BF_Edge;
Support DomainC_Mag ; Entity EdgesOf[All]; }
}
Constraint {
//condition ae = 0 on some boundaries
{ NameOfCoef ae; EntityType EdgesOf;
NameOfConstraint MagneticVectorPotential_3D; }
}
}
*/
{ Name Hgrad_vf_Mag ; Type Form0 ;
BasisFunction {
// v = v s + v s
// n n c,k c,k
{ Name sn ; NameOfCoef vn ; Function BF_Node ;
Support DomainC_Mag ; Entity NodesOf[ All, Not Electrode_Mag ] ; }
{ Name sck ;NameOfCoef vck ; Function BF_GroupOfNodes ;
Support DomainC_Mag ; Entity GroupsOfNodesOf[ Electrode_Mag ] ; }
}
SubSpace { // only for a PostOperation
{ Name vf ; NameOfBasisFunction sck ; }
}
GlobalQuantity {
{ Name V ; Type AliasOf ; NameOfCoef vck ; }
{ Name I ; Type AssociatedWith ; NameOfCoef vck ; }
}
Constraint {
{ NameOfCoef vn ;
EntityType NodesOf ; NameOfConstraint ElectricScalarPotential ; }
{ NameOfCoef V ;
EntityType GroupsOfNodesOf ; NameOfConstraint GlobalElectricPotential ; }
{ NameOfCoef I ;
EntityType GroupsOfNodesOf ; NameOfConstraint GlobalElectricCurrent ; }
}
}
}
Formulation {
{ Name MagDyn_avf; Type FemEquation ;
Quantity {
{ Name a ; Type Local ; NameOfSpace Hcurl_a_Gauge ; }
{ Name v ; Type Local ; NameOfSpace Hgrad_vf_Mag ; }
{ Name I ; Type Global ; NameOfSpace Hgrad_vf_Mag [I] ; }
{ Name V ; Type Global ; NameOfSpace Hgrad_vf_Mag [V] ; }
{ Name vf; Type Local ; NameOfSpace Hgrad_vf_Mag [vf] ; }
}
Equation {
Galerkin { [ nu[] * Dof{d a} , {d a} ]; In DomainC_Mag;
Jacobian Vol; Integration CurlCurl; }
Galerkin { DtDof [ sigma[] * Dof{a} , {a} ]; In DomainC_Mag;
Jacobian Vol; Integration CurlCurl; }
Galerkin { [ sigma[] * Dof{d v} , {a} ]; In DomainC_Mag;
Jacobian Vol; Integration CurlCurl; }
Galerkin { DtDof [ sigma[] * Dof{a} , {d v} ]; In DomainC_Mag;
Jacobian Vol; Integration CurlCurl; }
Galerkin { [ sigma[] * Dof{d v} , {d v} ]; In DomainC_Mag;
Jacobian Vol; Integration CurlCurl; }
GlobalTerm { [ Dof{I} , {V} ]; In Electrode_Mag; }
}
}
}
Resolution {
/* Frequency resolution */
{ Name MagDyn_vf;
System {
{ Name Sys_MagDyn; NameOfFormulation MagDyn_avf;
Type ComplexValue; Frequency {Freq1, Freq2, Freq3, Freq4, Freq5}; }
// Type ComplexValue; Frequency {Freq4}; }
}
Operation {
Generate[Sys_MagDyn]; Solve[Sys_MagDyn]; SaveSolution[Sys_MagDyn];
}
}
/* Time domain resolution */
{ Name MagDyn_avf_Time;
System {
{ Name Sys_MagDyn; NameOfFormulation MagDyn_avf; }
}
Operation {
InitSolution[Sys_MagDyn]; SaveSolution[Sys_MagDyn];
TimeLoopTheta[Mag_Time0, Mag_TimeMax, Mag_DTime[], Mag_Theta[]] {
Generate[Sys_MagDyn]; Solve[Sys_MagDyn]; SaveSolution[Sys_MagDyn];
}
}
}
}
PostProcessing {
{ Name MagDyn_avf ; NameOfFormulation MagDyn_avf ;
PostQuantity {
{ Name v ; Value { Term { [ {v} ] ; In DomainC_Mag ; } } }
{ Name h ; Value { Term { [ nu[] * {d a} ] ; In DomainC_Mag ; } } }
{ Name a ; Value { Term { [ {a} ] ; In DomainC_Mag ; } } }
// { Name j ; Value { Term { [ - sigma[]*(Dt[{a}]+ {d v}) ]; In DomainC_Mag ; } } }
{ Name j ;
Value {
Local { [ - sigma[]*(Dt[{a}]+ {d v}) ]; In DomainC_Mag ; Jacobian Vol; }
}
}
{ Name V ; Value { Term { [ {V} ] ; In Electrode_Mag ; } } }
{ Name I ; Value { Term { [ {I} ] ; In Electrode_Mag ; } } }
{ Name Z ; Value { Term { [ -{V}/{I} ] ; In Electrode_Mag ; } } }
{ Name Ipos ;
Value {
Integral {
// [ - sigma[] * {d v} * BF{d vf} ] ;
[ - sigma[]*(Dt[{a}]+ {d v}) * BF{d vf} ] ;
In DomainC_Mag ; Jacobian Vol ; Integration CurlCurl ;
}
}
}
{ Name roj2;
Value {
Integral { [ sigma[]*SquNorm[Dt[{a}]+{d v}] ]; In DomainC_Mag;
Jacobian Vol; Integration CurlCurl;}
}
}
}
}
}
/* --------------------------------------------------------------------------*/
/* --------------------------------------------------------------------------*/
-------------- next part --------------
/* -------------------------------------------------------------------
File "L3D.geo"
This file is the geometrical description used by GMSH to produce
the file "L3D.msh".
------------------------------------------------------------------- */
/* Definition of some parameters for geometrical dimensions, i.e.
w (width of 'Line')*/
w = 1.e-3 ;
h = 1.e-3 ;
/* Definition of parameters for local mesh dimensions */
s = 1. ;
p0 = w / 6. * s ;
/* Definition of gemetrical points */
Point(1) = { 0, 0, 0, p0} ;
Point(2) = { w, 0, 0, p0} ;
Point(3) = { 2*w, 0, 0, p0} ;
Point(4) = { 2*w, w, 0, p0} ;
Point(5) = { w, w, 0, p0} ;
Point(6) = { 0, w, 0, p0} ;
Point(7) = { 0, 2*w, 0, p0} ;
Point(8) = { w, 2*w, 0, p0} ;
Point(21) = { 0, 0, h, p0} ;
Point(22) = { w, 0, h, p0} ;
Point(23) = { 2*w, 0, h, p0} ;
Point(24) = { 2*w, w, h, p0} ;
Point(25) = { w, w, h, p0} ;
Point(26) = { 0, w, h, p0} ;
Point(27) = { 0, 2*w, h, p0} ;
Point(28) = { w, 2*w, h, p0} ;
/* Definition of gemetrical lines */
/* L in the plane z = 0*/
Line(1) = {1,2}; Line(2) = {2,5}; Line(3) = {5,6};
Line(4) = {6,1}; Line(5) = {2,3}; Line(6) = {3,4};
Line(7) = {4,5}; Line(8) = {5,8}; Line(9) = {6,7};
Line(10) = {7,8};
/* L in the plane z = h*/
Line(21) = {21,22}; Line(22) = {22,25}; Line(23) = {25,26};
Line(24) = {26,21}; Line(25) = {22,23}; Line(26) = {23,24};
Line(27) = {24,25}; Line(28) = {25,28}; Line(29) = {26,27};
Line(30) = {27,28};
/* Lines between L-s*/
Line(41) = {1, 21}; Line(42) = {2, 22};
Line(43) = {3, 23}; Line(44) = {4, 24};
Line(45) = {5, 25}; Line(46) = {6, 26};
Line(47) = {7, 27}; Line(48) = {8, 28};
/* Definition of geometrical surfaces */
/* Squares on L at z = 0 */
Line Loop(11) = {1,2,3,4}; Plane Surface(12) = {11};
Line Loop(13) = {5,6,7,-2}; Plane Surface(14) = {13};
Line Loop(15) = {-3,8,-10,-9}; Plane Surface(16) = {15};
/* Squares on L at z = h */
Line Loop(31) = {21,22,23,24}; Plane Surface(32) = {31};
Line Loop(33) = {25,26,27,-22}; Plane Surface(34) = {33};
Line Loop(35) = {-23,28,-30,-29}; Plane Surface(36) = {35};
/* Squares between L-s */
Line Loop(51) = {1, 42, -21, -41}; Plane Surface(52) = {51};
Line Loop(53) = {5, 43, -25, -42}; Plane Surface(54) = {53};
Line Loop(55) = {-3, 45, 23, -46}; Plane Surface(56) = {55};
Line Loop(57) = {-7, 44, 27, -45}; Plane Surface(58) = {57};
Line Loop(59) = {10, 48, -30, -47}; Plane Surface(60) = {59};
Line Loop(61) = {-9, 46, 29, -47}; Plane Surface(62) = {61};
Line Loop(63) = {4, 41, -24, -46}; Plane Surface(64) = {63};
Line Loop(65) = {-8, 45, 28, -48}; Plane Surface(66) = {65};
Line Loop(67) = {-2, 42, 22, -45}; Plane Surface(68) = {67};
Line Loop(69) = {-6, 43, 26, -44}; Plane Surface(70) = {69};
/* Definition of geometrical volumes */
Surface Loop(101) = {12, 32, 52, 68, 56, 64}; Volume(102) = {101};
Surface Loop(103) = {14, 34, 68, 54, 70, 58}; Volume(104) = {103};
Surface Loop(105) = {16, 36, 56, 66, 60, 62}; Volume(106) = {105};
/* Definition of Physical entities (volumes, surfaces). The Physical
entities tell GMSH the elements and their associated region numbers
to save in the file 'L3D.msh'. */
Physical Surface (110) = {60} ; /* Hot */
Physical Surface (111) = {70} ; /* Ground */
Physical Surface (112) = {12,14,16,32,34,36,62,64,52,54,58,66} ; /* Surf_Lat */
Physical Volume (100) = {102, 104, 106}; /* Mag + Conductor */