## How to output user defined variable to the vtu file?

Post processing utility for Elmer
Jenwel
Posts: 20
Joined: 22 Jun 2017, 07:07
Antispam: Yes

### Re: How to output user defined variable to the vtu file?

Hi Peter,

Thank you! The DG type of elemental field is necessary. And I still encounter the problem that I cannot define the similar variable Phys of Peddy for elemetal field.

I checked the source code "CalcFields .F90" and "FourierLoss.F90", I found I should define another dummy solver in an initialization subroutine? And maybe I should also define these variables in this dummy solver? As:

Code: Select all

! Create DG solver structures on-the-fly without actually solving the matrix equations.
ALLOCATE(Solvers(n+1))
Solvers(1:n) = Model % Solvers

Solvers(n+1) % Values => ListAllocate()
DGSolverParams => Solvers(n+1) % Values
CALL ListAddLogical( DGSolverParams, 'Discontinuous Galerkin', .TRUE. )
Solvers(n+1) % DG = .TRUE.
Solvers(n+1) % PROCEDURE = 0
Solvers(n+1) % ActiveElements => NULL()

CALL ListAddString( DGSolverParams, 'Exec Solver', 'never' )
CALL ListAddString( DGSolverParams, 'Equation', 'CoreLoss_Dummy' )
'CoreLoss CoreLossSolver_Dummy',.FALSE. )

CALL ListAddString( DGSolverParams, "Exported Variable 1 =", "Core Loss Phys e" )
CALL ListAddString( DGSolverParams, "Exported Variable 2 =", "Core Loss Peddy e")
CALL ListAddString( DGSolverParams, "Exported Variable 3 =", "CoreLoss_e" )

!PRINT *, "Initial Done", i
DEALLOCATE(Model % Solvers)
Model % Solvers => Solvers
Model % NumberOfSolvers = n+1

END SUBROUTINE CoreLossSolver_init0
An additional dummy solver is also defined:

Code: Select all

!------------------------------------------------------------------------------
SUBROUTINE CoreLossSolver_Dummy(Model,Solver,dt,Transient)
!------------------------------------------------------------------------------
USE MagnetoDynamicsUtils

IMPLICIT NONE
!------------------------------------------------------------------------------
TYPE(Solver_t) :: Solver
TYPE(Model_t) :: Model

REAL(KIND=dp) :: dt
LOGICAL :: Transient
!------------------------------------------------------------------------------
END SUBROUTINE CoreLossSolver_Dummy
!------------------------------------------------------------------------------
But when I want to link the pointer to the variable "Core Loss Phys e" in the actual subroutine CoreLossSolver, as:

Code: Select all

Ph_E => VariableGet( Mesh % Variables,'Core Loss Phys e' )
IF(.NOT. ASSOCIATED( Ph_E ))THEN
CALL Fatal('CoreLossSolver','Variable Core Loss Phys e does not exist')
END IF
The ERROR is reported:

Code: Select all

Model Input:  Unlisted keyword: [keddy] in section: [constants]
Model Input:  Unlisted keyword: [axiall] in section: [constants]
WARNING:: CheckTimer: Requesting time from non-existing timer: LoadMesh
WARNING:: CheckTimer: Requesting time from non-existing timer: MeshStabParams
MAIN:
MAIN: -------------------------------------
MAIN:  Time: 1/5   1.0000003999999999E-004
MAIN: -------------------------------------
MAIN:
WARNING:: CheckTimer: Requesting time from non-existing timer: PeriodicProjector
20 0.1011E+00
40 0.1733E-01
60 0.2544E-02
80 0.1113E-04
100 0.6558E-05
ComputeChange: NS (ITER=1) (NRM,RELC): (  2543.3687      1.0000000     ) :: mgdyn2d
20 0.7383E-03
40 0.2947E-03
60 0.2052E-03
80 0.1853E-04
100 0.1077E-06
ComputeChange: NS (ITER=2) (NRM,RELC): (  2556.8740     0.69178087E-01 ) :: mgdyn2d
20 0.6419E-05
40 0.5531E-05
60 0.4932E-05
80 0.9396E-07
ComputeChange: NS (ITER=3) (NRM,RELC): (  2550.2542     0.97899601E-02 ) :: mgdyn2d
ComputeChange: NS (ITER=4) (NRM,RELC): (  2550.8919     0.21344349E-03 ) :: mgdyn2d
ComputeChange: NS (ITER=5) (NRM,RELC): (  2550.8838     0.29937164E-05 ) :: mgdyn2d
ComputeChange: NS (ITER=6) (NRM,RELC): (  2550.8838     0.99764275E-07 ) :: mgdyn2d
WARNING:: GetPermittivity: Permittivity not defined in material, defaulting to that of vacuum
ComputeChange: NS (ITER=1) (NRM,RELC): ( 0.79328062E-01  2.0000000     ) :: calcfields
ComputeChange: NS (ITER=2) (NRM,RELC): ( 0.88136685E-01 0.10519973     ) :: calcfields
ComputeChange: NS (ITER=6) (NRM,RELC): ( 0.50781544E-03  2.0000000     ) :: calcfields
ComputeChange: NS (ITER=9) (NRM,RELC): (  545766.69      2.0000000     ) :: calcfields
Nsize =       12704
Nelem =       76602
CurrentTimeStep =           1
Ph % Values(1) =   0.0000000000000000
Nelem =       11362
ERROR:: CoreLossSolver: Variable Core Loss Phys e does not exist

It seems that the defined variable for elemental field is not linked, or defined.
The source code is uploaded. I do not know why did this problem occur?
Attachments
CoreLoss_Ver07.F90
Jenwel
Posts: 20
Joined: 22 Jun 2017, 07:07
Antispam: Yes

### Re: How to output user defined variable to the vtu file?

Hi,
Now the core loss of each node has been solved.
But, I do not know how to get the relationship between the element and its nodes, for the triangular type elements.
Is there any variable that denotes such mapping?
I can find some data in "mesh.elements", but I do not know how to invoke them.

If the relationship between the elements and nodes are determined, the average core loss of each element can be calculated.

BR
Jenwel
raback
Posts: 3744
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

### Re: How to output user defined variable to the vtu file?

Hi

Code: Select all

myave = MyVar % Values( MyVar % Perm(Element % NodeIndexes ) ) / Element % Type % NumberOfNodes

Note however, that when your field is discontinuous doing the galerkin fit to the nodes might not be the best thing. For the AV solver we usually look at elemental quantities that are computed first as elemental DG fields, and then averaged among elements in same body, if requested.

-Peter
Jenwel
Posts: 20
Joined: 22 Jun 2017, 07:07
Antispam: Yes

### Re: How to output user defined variable to the vtu file?

Hi, Peter

Thank you for the instruction. I agree that the fitting is not the best thing.
I found it difficult to programme a new solver, so I add the suboutine in the "CalcFields.F90". Now the elemtal DG field can be solved at first, as shown below:
CoreLoss_ElementField.jpg
CoreLoss_ElementalField
But, it is weird that the core loss also exists in the air gap.
I know it is caused by the principle to calculate the core loss. The core loss is calculated based on flux density, and actually the flux density also exists in the air gap region.
I do not know can I distinguish the material of the element and know it is the electrical steel or air? For the elements that belong to air, the core loss should be ZERO.

Thank you so much!

BR
Jenwel
raback
Posts: 3744
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

### Re: How to output user defined variable to the vtu file?

Hi

You should rather not use the B field computed to the nodes but use the vector potential at the integration points and evaluate the losses using them. This way your data is not smeared down.

-Peter
Jenwel
Posts: 20
Joined: 22 Jun 2017, 07:07
Antispam: Yes

### Re: How to output user defined variable to the vtu file?

Hi, Peter
Now I have calculated the elemental B field. I think maybe I can only calculate the core loss of the active elements.
In the solver "MagnetoDynamicsCalcFields", I use

Code: Select all

GetNOFActive()
to know the number of all active elements is 24518.
But the number of total elements is 76602.
Can I get the index of each active element in the total elements?

If this approach failed, I will try to use the A field at the integration points.
Thank you again!

BR
Jenwel
raback
Posts: 3744
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

### Re: How to output user defined variable to the vtu file?

Hi Jenwel,

I guess this is the physical model that you need implemented:

Code: Select all

		dBField(i) = BField(i) - PrevBField(i)
Ph % Values(i) = Hc * ABS(dBField(i)/dt) + Kh * ABS(BField(i)) * ABS(dBField(i)/dt)
Peddy % Values(i) = Keddy * ((dBField(i)/dt) ** 2)
Pcoreloss % Values(i) = Ph % Values(i) + Peddy % Values(i)

i.e. losses of style Hc*|B'|+Kh*|B||B'|+Keddy*|B'|^2 ?

It could well fit into CalcFields. It has already the possibility to compute losses for harmonic solution of style a*f^n*|B|+b*f^m*|B|^2. The FourierLoss solver can compute losses by making a fourier-transform on-the-fly and applying almost similar (more generic series) loss model. There is no loss model that directly uses the transient losses. I guess the background is rather similar since |B'| -> 2*\pi*f*|B| in harmonic case.

Would the loss model generalize to \sum_i c_i*|B|^n_i*|B'|^m_i?

-Peter
Jenwel
Posts: 20
Joined: 22 Jun 2017, 07:07
Antispam: Yes

### Re: How to output user defined variable to the vtu file?

Hi Peter,

That's the model I fitted into "CalcFields" now. Your codes basically agree with mine, except several variable names are different.
Now, the elemental loss can be solved. But the elements in air gap are also calculated by corresponding B field value, which is not expected.
Actually, the flux density exists in the air gap, but the air would not generate heat.
And I still do not know how to avoid calculating them.

So, I wonder if there is a variable that can denote the material or bodyID of the elements?
Before calculating the loss, the solver would know whether the element is air or iron by judging this variable.

Besides, I used the harmonic solution before. It can effectively solve the fourier loss. I have compared the results from the "FourierLoss solver" and the results from transient losses. I found that the results of "FourierLoss solver" is the mean (averaged) value of the transient losses in the integrate cycles.
The |B'| is equal to dB/dtheta*(2*pi*f) due to the fact that |B| varies with the angular position theta.

Honestly, I am pursuing the doctoral degree and the topic is the multi-objective optimization of the machine under the multi-physics fields.
So, at first, I want to solve the transient loss, then solve the thermal field of the machine. Maybe I will consider the air convection in the air gap and its impacts to the thermal field.
This is also the reason I begin to use Elmer. I believe it is the best choice to realize the FEA in multi-physics fields.
I also very appreciate you and other Elmer developers for creating such helpful Elmer.

Thank you so much!

BR
Jenwel
Jenwel
Posts: 20
Joined: 22 Jun 2017, 07:07
Antispam: Yes

### Re: How to output user defined variable to the vtu file?

Hi
I found the codes in "CalcFields"

Code: Select all

DO i = 1, GetNOFActive()
Element => GetActiveElement(i)
n = GetElementNOFNodes()
np = n*pSolver % Def_Dofs(GetElementFamily(Element),Element % BodyId,1)
nd = GetElementNOFDOFs(uSolver=pSolver)

CALL GetElementNodes( Nodes )
BodyId = GetBody()
Material => GetMaterial()
BodyForce => GetBodyForce()
I know the function

Code: Select all

Element => GetActiveElement(i)
can get the active element and link it with the pointer "Element ". Then the BodyID of this element can be somehow obtained via

Code: Select all

BodyId = GetBody()
But can I get the BodyID of all the elements? Including the non-active elements?

BR
Jenwel
Jenwel
Posts: 20
Joined: 22 Jun 2017, 07:07
Antispam: Yes

### Re: How to output user defined variable to the vtu file?

Hi !

OMG, after I read the data in the mesh document of my model, I realized that all the elements are divided into TWO main categories, the active elements and boundary elements.
In the file "mesh.elements" and "mesh.boundary", I found the number of active elements is 24518, the number of boundary elements is 1016.
And the size of elements is 25534, i.e. the sum of the number of active elements and boundary elements.

I found the code

Code: Select all

Element => Mesh % Elements(i)
can directly link the pointer Element.
And the code

Code: Select all

Element % BodyId
can obtain the BodyID of this element.

The number of flux density B for all elements, i.e. the variable "EL_MFD" in "CalcFields", is 229806.
It is equal to 3 multiplies the size of the defined variable Ph % Values, which is 76602.
I can understand that. It is because the B field "EL_MFD" has 3 Dofs in x- y- z-direction. And the defined variable Ph % Values is scalar value.

But I still do not understand. Why does the number of the elements is 25534, which is 1/3 of the size of the defined variable Ph % Values, i.e. 76602?
I am confused that, in my opinion, each element should have corresponding value of Ph % Values, so the size of the defined variable Ph % Values should also be equal to the number of elements?

Thank you!

BR
Jenwel