search

ElmerCSC / elmerfem

Official git repository of Elmer FEM software
http://www.elmerfem.org
Other
1.21k stars 319 forks source link

SEGFAULT from attempted adaptive remesh with discontinuity #597

Open NRJank opened 1 month ago

NRJank commented 1 month ago

Full original bug report in https://www.elmerfem.org/forum/viewtopic.php?t=6608 described a number of attempts to combine adaptive remesh with different discontinuity methods, which was triggering either odd results or segfaults in certain cases. Some changes were made to the code to catch the errors, halt operation, etc. Rechecking status of this on recent nightly build (ELMERFEM-nogui-nompi-Windows-AMD64.zip, 2024-10-23) it seems that Case 4 (combining heatgap with adaptive refinement) still produces a segfault,

adapt2.grd:

***** ElmerGrid input file for structured grid generation *****
Version = 210903

Coordinate System = Cartesian 2D
Subcell Divisions in 2D = 3 3 

Subcell Sizes 1 = 1e-3 1e-3 1e-3
Subcell Sizes 2 = 1e-3 1e-3 1e-3

Material Structure in 2D
  2 1 2
  2 2 2
  3 3 3
End

Materials Interval = 1 3

Boundary Definitions
! type     out      int      double   of the boundaries
  1        -1       3       1
  2        -3       1       1
  3        -1       2       2
End

Surface Elements = 10
Triangles = logical true  !needed for RGB adaptive meshing
Numbering = Horizontal
Element Degree = 1

adapt2.sif:

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "adapt2"
End

Simulation
  Max Output Level = 5
  Coordinate System = "Cartesian 2D"
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = "Steady State"
  Steady State Max Iterations = 10
  Post File = file "adapt.ep"
End

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

Body 2
  Target Bodies(1) = 3
  Equation = 1
  Material = 1
  Initial Condition = 1
End

Equation 1
  Heat Equation = True
  Active Solvers(1) = 1
End

Solver 1
  Exec Solver = "Always"
  Equation = "Heat Equation"
  Procedure = "HeatSolve" "HeatSolver"
  Variable = "Temperature"
  Variable Dofs = 1
  Optimize Bandwidth = True
  Linear System Solver = "Direct"
  Linear System Direct Method = "umfpack"
  Steady State Convergence Tolerance = 1.0e-5
  Stabilize = True
  Nonlinear System Convergence Tolerance = 1.0e-4
  Nonlinear System Max Iterations = 3
  Nonlinear System Newton After Iterations = 2
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1.0

  Adaptive Mesh Refinement = true
  Adaptive Remesh = false
  Adaptive Save Mesh = true
  Adaptive Error Limit = 1e-4
End

Material 1
   Density = 1
   Heat Conductivity = 1
   Heat Capacity = 1
End

Boundary Condition 1
  Target Boundaries(1) = 1
  Name = COLD
  Temperature = 0
End

Boundary Condition 2
  Target Boundaries(1) = 2
  Name = HOT
  Heat Flux BC = True
  Heat Flux = 100
End

Boundary Condition 3
  Target Boundaries(1) = 3
  Name = heatgap
  Heat Gap = True
  Heat Transfer Coefficient = 1000
End

segfault triggered by running:

elmergrid 1 2 adapt2.grd
elmersolver adapt2.sif

output:

ELMER SOLVER (v 9.0) STARTED AT: 2024/10/23 13:57:49
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-10-23)
MAIN:  Running one task without MPI parallelization.
MAIN:  Running with just one thread per task.
MAIN:  Lua interpreter linked in.
MAIN:  External optimization routines linked in.
MAIN: =============================================================
MAIN:
MAIN:
MAIN: -------------------------------------
MAIN: Reading Model: adapt2.sif
LoadInputFile: Scanning input file: adapt2.sif
LoadInputFile: Scanning only size info
LoadInputFile: First time visiting
LoadInputFile: Reading base load of sif file
LoadInputFile: Loading input file: adapt2.sif
LoadInputFile: Reading base load of sif file
LoadInputFile: Number of BCs: 3
LoadInputFile: Number of Body Forces: 0
LoadInputFile: Number of Initial Conditions: 0
LoadInputFile: Number of Materials: 1
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: ./adapt2
LoadMesh: Elapsed REAL time:     0.0020 (s)
MAIN: -------------------------------------
OptimizeBandwidth: Initial bandwidth for heat equation: 8
OptimizeBandwidth: Optimized bandwidth for heat equation: 5
MAIN: Number of timesteps to be saved: 10
MAIN:
MAIN: -------------------------------------
MAIN:  Steady state iteration:            1
MAIN: -------------------------------------
MAIN:
HeatSolver: Solving the energy equation for temperature
WARNING:: HeatSolver: The old way of dealing with HeatGap is obsolite!!
HeatSolve:
HeatSolve:
HeatSolve: -------------------------------------
HeatSolve:  TEMPERATURE ITERATION           1
HeatSolve: -------------------------------------
HeatSolve:
HeatSolve: Starting Assembly...
HeatSolve: Assembly:
WARNING:: CRS_AddToMatrixElement: Matrix element is to be added to a nonexistent position
WARNING:: CRS_AddToMatrixElement: Row: 11 Col: 16
WARNING:: CRS_AddToMatrixElement: Number of Matrix rows:24
WARNING:: CRS_AddToMatrixElement: Converting CRS to list
HeatSolve: Assembly done
ComputeChange: NS (ITER=1) (NRM,RELC): ( 0.73690407E-01  2.0000000     ) :: heat equation
HeatSolve: iter:    1 Assembly: (s)    0.00    0.00
HeatSolve: iter:    1 Solve:    (s)    0.00    0.00
HeatSolve:  Result Norm   :    7.3690406854581281E-002
HeatSolve:  Relative Change :    2.0000000000000000
HeatSolve:
HeatSolve:
HeatSolve: -------------------------------------
HeatSolve:  TEMPERATURE ITERATION           2
HeatSolve: -------------------------------------
HeatSolve:
HeatSolve: Starting Assembly...
HeatSolve: Assembly:
HeatSolve: Assembly done
ComputeChange: NS (ITER=2) (NRM,RELC): ( 0.73690407E-01  0.0000000     ) :: heat equation
HeatSolve: iter:    2 Assembly: (s)    0.00    0.01
HeatSolve: iter:    2 Solve:    (s)    0.01    0.01
HeatSolve:  Result Norm   :    7.3690406854581281E-002
HeatSolve:  Relative Change :    0.0000000000000000
RefineMesh:
RefineMesh: ----------- M E S H   R E F I N E M E N T --------------
RefineMesh: Initializing stuff on coarsest level!
RefineMesh: The new mesh consists of:
RefineMesh: Nodal points: 103
RefineMesh: Bulk elements: 163
RefineMesh: Boundary elements: 11
BuildQuandrantTree: Start
BuildQuandrantTree: Ready
OptimizeBandwidth: Initial bandwidth for heat equation: 82
OptimizeBandwidth: Optimized bandwidth for heat equation: 19
RefineMesh: ----------- E N D   M E S H   R E F I N E M E N T --------------
ComputeChange: SS (ITER=1) (NRM,RELC): ( 0.77382617E-01  2.0000000     ) :: heat equation
WritePostFile: Saving results in ElmerPost format to file adapt2/adapt.ep
MAIN:
MAIN: -------------------------------------
MAIN:  Steady state iteration:            2
MAIN: -------------------------------------
MAIN:
HeatSolver: Solving the energy equation for temperature
WARNING:: HeatSolver: The old way of dealing with HeatGap is obsolite!!
HeatSolve:
HeatSolve:
HeatSolve: -------------------------------------
HeatSolve:  TEMPERATURE ITERATION           1
HeatSolve: -------------------------------------
HeatSolve:
HeatSolve: Starting Assembly...
HeatSolve: Assembly:
WARNING:: CRS_AddToMatrixElement: Matrix element is to be added to a nonexistent position
WARNING:: CRS_AddToMatrixElement: Row: 55 Col: 20
WARNING:: CRS_AddToMatrixElement: Number of Matrix rows:103
WARNING:: CRS_AddToMatrixElement: Converting CRS to list

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0xebe1c3eb
#1  0xebd97034
#2  0x41672231
#3  0x69737ff7
#4  0x698928be
#5  0x69842553
#6  0x698913cd
#7  0x92eeebb2
#8  0x2ed0adb8
#9  0x2ed13566
#10  0x92cfa6a5
#11  0x92d13a43
#12  0x92d15b57
#13  0x92d17195
#14  0x9302a163
#15  0x9303380d
#16  0x416716ef
#17  0x41672ac3
#18  0x416713c0
#19  0x416714f5
#20  0x68bf7373
#21  0x6983cc90
#22  0xffffffff

case6 and case 7 catch similar attempts at combining adaptive refinement with newer discontinuity methods and issue: ERROR:: RefineMesh: Adaptive refinement not possible for discontinuous mesh!

case 8, however, seems to avoid the error being caught, runs to completion, but produces erroneous output where the adaptive solve refines around the discontinuity but the heatsolve ignores the heat gap as if it was 'leaky'. files for case 8:

adapt3.grd

***** ElmerGrid input file for structured grid generation *****
Version = 210903

Coordinate System = Cartesian 2D
Subcell Divisions in 2D = 3 3 

Subcell Sizes 1 = 1e-3 1e-3 1e-3
Subcell Sizes 2 = 1e-3 1e-3 1e-3

Material Structure in 2D
  2 1 2
  2 2 2
  3 3 3
End

Materials Interval = 1 3

Boundary Definitions
! type     out      int      double   of the boundaries
  1        -1       3       1
  2        -3       1       1
  3        -1       2       1
End

Surface Elements = 10
Triangles = logical true  !needed for RGB adaptive meshing
Numbering = Horizontal
Element Degree = 1

adapt3.sif

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "adapt3"
End

Simulation
  Max Output Level = 5
  Coordinate System = "Cartesian 2D"
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = "Steady State"
  Steady State Max Iterations = 10
  Post File = file "adapt.ep"
  Discontinuous Bulk Greedy = Logical false
End

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

Body 2
  Target Bodies(1) = 3
  Equation = 1
  Material = 1
  Initial Condition = 1
End

Equation 1
  Heat Equation = True
  Active Solvers(1) = 1
End

Solver 1
  Exec Solver = "Always"
  Equation = "Heat Equation"
  Procedure = "HeatSolve" "HeatSolver"
  Variable = "Temperature"
  Variable Dofs = 1
  Optimize Bandwidth = True
  Linear System Solver = "Direct"
  Linear System Direct Method = "umfpack"
  Steady State Convergence Tolerance = 1.0e-5
  Stabilize = True
  Nonlinear System Convergence Tolerance = 1.0e-4
  Nonlinear System Max Iterations = 3
  Nonlinear System Newton After Iterations = 2
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1.0

  Adaptive Mesh Refinement = true
  Adaptive Remesh = false
  Adaptive Save Mesh = true
  Adaptive Error Limit = 1e-4
End

Material 1
   Density = 1
   Heat Conductivity = 1
   Heat Capacity = 1
End

Boundary Condition 1
  Target Boundaries(1) = 1
  Name = COLD
  Temperature = 0
End

Boundary Condition 2
  Target Boundaries(1) = 2
  Name = HOT
  Heat Flux BC = True
  Heat Flux = 100
End

Boundary Condition 3
  Target Boundaries(1) = 3
  Name = heatgap
  Discontinuous Boundary  = Logical True
  Discontinuous Target Bodies (2) =  2 3
  Heat Gap = True
  Heat Transfer Coefficient = 1000
End

Note the only difference between Case 7 and Case 8 is changing Discontinuous Bulk Greedy in the Simulation block from true to false.

Would recommend in the short term both the segfault and the non-error case be caught and stopped with the same error message used in Case 6/7.

Would request in the long term being able to combine adaptive refinement and discontinuous boundaries.

full set of test cases from forum report attached for reference.

elmertest.zip

raback commented 1 month ago

There are some issues related to AMR & discontinuous meshes that hinder this. 1) Currently the error indicator is a nodal quantity and probably assumes that the primary field is nodal as well. There has been extension to p-elements within year or so, but not to DG-elements. This is probably a minor indexing thing. In parallel it could be somewhat more complicated since the parallel info should be different for nodal and DG fields. 2) How to consider DG jumps for the error indicator. Probably some alterations are needed. Not really my specialty...

NRJank commented 4 weeks ago

DG? (dual grid?)

I assumed this was a non-trivial item. Was bringing it up again mainly because i noticed the unhandled errors.

raback commented 4 weeks ago

DG = discontinous galerkin

(we have the standard + special version of DG where the discontinuity only happens between materials)