ListGetElementRealVec yields output of zeroes

[spacedout]

This function is to be found in file Lists.F90.

Its 2nd parameter, ngp, is the number of Gaussian integration points. ngp is a function of the element type. It will thus vary as I loop over all active elements in the mesh. Now when I call this function, setting this 2nd parameter to a constant equal to say the maximum of ngp over all element types then there is no problem and the result will be whatever the handle points to in the SIF. This is the desired behaviour. On the other hand, if I set it to the ngp of the current element being scanned, (the normal way to use the function obviously), then initially everything is OK until the ngp suddenly increases for a certain element , then the result will be all 0's totalling in number to ngp . This will happen again for all elements with that same greater value of ngp.

As this was never reported before, it may date back to the earliest version of ListGetElementRealVec.

It you need more details, I will gladly provide the simplest case that displays this issue.

[spacedout]

I managed to simplify the SIF and code tremendously and still get the erroneous output.

I ran on:

OS - RedHat or CentOS
Compiler - intel/2020.1.217

with the following commands (everything is attached in weirdbug.zip)

elmerf90 -o SoSimpleTest.so SoSimpleTest.F90
ElmerSolver timecoef.sif > printout.txt 2>&1

you get:

ELMER SOLVER (v 9.0) STARTED AT: 2023/11/20 17:31:57
ParCommInit: Initialize #PEs: 1
MAIN: OMP_NUM_THREADS not set. Using only 1 thread per task.
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: 58f9634a, Compiled: 2021-03-16)
MAIN: Running one task without MPI parallelization.
MAIN: Running with just one thread per task.
MAIN: =============================================================
LoadInputFile: Trying to read "Run Control" section only


.....................................................

: CUR ELEM 1099
: NUM NODES 3
: NUM DOFS 3
: NUM BDOFS 0
: TimeCOEF 1.0000000000000000
: TimeCOEF 1.0000000000000000
: TimeCOEF 1.0000000000000000
: CUR ELEM 1100
: NUM NODES 4
: NUM DOFS 4
: NUM BDOFS 0
: TimeCOEF 0.0000000000000000
: TimeCOEF 0.0000000000000000
: TimeCOEF 0.0000000000000000
: TimeCOEF 0.0000000000000000


The above is wrong and occurs with the following in SoSimpleTest.F90
TimeCoeff => ListGetElementRealVec( TimeCoeff_h, ngp, Basis, Element, Found )
! TimeCoeff => ListGetElementRealVec( TimeCoeff_h, 4, Basis, Element, Found )

but with
! TimeCoeff => ListGetElementRealVec( TimeCoeff_h, ngp, Basis, Element, Found )
TimeCoeff => ListGetElementRealVec( TimeCoeff_h, 4, Basis, Element, Found )

you get:


: CUR ELEM 1099
: NUM NODES 3
: NUM DOFS 3
: NUM BDOFS 0
: TimeCOEF 1.0000000000000000
: TimeCOEF 1.0000000000000000
: TimeCOEF 1.0000000000000000
: CUR ELEM 1100
: NUM NODES 4
: NUM DOFS 4
: NUM BDOFS 0
: TimeCOEF 1.0000000000000000
: TimeCOEF 1.0000000000000000
: TimeCOEF 1.0000000000000000
: TimeCOEF 1.0000000000000000

which is the desired result

[spacedout]

I have replaced

TimeCoeff => ListGetElementRealVec( TimeCoeff_h, 4, Basis, Element, Found )

by

IF (t >Active/2) THEN
TimeCoeff => ListGetElementRealVec( TimeCoeff_h, 5, Basis, Element, Found )
ELSE
TimeCoeff => ListGetElementRealVec( TimeCoeff_h, 4, Basis, Element, Found )
ENDIF


Now TimeCoeff points to a list of ones until we reach the first element number greater than Active/2. It then points to a list of zeroes for every element number greater than Active/2. I could be wrong but it seems to indicate there is a memory problem related to a DEALLOCATE followed by an ALLOCATE.

[raback]

Hi

The 2nd parameter is the number of gauss points and probably it is not a good idea to tamper it. This is the vectorized fetching of material parameters and it is assumed that the Basis vector is compatible with that. You could print out the "ngp" in your 1st test.

The default integration rule for a triangle is 3 nodes, and for quad its 4 nodes. I think there exists also a 5-point rule but you have to ask for this rule.

-Peter

[spacedout]

Hello

Well, fine I went back to the first test I uploaded to this forum and I printed ngp as well.

the minimum (maximum) of both ngp and nd is 3 (4)

I then commented out the DEALLOCATE code and altered the allocate code following it to

IF (.NOT. ALLOCATED(Basis)) THEN
ALLOCATE(Basis(4,4), dBasisdx(4,4,3), DetJ(4), STAT=allocstat)

IF (allocstat /= 0) THEN
CALL Fatal(Caller,'Local storage allocation failed')
END IF
END IF

but that did not change the weird results I already reported.

Ditto when I removed all code related to allocation and instead declare

REAL(KIND=dp) :: Basis(4,4), dBasisdx(4,4,3), DetJ(4)

So, I say the bug is within function ListGetElementRealVec and until it is resolved, I have no choice but to go back to the original modelPDE code which only uses non "Vec" routines like ListGetElementReal for example.

[raback]

Hi, maybe you can provide a full test case where this happens. -Peter

[spacedout]

But surely, no one can expect me to be able to debug the Elmer code itself as I never was involved in its design !

[raback]

Hi

I had now some time to look at it closer. Indeed, there was a bug! That should be fixed now. Thanx for reporting!

The bug manifested itself when there was different elementtypes within the same material such that the one with smaller number of gaussian integration points was first. If also the value of the keyword was constant in that section it was never initialized again after the number of integration points was increased. If the keyword was function of other fields then everything was ok, and also if the first element happened to have the max. number of integration points.

-Peter

[spacedout]

How do I integrate this update into my years old Elmer installation ? Or do I have to install the newest Elmer ? Anyhow I still don't see that update under Elmer GitHub. Thanks

[raback]

Hi,

You need to recompile / fetch new nightly build. This is in the library.

This is the commit:
https://github.com/ElmerCSC/elmerfem/co ... 3c42128635

-Peter