Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

ULARC

JohnNichols
Valued Contributor III
1,148 Views

I have been using this program since 1981.  It is a great program, but I have never used it with multiple load cases, really do not need to do that with coal loads - so big everything else is a nuisance.

I have been writing a paper over the last two days on beam elements with varying moments of inertia -- not a common element, but my interest in some old math stuff I have not touched since 2004 arose from a paper by Dr Lewis from England -- quite interesting in an Ivory Tower sort of way. I thought I would use the arch examples in the paper on ULARC, I had used it for a lot of arches.

Lewis gives an example of an arch with zero moments -- nifty really - so I was putting the data into ULARC and getting the output as 3 loads, self weight, applied load and total load. 

I was plotting the results and thought - they look strange.  Went got the manual -- seemed normal operation for using 3 loads one after the other -- but that is not the case.  Took me  awhile to track down the code location - total load at each node is summed after each analysis --

C
C     UPDATE MEMBER QUANTITIES
C
      DO 310 NE=1,NELS
C
              DO 270 I=1,3
                  STOT(NE,I)=STOT(NE,I)+RATIO*S(NE,I)
                  VP(NE,I)=VP(NE,I)*RATIO
  270         VPTOT(NE,I)=VPTOT(NE,I)+VP(NE,I)
C
C     ACCUMULATE MAXIMUM ROTATIONS
C
          DO 280 I=2,3
              IF (VPTOT(NE,I).GT.VPTP(NE,I)) VPTP(NE,I)=VPTOT(NE,I)
              IF (VPTOT(NE,I).LT.VPTN(NE,I)) VPTN(NE,I)=VPTOT(NE,I)
  280     CONTINUE
C
C     SET NEW HINGE CODES
C
          IF (IBUST.NE.0) GO TO 310
          IF (NE.NE.INEL) GO TO 300
          IHI(NE)=1
          write (1,290) NE,NODI(NE)
  290     FORMAT(//20H HINGE FORMS, MEMBER,I4,6H, NODE,I4)
  300     IF (NE.NE.JNEL) GO TO 310
          IHJ(NE)=1
          write (1,290) NE,NODJ(NE)
  310 CONTINUE

Weird

0 Kudos
20 Replies
Steven_L_Intel1
Employee
1,136 Views

Bumping this so people see it - it was "in limbo" for a while.

0 Kudos
jimdempseyatthecove
Honored Contributor III
1,136 Views

It would seem like you have missing code, before line 13:

          VPTP(NE,1)=VPTOT(NE,1) ! add this
          VPTN(NE,1)=VPTOT(NE,1) ! add this
          DO 280 I=2,3
              IF (VPTOT(NE,I).GT.VPTP(NE,I)) VPTP(NE,I)=VPTOT(NE,I)
              IF (VPTOT(NE,I).LT.VPTN(NE,I)) VPTN(NE,I)=VPTOT(NE,I)
  280     CONTINUE

The code without the statements may have worked by accident.
(IOW the largest/smallest were larger and smaller than 0.0, and were not the first element, and the compiler performed zero initialization)

Jim Dempsey

0 Kudos
Steven_L_Intel1
Employee
1,136 Views

jimdempseyatthecove wrote:

(IOW the largest/smallest were larger and smaller than 0.0, and were not the first element, and the compiler performed zero initialization)

There is no zero initialization unless you ask for it. You might happen to get zeroes anyway, but there's no guarantee.

0 Kudos
andrew_4619
Honored Contributor II
1,136 Views

That looks like a bang on the money on eagle-eyed spot Mr D.

Do we have some more spam attacks? It seems every post I make the last couple of days needs a captcha verification.

0 Kudos
JohnNichols
Valued Contributor III
1,136 Views

I enclose the complete program.

It has had some modernization since it was on punch cards, all variables declared and the output improved, but the core remains as written in about '72 by a Masters Student. I made a lot of money with this program in the early 80's.

When you run - input file is b.txt, any output name will do and it creates a file called R.CSV for doing plotting with EXCEL.

I really should update it -- but there is a lot to do and only so much time.

 

0 Kudos
JohnNichols
Valued Contributor III
1,136 Views
C
C     INITIALIZE IF NEW LOAD SEQUENCE
C
      IF (NLC.EQ.1) GO TO 40
      IF (LCODE.EQ.0) GO TO 90
   40 DO 50 NE=1,NELS
      IHI(NE)=0
   50 IHJ(NE)=0
      IBUST=0
      DO 60 M=1,NSJTS
          SPFRA(M)=0.
          SPFRB(M)=0.
   60 SPFRR(M)=0.
      DO 70 NJ=1,NJTS
          DO 70 I=1,3
   70 RRTOT(NJ,I)=0.
      DO 80 NE=1,NELS
          DO 80 I=1,3
              VPTOT(NE,I)=0.
              STOT(NE,I)=0.
              VPTP(NE,I)=0.
   80 VPTN(NE,I)=0.

I have to work out what they mean about a new load sequence. Also need to add FEAST - interesting playing with old code

0 Kudos
Steven_L_Intel1
Employee
1,136 Views

andrew_4619 wrote:

Do we have some more spam attacks? It seems every post I make the last couple of days needs a captcha verification.

Unfortunately, yes. We hope to turn this back off as soon as the spam is under control again. (Users with "Black Belt" status are exempt from this.) Sorry for the inconvenience.

0 Kudos
JohnNichols
Valued Contributor III
1,136 Views

I was black belt once -- sadness clouds my vision at my status loss

0 Kudos
jimdempseyatthecove
Honored Contributor III
1,136 Views

Your arrays were initialized to 0.0 (post #7), so then shouldn't your DO 280 loop in post #1 iterate from1,3?

*** Note, this would required 0.0 to be such that the (NE,I) index got overwritten when I==1 ***
It is much safer to simply unconditionally copy (NE,1), then conditionally copy (NE,2:3).
It is also more efficient and not subject to breaking code when VTOT(NE,1) is the highest number .OR. VTOT(NE,1) is the lowest number.

The original code (#1) never tested the (NE,1).

Jim Dempsey

0 Kudos
JohnNichols
Valued Contributor III
1,137 Views

Excellent question - I will look tonight

0 Kudos
JohnNichols
Valued Contributor III
1,137 Views

jimdempseyatthecove wrote:

Your arrays were initialized to 0.0 (post #7), so then shouldn't your DO 280 loop in post #1 iterate from1,3?

*** Note, this would required 0.0 to be such that the (NE,I) index got overwritten when I==1 ***
It is much safer to simply unconditionally copy (NE,1), then conditionally copy (NE,2:3).
It is also more efficient and not subject to breaking code when VTOT(NE,1) is the highest number .OR. VTOT(NE,1) is the lowest number.

The original code (#1) never tested the (NE,1).

Jim Dempsey

Your question reminded me of the Donald Sutherland line(s) in Kelly's Heroes.

"I cannot fix the tanks - I just drive them"

"Moriarity has made it so we can get out of battle faster than we got into it."

"The only way I have to keep them Tiger tanks busy is to let them blow holes in me."

"Its a Tiger tank and we have it by the STOP"

The end of another long day of coding and paper writing.

0 Kudos
John_Campbell
New Contributor II
1,137 Views

For loop 280, I=2,3 VP are the extreme end rotations or moments of the member, while I=1 is the axial force/extension in the member. The loop approach should only be for rotations 2,3.

I have looked back at my version of ULARC which was obtained in 1979 (and I last used in 1989). I also have a printed (quarto) manual somewhere. I have forgotten how multiple load cases were combined. If each load case goes independently to collapse, then it makes no sense to combine load cases, but if each load case is progressively applied to it's full load and then the last load case iterates to collapse, then the load cases could be accumulated. From memory, progressive loading is likely the case, as there is the capability for hinges to become refixed.

As John has mentioned previously, ULARC was a precursor to many programs where collapse was modelled. It is amazing for it's brevity.

 

0 Kudos
mecej4
Honored Contributor III
1,137 Views

The documentation, source code, sample input and output files of ULARC are available at https://nisee.berkeley.edu/elibrary/software.html . You have to sign up first.

0 Kudos
DavidWhite
Valued Contributor II
1,137 Views

Steve Lionel (Intel) wrote:

Quote:

 

Unfortunately, yes. We hope to turn this back off as soon as the spam is under control again. (Users with "Black Belt" status are exempt from this.) Sorry for the inconvenience.

We had better get the Black Belt Ninja's onto this :-)

0 Kudos
jimdempseyatthecove
Honored Contributor III
1,137 Views

>>The loop approach should only be for rotations 2,3.

OK, then is it proper to have the initial comparison against 0.0 for both the min and the max???

IOW for the max value, are the rotations always positive? (0.0 less than either VTOT(NE,2) or VTOT(NE,3))
IOW for the min value, are the rotations always negative? (0.0 greater than either VTOT(NE,2) or VTOT(NE,3))

The two above requirements are mutually exclusive. Do your not find that strange?

Bear in mind your post #6: It has had some modernization since it was on punch cards...

It was not unusual to drop a deck of cards, then pick them up in a different order (or missing a card or two that flew under equipment).

Jim Dempsey

0 Kudos
jimdempseyatthecove
Honored Contributor III
1,137 Views

>>The two above requirements are mutually exclusive. Do your not find that strange?

Let me answer my own question...

It would not be strange if element could express stress (compression force (+) or tension force (-)), and you do not know which, but you want the min and max about a 0.0 external force.

Jim Dempsey

0 Kudos
John_Campbell
New Contributor II
1,137 Views

Jim,

I am not sure how this relates to the reported problem. My recollection is that "VPTOT" is the accumulated member displacements for this iteration, while VPTP is the most positive value that occurs and VPTN is the most negative value that occurs over the range of iterations. These arrays were initialised in Subroutine LOAD for the first load case or if the load application code "LCODE" /= 0 (so loading could accumulate or restart). I don't see any problem with what is coded. ( Also this is structural analysis; tension is +ve !)

  IF (VPTOT(NE,I).GT.VPTP(NE,I)) VPTP(NE,I)=VPTOT(NE,I)
  IF (VPTOT(NE,I).LT.VPTN(NE,I)) VPTN(NE,I)=VPTOT(NE,I)

Actually, I wrote a program ULARCPLOT in about 1980 which ran on a Tektronix green screen. Lost now. I am not sure if we ever put hot spots on the hinges, but it would have looked good. Those screens made it difficult to do a moving plot. My recollection is there was not a lot of demand for ULARC back then as it was more for ultimate analysis and did not suit the then design standard.

John

0 Kudos
andrew_4619
Honored Contributor II
1,137 Views

John Campbell wrote:
I am not sure how this relates to the reported problem. 

To be honest having just re-read the OP what actually is the reported problem and what is the relevance of in initial code snipped to that problem?

0 Kudos
JohnNichols
Valued Contributor III
1,137 Views

I have provided the CODE - it has been updated a little from the stuff available at Berkeley so all the variables are declared.  I include a sample file, the sample file once you read the output there is not much missing. There are actually two manuals, one is a brief set of sheets which I have at home -- I will hunt them out -- the other is the masters report and I can pdf it.

ULAC is useful for research and for designing super-heavy civil things like Coal Bins, 4000 tonnes is not uncommon. Mix it with AXISHELL and you have a rapid analysis tool.  But you need to understand the underlying math well, if you have spent your life designing house beams this is not for you.  

The original comment related to the load cases and loop 270. When you are doing a large plastic analysis really you only use one load case, but I was looking at some interesting stuff by Lewis and I needed three load cases, dead, live and then combined load, but unfortunately I found out that the load cases are added as you proceed through -- not what I wanted -- I tracked the "feature" down to line 270 loop - I will fix it when I have a minute, but now I just do three runs.

The attached paper - submitted yesterday to Royal Society Proceedings A shows the use of the program. I added a quick CSV file so I could plot the displacement. Adding DISLIN to ULARC to plot the results is not hard -- when I have a bit of time.

Finally: I would rather see a 4th year structural student exposed to Intel Fortran, Gavin's stuff out of Duke Uni in terms of notes and shown a simple program like ULARC, AXISHELL and asked to write some code than a year using MAPLE or Mathematica to solve structures problems. Once you can set up a simple structures problem in a TEXT file - and you can compile the file yourself then the ABAQUS programs are trivial monsters.

John

 

 

0 Kudos
jimdempseyatthecove
Honored Contributor III
1,063 Views

John,

>>I was plotting the results and thought - they look strange.

In an iterative integration method round off errors are accumulative. These are especially prevalent when accumulating values that very greatly in magnitude. Using larger precision numbers can help, but better algorithms (or both) produce better results.

If I were to make an (uninformed) guess is your plots were produced by integrating the forces (in small steps) across the span of the arch. This often results in the integration as it progresses into a relatively large number as compared to each remaining step to be accumulated. IOW the effects of error tends to grow as the number of integrations steps increase. There are programming techniques you can use to reduce this effect.

real(8) :: TotalIntegral, TotalFraction, Temp
...
TotalIntegral = 0.0
TotalFraction = 0.0
DO I=1,LargeNumberOfIntegrationSteps
    ...
    Temp = ANINT(CurrentStepValue)          ! whole part
    TotalIntegral = TotalIntegral + Temp    ! accumulat whole part
    Temp = CurrentStepValue - Temp          ! fraction part
    TotalFraction = TotalFraction + Temp    ! accumulate fraction part (may produce results with whole part)
    Temp = ANINT(TotalFraction)             ! whole part of new TotalFraction
    TotalIntegral = TotalIntegral + Temp    ! propigate as carry
    TotalFraction = TotalFraction - Temp    ! reduce to fraction part
END DO

Essentially you are performing multi-precision via software.

*** Note, as you increase the number of integration steps, you should reduce the unit size, such that the unit size is on the order of the expected (average) CurrentStepValue. After integration, you can then convert back to the external units desired by your program.

Jim Dempsey

0 Kudos
Reply