smiths / vdisp

Software for 1D Settlement Analysis
2 stars 0 forks source link

Go through FORTRAN Code to find what each "chunk" does #14

Closed smiths closed 2 years ago

smiths commented 2 years ago

To create new code, we need to figure out what the old code does. We should go through the FORTRAN CODE and say what each "chunk" of code does. I think we can use the features of GitHub to make this relatively easy to do. For instance, we can do something like:

Code to initialize program, declare variables for PP, G and WC, and open files: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L1-L9

We should probably also keep a dictionary of variable names and their meaning. For instance, I don't know what PP, G and WC are. :smile: I feel like the W means water. G might be specific gravity? We'll figure it out.

EmilSoleymani commented 2 years ago

The following code opens input and output file. It reads the first line (the title) and copies it into the output file: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L6-L9

The following code reads the second line from the input file. It parses NPROB, NOPT, NBPRES, NNP, NBX, NMAT, DX. It first prints NPROB, NNP, NBX, NMAT, DX. NOPT and NBPRES are used to determine soil type and foundation type. The if statements select what FORMAT string to print from lines 104-114: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L15-L24 https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L104-L114

The following code calculates DEPF and DEPPR and prints them with format strings on lines 115 and 116 : https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L26-L29 https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L115-L116

Here is a dictionary I have come up with so far for these code snippets, based on the format strings:

Having a hard time following the specifics of this part, however I understand it prints out the first table which is Number of Soil vs. Element. It prints each entry with format 32 on line 118. At the end it prints the header for the next table, which shows the specific gravity, water content and void ratio for each material: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L30-L49

Line 32 just writes the header of the table that is about to be printed. My problem in understanding comes from IE(N,1) on line 33. I understand that IE(NQ,1) is defined in a COMMON block on line 2 in order to be accessed in other subroutines, however I don't understand what it is - a function? a variable? what are we passing into it? how are we reading from it if I can't find its definition anywhere in the document? However, I do understand what is happening, just now HOW it's happening. In VDIN.DAT we have:

1 1 

12 2

16 2

This is read by the code snippet, and in a table fills in all the intermediate values ( > 1 and < 12 is 1, > 12 and < 16 is 2) like so:

ELEMENT NUMBER OF SOIL

    1            1
    2            1
    3            1
    4            1
    5            1
    6            1
    7            1
    8            1
    9            1
   10            1
   11            1
   12            2
   13            2
   14            2
   15            2
   16            2

Now after printing this, and the header for the next table ( MATERIAL SPECIFIC GRAVITY WATER CONTENT,% VOID RATIO ), the following code enumerates the materials (number of different soil layers) and lists each ones Specific gravity, water content and void ratio: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L51-L53

Thus, we can deduce the following naming:

EmilSoleymani commented 2 years ago

The following code snippet reads in a bunch of values from VDIN.DAT. It prints out some info depending on the values (specifics listed below): https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L55-L64

EmilSoleymani commented 2 years ago

The following code calculates effective overburden pressure (According to this link it could also be called effective stress which is what our geotechnical engineering notes call it). There are two arrays, P and PP, which seem to always both have the same value within this block of code (thanks to PP(I)=P(I) on line 77). It seems to loop through all the nodal points (I=2,NNP), retrieving the value from the table we calculated before here, which tell us the material at each nodal point? We use this material's Specific Gravity, G(MTYP), water content, WC(MTYP) (divided by 100 probably for unit conversions - UPDATE: DIVIDE BY 100 SINCE INPUT IS GIVEN AS PERCENT BUT WE NEED IT AS DECIMAL), and void ratio, EO(MTYP) to make the calculations for some value GAMM (Gamma?). There is a constant GAW=0.03125 that I can't figure out the name for that is also used in these calculations. The pressure of each nodal point is the pressure of the previous one + DX*GAMM. Calculations seem similar to what is happening in this Chegg article about calculating effective overburden pressure. Line 81 skips directly to calling of SLAB subroutine if IOPTION=1 or NOPT != 0, else there are a bunch of calculations with variables who's names or purpose I can not currently figure out prior to calling the subroutine: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L66-L90

Note: There is a typo on line 85, it should say DO I=1,MO

smiths commented 2 years ago

Great work so far @EmilSoleymani. I have some comments on the parts of your review that I can clarify:

EmilSoleymani commented 2 years ago

Slab Subroutine

The SLAB subroutine seems to use some of the previously seen data to find the respective surcharge pressure values and store them in the array P. https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L221-L263

Recall that the variables NBPRES and MRECT help determine the type of foundation and whether the calculation is at the edge of the slab or at the centre.

We can use Julia's @enum macro for NBPRES and a simple Boolean for MRECT.

The code is just a series of calculations. It is hard to understand the meaning of all the intermediate variables used in the calculations, assuming they have one at all, however it will be trivial to replicate the same process.

smiths commented 2 years ago

We should try to find the SLAB equations in the settlement book (pdf). I agree that they are hard to deduce from reading the code. We can just translate the code to Julia, but I would like to also figure out what theory this code is based on. We'll need this information for our documentation.

EmilSoleymani commented 2 years ago

Appendix F on Page 97 of the Settlments Book will be very useful, I was not aware there was an appendix about the code. It has helped me add to our dictionary of variables:

The remaining five subroutines are explained in the appendix as well and I will be looking over those next. Only one of them gets called depending on NOPT.

smiths commented 2 years ago

Yes, this is very valuable information. I'm glad to see that our guesses were on the right track, but this is a quicker way to understand the code.

EmilSoleymani commented 2 years ago

Subroutine MECH

Calculates heave or consolidation from result of one-dimensional consolidometer swell tests performed on cohesive soil for each layer in soil profile. After examining the settlement book and geotechnical engineering notes, I seem to have matched lines 159-165 to some formulas/theory! https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L159-L165

Examining line 164, we see a calculation for final void ratio. This calculation includes taking log_10 of a variable called CA. Examining line 161 indicates that CA is the ratio of two swell pressures. image This matches formula (3-21) on page 28 of the settlements book perfectly. In the picture they are finding change in void ratio, but changing that to (final void ratio - initial void ratio) and moving the initial to RHS produces the calculations on line 164 exactly

Examining line 165, we see that if PR > PM(MTYP) we do a similar calculation with an extra log term. This can be matched to equation (3-23) on page 30 of the settlements book: image The inequality on the lines below the equation explain the IF statement in the FORTRAN code. We can now deduce the following:

EmilSoleymani commented 2 years ago

Subroutine Mech (continued)

On page 30 of the settlements textbook, there is a step-by-step algorithm outlining the procedure of primary consolidation calucation: image In the previous comment ( lines 159-165 ) we performed the calculations for Step 6 of this procedure. Examining the lines that come right after, lines 167-173, we see a similar structure to the formula in Step 7: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L167-L173 The equation is equation (3-20) on page 28 of the settlement analysis book: image Note that the numerator here is what we calculate in Step 6 which was covered in the previous comment.

It seems Step 7, which is just summing up all the individual settlements, is incorporated into the code by just adding current DEL to DELH1, thus. accumulating all settlements into DELH1. Lines 174-194 seem to be a repeat of the same procedure with a different variable DELH2, just with a different DXX value. This code seems to execute when NOUT=0. Regardless, DELH2 is added is added onto DELH1. From my understanding, since NOUT=0 means we only output heave computations, DELH2 is the result of heave computations. Following that logic, since NOUT=1 means we have heave computations AND pressure at each depth increment, DELH1 stores the value of sum of pressures at each depth after execution of lines 159-173 as discussed previously, then when DELH2 is added to it it stores the final value that the user wants.

EmilSoleymani commented 2 years ago

While reading about the LEON subroutine and its (many) resulting formulas, I think I might have stumbled upon what the constant GAW is on page 4 of the settlements textbook: image This makes sense since gamma_w could be abbreviated to GAW. The textbook gives a version accurate to 3 decimal places, while the code uses the slightly more accurate 0.03125.

smiths commented 2 years ago

@EmilSoleymani, I'm sure you are right about GAW being the specific weight of water. However, the 0.031 tsf(?) is weird. The specific weight of water is around 9.807 kN/m^3 in SI. I thought that tsf might be the imperial equivalent, but my google search suggests that tsf is tonnes per square foot. tsf is a pressure unit. It applies to an area, while specific weight applies to a volume. I'm not sure how to rationalize this. Your text is cut off though. Is there something after tsf that could be a clue?

smiths commented 2 years ago

Your work on deducing the meaning of the code fragments looks great. This will take some iteration, but you are making great progress.

With respect to your question about math in GitHub markdown, I googled it. There appears to be a hack that works. Here is an example:

If you need them for some reason, here is a list of GitHub markdown emoji markup:

https://gist.github.com/rxaviers/7360908

EmilSoleymani commented 2 years ago

Subroutine Leon

Subroutine Leon calculates immediate settlement of shallow foundations on granular soils (Leonard's and Frost Model). It's formula can be found on page 97: image The book gives a complicated breakdown of each variable and I had to go all over the book looking for intermediate formulas. I took notes on everything I found and they are available for download in the following pdf: Leonard’s and Frost Model.pdf I will now be comparing the code to my findings in my notes. First of all, the delta p is referred to as QNET in the code. It is calculated and printed in lines 297-299 of the code: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L297-L299 Constant c1 is calculated on lines 301-302. The formula matches perfectly with my notes, and as stated in book this constant must be greater than 0.5 thus it is adjusted if needed on line 302. https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L301-L302 UW is calculated in the textbook by multiplying z_w and gamma_w (GAW), and that is shown exactly in the code. image https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L310 The value of K_d, the horizontal stress index, is calculated on line 311 as the variable AKD. Its formula matches my notes. https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L311 I am not sure what is being calculated on line 312. I have searched all over the book for the calculation (P_1 - P_0)/(P_0 -u_w) with no success. It might be the result of dividing two formulas by each other. I will keep looking for this variables meaning. The next line calculates the value of E_D just as expected: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L313 The next line calculates q_c/(sigma_v prime) which is part of the last term of calculating the coefficient of earth pressure, K_oc. https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L314 The next line calculates the coefficient of earth pressure and calls it AK0, with the exact formula as in my notes: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L315 Lines 316-322 deal with calculating constant phi_ps, the plane strain angle, which as in my notes is calculated using the BICUBE subroutine. At the moment I have no clue what is going on in that subroutine and the book provides no help with regards to this, so I will be treating this as a blackbox for now that just outputs us phi_ps. Line 322 converts PHIPS from degrees to radians. https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L316-L322 Lines 323-333 deal with calculating phi_ax, the axial friction angle, from phi_s. Once again, the textbook makes no effort in explaining this. They give a formula which is in my notes that is completely different to what is going on in the code. Nevertheless, the correct value of phi_ax is stored in variable PHIAX, then converted to radians and stored in PHI in line 333. https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L323-L333 Line 334 calculates the over consolidation ratio, OCR, with the same formula I have in my notes: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L334 The next line calculates sigma_p prime, the preconsolidation/max past pressure, by multiiplying OCR and sigma_v prime. https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L335 Lines 336-337 calculate R_izoc, the ratio of stress increment of overconsolidated soil in layer to total stress increment in layer. Formula is same as in my notes. Value stored in variable ROC. It is adjusted to 0 if it was negative. https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L336-L337 Lines 338-339 calculate R_iznc, the ratio of stress increment of normally consolidated soil in layer to total stress increment in layer. Formula is same as in my notes. Value stored in variable RNC. It is adjusted to 0 if it was negative. https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L338-L339

EmilSoleymani commented 2 years ago

image image Here are the paragraphs on page 4 of the settlements textbook which refer to the unit weight of water with the units tsf. I too was confused upon reading this.

smiths commented 2 years ago

Great analysis @EmilSoleymani. I'm not surprised there are things that are not yet clear. There is a lot going on in this code. We will be able to ask Dr. Stolle questions too. The variables with K in them, like AKD, may be related to the bulk modulus. The A is probably just there as a prefix to show the variable is a float.

smiths commented 2 years ago

I wonder if the units of tsf work with gamma because they are talking about a unit width of soil. Gamma is in N/m^3 in general, but for a unit width of soil, it is in N/m^2.

We should be careful with the units. We'll want the code to also work for SI units.

EmilSoleymani commented 2 years ago

Subroutine Leon (continued)

Last time I stated that the phi_ax calculation was not matching the formula in my notes, however it actually was (almost). Lines 331-333 seem to be making the phi_ax calculation based on the formula in my notes where the value of phi_ps is stored in variable UC. There might be a typo in the textbook formula since it clearly says phi_ax = phi_px - [ (phi_ps - 32) / 3] while the calculations treat phi_px as phi_ps. This wouldn't be the first time I've noticed an error in the subscript of formulas in the textbook (looking at the picture for the main formula in my previous post about this subroutine, the constant C_1 is labeled as C_i? In the K_0 calculations in the book, K_D was written as K_p but explained as K_D?)

image

https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L331-L333 Note: I just found an online version of the book here that had the correct phi_ax formula:

image

Note: From another book I found out why there was a conversion from degrees to radians. The formula labelled as (F-4) above uses degrees in its calculations. This will be important for us when we reproduce this code so we don't perform conversions in the wrong steps:

image

Lines 343-348 deal with the calculation of the influence factor, I_zp, which is better explained in the Schmertmann model (from my understanding, Leon and Frosts model is a slight improvement upon the Schmertman model - they share many similarities):

image

https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L343-L348 Lines 350-357 are a repeat of 341-349, but with some slight differences. The variable NBPRES decides which of the two blocks executes: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L340-L357 There seems to be another typo in the settlements book, this time relating to the elastic soil moduli, E_izoc and E_iznc. They are both equal to some constant times E_D, while the book writes E_p:

image

Those formulas, paired with the following special case give us the calculations for the elastic soil moduli which are reflected in lines 358-360 of the code:

image

https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L358-L360 As seen above, line 360 does the final calculation which is added onto variable DELH, which is the value this subroutine writes to the output file.

smiths commented 2 years ago

This is a great review of the code and the settlement book. I know it is frustrating when the book and the code don't agree. This is the motivation for having one source for the documentation and the code, as done in literate programming and even more dramatically with Drasil.

The examples where the documentation and the code disagree, or where there are mistakes in the documentation, will be very valuable for future papers. I'm going to create an issue for us to keep track of those.

EmilSoleymani commented 2 years ago

Subroutine Schmert

Subroutine Schmert calculates immediate settlement of shallow foundations on granular soils (Schmertmann model). As mentioned in subroutine Leon explanation, the Leonards and Frost model is very similar to the Schmertmann model. Consequently, most of the calculations remain the same with some slight differences being present. The formula of the Schmertmann model is as follows, found on p.104 of the settlments textbook:

image

Note: This subroutine is much simpler than the LEON subroutine. The calculations for R_izoc and R_iznc (and all their subsequent formulas like OCR, phi_ax, phi_ps) were the reason for many lines of code. Since these variables don't exist, the remaining calculations are quite simple.

The first big difference of this model is the presence of a new constant, C_t. This is the correction for time-dependant increase in settlement. C_t can be found by the following formula where t is the time in years following construction when settlement is to be calculated for:

image

The second difference is the calculation of the soil modulus, E_si:

image

An interesting point here is that when NOPT=2 we call SCHMERT and expect to read q_c data from the input file, but when NOPT=4, E_si data for each layer of soil is passed in directly (useful for when q_c data is not availabl, but E_si can be determined through tests). Lines 623-628 are same as in the LEON subroutine - they initialize some variables and calculate delta_p in variable QNET. Then the constant c_1 is calculated and stored in variable C1, which will have a value no less than 0.5: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L623-L628 Lines 629 and 630 calculate constant C_t based on formula given above: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L629-L630 Lines 636-638 calculate E_si as described in textbook: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L636-L638 Lines 639-656 calculate the influence factor, I_zp: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L639-L656 Finally, lines 657-658 put it all together, calculating the final value of the formula which is what is written to the output file: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L657-L658

smiths commented 2 years ago

Another great summary @EmilSoleymani.

EmilSoleymani commented 2 years ago

Just to help me remember what I have left to do:

EmilSoleymani commented 2 years ago

Update:

In subroutine LEON, I previously had no idea where the calculations for variables AKA and AKP came from on lines 323 and 324. I found their formulas in the Geotechnical Engineering notes: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L323-L324

image image

Still not sure why they are used. They seem to be used as some sort of "bounds". When the calculated K_o value is less than K_a or greater than K_p, it is adjusted back to within bounds. It seems when the changes exceed a difference of 0.01, we redo the calculations by calling the BICUBE subroutine: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L323-L329

smiths commented 2 years ago

Your explanation (about bounds) sounds reasonable @EmilSoleymani. Once we have a set of questions along these lines, we can talk to Dr. Stolle about them.

EmilSoleymani commented 2 years ago

Block Data

In order to figure out what is going on in the BLOCK DATA code chunk, I wrote a quick script which replicated the initialization of XX, YY, and U arrays and wrote their outputs to console:

      PROGRAM DATATEST
      REAL XX(100), YY(100), U(100)
      DATA (XX(I),I=1,99,11)/9*10./,(XX(I),I=2,99,11)/9*20./,(XX(I),I=3,
     199,11)/9*30./,(XX(I),I=4,99,11)/9*50./,(XX(I),I=5,99,11)/9*100./,
     2(XX(I),I=6,99,11)/9*200./,(XX(I),I=7,99,11)/9*300./,(XX(I),I=8,99,
     311)/9*500./,(XX(I),I=9,99,11)/9*1000./,(XX(I),I=10,99,11)/9*2000./
     4,(XX(I),I=11,99,11)/9*3000./
      DATA (YY(I),I=1,11)/11*0.16/,(YY(I),I=12,22)/11*0.20/,(YY(I),I=23,
     133)/11*0.4/,(YY(I),I=34,44)/11*0.6/,(YY(I),I=45,55)/11*0.8/,(YY(I)
     2,I=56,66)/11*1.0/,(YY(I),I=67,77)/11*2.0/,(YY(I),I=78,88)/11*4./,
     3(YY(I),I=89,99)/11*6./
      DATA (U(I),I=1,99)/25.,30.1,33.2,36.4,39.9,42.8,44.4,46.,48.5,50.5
     1,52.,24.8,30.,33.,36.2,39.7,42.6,44.2,45.8,48.2,50.2,51.5,24.5,29.7
     2,32.6,35.6,39.3,42.1,43.7,45.4,47.5,49.7,51.,24.2,29.5,32.2,35.1
     3,38.8,41.7,43.3,45.,47.2,49.,50.,24.,29.2,31.7,34.7,38.4,41.4,42.9
     4,44.6,46.8,48.6,49.7,23.8,28.8,31.5,34.4,38.,41.,42.5,44.3,46.5,48.
     5 4,49.5,23.,27.5,30.,33.,36.6,39.6,41.2,43.,45.4,47.2,48.4,22.,26.
     6,28.3,31.2,34.5,37.7,39.7,41.5,43.7,45.7,47.,21.,25.,27.2,30.,33.8
     7,36.,38.2,40.3,42.7,44.8,46.1/

      WRITE (*,*) "Testing Block Data from VDispl.for"

      WRITE (*,*) "Printing values in XX:"
      DO I=1,99
      WRITE (*,10) I, ": ", XX(I)
   10 FORMAT(I4, A, F8.1)
      END DO

      WRITE (*,*) "Printing values in YY:"
      DO I=1,99
      WRITE (*,20) I, ": ", YY(I)
   20 FORMAT(I4, A, F8.3)
      END DO

      WRITE (*,*) "Printing values in U:"
      DO I=1,99
      WRITE (*,30) I, ": ", U(I)
   30 FORMAT(I4, A, F8.3)
      END DO

      END PROGRAM DATATEST

Output:

 Testing Block Data from VDispl.for
 Printing values in XX:
   1:     10.0
   2:     20.0
   3:     30.0
   4:     50.0
   5:    100.0
   6:    200.0
   7:    300.0
   8:    500.0
   9:   1000.0
  10:   2000.0
  11:   3000.0
  12:     10.0
  13:     20.0
  14:     30.0
  15:     50.0
  16:    100.0
  17:    200.0
  18:    300.0
  19:    500.0
  20:   1000.0
  21:   2000.0
  22:   3000.0
  23:     10.0
  24:     20.0
  25:     30.0
  26:     50.0
  27:    100.0
  28:    200.0
  29:    300.0
  30:    500.0
  31:   1000.0
  32:   2000.0
  33:   3000.0
  34:     10.0
  35:     20.0
  36:     30.0
  37:     50.0
  38:    100.0
  39:    200.0
  40:    300.0
  41:    500.0
  42:   1000.0
  43:   2000.0
  44:   3000.0
  45:     10.0
  46:     20.0
  47:     30.0
  48:     50.0
  49:    100.0
  50:    200.0
  51:    300.0
  52:    500.0
  53:   1000.0
  54:   2000.0
  55:   3000.0
  56:     10.0
  57:     20.0
  58:     30.0
  59:     50.0
  60:    100.0
  61:    200.0
  62:    300.0
  63:    500.0
  64:   1000.0
  65:   2000.0
  66:   3000.0
  67:     10.0
  68:     20.0
  69:     30.0
  70:     50.0
  71:    100.0
  72:    200.0
  73:    300.0
  74:    500.0
  75:   1000.0
  76:   2000.0
  77:   3000.0
  78:     10.0
  79:     20.0
  80:     30.0
  81:     50.0
  82:    100.0
  83:    200.0
  84:    300.0
  85:    500.0
  86:   1000.0
  87:   2000.0
  88:   3000.0
  89:     10.0
  90:     20.0
  91:     30.0
  92:     50.0
  93:    100.0
  94:    200.0
  95:    300.0
  96:    500.0
  97:   1000.0
  98:   2000.0
  99:   3000.0
 Printing values in YY:
   1:    0.160
   2:    0.160
   3:    0.160
   4:    0.160
   5:    0.160
   6:    0.160
   7:    0.160
   8:    0.160
   9:    0.160
  10:    0.160
  11:    0.160
  12:    0.200
  13:    0.200
  14:    0.200
  15:    0.200
  16:    0.200
  17:    0.200
  18:    0.200
  19:    0.200
  20:    0.200
  21:    0.200
  22:    0.200
  23:    0.400
  24:    0.400
  25:    0.400
  26:    0.400
  27:    0.400
  28:    0.400
  29:    0.400
  30:    0.400
  31:    0.400
  32:    0.400
  33:    0.400
  34:    0.600
  35:    0.600
  36:    0.600
  37:    0.600
  38:    0.600
  39:    0.600
  40:    0.600
  41:    0.600
  42:    0.600
  43:    0.600
  44:    0.600
  45:    0.800
  46:    0.800
  47:    0.800
  48:    0.800
  49:    0.800
  50:    0.800
  51:    0.800
  52:    0.800
  53:    0.800
  54:    0.800
  55:    0.800
  56:    1.000
  57:    1.000
  58:    1.000
  59:    1.000
  60:    1.000
  61:    1.000
  62:    1.000
  63:    1.000
  64:    1.000
  65:    1.000
  66:    1.000
  67:    2.000
  68:    2.000
  69:    2.000
  70:    2.000
  71:    2.000
  72:    2.000
  73:    2.000
  74:    2.000
  75:    2.000
  76:    2.000
  77:    2.000
  78:    4.000
  79:    4.000
  80:    4.000
  81:    4.000
  82:    4.000
  83:    4.000
  84:    4.000
  85:    4.000
  86:    4.000
  87:    4.000
  88:    4.000
  89:    6.000
  90:    6.000
  91:    6.000
  92:    6.000
  93:    6.000
  94:    6.000
  95:    6.000
  96:    6.000
  97:    6.000
  98:    6.000
  99:    6.000
 Printing values in U:
   1:   25.000
   2:   30.100
   3:   33.200
   4:   36.400
   5:   39.900
   6:   42.800
   7:   44.400
   8:   46.000
   9:   48.500
  10:   50.500
  11:   52.000
  12:   24.800
  13:   30.000
  14:   33.000
  15:   36.200
  16:   39.700
  17:   42.600
  18:   44.200
  19:   45.800
  20:   48.200
  21:   50.200
  22:   51.500
  23:   24.500
  24:   29.000
  25:   32.600
  26:   35.600
  27:   39.300
  28:   42.100
  29:   43.700
  30:   45.400
  31:   47.500
  32:   49.700
  33:   51.000
  34:   24.200
  35:   29.500
  36:   32.200
  37:   35.100
  38:   38.800
  39:   41.700
  40:   43.300
  41:   45.000
  42:   47.200
  43:   49.000
  44:   50.000
  45:   24.000
  46:   29.200
  47:   31.700
  48:   34.700
  49:   38.400
  50:   41.400
  51:   42.900
  52:   44.600
  53:   46.800
  54:   48.600
  55:   49.700
  56:   23.800
  57:   28.800
  58:   31.500
  59:   34.400
  60:   38.000
  61:   41.000
  62:   42.500
  63:   44.300
  64:   46.500
  65:  484.000
  66:   49.500
  67:   23.000
  68:   27.500
  69:   30.000
  70:   33.000
  71:   36.600
  72:   39.600
  73:   41.200
  74:   43.000
  75:   45.400
  76:   47.200
  77:   48.400
  78:   22.000
  79:   26.000
  80:   28.300
  81:   31.200
  82:   34.500
  83:   37.700
  84:   39.700
  85:   41.500
  86:   43.700
  87:   45.700
  88:   47.000
  89:   21.000
  90:   25.000
  91:   27.200
  92:   30.000
  93:   33.800
  94:   36.000
  95:   38.200
  96:   40.300
  97:   42.700
  98:   44.800
  99:   46.100

At the moment I have no idea where these values are coming from. In the first array, XX, all indices i congruent modulo 11 seem to have the same value but I am not sure where the values 10,20,30,50,100,200,300,500,1000,2000,3000 come from. The other two arrays seem ever more random.

smiths commented 2 years ago

Good idea to print out the values @EmilSoleymani. It looks to me like the XX and YY values form a grid in 2D space. I'm guessing that it is a grid underneath the footing. I'm assuming that the x direction is horizontal and the y direction is vertical. The numbers in the x direction are close to the centre of the footing to start with, but then a good distance away at the end. I'm guessing u is the pore water pressure? I would think that this would be something that would be read in for each problem, but maybe it is hard-coded for some reason. Maybe u is something else? I'm fairly sure that u varies over the x-y grid though.

EmilSoleymani commented 2 years ago

Subroutine BICUBE Data

Similarly, the BICUBE subroutine has long hardcoded data. There is a 100x4 array of ints called KE. I replicated its initialization and checked what the values are. This array is accessed on lines 556-559 through a DO loop: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L556-L559 To make things even harder to wrap your head around, these variables are used to access the XX and YY arrays seen before: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L565-L568 Here is what the KE array looks like:

 Testing BICUBE Data Blocks
   M | I | J | K | L
   -----------------
   1   1   2  13  12
   2   2   3  14  13
   3   3   4  15  14
   4   4   5  16  15
   5   5   6  17  16
   6   6   7  18  17
   7   7   8  19  18
   8   8   9  20  19
   9   9  10  21  20
  10  10  11  22  21
  11  12  13  24  23
  12  13  14  25  24
  13  14  15  26  25
  14  15  16  27  26
  15  16  17  28  27
  16  17  18  29  28
  17  18  19  30  29
  18  19  20  31  30
  19  20  21  32  31
  20  21  22  33  32
  21  23  24  35  34
  22  24  25  36  35
  23  25  26  37  36
  24  26  27  38  37
  25  27  28  39  38
  26  28  29  40  39
  27  29  30  41  40
  28  30  31  42  41
  29  31  32  43  42
  30  32  33  44  43
  31  34  35  46  45
  32  35  36  47  46
  33  36  37  48  47
  34  37  38  49  48
  35  38  39  50  49
  36  39  40  51  50
  37  40  41  52  51
  38  41  42  53  52
  39  42  43  54  53
  40  43  44  55  54
  41  45  46  57  56
  42  46  47  58  57
  43  47  48  59  58
  44  48  49  60  59
  45  49  50  61  60
  46  50  51  62  61
  47  51  52  63  62
  48  52  53  64  63
  49  53  54  65  64
  50  54  55  66  65
  51  56  57  68  67
  52  57  58  69  68
  53  58  59  70  69
  54  59  60  71  70
  55  60  61  72  71
  56  61  62  73  72
  57  62  63  74  73
  58  63  64  75  74
  59  64  65  76  75
  60  65  66  77  76
  61  67  68  79  78
  62  68  69  80  79
  63  69  70  81  80
  64  70  71  82  81
  65  71  72  83  82
  66  72  73  84  83
  67  73  74  85  84
  68  74  75  86  85
  69  75  76  87  86
  70  76  77  88  87
  71  78  79  90  89
  72  79  80  91  90
  73  80  81  92  91
  74  81  82  93  92
  75  82  83  94  93
  76  83  84  95  94
  77  84  85  96  95
  78  85  86  97  96
  79  86  87  98  97
  80  87  88  99  98

There seems to be a pattern where every 11th number is skipped. Each column just starts at a different point. What is the significance of the number 11? XX has same values every 11 numbers, YY has the same values for each group of 11 numbers, and now KE entries are skipping each 11th number

smiths commented 2 years ago

I'm fairly sure that the data blocks represent a sequence of square regions/elements. The first number (M) is the element number, the next 4 numbers are the connectivity of its vertices (nodes), as shown in this figure.

ElementConnectivity

Element 1 is made up of nodes 1, 2, 13, 12. Element 2 is made up of nodes 2, 3, 14, 13, etc.

EmilSoleymani commented 2 years ago

That makes sense! Thus lines 565 and 566 are calculating the length/width of one of these squares assuming XX and YY are their coordinates in 2D space: https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L565-L566

EmilSoleymani commented 2 years ago

Having began the code for calculating effective stresses, I have matched another equation to the theory:

https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L74

image

The process follows exactly what Dr. Stolle presented to us:

image

You can see the code subtracts GAW from GAMM once the depth is past the depth of the ground water table:

https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L75

This matches the slide, where effective stress = (gamma_sat - gamma_w)z

smiths commented 2 years ago

Great analysis @EmilSoleymani. At some point we should document the effective stress concept in our SRS. This will definitely be part of our SRS.

We can talk about it more as the summer goes along, but I think the best approach for building the SRS will be to focus on the theoretical models first. We shouldn't try to classify them, just dump likely equations into the document. We can then start to organize as we move forward.

EmilSoleymani commented 2 years ago

The following lines seem to make an adjustment to the effective stress calculation:

https://github.com/smiths/vdisp/blob/0ee717473b63e17b6b0152c12a2181bdb9aff6e5/fortran_src/vdisp.for#L81-L88

These adjustments follow method 1 from Figure 5-1 on p.99 of the Settlements Analysis book.

image

In our code, the if statement NOPT.NE.0.OR.IOPTION.EQ.0 becomes model != ConsolidationSwell || !equilibriumMoistureProfile. So in this case it appears we add gamma_w * (z - z_a) to each entry, P(i). I think line 86, DGWT/DX - FLOAT(I-1) is performing the (z - z_a) calculation, and line 87 is altering the P array.

Note: The textbook correctly gives gamma_w units of ton/ft^3 here!

smiths commented 2 years ago

@EmilSoleymani do you know the definition of DGWT in the code? Is it the same as z? Since DGWT/DX - FLOAT(I-1) is multiplied by DX, we get DGWT - FLOAT(I-1)*DX. Assuming that DX is constant, this does look like z_a, since I-1 seems to be the number of increments.

I remember we talked about keeping a file of inconsistencies that were found in the code/documentation. Did you start that file? We should be sure to have the GAW example in there. (It is worth noting that in Drasil we couldn't have this kind of inconsistency - it would either be right everywhere, or wrong everywhere.) :smile:

EmilSoleymani commented 2 years ago

@smiths DGWT = depth to ground water table

EmilSoleymani commented 2 years ago

I think we can confidently state that we have completed our analysis of the FORTRAN code.

smiths commented 2 years ago

Agreed