Improved visualization of discontinuous fields

Numerical methods and mathematical models of Elmer
Post Reply
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Improved visualization of discontinuous fields

Post by raback »

Hi All,

The visualization of the results obtained from MagnetoDynamicsCalcFields, and other discontinuous fields, have been somewhat suboptimal until recently. Particularly with the adoption of quadratic edge elements the old options were no longer satisfactory. This mail summarizes the recent changes related to the optimal visualization. I hope that with these new features the excellent 2nd order edge elements can show their full potential, but also the 1st order edge elements will benefit from these. The text is related to the current devel version.

The CalcFields solver can basically compute nodal quantities and elemental quantities. For fields known to be continuous the nodal quantities are good as such. However, for quantities having potential jumps between element interfaces (such as current density) the elemental fields are much more suitable. The elemental fields are actually standard nodal fields with the degrees of freedom being independent in each element (as in the Discontinuous Galerkin method). This makes the amount of data potentially quite huge. One may choose between three different ways to save the results in VTU output:

1) Constant (averaged Cell data) within each element (default)
2) "Discontinuous Galerkin" data where all nodes are different in each elements
3) "Discontinuous Bodies" data where nodes are made different only at the interfaces.

The preferred method of choice for visualization is usually 3). It also automatically averages the different shared instances of the DG fields to one single value. Unfortunately for quadratic elements there seems to be a problem. It seems that the rendering in Paraview cannot handle discontinuities well in conjunction with quadratic elements. To overcome that problem there is a possibility to save the result in linear basis even though having quadratic mesh. Actually this may make sense since the derived fields are still often just linear in accuracy.

Below is an example of a sif file that seems to provide nice visualizations for a number of cases both in serial and parallel

Code: Select all

Solver 3
  Exec Solver = after timestep
  Equation = "ResultOutput"
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Output File Name = f
  Vtu format = True
  Discontinuous Bodies = True
! Discontinuous Galerkin = True ! bloated alternative for the above
  Save Geometry Ids = True
  Save Linear Elements = True  ! use this only in conjunction with quadratic mesh !  Save Bulk Only = True ! Save Boundaries Only = True End
Unfortunately there is some loss of information when performing these
operations and line plotting in Paraview is not optimal. The alternative
is to save the data directly using SaveLine. Previously this routine had
problems with discontinuous fields also since it only saves data at the
intersection of faces, and assumes value in between to be linear. If the
face values are not unique this does not work. Now there is a remedy to
give number of divisions for each line. When using the "Polyline
Divisions" to specify the points also the DG fields are considered
correctly, for example

Code: Select all

Solver 4
  Exec Solver = after timestep
  Equation = SaveLine
  Procedure = "SaveData" "SaveLine"

  Filename = g.dat
  
  Polyline Coordinates(8,3) =  0   -0.05  -0.07  0    -0.05 0.115  \
                               0.0  0.05  -0.07  0.0   0.05 0.115 \
                              -0.05 0.0   -0.07 -0.05  0.0  0.115 \
                               0.05 0.0   -0.07  0.05  0.0  0.115 
  Polyline Divisions(4) = 100 100 100 100 
End
In parallel this creates independent files that can be easily joined
afterwords, e.g.:

Code: Select all

cat g*.dat > g.dat
Note that if one chooses "Discontinuous Bodies" in the VTU output it
averages the data. Hence if you want the original discontinuous (linear or
quadratic) data give this solver a smaller Id number than for the
ResultOutputSolver.


Finally there is a possibility to also use SaveGridData solver would you
want to save the data in a uniform grid. This also now works with DG
fields and can therefore capture the true results from CalcFields. This is
probably not that useful, but an example is provided here for
completeness.

Code: Select all

Solver 5
  Equation = SaveGridData
  Procedure = "SaveGridData" "SaveGridData"

  Grid dx = Real 0.05
  Grid dy = Real 0.05
  Grid dz = Real 0.00005

! Reduce the box where the points are saved in 
  Min Coordinate 1 = Real -0.07
  Max Coordinate 1 = Real 0.07
  Min Coordinate 2 = Real -0.07
  Max Coordinate 2 = Real 0.07

! If nodes are exactly at the interface eliminate redundant ones
  Check For Duplicates = Logical True
  Binary Output = Logical True
  Single Precision = Logical True

  Filename Prefix = h

! save in vtu format in paraview (use glyphs)
  Vtu Format = Logical True

! save in ascii format for Matlab etc. 
  Table Format = Logical True
End
I hope people using the great MagnetodDynamics solver will find these instructions useful.

-Peter
spanda
Posts: 11
Joined: 19 Feb 2018, 14:24
Antispam: Yes

Re: Improved visualization of discontinuous fields

Post by spanda »

Hi Peter,

I am trying to simulate test case circuits_harmonic_massive (1962.sif in your github repository) with quadratic approximation and results drastically differ (in Paraview) from the expected results that I have been gotten after simulation with first order edge elements.

Here is the sif file I 'm using. It is just a little udjusted Takala's sif file from github (Quadratic Approximation = Logical True
Use Piola Transform = Logical True).

Code: Select all

! A simple testcase for testing CircuitsAndDynamics module.
! Author: Eelis Takala, Trafotek Oy
! Original date: October 2015
! eelis.takala@trafotek.fi

Check Keywords "Warn"
INCLUDE sif/1962-circuits.definitions
INCLUDE 1962/mesh.names

Header 1
   Mesh DB "." "1962"
End

Constants 1
   Permittivity of Vacuum = 8.8542e-12
End

Simulation 1
   Max Output Level = 3
   Coordinate System = Cartesian 3D
   Coordinate Mapping(3) = 1 2 3
   Coordinate Scaling = 0.001
   Simulation Type = Steady
   Steady State Max Iterations = 1
   BDF Order = 1
   Output Intervals = 1
   Post File = rezultat_2th.vtu
   Angular Frequency = 314.159265359
!   output file = "runoutput.dat"
End

Solver 1
   Exec Solver = Before All
   Procedure = "WPotentialSolver" "Wsolve"
   Equation = "Wire direction"
   Variable = W
   Linear System Solver = Iterative
   Linear System Iterative Method = CG
   Linear System Max Iterations = 10000
   Linear System Convergence Tolerance = 1.0e-10
   Linear System Preconditioning = ILU0
   Linear System Abort Not Converged = True
   Linear System Residual Output = 1000
End

Solver 2
   Exec Solver = Always
   Equation = Circuits
   Variable = X
   No Matrix = Logical True
   Procedure = "CircuitsAndDynamics" "CircuitsAndDynamicsHarmonic"
End

Solver 3
   Equation = "MGDynamics"
   Variable = "A[A re:1 A im:1]"
   Procedure = "MagnetoDynamics" "WhitneyAVHarmonicSolver"
   Angular Frequency = 314.159265359
   Export Lagrange Multiplier = Logical True
   Linear System Complex = Logical True
   Quadratic Approximation = Logical True
   Use Piola Transform = Logical True
   Linear System Solver = Iterative
   Linear System Iterative Method = BicgstabL
   Linear System preconditioning = Circuit
   Linear System Convergence Tolerance = 1.e-7
   Linear System Max Iterations = 3000
   Linear System Residual Output = 1000
   BicgStabL Polynomial Degree = 4
   Linear System Abort not Converged = Logical False
   Steady State Convergence Tolerance = 1e-06
End

Solver 4
   Equation = "MGDynamicsCalc"
   Procedure = "MagnetoDynamics" "MagnetoDynamicsCalcFields"
   Linear System Symmetric = True
   Potential Variable = String "A"
   Calculate Current Density = Logical True
   Loss Estimation = Logical True
   Calculate Joule Heating = Logical True
   Steady State Convergence Tolerance = 0
   Linear System Solver = "Iterative"
   Linear System Preconditioning = None
   Linear System Residual Output = 1000
   Linear System Max Iterations = 5000
   Linear System Iterative Method = CG
   Steady State Convergence Tolerance = 1e-6
   Linear System Convergence Tolerance = 1.0e-8
End

Solver 5
   Exec Solver = after timestep
   Equation = "ResultOutput"
   Procedure = "ResultOutputSolve" "ResultOutputSolver"
   Output File Name = 2ndNedelec
   Vtu format = Logical True
   Discontinuous Bodies = True 
   Save Geometry Ids = Logical True
   Save Linear Elements = True
End

Solver 6
   Exec Solver = After timestep
   Equation = Circuits Output
   Procedure = "CircuitsAndDynamics" "CircuitsOutput"
End

Solver 7
   Exec Solver = Never
   Equation = "sv"
   Procedure = "SaveData" "SaveScalars"
   Filename = 1962/dat/1962.dat
End

Equation 1
   Active Solvers(2) = 3 4
End

Equation 2
   Active Solvers(4) = 1 2 3 4
End

Material 1
   Name = perm1e1
   Relative Permittivity = Real 1
   Electric Conductivity = Real 1000000
   Relative Permeability = Real 1
End

Material 2
   Name = air
   Relative Permittivity = Real 1
   Electric Conductivity = Real 0
   Relative Permeability = Real 1
End

Material 3
   Name = aluminium
   Relative Permittivity = Real 1
   Relative Permeability = Real 1
   Electric Conductivity = 3.07e7
End

Body 1
   Name = core
   Target Bodies(1) = $ core
   Equation = 1
   Material = 1
   Initial Condition = 1
End

Body 2
   Name = air
   Target Bodies(1) = $ air
   Equation = 1
   Material = 2
   Initial Condition = 1
End

Body 3
   Name = L1
   Target Bodies(1) = $ L1
   Equation = 2
   Material = 3
   Initial Condition = 1
   Body Force = 1
End

Component 1
   Name = String L1
   Master Bodies = Integer 3
   Coil Type = String massive
   Foil Winding Voltage Polynomial Order = Integer 1
End

Body Force 1
   Name = "Circuit"
   testsource Re = Real $1.41421356237*cos(0.0*pi/3)
   testsource Im = Real $1.41421356237*sin(0.0*pi/3)
End

Boundary Condition 1
   Name = BCn Flux Parallel
   Target Boundaries(2) = $ coreface_xy xy0
   A re {e} = Real 0
   A im {e} = Real 0
End

Boundary Condition 2
   Name = ground
   Target Boundaries = $ L1_gamma1
   W = Real 0
   A re {e} = Real 0
   A im {e} = Real 0
End

Boundary Condition 3
   Name = current in foil winding
   Target Boundaries = $ L1_gamma0
   W = Real 1
   A re {e} = Real 0
   A im {e} = Real 0
End

Solver 3 :: Reference Norm = Real 3.80171554E-05
Solver 3 :: Reference Norm Tolerance = Real 1e-03 
$fprintf( stderr, "TEST CASE 1\n");
RUN

Is this the right way for setting options for second order edge element simulation and visualization.

Here is Takala's original sif file from github:

Code: Select all

! A simple testcase for testing CircuitsAndDynamics module.
! Author: Eelis Takala, Trafotek Oy
! Original date: October 2015
! eelis.takala@trafotek.fi

Check Keywords "Warn"
INCLUDE sif/1962-circuits.definitions
INCLUDE 1962/mesh.names

Header 1
   Mesh DB "." "1962"
End

Constants 1
   Permittivity of Vacuum = 8.8542e-12
End

Simulation 1
   Max Output Level = 3
   Coordinate System = Cartesian 3D
   Coordinate Mapping(3) = 1 2 3
   Coordinate Scaling = 0.001
   Simulation Type = Steady
   Steady State Max Iterations = 1
   BDF Order = 1
   Output Intervals = 1
   Angular Frequency = 314.159265359
!   output file = "runoutput.dat"
End

Solver 1
   Exec Solver = Before All
   Procedure = "WPotentialSolver" "Wsolve"
   Equation = "Wire direction"
   Variable = W
   Linear System Solver = Iterative
   Linear System Iterative Method = CG
   Linear System Max Iterations = 10000
   Linear System Convergence Tolerance = 1.0e-10
   Linear System Preconditioning = ILU0
   Linear System Abort Not Converged = True
   Linear System Residual Output = 1000
End

Solver 2
   Exec Solver = Always
   Equation = Circuits
   Variable = X
   No Matrix = Logical True
   Procedure = "CircuitsAndDynamics" "CircuitsAndDynamicsHarmonic"
End

Solver 3
   Equation = "MGDynamics"
   Variable = "A[A re:1 A im:1]"
   Procedure = "MagnetoDynamics" "WhitneyAVHarmonicSolver"
   Angular Frequency = 314.159265359
   Export Lagrange Multiplier = Logical True
   Linear System Symmetric = Logical True
   Linear System Complex = Logical True
   Linear System Solver = Iterative
   Linear System Iterative Method = BicgstabL
   Linear System preconditioning = Circuit
   Linear System Convergence Tolerance = 1.e-7
   Linear System Max Iterations = 3000
   Linear System Residual Output = 1000
   BicgStabL Polynomial Degree = 4
   Linear System Abort not Converged = True
   Steady State Convergence Tolerance = 1e-06
End

Solver 4
   Equation = "MGDynamicsCalc"
   Procedure = "MagnetoDynamics" "MagnetoDynamicsCalcFields"
   Linear System Symmetric = True
   Potential Variable = String "A"
   Calculate Current Density = Logical True
   Loss Estimation = Logical True
   Steady State Convergence Tolerance = 0
   Linear System Solver = "Iterative"
   Linear System Preconditioning = None
   Linear System Residual Output = 1000
   Linear System Max Iterations = 5000
   Linear System Iterative Method = CG
   Steady State Convergence Tolerance = 1e-6
   Linear System Convergence Tolerance = 1.0e-8
End

Solver 5
   Exec Solver = Never
   Equation = "ResultOutput"
   Procedure = "ResultOutputSolve" "ResultOutputSolver"
   Output File Name = 1962-results
   Vtu format = Logical True
   Save Geometry Ids = Logical True
End

Solver 6
   Exec Solver = After timestep
   Equation = Circuits Output
   Procedure = "CircuitsAndDynamics" "CircuitsOutput"
End

Solver 7
   Exec Solver = Never
   Equation = "sv"
   Procedure = "SaveData" "SaveScalars"
   Filename = 1962/dat/1962.dat
End

Equation 1
   Active Solvers(2) = 3 4
End

Equation 2
   Active Solvers(4) = 1 2 3 4
End

Material 1
   Name = perm1e1
   Relative Permittivity = Real 1
   Electric Conductivity = Real 0
   Relative Permeability = Real 1e1
End

Material 2
   Name = air
   Relative Permittivity = Real 1
   Electric Conductivity = Real 0
   Relative Permeability = Real 1
End

Material 3
   Name = aluminium
   Relative Permittivity = Real 1
   Relative Permeability = Real 1
   Electric Conductivity = 3.07e7
End

Body 1
   Name = core
   Target Bodies(1) = $ core
   Equation = 1
   Material = 1
   Initial Condition = 1
End

Body 2
   Name = air
   Target Bodies(1) = $ air
   Equation = 1
   Material = 2
   Initial Condition = 1
End

Body 3
   Name = L1
   Target Bodies(1) = $ L1
   Equation = 2
   Material = 3
   Initial Condition = 1
   Body Force = 1
End

Component 1
   Name = String L1
   Master Bodies = Integer 3
   Coil Type = String massive
   Foil Winding Voltage Polynomial Order = Integer 1
End

Body Force 1
   Name = "Circuit"
   testsource Re = Real $1414.21356237*cos(0.0*pi/3)
   testsource Im = Real $1414.21356237*sin(0.0*pi/3)
End

Boundary Condition 1
   Name = BCn Flux Parallel
   Target Boundaries(2) = $ coreface_xy xy0
   A re {e} = Real 0
   A im {e} = Real 0
End

Boundary Condition 2
   Name = ground
   Target Boundaries = $ L1_gamma1
   W = Real 0
   A re {e} = Real 0
   A im {e} = Real 0
End

Boundary Condition 3
   Name = current in foil winding
   Target Boundaries = $ L1_gamma0
   W = Real 1
   A re {e} = Real 0
   A im {e} = Real 0
End

Solver 3 :: Reference Norm = Real 3.80171554E-05
Solver 3 :: Reference Norm Tolerance = Real 1e-03 
$fprintf( stderr, "TEST CASE 1\n");
RUN
I hope it is OK question :)

Thank you
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Improved visualization of discontinuous fields

Post by raback »

Hi spanda,

How fresh is your Elmer version?

Until a month or so the CalcFields was expecting the quadratic/piola settings to be consistent. Currently they should follow the primary solver.

There could be also some other reason. I don't know how extensively the circuit simulator has been tested in conjuntion with the quadratic Hcurl elements.

-Peter
spanda
Posts: 11
Joined: 19 Feb 2018, 14:24
Antispam: Yes

Re: Improved visualization of discontinuous fields

Post by spanda »

Hi Peter,
Thanks for response :)

I installed binary one year ago. So, it is offical version from May, 2017.

I have tried to simulate another test case without circuit simulator and it worked fine.
I thought maybe the problem is in Paraview visualization of simulation that uses second-order edge elements. Now, I don't know if I am missing something in the SIF file or something else or circuit simulator is really the problem. I will try more things.

Sorry for my bad English

Thanks Peter

-Spanda
Post Reply