Go through FORTRAN Code to find what each "chunk" does #14
Comments
|
The following code opens input and output file. It reads the first line (the title) and copies it into the output file: Lines 6 to 9 in 0ee7174 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 Lines 15 to 24 in 0ee7174 Lines 104 to 114 in 0ee7174 The following code calculates DEPF and DEPPR and prints them with format strings on lines 115 and 116 : Lines 26 to 29 in 0ee7174 Lines 115 to 116 in 0ee7174 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: Lines 30 to 49 in 0ee7174 Line 32 just writes the header of the table that is about to be printed. My problem in understanding comes from 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: 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: Lines 51 to 53 in 0ee7174 Thus, we can deduce the following naming:
|
|
The following code snippet reads in a bunch of values from Lines 55 to 64 in 0ee7174
|
|
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 Lines 66 to 90 in 0ee7174
|
|
Great work so far @EmilSoleymani. I have some comments on the parts of your review that I can clarify:
|
Slab SubroutineThe Lines 221 to 263 in 0ee7174 Recall that the variables
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. |
|
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. |
|
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:
As we have figured out so far, the program assigns a layer number to each soil element. It then calculates the effective overburden pressure which is prior to the placement of the foundation. The subroutine 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 |
|
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. |
Subroutine MECHCalculates 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! Lines 159 to 165 in 0ee7174 Examining line 164, we see a calculation for final void ratio. This calculation includes taking log_10 of a variable called Examining line 165, we see that if
|
Subroutine Mech (continued)On page 30 of the settlements textbook, there is a step-by-step algorithm outlining the procedure of primary consolidation calucation: Lines 167 to 173 in 0ee7174 The equation is equation (3-20) on page 28 of the settlement analysis book:  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 |
|
While reading about the |
|
@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? |
|
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: |
Subroutine LeonSubroutine Leon calculates immediate settlement of shallow foundations on granular soils (Leonard's and Frost Model). It's formula can be found on page 97: Lines 297 to 299 in 0ee7174 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. Lines 301 to 302 in 0ee7174 UW is calculated in the textbook by multiplying z_w and gamma_w (GAW), and that is shown exactly in the code. Line 310 in 0ee7174 The value of K_d, the horizontal stress index, is calculated on line 311 as the variable AKD. Its formula matches my notes.Line 311 in 0ee7174 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: Line 313 in 0ee7174 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. Line 314 in 0ee7174 The next line calculates the coefficient of earth pressure and calls it AK0, with the exact formula as in my notes:Line 315 in 0ee7174 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.Lines 316 to 322 in 0ee7174 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.Lines 323 to 333 in 0ee7174 Line 334 calculates the over consolidation ratio, OCR, with the same formula I have in my notes: Line 334 in 0ee7174 The next line calculates sigma_p prime, the preconsolidation/max past pressure, by multiiplying OCR and sigma_v prime. Line 335 in 0ee7174 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.Lines 336 to 337 in 0ee7174 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.Lines 338 to 339 in 0ee7174 |
|
|
|
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. |
|
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. |
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 Lines 331 to 333 in 0ee7174 Note: I just found an online version of the book here that had the correct phi_ax formula:  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:  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):  Lines 343 to 348 in 0ee7174 Lines 350-357 are a repeat of 341-349, but with some slight differences. The variable NBPRES decides which of the two blocks executes:Lines 340 to 357 in 0ee7174 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:  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:  Lines 358 to 360 in 0ee7174 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.
|
|
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. |
Subroutine SchmertSubroutine 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: 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: Lines 623 to 628 in 0ee7174 Lines 629 and 630 calculate constant C_t based on formula given above: Lines 629 to 630 in 0ee7174 Lines 636-638 calculate E_si as described in textbook: Lines 636 to 638 in 0ee7174 Lines 639-656 calculate the influence factor, I_zp: Lines 639 to 656 in 0ee7174 Finally, lines 657-658 put it all together, calculating the final value of the formula which is what is written to the output file: Lines 657 to 658 in 0ee7174 |
|
Another great summary @EmilSoleymani. |
|
Just to help me remember what I have left to do:
|
Update:In subroutine Lines 323 to 324 in 0ee7174   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:Lines 323 to 329 in 0ee7174 |
|
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. |
Block DataIn order to figure out what is going on in the 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 DATATESTOutput: 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.100At 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. |
|
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. |
Subroutine BICUBE DataSimilarly, the Lines 556 to 559 in 0ee7174 To make things even harder to wrap your head around, these variables are used to access the XX and YY arrays seen before:Lines 565 to 568 in 0ee7174 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 98There 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 |
|
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.
Element 1 is made up of nodes 1, 2, 13, 12. Element 2 is made up of nodes 2, 3, 14, 13, etc. |
|
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: Lines 565 to 566 in 0ee7174 |
|
Having began the code for calculating effective stresses, I have matched another equation to the theory: Line 74 in 0ee7174
The process follows exactly what Dr. Stolle presented to us:
You can see the code subtracts Line 75 in 0ee7174 This matches the slide, where effective stress = (gamma_sat - gamma_w)z |
|
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. |
|
The following lines seem to make an adjustment to the effective stress calculation: Lines 81 to 88 in 0ee7174 These adjustments follow method 1 from Figure 5-1 on p.99 of the Settlements Analysis book.
In our code, the if statement
|
|
@EmilSoleymani do you know the definition of 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.) 😄 |
|
@smiths DGWT = depth to ground water table |
|
I think we can confidently state that we have completed our analysis of the FORTRAN code. |
|
Agreed |
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:
vdisp/fortran_src/vdisp.for
Lines 1 to 9 in 0ee7174
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. 😄 I feel like the W means water. G might be specific gravity? We'll figure it out.
The text was updated successfully, but these errors were encountered: