Solver stops after Gebhard factors

[Lhoeltke]

Total Elmer and FEM newbie here. I want to simulate the thermal equilibrium between a small heated crystal and a heated metallic sphere through radiation, The model represents a small crystal in a ultra-high vacuum chamber. But I got a problem which I could not resolve.
The solver just stops after calculating the gebhard factors in this case as well as for the crystal alone (radiating into the vacuum), but shows no error messages, except for some Matrix access errors way earlier, but these were often present in simulations that ran through.

Code: Select all

ELMER SOLVER (v 9.0) STARTED AT: 2024/02/27 20:56:43

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: 9.0 (Rev: Release, Compiled: 2024-02-22)
MAIN:  Running one task without MPI parallelization.
MAIN:  Running with just one thread per task.
MAIN:  Lua interpreter linked in.
MAIN: =============================================================

MAIN: 
MAIN: 
MAIN: -------------------------------------
MAIN: Reading Model: case.sif

LoadInputFile: Scanning input file: case.sif
LoadInputFile: Scanning only size info
LoadInputFile: First time visiting
LoadInputFile: Reading base load of sif file

LoadInputFile: Loading input file: case.sif
LoadInputFile: Reading base load of sif file

LoadInputFile: Number of BCs: 4
LoadInputFile: Number of Body Forces: 0
LoadInputFile: Number of Initial Conditions: 1
LoadInputFile: Number of Materials: 2
LoadInputFile: Number of Equations: 1
LoadInputFile: Number of Solvers: 1
LoadInputFile: Number of Bodies: 2

ListTagKeywords: Setting weight for keywords!
ListTagKeywords: No parameters width suffix: normalize by area
ListTagKeywords: Setting weight for keywords!
ListTagKeywords: No parameters width suffix: normalize by volume
ElmerAsciiMesh: Base mesh name: ./.

MapCoordinates: Scaling coordinates: 1.000E+00 1.000E+00 1.000E+00

LoadMesh: Elapsed REAL time:     0.1270 (s)

MAIN: -------------------------------------
AddVtuOutputSolverHack: Adding ResultOutputSolver to write VTU output in file: case

RadiationFactors: ----------------------------------------------------
RadiationFactors: Computing radiation factors for heat transfer
RadiationFactors: ----------------------------------------------------

RadiationFactors: Total number of Radiation Surfaces 4480 out of 9728

RadiationFactors: Temporarily updating the mesh.nodes file!

ViewFactors: 
ViewFactors: ==================================================
ViewFactors:  E L M E R  V I E W F A C T O R S,  W E L C O M E
ViewFactors: ==================================================

ViewFactors: 
ViewFactors: 
ViewFactors: Reading Model...
ViewFactors: Computing view factors as defined in file: case.sif

LoadInputFile: Scanning input file: case.sif
LoadInputFile: Scanning only size info
LoadInputFile: First time visiting
LoadInputFile: Reading base load of sif file

LoadInputFile: Loading input file: case.sif
LoadInputFile: Reading base load of sif file

LoadInputFile: Number of BCs: 4
LoadInputFile: Number of Body Forces: 0
LoadInputFile: Number of Initial Conditions: 1
LoadInputFile: Number of Materials: 2
LoadInputFile: Number of Equations: 1
LoadInputFile: Number of Solvers: 1
LoadInputFile: Number of Bodies: 2

ListTagKeywords: Setting weight for keywords!
ListTagKeywords: No parameters width suffix: normalize by area
ListTagKeywords: Setting weight for keywords!
ListTagKeywords: No parameters width suffix: normalize by volume

ElmerAsciiMesh: Base mesh name: ./.

MapCoordinates: Scaling coordinates: 1.000E+00 1.000E+00 1.000E+00

LoadMesh: Elapsed REAL time:     0.1690 (s)

ViewFactors: Computing view factors for radiation body1

ViewFactors: Number of surfaces participating in radiation: 4480

ViewFactors: Computing viewfactors...

surfs: 4480, min(353)=0.00, max(29)=1.00, avg=0.99
ViewFactors: View factors computed in time (s):   10.40   10.40

ViewFactors: 
ViewFactors: Viewfactors before manipulation:
ViewFactors: Minimum row sum:    0.000E+00
ViewFactors: Maximum row sum:    9.999E-01
ViewFactors: Symmetrizing Factors...

ViewFactors: View factors manipulated in time (s):    0.12
ViewFactors: 
ViewFactors: Viewfactors after manipulation:
ViewFactors: Minimum row sum:    0.000E+00
ViewFactors: Maximum row sum:    9.999E-01

ViewFactors: Saving view factors in ascii mode

ViewFactors: *** ALL DONE ***

RadiationFactors: Loading view factors from ascii file: ././ViewFactors.dat

RadiationFactors: Computing factors...

ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      932044         353
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      940739         365
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      940739         365
ERROR:: MakeMatrixIndex:  Trying to access non-existent column:      966848         373

RadiationFactors: View factors filling (%)               9.2932E+01

RadiationFactors:    Solution:   0 % done

RadiationFactors:    Solution:   7 % done

RadiationFactors:    Solution:  16 % done

RadiationFactors:    Solution:  24 % done

RadiationFactors:    Solution:  33 % done

RadiationFactors:    Solution:  42 % done

RadiationFactors:    Solution:  50 % done

RadiationFactors:    Solution:  60 % done

RadiationFactors:    Solution:  69 % done

RadiationFactors:    Solution:  79 % done

RadiationFactors:    Solution:  89 % done

RadiationFactors:    Solution:  98 % done

RadiationFactors: Minimum Gebhart factors sum          0.000000E+00
RadiationFactors: Maximum Gebhart factors sum          8.910557E-04
RadiationFactors: Maximum share of omitted factors              NaN
RadiationFactors: Gebhart factors filling (%)          1.992985E-05
RadiationFactors: Gebhart factors count                4

Can you help me resolve this problem, so that it runs the thermal simulation?
When running the radiation tutorial file as well as the radiation equilibrium between two crystals everything works just fine, as long as i dont define a surface radiating onto nothing.
I also ask myself if Elmer even can simulate heat loss through radiation without another body receiving the radiation, because i couldnt really find anything in the manuals on it.
The project folder is attached in a .zip file, but In have ommited the view factors file, because it was giant.
Thank you for your help in advance.

[kevinarden]

I ran it in command window and increased the message level to get more out put. These are the errors, the program is halting because of an invalid memory reference.

out.txt

(2.67 KiB) Downloaded 11 times

[raback]

Hi

This is a strange one. The computation of view factors is one of the older pieces of Elmer and still there seems to be some strange bug. This piece of code has been used alot!

I used quite a bit of time today to find the bug but I could not really find the main culprint. What happens is that ~25 of the boundary elements in the inner circle have a normal pointing to the wrong direction (out of >4000). Because of this they do not get any hits and the sum of view factors becomes zero.

There are basically tests to ensure the correct direction. The way used in this code is to identify one node of the parent that is not on the boundary element and use the coordinates of this node and triangle center to check the direction. I can see that this could fail if the tetrahedron is broken so that the 4th node also lies on the same plane.

When the direction of these normals is altered everthing goes fine.

I added some verbosity to the code that is responsible for this. This was done in a bit of hurry but all 18 tests related to radiation still pass...
https://github.com/ElmerCSC/elmerfem/co ... fc353cfd68

Maybe you can study the mesh for possible issues.

Just run "ViewFactors" to see the issues.

-Peter

[raback]

This piece of debug code in ViewFactors.F90 let's this case pass. Just checks the normal direction when R>20.

Code: Select all

       ! This fixes the case by "Lhoeltke" on forum. Not a generic one!
       ! In line ~530 of ViewFactors. 
#if 1       
       BLOCK
         REAL(KIND=dp) :: r0(3), n0(3), rad,dotpr
         j = 0
         k = 0
         
         DO i=1,n         
           Element => RadElements(i)
           r0(1) = SUM(model % nodes % x(element % nodeindexes)) / element % Type % NumberofNodes
           r0(2) = SUM(model % nodes % y(element % nodeindexes)) / element % Type % NumberofNodes
           r0(3) = SUM(model % nodes % z(element % nodeindexes)) / element % Type % NumberofNodes

           rad = SQRT(SUM(r0**2))
           IF(rad < 20.0) CYCLE
           
           n0 = normals(3*(i-1)+1:3*(i-1)+3)           
           dotpr = SUM(n0*r0)

           ! We should not have normals pointing outwards.
           IF( dotpr > 0.0_dp ) THEN
             j =j+1
             PRINT *,'bc elem:',i,r0, n0 
             normals(3*(i-1)+1:3*(i-1)+3) = -n0          
           ELSE
             k = k+1
           END IF           
         END DO
         PRINT *,'Wrong normal directions:',j
       END BLOCK
#endif

[kevinarden]

I tried a simpler 2D case same error

RadiationSphere.zip

(36.6 KiB) Downloaded 10 times

[kevinarden]

An axis symmetric case works

2daxis.zip

(97.32 KiB) Downloaded 9 times

[Lhoeltke]

Thanks for all the replys, they were really helpfull. I can run both of Kevinardens cases without problems. The Problem in my case was indeed the normals pointing into the wrong direction, after checking a new mesh for the normal directions, everything worked fine.

Two smaller problems remain however:
- The solver sometimes exit with "Error occurred in umf4num: -1.0000000", another thread (viewtopic.php?t=3110) suggests this is because of insufficient memory, so I will just throw more RAM at it.
- I sadly dont know what Peter means with

Just run "ViewFactors" to see the issues.

Where do i run this to check if my meshes are correct? This would be very usefull as my meshes will only get more complex.

If the view factor sum cannot be Zero, then simulations of meshes/surfaces just radiating outward into/onto nothing are not possible, is that correct? This seems odd to me because why should the calculation not proceed?

[raback]

Hi

The "ViewFactors" binary is called by "ElmerSolver" automatically. However, if the problem is the view factors then just compute those. If the view factors exist and the nodes are the same they are not recomputed unless requested.

If you decleare your boundaries open then the rowsum can be less than one. However, if the system is closed it makes sense to use this info since it will normalize the rowsum to one thereby reducing the errors. Your case was fully closed. The errors were found earlier then.

What was broken in your mesh?

-Peter

[Lhoeltke]

I had to get home to my better PC first to check the mesh using gmesh, and it seems there are no normals pointing in the wrong direction as I expected. So I am back to square one on that issue . The old simulation also still stops at the same point with the new code you implemented.

In the new meshing procedure in gmsh I only use 3D meshing once and before that I only split the 2D meshes, which leads to meshes that (mostly) work. But only with the new code you implemented.

Also the following error often comes up in the simulations with the new meshes as soon as they get too big:

Code: Select all

ELMER SOLVER (v 9.0) STARTED AT: 2024/03/03 15:07:20

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: 9.0 (Rev: Release, Compiled: 2024-02-22)
MAIN:  Running one task without MPI parallelization.

MAIN:  Running with just one thread per task.
MAIN:  Lua interpreter linked in.
MAIN: =============================================================

MAIN: 
MAIN: 
MAIN: -------------------------------------
MAIN: Reading Model: case.sif

LoadInputFile: Scanning input file: case.sif
LoadInputFile: Scanning only size info
LoadInputFile: First time visiting
LoadInputFile: Reading base load of sif file

LoadInputFile: Loading input file: case.sif
LoadInputFile: Reading base load of sif file

LoadInputFile: Number of BCs: 4
LoadInputFile: Number of Body Forces: 0
LoadInputFile: Number of Initial Conditions: 1
LoadInputFile: Number of Materials: 2
LoadInputFile: Number of Equations: 1
LoadInputFile: Number of Solvers: 1
LoadInputFile: Number of Bodies: 2

ListTagKeywords: Setting weight for keywords!
ListTagKeywords: No parameters width suffix: normalize by area
ListTagKeywords: Setting weight for keywords!
ListTagKeywords: No parameters width suffix: normalize by volume
Loading user function library: [HeatSolve]...[HeatSolver_Init0]
LoadMesh: Loading serial mesh!
ElmerAsciiMesh: Performing step: 1
ElmerAsciiMesh: Base mesh name: ./.
ReadHeaderFile: Reading header info from file: ././mesh.header
InitializeMesh: Number of nodes in mesh: 4877
InitializeMesh: Number of bulk elements in mesh: 14899
InitializeMesh: Number of boundary elements in mesh: 9728
InitializeMesh: Initial number of max element nodes: 4
ElmerAsciiMesh: Performing step: 2

ReadNodesFile: Reading nodes from file: ././mesh.nodes

MapCoordinates: Performing coordinate mapping
SetMeshDimension: Dimension of mesh is: 3
SetMeshDimension: Max dimension of mesh is: 3

MapCoordinates: Scaling coordinates: 1.000E+00 1.000E+00 1.000E+00
ElmerAsciiMesh: Performing step: 3

ReadElementsFile: Reading bulk elements from file: ././mesh.elements

ElmerAsciiMesh: Performing step: 4

ReadBoundaryFile: Reading boundary elements from file: ././mesh.boundary

PermuteNodeNumbering: Performing node mapping

MapBodiesAndBCs: Remapping bodies
MapBodiesAndBCs: Minimum initial body index: 1
MapBodiesAndBCs: Maximum initial body index: 2
MapBodiesAndBCs: Remapping boundaries
MapBodiesAndBCs: Minimum initial boundary index: 1
MapBodiesAndBCs: Maximum initial boundary index: 8

ElmerAsciiMesh: Performing step: 5
ElmerAsciiMesh: Performing step: 6
LoadMesh: Loading mesh done
LoadMesh: Preparing mesh done
LoadMesh: Elapsed REAL time:     0.2010 (s)
MeshStabParams: Computing stabilization parameters

MeshStabParams: Elapsed REAL time:     0.2010 (s)
LoadModel: Defined radition solver by Equation name "heat equation"
MAIN: -------------------------------------
AddVtuOutputSolverHack: Adding ResultOutputSolver to write VTU output in file: case
AppendNewSolver: Increasing number of solvers to: 2
AddMeshCoordinates: Setting mesh coordinates and time
AddSolvers: Setting up 2 solvers
AddSolvers: Setting up solver 1: heat equation
AddEquationBasics: Setting up keywords internally for legacy solver: heat equation
AddEquationBasics: Using procedure: HeatSolve HeatSolver
AddEquationBasics: Setting up solver: heat equation

Loading user function library: [HeatSolve]...[HeatSolver_Init]
Loading user function library: [HeatSolve]...[HeatSolver_bulk]
Loading user function library: [HeatSolve]...[HeatSolver]
AddEquationBasics: Creating standard variable: temperature
RadiationFactors: ----------------------------------------------------
RadiationFactors: Computing radiation factors for heat transfer
RadiationFactors: ----------------------------------------------------
RadiationFactors: Using sparse matrix format for factor computations.
RadiationFactors: Using direct solver for radiation factors

RadiationFactors: Total number of Radiation Surfaces 4480 out of 9728

RadiationFactors: Temporarily updating the mesh.nodes file!

ViewFactors: 
ViewFactors: ==================================================
ViewFactors:  E L M E R  V I E W F A C T O R S,  W E L C O M E
ViewFactors: ==================================================

ViewFactors: 
ViewFactors: 
ViewFactors: Reading Model...
ViewFactors: Computing view factors as defined in file: case.sif

LoadInputFile: Scanning input file: case.sif
LoadInputFile: Scanning only size info
LoadInputFile: First time visiting
LoadInputFile: Reading base load of sif file

LoadInputFile: Loading input file: case.sif
LoadInputFile: Reading base load of sif file

LoadInputFile: Number of BCs: 4
LoadInputFile: Number of Body Forces: 0
LoadInputFile: Number of Initial Conditions: 1
LoadInputFile: Number of Materials: 2
LoadInputFile: Number of Equations: 1
LoadInputFile: Number of Solvers: 1
LoadInputFile: Number of Bodies: 2

ListTagKeywords: Setting weight for keywords!
ListTagKeywords: No parameters width suffix: normalize by area
ListTagKeywords: Setting weight for keywords!
ListTagKeywords: No parameters width suffix: normalize by volume
Loading user function library: [HeatSolve]...[HeatSolver_Init0]
LoadMesh: Loading serial mesh!
ElmerAsciiMesh: Performing step: 1
ElmerAsciiMesh: Base mesh name: ./.
ReadHeaderFile: Reading header info from file: ././mesh.header
InitializeMesh: Number of nodes in mesh: 4877
InitializeMesh: Number of bulk elements in mesh: 14899
InitializeMesh: Number of boundary elements in mesh: 9728
InitializeMesh: Initial number of max element nodes: 4

ElmerAsciiMesh: Performing step: 2

ReadNodesFile: Reading nodes from file: ././mesh.nodes

MapCoordinates: Performing coordinate mapping
SetMeshDimension: Dimension of mesh is: 3
SetMeshDimension: Max dimension of mesh is: 3
MapCoordinates: Scaling coordinates: 1.000E+00 1.000E+00 1.000E+00
ElmerAsciiMesh: Performing step: 3

ReadElementsFile: Reading bulk elements from file: ././mesh.elements

ElmerAsciiMesh: Performing step: 4

ReadBoundaryFile: Reading boundary elements from file: ././mesh.boundary

PermuteNodeNumbering: Performing node mapping

MapBodiesAndBCs: Remapping bodies
MapBodiesAndBCs: Minimum initial body index: 1
MapBodiesAndBCs: Maximum initial body index: 2
MapBodiesAndBCs: Remapping boundaries

MapBodiesAndBCs: Minimum initial boundary index: 1
MapBodiesAndBCs: Maximum initial boundary index: 8
ElmerAsciiMesh: Performing step: 5
ElmerAsciiMesh: Performing step: 6
LoadMesh: Loading mesh done

LoadMesh: Preparing mesh done
LoadMesh: Elapsed REAL time:     0.0580 (s)
MeshStabParams: Computing stabilization parameters

MeshStabParams: Elapsed REAL time:     0.0100 (s)
LoadModel: Defined radition solver by Equation name "heat equation"
ViewFactors: Number of nodes in mesh: 4877

ViewFactors: Computing view factors for radiation body1

ViewFactors: Number of surfaces participating in radiation: 4480

ViewFactors: Computing viewfactors...

surfs: 4480, min(3430)=1.00, max(3762)=1.00, avg=1.00
ViewFactors: View factors computed in time (s):   10.68   10.68

ViewFactors: 
ViewFactors: Viewfactors before manipulation:
ViewFactors: Minimum row sum:    9.981E-01
ViewFactors: Maximum row sum:    1.001E+00
ViewFactors: Symmetrizing Factors...

ViewFactors: View factors manipulated in time (s):    0.11
ViewFactors: 
ViewFactors: Viewfactors after manipulation:
ViewFactors: Minimum row sum:    9.981E-01
ViewFactors: Maximum row sum:    1.001E+00
WARNING:: ViewFactors: Rowsum of view factors should not be larger than one!
ViewFactors: Printing some lumped information of view factors

ViewFactors: BC 1 with 4224 elems perm: 1
ViewFactors: BC 3 with 256 elems perm: 2
ViewFactors: Lumped areas:
ViewFactors:    125474.77915110193        128.00000000000000

ViewFactors: Lumped View Factor Matrix:
ViewFactors:   0.99920811610221894        1.0198139438674849E-003
ViewFactors:   0.99969476079676289        0.0000000000000000

ViewFactors: Saving view factors in ascii mode

ViewFactors: *** ALL DONE ***

RadiationFactors: Loading view factors!

RadiationFactors: Loading view factors from ascii file: ././ViewFactors.dat

RadiationFactors: Computing factors...

CreateRadiationMatrix: Number of entries in matrix: 18892474

RadiationFactors: View factors filling (%)               9.4131E+01

TabulateEmissivity: We used real emissivity for all elements!

DirectSolver: Using direct method: umfpack

 Error occurred in umf4num:   -1.0000000000000000     

STOP 1

even though this machine has more memory than my Laptop (8GB vss. 16 GB). It also only allocates 3-4 GB before throwing the error. The ttterative solver always throws convergence errors.

[kevinarden]

I can control the element normal in salome, do not know how to in gmsh.
This fixed the 3d issue with view factors, however I also get


DirectSolver: Using direct method: umfpack

Error occurred in umf4num: -1.0000000000000000

STOP 1

with 32GB Ram