[Solved!] Mgdyn simulation problem (is it boundary conditions?)

[Nick_99]

Hi all,

Some background
I'm doing a masters degree and have been playing around with Elmer, for the past few weeks, so that I can get acquainted with it for my research. I'm particularly interested in the mgdyn and mgdyn2D solvers. I followed the tutorial to simulate the 2D horseshoe magnet and was able to do that successfully. I then replicated that tutorial from scratch (so modeled the horseshoe etc) and was also able to do that perfectly and got results that agree with the tutorial.

Edit: I'm using ElmerGUI and I'm running everything on Ubuntu 22.04.2 LTS.

My issue
Two weeks ago, I started trying to extend the 2D horseshoe example into a 3D simulation. I've managed to make a model in gmsh but have run into nonsensical results with Elmer. Recently, I have simplified the 3D horseshoe down to just two bar magnets but am still running into similar issues. When I load up the results of the sim into ParaView, I can see that there is zero flux on any of the bodies and the magnetic field strength (which appears to be wrong as well) is only non-zero on the two bar magnets (the air between them is also zero).

My boundary condition is set to

Code: Select all

AV {e} = Real 0

and is applied to the region surrounding the bar magnets (i.e. the air around them). Is this perhaps the issue? Is this BC causing the flux and field strength to be zero outside of the magnets?

I also see the solver outputting 'SolveLinearSystem: Solution trivially zero!' when I run the sim. Does anyone know what this means? I did find this post but the solution provided by raback didn't seem to help me

I've reduced the mesh size quite a bit because I was fiddling with the BCs and didn't want to wait 5 minutes per simulation with a high resolution mesh. I'm not too concerned with accuracy at the moment, I just want to see the field flowing between the magnets

My SIF file

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 5
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Steady state
  Steady State Max Iterations = 1
  Output Intervals(1) = 1
  Coordinate Scaling = 1e-3
  Solver Input File = case.sif
  Post File = case.vtu
End

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.670374419e-08
  Permittivity of Vacuum = 8.85418781e-12
  Permeability of Vacuum = 1.25663706e-6
  Boltzmann Constant = 1.380649e-23
  Unit Charge = 1.6021766e-19
End

Body 1
  Target Bodies(1) = 1
  Name = "Body 1"
  Equation = 1
  Material = 3
End

Body 2
  Target Bodies(1) = 2
  Name = "Body 2"
  Equation = 1
  Material = 2
End

Body 3
  Target Bodies(1) = 3
  Name = "Body 3"
  Equation = 1
  Material = 1
End

Solver 1
  Equation = MgDynPost
  Procedure = "MagnetoDynamics" "MagnetoDynamicsCalcFields"
  Calculate Magnetic Field Strength = True
  Exec Solver = Always
  Stabilize = True
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-10
  BiCGstabl polynomial degree = 2
  Linear System Preconditioning = ILU0
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 10
  Linear System Precondition Recompute = 1
End

Solver 2
  Equation = MgDyn
  Procedure = "MagnetoDynamics" "WhitneyAVSolver"
  Exec Solver = Always
  Stabilize = True
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-10
  BiCGstabl polynomial degree = 2
  Linear System Preconditioning = ILU0
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 10
  Linear System Precondition Recompute = 1
End

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

Material 1
  Name = "Air"
  Relative Permeability = 1
End

Material 2
  Name = "Mag_S"
  Magnetization 1 = Real -780e3
  Relative Permeability = 5000
End

Material 3
  Name = "Mag N"
  Relative Permeability = 5000
  Magnetization 1 = Real 780e3
End

Boundary Condition 1
  Target Boundaries(6) = 1 2 3 4 5 6 
  Name = "BoundaryCondition 1"
  AV {e} = Real 0
End

Some images

Flux density

Screenshot 2023-07-12 091955.png

(42.9 KiB) Not downloaded yet

Magnetic field strength

Screenshot 2023-07-12 092111.png (117.51 KiB) Viewed 1121 times

As an aside, in ParaView, when I'm selecting a quantity to visualize (say, magnetic field strength), some quantities will be accompanied by another quantity of the same name but ending in an 'e' (for example, magnetic field strength e). What does this mean? Is it the error in the computation?

Thanks,
Nick

[kevinarden]

Solution trivially zero means there was no solution, typically because there is no right hand side to the equation defined.

The problem in your sif file is the post processor is going first and then the whitney solver
Since post goes first it reports no solution because there is not a solution to post process
Switch the order of the solvers and it works

e after the variable means that it is elemental data not nodal data, many solver save nodal data by default but the magdyn will typically save nodal and elemental data

[Nick_99]

Oh my word I can't believe it was that simple. I've literally been pulling my hair out for days now. Thank you so so much!!!!

[Nick_99]

Hi Kevin (or anyone else ),

I've run into a smaller issue and was wondering if you could give me some advice.

After you helped me, I went and modified my example, attached to my first post, so that the magnets were more rectangular (rather than cuboidal) in shape. I then proceeded to fiddle with the mesh size. Apart from that, the only other thing I've changed from the original case.sif is the maximum number of iterations for linear system (increased to 500 because the simulation was ending before reaching the convergence criterion).

If I keep the mesh size low, then this issue does not appear. However, with a low mesh size, my results are useless because the model has a low accuracy.

When I look at the results in ParaView, I see this oddity:

issue.png

(293.9 KiB) Not downloaded yet

I'm not sure why the results looks like this

My case.sif is:

case.sif

(2.92 KiB) Downloaded 45 times

My geometry file is:

Code: Select all

SetFactory("OpenCASCADE");

// lm corresponds to the mesh size on and inside the magnets' geometry.
// lr specifies the mesh size on and inside the region's geometry.
lm = 1;
lr = 5;

Point(1) = {-50, 7, -5, lm};
Point(2) = {-50, -7, -5, lm};
Point(3) = {-10, -7, -5, lm};
Point(4) = {-10, 7, -5, lm};
Point(5) = {-50, 7, 5, lm};
Point(6) = {-50, -7, 5, lm};
Point(7) = {-10, -7, 5, lm};
Point(8) = {-10, 7, 5, lm};

Line(1) = {5, 8};
Line(2) = {5, 6};
Line(3) = {6, 2};
Line(4) = {2, 1};
Line(5) = {1, 5};
Line(6) = {1, 4};
Line(7) = {4, 8};
Line(8) = {3, 4};
Line(9) = {3, 2};
Line(10) = {3, 7};
Line(11) = {7, 8};
Line(12) = {7, 6};

Curve Loop(1) = {6, -8, 9, 4};
Plane Surface(1) = {1};

Curve Loop(2) = {6, 7, -1, -5};
Plane Surface(2) = {2};

Curve Loop(3) = {7, -11, -10, 8};
Plane Surface(3) = {3};

Curve Loop(4) = {4, 5, 2, 3};
Plane Surface(4) = {4};

Curve Loop(5) = {9, -3, -12, -10};
Plane Surface(5) = {5};

Curve Loop(6) = {1, -11, 12, -2};
Plane Surface(6) = {6};

Surface Loop(1) = {6, 2, 1, 3, 5, 4};
Volume(1) = {1};

Point(9) = {10, 7, 5, lm};
Point(10) = {10, -7, 5, lm};
Point(11) = {10, -7, -5, lm};
Point(12) = {10, 7, -5, lm};
Point(13) = {50, 7, -5, lm};
Point(14) = {50, 7, 5, lm};
Point(15) = {50, -7, 5, lm};
Point(16) = {50, -7, -5, lm};

Line(13) = {11, 16};
Line(14) = {16, 13};
Line(15) = {13, 12};
Line(16) = {12, 11};
Line(17) = {11, 10};
Line(18) = {12, 9};
Line(19) = {13, 14};
Line(20) = {16, 15};
Line(21) = {15, 14};
Line(22) = {10, 9};
Line(23) = {9, 14};
Line(24) = {15, 10};

Curve Loop(7) = {13, 14, 15, 16};
Plane Surface(7) = {7};
Curve Loop(8) = {15, 18, 23, -19};
Plane Surface(8) = {8};
Curve Loop(9) = {19, -21, -20, 14};
Plane Surface(9) = {9};
Curve Loop(10) = {17, -24, -20, -13};
Plane Surface(10) = {10};
Curve Loop(11) = {16, 17, 22, -18};
Plane Surface(11) = {11};
Curve Loop(12) = {22, 23, -21, 24};
Plane Surface(12) = {12};

Surface Loop(2) = {8, 7, 10, 11, 12, 9};
Volume(2) = {2};

Point(17) = {-80, -50, -40, lr};
Point(18) = {-80, 50, -40, lr};
Point(19) = {-80, 50, 40, lr};
Point(20) = {-80, -50, 40, lr};
Point(21) = {80, -50, -40, lr};
Point(22) = {80, -50, 40, lr};
Point(23) = {80, 50, 40, lr};
Point(24) = {80, 50, -40, lr};

Line(25) = {18, 24};
Line(26) = {24, 23};
Line(27) = {23, 19};
Line(28) = {19, 18};
Line(29) = {18, 17};
Line(30) = {17, 21};
Line(31) = {21, 24};
Line(32) = {21, 22};
Line(33) = {22, 23};
Line(34) = {19, 20};
Line(35) = {20, 17};
Line(36) = {20, 22};

Curve Loop(13) = {29, 30, 31, -25};
Plane Surface(13) = {13};
Curve Loop(14) = {25, 26, 27, 28};
Plane Surface(14) = {14};
Curve Loop(15) = {29, -35, -34, 28};
Plane Surface(15) = {15};
Curve Loop(16) = {31, 26, -33, -32};
Plane Surface(16) = {16};
Curve Loop(17) = {30, 32, -36, 35};
Plane Surface(17) = {17};
Curve Loop(18) = {34, 36, 33, 27};
Plane Surface(18) = {18};

Surface Loop(3) = {18, 15, 13, 17, 16, 14};

Surface Loop(4) = {6, 2, 1, 3, 5, 4};
Surface Loop(5) = {12, 11, 7, 10, 9, 8};

Volume(3) = {3, 4, 5};

Here is also the log for the solver. This was run with MPI using 8 threads.

latest.log

(40.48 KiB) Downloaded 45 times

EDIT: So really my question boils down to me asking if this is an issue with how I've setup Elmer or is it an issue with ParaView?

Kind regards,
Nick

[kevinarden]

I think the issue is that it is not a actually a converged solution. Not sure why it is not converging.

[Nick_99]

Hmm ok thank you!! I'll look into it more. I have a hunch it may be my boundary condition.

[Rich_B]

Hello,

Take a look at the test, elmerfem\fem\tests\mgdyn_steady_2ndorder. The test uses different solver settings. See the attached modified sif file, which runs to completion with your geometry, properly converged.

mag-flux.png (88.29 KiB) Viewed 1021 times

Contents of the mesh.header file, 40,517 nodes, 241,964 elements:

40517 241964 17266
2
303 17266
504 241964

Rich.

[Nick_99]

Hi Rich,

Thank you so so much!! You're a lifesaver
I will have a look at the test you linked to.

I appreciate your (likewise for Kevin) time and effort to assist me.

Kind regards,
Nick

[kevinarden]

I took you magnet geometry and instead of wrapping it in a rectangle of air I surrounded them in a sphere of air. It then converged very rapidly using the default solver settings. I think the solver struggles with the infinity boundary being a rectangle.

[Nick_99]

Ok I will give that a try. I used a rectangle because I was struggling with the sphere's mesh intersecting with the geometries it encompasses.

Thanks a ton