FEAP User Forum
FEAP => Programming => Topic started by: kiakarim on August 30, 2021, 07:17:57 AM
-
Hello everyone,
is it possible that dt of tdata.h is wrongly updated? I have to type in the dt,,... command two times in order to receive the correct timestep in my user element.
feap 8.6.1j
-
That should not be the case. The code for setting dt is in program/pmacr2.f and it is exceedingly simple:
case (2)
dtold = dt
dt = ct(1,l)
-
Hello, I'm sorry to bother you again.
The previous question for the timestep was due to me having problems transferring an element from ver8.5.2i to ver8.6.1j, which uses an implicit time discretization.
Attached you find a dummy element, where I added the line p(1,1) = 1d0/dt, which should simulate a time discretizing term.
As well you can find a simple input file which calls the dummy element.
The problem now is, that in 8.5 the computation continues if i issue the TIME command, whereas version 8.6 immediately breaks.
Also if add the command "write(*,*) 1/dt " ver85 outputs infinity, while ver86 breaks.
Is this a known problem to you?
-
Strange, my 8.6 reveals the correct value. However you should check dt before dividing by it in case you ever want to do dt = 0.0
Maybe you need to recompile the program?
-
Actually I did a plain reinstall of the patch you uploaded on ubuntu 20.04.02 LTS with a few modifications:
- $FEAPHOME8_6 : mv makefile.in_gfortran makefile.in
- $FEAPHOME8_6/main : mv makefile_default makefile
- $FEAPHOME8_6/main/makefile : uncomment line with dsymutil
- cd $FEAPHOME8_6; make archive
Auxiliar packages ::
gfortran-10 // gcc-10 : 10.3.0
liblapack-dev// 3.9.0-1build1
libx11-dev // 2:1.6.9-2ubuntu1.2
Additional note:
The time step is set in the Inputfile Ielmt and is nonzero.
Do you receive "infinity" as output from the user element, when write command is uncommented?
2image.jpg is the output of version 8.5
-
Can you uncomment the save line in your element.
Also place write(*,*) 'Elmt01: ',dt before the if statement in your element
and place write(*,*) 'Pmacr2: ',dt on (about) line 209 in program/pmacr2.f,
just after the dt=ct(1,l).
run make archive in the program folder, and then recompile and link in your element, and post
what gets printed out.
-
The code changes and the program output are attached.
-
Sorry one other thing.
under isw.eq.6 .or. isw.eq.3
comment out the setting of p(1,1) and uncomment the print statement and change it to print dt instead of p(1,1).
-
The changes and the new output are attached.
-
I see the problem. This is actually a bug in version 8.6 -- though it would be pretty rare for this to affect anyone; I guess you are an exception :)
Here is a quick fix, that should resolve the issue:
Edit program/pmacr2.f. At (about) line 598 insert
if(nstep.gt.0) thenAt (about) line 623 insert
endif
-
Okay, well I thought this problem should be quite common, considering the fact, that transient heat transfer requires the division by dt?!?
Nevertheless, the problem seems to be resolved, although I cannot get why. Attached you can find the code change in pmacr2.f.
Thank You!
-
Good the problem is resolved for you.
I suggest you check the value of dt before you divide (maybe define a variable for reciprocal time to avoid all the divides - since divide is the most expensive of basic arithmetic operations, +/-/*)
Also I suggest you read the programmer manual about using the feap time integration algorithms and the use of "ctan(*)" instead of the 1/dt.
-
I just finished implementing your suggestion and now it works like charms and without any compability problems between 8.5 and 8.6.
I'm further interested in your opinion on the following topic
What do you think is the best way to solve a rather big system of linear equations inside a elmt routine. Until now I used the lapack utility dgesv(*).
Thank you! ;D
-
Depends on the structure of the equations. If dense, then any good solver, if it has structure you may be able to find a more efficient ordering.