Elmer Magnetics 3D comparison with 2D Axisymmetric Lumpy

[crobar]

Hello,

I've been continuing to work with the 3D WhitneyAV Magnetodynamics solver, and attempting to benchmark it against other programs. I do have access to a commercial 3D solver, but first I am running a test case using FEMM a high quality 2D solver. FEMM can handle axisymmetric problems in 2D which allows us to compare the results of Elmer in 3D to the same axisymmetric problem. I've made a test case which is a cylindrical magnet, magnetised in the axial direction, near a cylindrical coil wrapped around a cylinder of Iron (with a nonlinear B-H curve). These are all then enclosed in two boxes of air, an inner box and an outer box to allow a fine mesh close to the problem and looser further out. Some images of the geometry are shown below. The air boxes are only shown on the left, the transparent cylinders.

Images of problem geometry

geom.png (48.58 KiB) Viewed 7848 times

After solving the problem in Elmer, and FEMM, I took samples of the magnetic flux density and magnetic vector potential along a line going parallel to the z axis, with y=0 and x at a position which put the line 1mm outside the magnet and core surface,(and 1mm inside the coil inner diameter). I took 1000 points and used natrual neighbour interpolation to get the Elmer results. The results of this sample are shown in the figures below:

field_samples_elmer_vs_femm.png

(50.97 KiB) Not downloaded yet

In my post preview the image isn't appearing, but let's hope it works for the real thing.

Either way, in this image you can see that Elmer has much 'lumpier' magnetic vector potential results in FEMM (which I believe my do some smoothing), it also has very lumpy flux density in the z direction (parallel to the sample line). Is there anything that can be done to smooth the Elmer results?

There is also a significant difference in the B field in the tangential direction, even if we took a smoothed curve. Any comments? I would be pretty confident about FEMM's result, it has years of testing against real-life things.

These were produced with a pretty fine mesh around the region of interest (over 2 million elements). The entire case include FEMM input file is attached. I didn't include the mesh as it would be too big, so you would have to generate it using Salome from the hdf file, export to unv, and convert to elmer format using ElmerGrid to actually run the case.

[crobar]

In the previous post I mis-spoke, I was actually doing linear interpolation on the vector field, see the following results for natural interpolation, but not much has changed:

B_and_A_sample_comparisons_natural_interp.png

(52.34 KiB) Not downloaded yet

[raback]

Hi

There are two things you could do:

1) Use quadratic edge elements for AV solver:
In Solver section set "Quadratic Approximation = .TRUE." and use a quadratic nodal mesh for basis.

2) Use optimal interpolation for visualization within Elmer, see
viewtopic.phpf=3&t=4227&p=14936&hilit=q ... 5fc#p14936

It may be noted that in electromagnetics the Hcurl conforming edge elements in 3D results to a huge increase in complexity and computational effort. If many application areas (such as HeatSolver) the same machinery can deal with both 2D and 3D problems but the \curl\curl operator is quite a different creature in 2D and 3D.

-Peter

[crobar]

Hi Peter,

Thanks, I will try the quadratic mesh, although googling earlier indicated there might be an issue converting unv quadratic meshes from salome to Elmer format. I've started generating one now anyway to see how it goes.

Your second link leads to a Not Found (404), because it's missing the start of the url, I think you mean

viewtopic.php?f=3&t=4227&p=14936&hilit= ... b41#p14936

which does look useful, thanks!

I will try and report back on the different methods.

[crobar]

1) Use quadratic edge elements for AV solver:
In Solver section set "Quadratic Approximation = .TRUE." and use a quadratic nodal mesh for basis.

I've just had a go at this, but I'm getting the error:

Code: Select all

ERROR:: GetEdgeBasis: Can't handle but linear elements, sorry.

Here's the full sif file:

Code: Select all

$filename  = "elmer_vs_femm_quad"  

! ======================================================================!
!	HEADER								!
! ======================================================================!

Header
   CHECK KEYWORDS "Warn"
   Mesh DB "/home/rcrozier/build/fea_temp_files/elmer_vs_femm" $filename
End

! ======================================================================!
!		Simulation						!
! ======================================================================!

Simulation
   Max Output Level = 7
   Coordinate System = "Cartesian 3D"
   Coordinate Mapping(3) = 1 2 3 
   Simulation Type = Steady state
   Steady State Max Iterations = 1
   Output Intervals = 1
   Use Mesh Names = Logical True
End

! ======================================================================!
!		Constants						!
! ======================================================================!

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.67e-08
  Permittivity of Vacuum = 8.8542e-12
  Boltzmann Constant = 1.3807e-23
  Unit Charge = 1.602e-19
End

! ======================================================================!
!		Body, Material,	Body Force, Initial Cond.		!
! ======================================================================!

Material 1
  Name = "Air"
  Density = 1000
  Electric Conductivity = 0.0
  Relative Permeability = 1.0
End

Material 2
  Name = "Magnet"
  Density = 1000
  Electric Conductivity = 0.0
  Relative Permeability = 1.05
End

Material 3
  Name = "Iron"
  Density = 7500
  Electric Conductivity = 0.0
  Relative Permeability = 1000.0
  
  H-B Curve = Variable "dummy"
   Real  Cubic Monotone
     INCLUDE HB_Iron
   End
End

Body 1
   Name = "magnet"
   Equation = 1
   Material = 2
   Body Force = 1
End

Body 2
   Name = "core"
   Equation = 1
   Material = 3
End

Body 3
   Name = "coil"
   Equation = 1
   Material = 1
End

Body 4
   Name = "inner_air_box"
   Equation = 1
   Material = 1
End

Body 5
   Name = "outer_air_box"
   Equation = 1
   Material = 1
End

Body Force 1
   Name = "Magnetization"
   Magnetization 1 = Real 0.0
   Magnetization 2 = Real 0.0
   Magnetization 3 = Real 979000
End

! ======================================================================!
!			Equations and Solvers				!
! ======================================================================!

Solver 1
   Equation = "MGDynamics"

   Variable = "A"

   Procedure = "MagnetoDynamics" "WhitneyAVSolver"
   Fix Input Current Density = Logical True

   Newton-Raphson Iteration = Logical True
   Nonlinear System Max Iterations = 30
   Nonlinear System Convergence Tolerance = 1e-6

   Linear System Symmetric = Logical True
   Linear System Solver = "Iterative"
   Linear System Preconditioning = None
   Linear System Convergence Tolerance = 1e-8
   Linear System Residual Output = 100
   Linear System Max Iterations = 5000
   Linear System Iterative Method = CG
   Steady State Convergence Tolerance = 1e-8

   Quadratic Approximation = Logical True

End

Solver 2
   Equation = "MGDynamicsCalc"
   Procedure = "MagnetoDynamics" "MagnetoDynamicsCalcFields"
   Linear System Symmetric = True
   Potential Variable = String "A"

   Calculate Magnetic Field Strength = Logical True
   Calculate Magnetic Vector Potential = Logical True
   Steady State Convergence Tolerance = 0
   Linear System Solver = "Iterative"
   Linear System Preconditioning = None
   Linear System Residual Output = 0
   Linear System Max Iterations = 5000
   Linear System Iterative Method = CG
   Steady State Convergence Tolerance = 1e-8
   Linear System Convergence Tolerance = 1.0e-8

   Quadratic Approximation = Logical True

End

Solver 3
  Exec Solver = after all
  Equation = "ResultOutput"
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Output File Name = "/home/rcrozier/build/fea_temp_files/elmer_vs_femm/"$filename
  Save Geometry Ids = Logical True
  Vector Field 1 = String Magnetic Field Strength
  Vector Field 2 = String Magnetic Flux Density
  Vector Field 3 = String Magnetic Vector Potential
  Potential Variable = String av
  Show Variables = Logical True
  Vtu format = Logical True

  !Quadratic Approximation = Logical True
End


Equation 1
  Name = "Coupled Equations"
  Active Solvers(2) = 1 2
End


! ======================================================================!
!		Boundary Condition Section				!
! ======================================================================!

! ----- names for boundaries -----
!$ outer_bound = 1

Boundary Condition 1
  Name = "outer_boundary"
  A = Real 0.0
  A {e} = Real 0.0
End

Here's all the output:

Code: Select all

ELMER SOLVER (v 8.2) STARTED AT: 2016/08/10 15:49:05
ParCommInit:  Initialize #PEs:            1
MAIN: 
MAIN: =============================================================
MAIN: ElmerSolver finite element software, Welcome!
MAIN: This program is free software licensed under (L)GPL
MAIN: Copyright 1st April 1995 - , CSC - IT Center for Science Ltd.
MAIN: Webpage http://www.csc.fi/elmer, Email elmeradm@csc.fi
MAIN: Version: 8.2, Compiled: 2016-08-03)
MAIN: =============================================================
MAIN: 
MAIN: 
MAIN: -------------------------------------
MAIN: Reading Model: case_quad.sif
Model Input:  Unlisted keyword: [a] in section: [boundary condition 1]
Model Input:  Unlisted keyword: [a {e}] in section: [boundary condition 1]
Loading user function library: [MagnetoDynamics]...[WhitneyAVSolver_Init0]
Loading user function library: [MagnetoDynamics]...[MagnetoDynamicsCalcFields_Init0]
Loading user function library: [ResultOutputSolve]...[ResultOutputSolver_Init0]
Loading user function library: [MagnetoDynamics]...[MagnetoDynamics_Dummy_Init0]
LoadMesh: Base mesh name: /home/rcrozier/build/fea_temp_files/elmer_vs_femm/elmer_vs_femm_quad
LoadMesh: Performing node mapping
ReadTargetNames: Reading names info from file: /home/rcrozier/build/fea_temp_files/elmer_vs_femm/elmer_vs_femm_quad/mesh.names
WARNING:: ReadTargetNames: Could not map name to Body nor BC: bnry7
LoadMesh: Minimum initial body index: 2
LoadMesh: Maximum initial body index: 6
LoadMesh: Minimum initial boundary index: 1
LoadMesh: Maximum initial boundary index: 7
LoadMesh: Elapsed time (CPU,REAL):    16.4160   21.6721 (s)
MeshStabParams: Computing stabilization parameters
MeshStabParams: Elapsed time (CPU,REAL):     0.7560    0.7599 (s)
MAIN: -------------------------------------
Loading user function library: [MagnetoDynamics]...[WhitneyAVSolver_Init]
Loading user function library: [MagnetoDynamics]...[WhitneyAVSolver_bulk]
Loading user function library: [MagnetoDynamics]...[WhitneyAVSolver]
OptimizeBandwidth: ---------------------------------------------------------
OptimizeBandwidth: Computing matrix structure for: mgdynamics...done.
OptimizeBandwidth: Half bandwidth without optimization: 1954962
OptimizeBandwidth: 
OptimizeBandwidth: Bandwidth Optimization ...done.
OptimizeBandwidth: Half bandwidth after optimization: 124651
OptimizeBandwidth: ---------------------------------------------------------
Loading user function library: [MagnetoDynamics]...[MagnetoDynamicsCalcFields_Init]
Loading user function library: [MagnetoDynamics]...[MagnetoDynamicsCalcFields_bulk]
Loading user function library: [MagnetoDynamics]...[MagnetoDynamicsCalcFields]
OptimizeBandwidth: ---------------------------------------------------------
OptimizeBandwidth: Computing matrix structure for: mgdynamicscalc...done.
OptimizeBandwidth: Half bandwidth without optimization: 2241616
OptimizeBandwidth: 
OptimizeBandwidth: Bandwidth Optimization ...done.
OptimizeBandwidth: Half bandwidth after optimization: 198555
OptimizeBandwidth: ---------------------------------------------------------
Loading user function library: [ResultOutputSolve]...[ResultOutputSolver_Init]
Loading user function library: [ResultOutputSolve]...[ResultOutputSolver_bulk]
Loading user function library: [ResultOutputSolve]...[ResultOutputSolver]
Loading user function library: [MagnetoDynamics]...[MagnetoDynamics_Dummy_Init]
Loading user function library: [MagnetoDynamics]...[MagnetoDynamics_Dummy_bulk]
Loading user function library: [MagnetoDynamics]...[MagnetoDynamics_Dummy]
MAIN: 
MAIN: -------------------------------------
MAIN:  Steady state iteration:            1
MAIN: -------------------------------------
MAIN: 
SingleSolver: Attempting to call solver
SingleSolver: Solver Equation string is: mgdynamics
WhitneyAVSolver: Solving the AV equations with edge elements
OptimizeBandwidth: ---------------------------------------------------------
OptimizeBandwidth: Computing matrix structure for: mgdynamics...done.
OptimizeBandwidth: Half bandwidth without optimization: 2241616
OptimizeBandwidth: 
OptimizeBandwidth: Bandwidth Optimization ...done.
OptimizeBandwidth: Half bandwidth after optimization: 198555
OptimizeBandwidth: ---------------------------------------------------------
DefUtils::DefaultDirichletBCs: Setting Dirichlet boundary conditions
SetDirichletBoundaries: Number of dofs set: 0
DefUtils::DefaultDirichletBCs: Dirichlet boundary conditions set
SolveSystem: Solution trivially zero!
ERROR:: GetEdgeBasis: Can't handle but linear elements, sorry.

Which solver should I be putting the keyword into?

[raback]

Hi

I think you need also to set "Use Piola Transform = Logical True". This is quite novel and the keyword combinations are not harmonized. Quadratic elements always require this Piola transformed family of elements.

-Peter

[crobar]

Results without quadatic elements, but with Discontinuous Bodies = True. There's a slight difference but not much really.

B_and_A_sample_comparisons_noquad_discontinuous_bodies.png

(51.77 KiB) Not downloaded yet

[crobar]

Hi

I think you need also to set "Use Piola Transform = Logical True". This is quite novel and the keyword combinations are not harmonized. Quadratic elements always require this Piola transformed family of elements.

-Peter

Thanks, I will try

[crobar]

Results without quadratic but with Discontinuous Galerkin = True, v similar to Discontinuous Bodies = True

B_and_A_sample_comparisons_noquad_discontinuous_galerkin.png

(50.38 KiB) Not downloaded yet

[crobar]

Ok, so I tried with Piola Transform = Logical True, just put it in all solvers along with Quadratic Approximation = Logical True but Elmer works for a while then just stops with the final message "Kiiled". I've posted all the output below:

Code: Select all

ELMER SOLVER (v 8.2) STARTED AT: 2016/08/11 08:39:31
ParCommInit:  Initialize #PEs:            1
MAIN: 
MAIN: =============================================================
MAIN: ElmerSolver finite element software, Welcome!
MAIN: This program is free software licensed under (L)GPL
MAIN: Copyright 1st April 1995 - , CSC - IT Center for Science Ltd.
MAIN: Webpage http://www.csc.fi/elmer, Email elmeradm@csc.fi
MAIN: Version: 8.2, Compiled: 2016-08-09)
MAIN: =============================================================
MAIN: 
MAIN: 
MAIN: -------------------------------------
MAIN: Reading Model: case_quad.sif
Model Input:  Unlisted keyword: [a] in section: [boundary condition 1]
Model Input:  Unlisted keyword: [a {e}] in section: [boundary condition 1]
Loading user function library: [MagnetoDynamics]...[WhitneyAVSolver_Init0]
Loading user function library: [MagnetoDynamics]...[MagnetoDynamicsCalcFields_Init0]
Loading user function library: [ResultOutputSolve]...[ResultOutputSolver_Init0]
Loading user function library: [MagnetoDynamics]...[MagnetoDynamics_Dummy_Init0]
LoadMesh: Base mesh name: /home/rcrozier/build/fea_temp_files/elmer_vs_femm/elmer_vs_femm_quad
LoadMesh: Performing node mapping
ReadTargetNames: Reading names info from file: /home/rcrozier/build/fea_temp_files/elmer_vs_femm/elmer_vs_femm_quad/mesh.names
WARNING:: ReadTargetNames: Could not map name to Body nor BC: bnry7
LoadMesh: Minimum initial body index: 2
LoadMesh: Maximum initial body index: 6
LoadMesh: Minimum initial boundary index: 1
LoadMesh: Maximum initial boundary index: 7
LoadMesh: Elapsed time (CPU,REAL):    16.5280   25.0838 (s)
MeshStabParams: Computing stabilization parameters
MeshStabParams: Elapsed time (CPU,REAL):     0.7240    0.7235 (s)
MAIN: -------------------------------------
Loading user function library: [MagnetoDynamics]...[WhitneyAVSolver_Init]
Loading user function library: [MagnetoDynamics]...[WhitneyAVSolver_bulk]
Loading user function library: [MagnetoDynamics]...[WhitneyAVSolver]
Killed

And here is the full sif file:

Code: Select all

$filename  = "elmer_vs_femm_quad"  

! ======================================================================!
!	HEADER								!
! ======================================================================!

Header
   CHECK KEYWORDS "Warn"
   Mesh DB "/home/rcrozier/build/fea_temp_files/elmer_vs_femm" $filename
End

! ======================================================================!
!		Simulation						!
! ======================================================================!

Simulation
   Max Output Level = 7
   Coordinate System = "Cartesian 3D"
   Coordinate Mapping(3) = 1 2 3 
   Simulation Type = Steady state
   Steady State Max Iterations = 1
   Output Intervals = 1
   Use Mesh Names = Logical True
End

! ======================================================================!
!		Constants						!
! ======================================================================!

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.67e-08
  Permittivity of Vacuum = 8.8542e-12
  Boltzmann Constant = 1.3807e-23
  Unit Charge = 1.602e-19
End

! ======================================================================!
!		Body, Material,	Body Force, Initial Cond.		!
! ======================================================================!

Material 1
  Name = "Air"
  Density = 1000
  Electric Conductivity = 0.0
  Relative Permeability = 1.0
End

Material 2
  Name = "Magnet"
  Density = 1000
  Electric Conductivity = 0.0
  Relative Permeability = 1.05
End

Material 3
  Name = "Iron"
  Density = 7500
  Electric Conductivity = 0.0
  Relative Permeability = 1000.0
  
  H-B Curve = Variable "dummy"
   Real  Cubic Monotone
     INCLUDE HB_Iron
   End
End

Body 1
   Name = "magnet"
   Equation = 1
   Material = 2
   Body Force = 1
End

Body 2
   Name = "core"
   Equation = 1
   Material = 3
End

Body 3
   Name = "coil"
   Equation = 1
   Material = 1
End

Body 4
   Name = "inner_air_box"
   Equation = 1
   Material = 1
End

Body 5
   Name = "outer_air_box"
   Equation = 1
   Material = 1
End

Body Force 1
   Name = "Magnetization"
   Magnetization 1 = Real 0.0
   Magnetization 2 = Real 0.0
   Magnetization 3 = Real 979000
End

! ======================================================================!
!			Equations and Solvers				!
! ======================================================================!

Solver 1
   Equation = "MGDynamics"

   Variable = "A"

   Procedure = "MagnetoDynamics" "WhitneyAVSolver"
   Fix Input Current Density = Logical True

   Newton-Raphson Iteration = Logical True
   Nonlinear System Max Iterations = 30
   Nonlinear System Convergence Tolerance = 1e-6

   Linear System Symmetric = Logical True
   Linear System Solver = "Iterative"
   Linear System Preconditioning = None
   Linear System Convergence Tolerance = 1e-8
   Linear System Residual Output = 100
   Linear System Max Iterations = 5000
   Linear System Iterative Method = CG
   Steady State Convergence Tolerance = 1e-8

   Quadratic Approximation = Logical True
   Use Piola Transform = Logical True

End

Solver 2
   Equation = "MGDynamicsCalc"
   Procedure = "MagnetoDynamics" "MagnetoDynamicsCalcFields"
   Linear System Symmetric = True
   Potential Variable = String "A"

   Calculate Magnetic Field Strength = Logical True
   Calculate Magnetic Vector Potential = Logical True
   Steady State Convergence Tolerance = 0
   Linear System Solver = "Iterative"
   Linear System Preconditioning = None
   Linear System Residual Output = 0
   Linear System Max Iterations = 5000
   Linear System Iterative Method = CG
   Steady State Convergence Tolerance = 1e-8
   Linear System Convergence Tolerance = 1.0e-8

   Quadratic Approximation = Logical True
   Use Piola Transform = Logical True

End

Solver 3
  Exec Solver = after all
  Equation = "ResultOutput"
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Output File Name = "/home/rcrozier/build/fea_temp_files/elmer_vs_femm/"$filename
  Save Geometry Ids = Logical True
  Vector Field 1 = String Magnetic Field Strength
  Vector Field 2 = String Magnetic Flux Density
  Vector Field 3 = String Magnetic Vector Potential
  Potential Variable = String av
  Show Variables = Logical True
  Vtu format = Logical True

  Quadratic Approximation = Logical True
  Use Piola Transform = Logical True
End


Equation 1
  Name = "Coupled Equations"
  Active Solvers(2) = 1 2
End


! ======================================================================!
!		Boundary Condition Section				!
! ======================================================================!

! ----- names for boundaries -----
!$ outer_bound = 1

Boundary Condition 1
  Name = "outer_boundary"
  A = Real 0.0
  A {e} = Real 0.0
End