FLOW3.DOC 24 July 1995 Summary of results. NX=41, NY=13, element size =1/4 Consider maximum discrepancy between sensitivity dU/dAlpha2 and finite difference Del U/Del Alpha2. (These are runs x12, x13, x14, x18 and x25). Smoothing options: 0, no smoothing 1, smooth at P nodes using ZZ. 2, smooth at P nodes using serendipity. 3, smooth at all nodes, average the multiple values at non-P nodes. 4, same as 3, but only smooth dUdY and dVdY, not dPdY. Max Discrepancy Smooth U V p 0 0.045 0.041 0.459 1 0.050 0.043 0.640 2 0.057 0.049 0.815 3 0.043 0.033 0.641 4 0.043 0.033 0.465 Smoothing option 4 seems the best, but not by much. I'll try a finer mesh before despairing, though. Runs x26, x27: NX=81, NY=25, element size=1/8 Max Discrepancy Smooth U V p 0 0.0182 0.0191 0.265 4 0.0167 0.0144 0.276 Note that we're also still seeing the roughly O(h) behavior in the decrease of the discrepancy from h=1/4 to h=1/8 when O(h*h) would be so much nicer. 23 July 1995 I programmed the Zienkiewicz-Zhou method for retrieving improved values of dUdY. And I don't get good values. I've tried a number of tricks, and now I'm going to try to solve a Stokes problem. (Both Max and Janet think it possible that the "ultraconvergence" mentioned by Z&Z occurs only in linear problems). OK, I can solve a simple Stokes problem. Now let me look at the behavior of Z&Z on a more interesting Stokes problem. 07 June 1995 Added FILEB so that I can run multiple jobs at one time. 06 June 1995 On the Cray, the finite difference calculation is giving very bad values. So first I'm making EPSDIF an input quantity. GBSOL: the ROW subroutine could be suitable as long as I had set up suitable sparse storage for the nonzero elements of the matrix. I wrote a routine SETNAB which is a start on this. Second question: fast I/O for reads and writes. 05 Jun 1995 Just submitted an NX=161 job on the Cray. I'm hoping it runs eventually. I would use this one (ALPHA=0) and two others (ALPHA=0.5, ALPHA=1.0) for runs on the ALPHA using IBC=3. Also, how about grading the mesh, in both X and Y, rather than just the closest row? The Cray runs don't seem correct, so I've copied the latest source file over, recompiling and rerunning shortly. 04 June 1995 I just ran an NX=81 job on the Cray. I'm going to try to bump it up to the next size. 27 May 1995 I got the FLOW3 program running on the PK Cray again. The output doesn't look sensible, though I haven't run INPUT.901 in a while. I'll check through it later. 24 May 1995 I am now ready to run, zero bump, with IBC=0, 1, 2, 3. I should also try using a refined grid near the bump. (Using IGRID=1, and XBORD, YBORD inputs...) Also, deleted MX, MY. 23 May 1995 After exiting GETGRD, you have restored the solution, but NOT the bump position, or whatever data was perturbed in the last finite difference calculation. So, I modified GETGRD so that it calls FLOSOL for the original data one last time. One more damn complication! Once I get some satisfactory answers for a zero bump, I will move on to use a bump parameter of 0.5 or 1.0. 22 May 1995 I implemented a three point difference formula for dUdY and got data about as good as using the finite element data, and much better than the two point finite difference. 21 May 1995 Merged SIMMAT into FPRIME. Dropped ICORR. Added "IBC", 0 means use finite element estimate for dUdY, 1 means use finite difference. Then add 2 to use finer mesh estimate. 15 May 1995 Max wants to try to produce some papers out of my thesis. The first topic to investigate is the influence of the boundary condition on the accuracy of the sensitivities. He suggests trying to get a more accurate boundary condition by: 1) Using a finer grid for the flow solution, to get good values of dudx. (Hard to program) 2) Refining the grid near the bump to get better values of dudx. Also, looking into getting access to the old Crays again. Alas, all I managed to do today was to get the program running, and put in a common output editor for QSOLVE, and the other routines. How about a parameter JGRID? But first, print out the sensitivities. 04 May 1995 I'm trying to replace the "-2" conditions by actual equations for the velocities on the bump. This is not easy, particularly since the quantity "du/dy" is not well defined at nodes, since du/dy is only well defined inside a particular element. I have replaced the -2 and -1 conditions. Also, I have made the V velocities along the outflow correspond to unknowns as well. 06 April 1995 The linear/quadratic bump problem has been rectified. Somehow, it had to do with the fact that I had made a copy of IOPT, called IOPT2, so that I could vary the Reynolds number if necessary, during a continuation. At the same time, I made it possible NOT to vary certain other things, and somehow the bump parameter was one of them. But only for linear or quadratic? 05 April 1995 Modified NEWTON to restart at best point, as SIMSOL already does. I really need to modify SOLCON to start from RE=1. I submitted INPUT.826, the first of a sequence for Chapter 14, solving at bigger and bigger Reynolds numbers. Something is wrong with the linear and quadratic bump descriptions. If I run an optimization with ishapb=1 or 2, I get zero bump, and zero bump sensitivities. I take this to mean that I'm accidentally computing a zero basis function... 01 April 1995 I really want a few more runs to fill up some tables and chapters. I left these jobs running: 816 - Trying unfeasible target at RE=100. Re=1 was a disaster. Could not some points. 821 - Trying NX=61, NY=19, hoping bump will simply plot properly. Seems to be fixed now... 823 - Trying TOLNEW=1E-7. TOLNEW=1.E-5 didn't get past RE=120. Could not compute target point. What seems to be different about the other runs? First, retry 823 with ISHAPB=ISHAPF=2. Retry 816 with ICORR=1, IGRAD=3. Ridiculous negative bumps. Replaced old 816 by new one: 3 varying inflow parameters, bump fixed at 0.5, RE=100. Target flow is parabolic. I might have found a flaw in 816: XBRTAR was not set, and so was at its default value of 3, while the feasible space was at 4. This put the variable count off, and made the cost functional evaluation ridiculously wrong. It's 10:30 on Saturday night, and I'm going home. 29 March 1995 Got rid of stupid JGRID variable. Max suggested trying Re=1000 jobs on a finer mesh on the Alpha3. I'm trying to get a new target data calculation off the ground. 22 March 1995 Presumably, the good combinations are: IPRED=IFDS, IDS+CORR IGRAD=IDS, IFDS-CORR, IDFD 21 March 1995 In order to make the argument that discretized sensitivities or finite difference sensitivities are needed, I restored OSOLVE, the functional value only version. In particular, note that OSOLVE would not necessarily know that a solution was being done at a finite difference point (though I fear I told OSOLVE that...) You need to clean up your input. You should have one set of input that tells you what things you compute: IDS, IFDS, ICORR, IDFD and another which tells you what you USE: IPRED=0, 1=IDS, 2=IFDS, 3=IFDS-CORR, 4=IDS+CORR IGRAD=0, 1=IDS, 2=IFDS, 3=IFDS-CORR, 4=IDFD I haven't tried the "corrected sensitivities". This would be useful for Newton predictions, where IDS+CORR is presumably more accurate than IDS. Now that I've split out GETFIX, it is easier to start treating the movement correction as a separate entity. 20 March 1995 1) I dropped IFIX. I want to be able to use uncorrected finite difference estimates for Newton prediction, but corrected finite differences for gradient estimation. I haven't tried this out yet. 2) I am trying to figure out why problem 568 gives a negative bump. So I moved down to REYNLD=1, and I am doing a march. Even with V and P discrepancies, I can't seem to avoid a negative bump. So I'm going to try a problem where the limits are the same, but the target and feasible shapes are different. Corrected an error in GETDU. NPARF was zero. 17 March 1995 I added a simple scaling scheme to FLOW3 today. I divide the momentum equations by REMAX, leaving the continuity and boundary conditions alone. 13 March 1995 It's so dissatisfying to see my cheapo Newton scheme being wiped out by calls to simple iteration, which force me to evaluate a different matrix, that I've considered the following two options: A) make two different matrix arrays. B) only call simple iteration if Newton fails. In that case, restart from original point, do simple, then Newton. Also, try SIMSOL's trick of saving the best point, rather than the first one. Also, get the damn condition number estimate, or the determinant. I'm worried about high Reynolds number flow, and about situations where the inflow is zero or very close to zero (finite difference). 11 March 1995 I also though, stupidly, that the simple iteration could be made more efficient by reusing the system matrix. But it's a linear system. So if you resolve it, you get the same thing. At least I know, but now I have to redo a bunch of initial results... One possible improvement to SIMSOL: save the point with the lowest F norm, and return that, whatever it is. This could be done simply by overwriting G. Yes, I got that to work, and I like it. 10 March 1995 I thought that, with a bit of work, the simple iteration and Newton iteration could share a matrix. But this cannot be done. 09 March 1995 Don't know what I can accomplish in the time between now and 21 April. Just restored the simple iteration code to FLOW3, to try to get some results at higher Reynolds numbers. The behavior, at least at RE=100, is vastly improved. 17 February 1995 I correctly deduced that I was not computing the sensitivity of the discrete equations, but rather the discretization of the sensitivity equations. I foolishly tried to compute the former; getting poor results, I realized that I needed not just the boundary integral (which I had), but the dependence of the basis functions on the nodes which in turn depend on alpha. At that point, I raised my hands and said, "I don't have time to pursue this now." Now I must be careful to distinguish the discretized sensitivities I am using from the sensitivities for the discrete problem. I tempted fate: I changed the code so that nodes on the fixed wall now generate equations for U=0 and V=0. The code seems to be working. Now all I have to do is handle the U, V inflow, the U, V outflow, and the bump, and every element will have a full set of coefficients. 13 February 1995 I want to go back to computing DelJ/DelBeta like I used to, as a good check of the partial derivative. The new variable JGRAD controls whether I use chain-rule-sensitivities, chain-rule-differences, or finite differences to compute the cost gradients. I'm glad I did that! I now have independent confirmation that my sensitivitity calculation is inaccurate. Running INPUT.614, I get the following: PrCSen: Cost gradients from the chain rule and sensitivities: -0.86784 -0.34021E-01 -0.10527E-03 PrCSen: Cost gradients from the chain rule and finite differences: -0.86784 -0.31578E-01 -0.10501E-03 PrCSen: Cost gradients from finite differences: -0.86784 -0.31578E-01 -0.10501E-03 10 February 1995 PLTWRT was not writing out the value of the sensitivity at boundary nodes, since it's not an unknown there. But we can compute it, it's just dUdAlpha = -dUdY*Shape. Of course, dUdY is not well defined at a node, since U is not continuously differentiable there. But, rather than not print anything, I compute DUDYN(I)=Value of dUdY at node I, averaged over all elements in which node I appears. 08 February 1995 So what's the result of all this? * Write something about how the finite element solution is not space differentiable. * Write something about how this affects the usefulness of the formula relating partials and totals. * Write up your results about the interpretation of FD and SEN. I discovered something I should have known. The finite element solution is not differentiable at the boundaries of an element. That is, quantities like dUdY have two or six values, depending on where you look. For instance, at a corner node, I got the following values for dUdY: GetDU DuDy(81)= 0.708792324314766 GetDU DuDy(81)= 0.673444964504792 GetDU DuDy(81)= 0.601876028295308 GetDU DuDy(81)= 0.673444964504792 GetDU DuDy(81)= 0.601876028295308 GetDU DuDy(81)= 0.678665874508058 What this means is that I had best do my work in element interiors, although I think this invalidates the formula du/dalpha = partial u/partial y partial y/partial alpha + partial u/partial alpha Well, no, but it means I can't compare these quantities at a node, say, but in the interior. That just eliminates the ambiguity in the meaning of dUdY but doesn't clear up the accuracy, I would say. I really want to drop the special treatment of variables at the boundary. As a signficant step to doing so, I will set up a character variable, of length NEQN, containing 'U', 'V', 'P', and later I can add 'UI' for U inflow, 'UW' for U on wall, 'UO' for U on outflow, and 'UB' for U on the bump. CHARACTER EQN(NEQN)*2 To try to convince myself that SEN is OK, (in which case maybe I'm just not getting an accurate value of GRADF), I am printing out some solution values at ALPHA=0.5, 0.51. To compare them at a fixed point, I have written some routines UPVAL and UPVALQ. 07 February 1995 OK, I have some confidence that GRAD correctly computes the expected change in a solution component G(I), which is presumed to be attached to a node that moves as the shape parameter changes. I have some evidence that the correction term is computed correctly, or at least corresponds to the values of GRADF=DUDY*DYDP. I have evidence that the values of SEN and GRAD-GRADF are still far off. Now I want to see whether SEN correctly estimates the change in U at a fixed point. OK, here's a simple discrepancy to explain. I'm looking at the finite differences and sensitivities computed at the target point, where ALPHA=0.5. (INPUT.710, NX=41, ISHAPB=1) IGRID, IFIX, MAX(Sen-FD) 0, 0 1.12, 0.16, 4.7 0, 1 0.213, 0.17, 5.5 1, 0 0.102, 0.13, 4.7 1, 1 0.102, 0.13, 5.5 So, my question: Which choice of IGRID and IFIX is appropriate? Why can't you compute the shape derivatives to great accuracy? First of all, when computing finite differences, you are always comparing at least some points incorrectly. Either all the points above the bump move, or at least the points ON the bump move. Your finite difference estimate subtracts these two values as though they applied to the same physical location, when they don't. FIXGRD attempts to correct this effect, by subtracting off the term dudy*dydp from the computation, for those nodes that move. Increasing the size of the finite difference perturbation EPSDIF by a factor of 100 didn't make any difference to the outcome of INPUT.710. So that's not the problem. I will also try ISHAPB=1 again. (INPUT.710, NX=21, IGRID=0, IFIX=1) ISHAPB=1 produced a much better agreement between SEN and FD at ALPHA=0: PrGS2, compute Maximum (Sensitivity-Finite Difference) Par Var Sen Index Node X Y Row Column 2 U 0.2549 184 106 2.000 0.2500 2 9 2 V 0.1250 136 81 1.500 0.5000 3 7 2 P 1.910 132 79 1.500 0.5268E-07 1 7 Par Var FD Index Node X Y Row Column 2 U 0.2545 184 106 2.000 0.2500 2 9 2 V 0.1249 136 81 1.500 0.5000 3 7 2 P 1.908 132 79 1.500 0.5268E-07 1 7 Par Var Max(Sen-FD) Index Node X Y Row Column 2 U 0.6181E-03 191 109 2.000 1.000 5 9 2 V 0.1033E-03 261 142 2.500 2.750 12 11 2 P 0.1960E-02 262 143 2.500 3.000 13 11 but somewhat poorer agreement at ALPHA=0.5: PrGS2, compute Maximum (Sensitivity-Finite Difference) Par Var Sen Index Node X Y Row Column 2 U 0.5662 184 106 2.000 0.7083 2 9 2 V 0.2522 139 82 1.500 0.9375 4 7 2 P 5.397 132 79 1.500 0.2500 1 7 Par Var FD Index Node X Y Row Column 2 U 0.7371 184 106 2.000 0.7083 2 9 2 V 0.2183 117 70 1.250 1.083 5 6 2 P 3.320 132 79 1.500 0.2500 1 7 Par Var Max(Sen-FD) Index Node X Y Row Column 2 U 0.1709 184 106 2.000 0.7083 2 9 2 V 0.1283 162 93 1.750 0.5938 2 8 2 P 2.077 132 79 1.500 0.2500 1 7 06 February 1995 Notice that the good agreement in problem 710 occurs for ALPHA=0, where the bump is flat. Is this something I need to factor in? I will run 710 with a finer mesh, and if the agreement at ALPHA=0 improves correctly, then I will work on the hypothesis that I need to deal with nonzero bumps... PROBLEM 710 RESULTS FOR NX=21 Taylor - Get base point solution. Inflow 0.500000 Bump 0.000000E+00 REYNLD 1.00000 Cost at Newton starting point 0.430 Cost functional 0.921E-02 PrGS2, compute Maximum (Sensitivity-Finite Difference) Par Var Sen Index Node X Y Row Column 2 U 0.3034 184 106 2.000 0.2500 2 9 2 V 0.1539 113 68 1.250 0.5000 3 6 2 P 1.929 81 53 1.000 0.0000E+00 1 5 Par Var FD Index Node X Y Row Column 2 U 0.3033 184 106 2.000 0.2500 2 9 2 V 0.1529 113 68 1.250 0.5000 3 6 2 P 1.907 81 53 1.000 0.0000E+00 1 5 Par Var Max(Sen-FD) Index Node X Y Row Column 2 U 0.7482E-03 286 158 3.000 0.2500 2 13 2 V 0.1314E-02 83 54 1.000 0.2500 2 5 2 P 0.2213E-01 81 53 1.000 0.0000E+00 1 5 PROBLEM 710 RESULTS FOR NX=41 Taylor - Get base point solution. Inflow 0.500000 Bump 0.000000E+00 REYNLD 1.00000 Cost at Newton starting point 0.428 Cost functional 0.974E-02 PrGS2, compute Maximum (Sensitivity-Finite Difference) Par Var Sen Index Node X Y Row Column 2 U 0.4604 796 402 2.000 0.1250 2 17 2 V 0.1571 489 255 1.250 0.5000 5 11 2 P 2.494 480 251 1.250 0.4610E-07 1 11 Par Var FD Index Node X Y Row Column 2 U 0.4604 796 402 2.000 0.1250 2 17 2 V 0.1571 489 255 1.250 0.5000 5 11 2 P 2.495 480 251 1.250 0.4610E-07 1 11 Par Var Max(Sen-FD) Index Node X Y Row Column 2 U 0.9700E-04 1216 602 3.000 0.1250 2 25 2 V 0.1640E-03 377 202 1.000 0.1250 2 9 2 P 0.5458E-02 375 201 1.000 0.0000E+00 1 9 05 February 1995 I am back to despair about the poor agreement between the finite difference estimates and the sensitivities for shape parameters. I ran a sequence of problems with NX=21, 41, 61, and did not get reasonable behavior of the discrepancies. In particular, for IGRID=0: NX=21 Sen FD Diff U 1.047 1.252 0.7690 V 0.337 0.286 0.3162 P 1.077 0.611 0.9171 NX=41 Sen FD Diff U 1.321 1.412 0.8278 V 0.384 0.314 0.4043 P 2.159 0.622 1.971 NX=61 Sen FD Diff U ? V ? P ? for IGRID=1 NX=21 Sen FD Diff U 1.047 1.061 0.3867 V 0.337 0.249 0.1580 P 1.077 0.681 0.3960 NX=41 Sen FD Diff U 1.321 1.332 0.3380 V 0.384 0.310 0.1734 P 2.159 1.306 0.8531 NX=61 Sen FD Diff U 1.457 1.459 0.1976 V 0.351 0.352 0.1546 P 2.293 1.396 0.8977 I should have used IGRID=1 the whole time. Nonetheless, I see that I'm doing a terrible job, assuming that Sen and FD must agree. Strange about the big difference in the pressures! Here's my latest clue. Switching to IBUMP=0, using non-isoparametric elements, I get exactly zero sensitivities. This shouldn't be. 04 February 1995 1) You should compute partial J/partial beta(i) in GETGRD as a check. 2) You should have more grid near the bump. Things happen there. Simplifications: dropped routines BLEND, GRID2N, GRID4N, and variables SXY, TXY and LONG. 30 January 1995 I am taking a look at a crude one shot optimization which replaces Problem I: Find PARA which determines U,V,P satisfying F(X)=0 which minimize COST. Solution: Unconstrained minimization of COST as a function of PARA. by Problem II. Find X=(U,V,P,PARA) so that F(U,V,P)=0 and dCost/dPara=0 Solution: Nonlinear equation solver. The main difficulty seems to be computing the jacobian of the enlarged system. I have started a mockup of this routine, which will correspond to option ITYPE=7. 25 January 1995 I got into a little problem with the varying Newton tolerance. When it's set real loose, then the gradient points were getting accepted immediately. But that made the gradient identically zero. So I had to go into NEWTON and force at least one Newton step for gradient calculations. 24 January 1995 Still studying efficiency. Setting up ITOL, TOLMIN, TOLMAX, TOLFAC, so that Newton tolerance is set by TOLNEW=TOLFAC*CGRAD where CGRAD is the max norm of the cost gradient. 23 January 1995 I want to add another jacobian option, an extremely miserly option that waits until Newton fails before evaluating the jacobian. I have set this up as JJAC. I am rerunning INPUT.717 to see if IJAC=1 corresponds to my previous results. If so, I will then run INPUT.720, with IJAC=2 (miserly option). 22 January 1995 I'm trying to come up with a reasonable plan of action. My first idea is to look at questions of efficiency. I want to compare using gradients versus sensitivities on the standard problem. So I added IPRED=3, 4, so that ONLY one of the finite difference and sensitivity estimates will be made, so that I can monitor the work involved. Runs 717, 718. My calculations are coming out OK, except now I want to pass IOPT to GETGRD, so that variables we aren't varying won't add to our cost. Got that OK. Now I noticed that I don't evaluate the jacobian if I'm doing a gradient. For a fair first comparison, I want to do that, so I'm temporarily turning that off. 19 January 1995 A first good sign. I finally got the correction to the gradients working properly. Now, for the base problem, I had a large difference in the gradients for U: Sen FD Max(Sen-FD) U 0.3169 0.2661 0.5087 U 0.3169 0.3154 0.54E-2 The V and P gradients weren't changed, but that's understandable. Their biggest differences occurred at X=3, just outside the range of influence of the parameter. But the important thing is that I got a correction to the finite difference gradients, and it seems to improve the agreement with the sensitivities. I found a possible problem in REFBSP, which might improve the matching of pressure sensitivities. I had the roles of XSI and ETA reversed, which affected the returned values of the basis function Q and its derivatives DQDX and DQDY. I'm getting much improved data for problem 710: PrCost: 0.921E-02 PrGS, compute Maximum (Sensitivity-Finite Difference) PAR Var Sen FD Max(Sen-FD) Index Node X Y 2 U 0.3034 0.3033 0.7483E-03 286 158 3.0 0.2500 2 V 0.1539 0.1529 0.1314E-02 83 54 1.0 0.2500 2 P 1.929 1.907 0.2213E-01 81 53 1.0 0.0 PrTayl: Taylor analysis Varying parameter 2 Stepsize h= 1.000000000000000E-002 Par Var G Gnew-G G+h*FD-Gnew G+h*Sen-Gnew 2 U 0.5000 0.2670E-02 0.5088E-02 0.5090E-02 2 V 0.3613E-14 0.1538E-02 0.1461E-04 0.1335E-04 2 P 4.444 0.1923E-01 0.1603E-03 0.1892E-03 2 J 0.9211E-02 0.2208E-03 0.1861E-05 0.2085E-05 2 U/h 0.5000 0.2670 0.5088 0.5090 2 V/h 0.3613E-14 0.1538 0.1461E-02 0.1335E-02 2 P/h 4.444 1.923 0.1603E-01 0.1892E-01 2 J/h 0.9211E-02 0.2208E-01 0.1861E-03 0.2085E-03 Just to summarize, here's how the IGRID, IFIX situation stands: Prob IGRID IFIX Udif Vdif Pdif #709 0 0 0.509 0.001 0.02 #706 1 0 0.005 0.002 0.03 #710 0 1 0.0007 0.001 0.2 #711 1 1 0.005 0.002 0.3 The conclusion seems to be to use IGRID=0, IFIX=1. 18 January 1995 Made a stab at correcting the gradients, by computing terms like Grad(inode,ipar) = Grad(inode,ipar) - du(i)/dy * dy/dpar(ipar) However, I don't see how it can help in the case IGRID=1, which is where I had my most success in bringing down the errors. And I don't even expect it to work for IGRID=0. IGRID IFIX 706: 1 0 709 0 0 710 0 1 711 1 1 Printing out the values of DUDY, DVDY, DPDY, I believe them. The last two are essentially zero, and DUDY is: y 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 dudy 0.55 0.44 0.33 0.22 0.11 0.00 -0.11 -0.22 -0.33 -0.44 -0.55 which accords with the exact solution: u(x,y) = 0.5 * (4/9) * y*(3-y) dudy(x,y) = 0.5 * (4/9) * (3-2y) Now I'll look at dy/d alpha. At least along the bump, I'm getting good results: x 1.25 1.50 1.75 2.00 2.25 2.50 2.75 dydp 0.43 0.75 0.93 1.00 0.93 0.75 0.43 Distressingly, the four runs I made seem to show that for both IGRID=0 and IGRID=1, IFIX=1 is the same as IFIX=0. That surely can't be. I will check it out. 17 January 1995 After a fair amount of work, I have gotten INPUT.706 to the point where, for a bump parameter of zero, the sensitivities and gradients agree something like this: PAR, Var, Sen, FD, Max(Sen-FD), Index, Node, X, Y 2 U 0.3169 0.3141 0.7391E-02 265 146 2.75 0.5000 2 V 0.1828 0.1848 0.1526E-01 287 158 3.00 0.2500 2 P 2.409 2.468 0.1654 285 157 3.00 0.0000E+00 which is perhaps the best they've ever been. So I think I was on the right track with the correction in BMPSEN. Here is the Taylor analysis for the same solution: Varying parameter 2 Stepsize h= 1.000000000000000E-002 Var G Gnew-G G+h*FD-Gnew G+h*Sen-Gnew 2 U 0.5000 0.2657E-02 0.5077E-02 0.5080E-02 2 V 0.4532E-14 0.1912E-02 0.7294E-04 0.1285E-03 2 P 4.444 0.2508E-01 0.7535E-03 0.1219E-02 2 J 0.9316E-02 0.2309E-03 0.6202E-06 0.2178E-05 2 U/h 0.5000 0.2657 0.5077 0.5080 2 V/h 0.4532E-14 0.1912 0.7294E-02 0.1285E-01 2 P/h 4.444 2.508 0.7535E-01 0.1219 2 J/h 0.9316E-02 0.2309E-01 0.6202E-04 0.2178E-03 Which shows that U is great, though V and P are not. Meanwhile, here's the behavior of the finite difference gradients when plugged into the equations defining the sensitivities: PrFx2: Base Point Bump Sensitivity Residuals using Gradients MaxNorm of residual is 3.223062372448425E-002 at index 286 Var Eqn Node Solution Residual *U 133 80 -0.240341 0.233246E-01 *U 135 81 -0.825379E-01 -0.115740E-01 *U 184 106 -0.311250 0.480582E-02 *U 186 107 -0.116932 -0.231460E-02 *U 235 132 -0.234240 -0.137126E-01 *U 237 133 -0.118208 0.694462E-02 *U 286 158 -0.169894 -0.322306E-01 *U 288 159 -0.105319 0.162037E-01 However, in the same problem, the sensitivities and finite differences begin to pull apart for other points, such as the target point itself: PAR, Var, Sen, FD, Max(Sen-FD), Index, Node, X, Y 2 U 0.6172 0.6498 0.1279 267 147 2.750 0.9141 2 V 0.2575 0.2420 0.9355E-01 268 147 2.750 0.9141 2 P 4.582 3.318 1.648 239 133 2.500 0.8125 Because I believe that the shape sensitivity equations can be solved exactly for the Poiseuille problem, I am going to sit down and work out the form of the right hand side. I may be wrong in my assumption, but I think the easiest way of understanding what's going on is to try to explain why, for this problem, the finite difference gradients don't satisfy the sensitivity equations. It's always possible that I'm forgetting the fact the the nodes move. 16 January 1995 The results for IGRID=1 are much better. I still, however, have a discrepancy. I am going to look at a similar problem, but with the bump parameter equal to zero, which should make analysis easier. This will be problem INPUT.706. Also, I have added IGRID=2, which absolutely refuses to update internal node coordinates after they are set for the target or first point. I'm glad I turned to looking at the problem with the bump parameter equal to zero. I discovered that, at the point X=1.25, Y=0.25, which is a quadrature point, as well as a grid point, UVAL2 returns the value U=0.0763, P=1.94, whereas PRSOL shows that U=0.152778 and we would expect by interpolation that P=3.9 or so. I suspect an error in the routine UVAL2, or else in the evaluation of the PHI's, so I wrote PRBAS to find out what's happening. What a truckload of worms! I have verified that, for element 26, at quadrature point 1, which is the average of nodes 1 and 2, hence equal to node 4, the quadratic basis functions are all zero, except for the fourth one, which is exactly 1. Now I am rerunning the problem, and having UVAL2 print out the individual pieces of the computation of U. I am hoping that I set it up correctly, so that it analyzes exactly the call to UVAL2 that I have been studying so far. 13 January 1995 The first bit of progress: I am now calculating in XQUAD, YQUAD the physical locations of the quadrature points, even for isoparametric elements. I got my first confirmation that I am doing this correctly today. HOwever, my basic error remains. If I compute DUDY two different ways, I get different results. One way uses the fact that I store the values of DUDY at all quadrature points. For NQUAD=3, the midside nodes are also quadrature points, so I have one value this way. The other method is to evaluate from scratch. Something is wrong because these two methods should agree. Only then can I expect my sensitivity results to improve. I have reached a puzzling point where I have made some corrections, and now, oddly, I get very good agreement in the dJ/dX term, comparing sensitivities to finite differences, but the dU/dX term is really bad. I shall have to think some more. 12 January 1995 I am convinced that I have found the problem in BMPSEN. The whole calculation boils down to evaluating the term -dudy*shape but I was evaluating dudy at a quadrature point, rather than at the appropriate boundary point. I have tried to replace this code by a corrected version, but I am not getting better answers. Now I want to make sure that the dudy that I compute is correct, or at least, agrees with the value computed by UVAL2. Noticed that arguments in REFBSP were ...,ETA,XSI) but REFBSP was being called with ...,XQ,YQ) which is incorrect. Turned them around. Added ETAQ, XSIQ to SETQUD and SETBAS, in attempt to cut down on the ambiguity of global versus reference coordinates. Added ETAN, XSIN which are the XSI and ETA coordinates of the nodes. 11 January 1995 I am getting decent results from TAYLOR, although there is one discrepancy I can't explain. For the inflow and RE parameters, the sensitivities and finite differences give the same results. For the shape parameter, however, the finite differences predict the value of the state variables better than the sensitivities do. But for the cost function, it seems as though the sensitivities do a better job. I am hoping this is a programming error, since otherwise it flies in the face of common sense. If I still can't dispel it, I will move the test point around a bit and see if I can fix it that way. I ran INPUT.704, to concentrate on a Taylor analysis of the shape parameters. Interestingly enough, I am getting almost the exact same residual for h=0.1, 0.01, and 0.0001: Max( (Unew-Uold)/h - Usen(old) ) = 1.011 This suggests, as usual, that I am missing some term in the definition of the sensitivities. Check #1: Repeat the analysis for ISHAPB=2, so that any problems caused by incorrect spline evaluation will be avoided. This will be INPUT.705. Check #2: Rewrite from scratch the routine BMPSEN. 10 January 1995 I am looking at INPUT.703. I want to see if I can get some nice Taylor results, showing that I can predict U, V, P and J based on the derivatives. Trying to make DPARA1 and DPARA2 variables, the cost derivatives using sensitivities and finite differences, rather than using DPARA. This necessitates some more monkey business in DISDER, and passing DPARA1 and DPARA2 to QSOLVE, GETTAR, and so on, which I haven't done yet. THis is mainly so that I can compare the two methods of getting cost derivatives in TAYLOR. 16 December 1994 I want to be prepared to start using Taylor approximations to the flow, that is, to use the sensitivities to compute an estimate like: U(Para+DPara) = U(Para) + Sum(i) dU/dPara(i) * dPara(i) Of course, there is a whole underlying problem of what to do when the parameter being varied is a shape parameter. Particularly for large parameter increments, I have to think this through. I have written TAYLOR and PRTAYL, which compute and display the results of using sensitivites and finite differences. Problem INPUT.703 is my test input for this routine. I will have to continue this work in the new year. There's something wrong with the continuation option. If I try to compute the target point partar(1)=0.5 partar(2)=0.375 partar(3)=0.5 partar(4)=0.375 partar(5)=100.0 I fail, but if I set partar(5)=50.0, I can do it. Moreover, I can then recompute with partar(5)=100.0. And yet, SOLCON cannot manage to get to partar(5)=100. This suggests that you could get better performance from SOLCON if you could track down the problem. Hey, what's worse, INPUT.531 now fails, and it ran a few days ago. I will say that, for the benefit of ITYPE=6, I did just force the code to do both gradients and sensitivities, and I wasn't careful about making that change... 15 December 1994 Working on some Taylor routines for FLOW3. 13 December 1994 I got the results back from INPUT.641, which was a 2D marching run for a 41 by 13 problem, to compare to INPUT.531. The gradient field was much nicer. I made some nice plots for the thesis. However, it would be nice to refer to the theorem of Carter that Jeff mentioned, where inaccurate gradient data is acceptable as long as it stays roughly in the right direction. 12 December 1994 I haven't been running the program much, lately, because I've been concentrating on writing up the paper for Bozeman, and writing up chapters of the thesis. I have just come back to the program because I'm trying to specify exactly which input files were used to create the graphs. I'm at a momentary standstill with the 2D sensitivity versus finite difference comparison graphs; I can't seem to get the optimizer to halt at the place I claim. Also, I'm having the usual difficulties at REYNLD=100. I've stuck in a little damping, and I'll also see if I can't do some continuation on RE... I am trying to verify that refining the mesh improves the approximating power of the sensitivities, so I'm running a 41 by 13 problem (at Reynld=100), and waiting a long time for a response from the Alpha. 11 November 1994 I roughed out two routines, BMPDIF and SETBAS3, which can be used to estimate the derivatives like D2XSI/DXDALPHA. 10 November 1994 Maybe I need to think about what I want to do. I've lost track of the direction Max sent me in. Also, why is FPRIME called 60 times when GETDER and others are called 30? I fixed it. It's because IJAC=1, and we're taking 2 Newton steps, on average. Modify ISOTRI: 0, triangle is not isoparametric, and nodes never move. 1, triangle is not isoparametric, but nodes may move. 2, triangle is isoparametric. Then modify SETBAS so that only triangles in classes 1 and 2 are recomputed. Well, no thanks to the stupid debuggers on the ALPHA and SUN, I think I tracked down a problem that has been causing errors for some time. The main program did not dimension PHI properly. I discovered this from some bedeviled basis calculations. I was sure that I only had to update the basis functions for the moving elements, but when I tried to do this, I was getting nonsensical results. The modifications to SETBAS are working. I evaluate all the basis functions once, and then update the basis functions for elements of type 1 or 2. The resulting code only calls TRANS 27360 times, versus 133,920. The resulting code ran in 29 CPU seconds, versus 40. Even better results could be expected if we could run IBUMP=1. 09 November 1994 Some change that I made to get 640 to run (shifting the solution from one grid to another), has resulted in 639 no longer running. I'm getting a singular jacobian now, on the Sun. But I should be glad, because on the ALPHA I just get a floating divide by zero, and no traceback no matter what I try to do to get one. The SGI's not available, or I could try it there. I give up on my chances of tracking down the error. I'm going to cache the current copy, retrieve flow3.11.07 and try rebuilding what I have now. In particular, this means I have to delete NPARR, FLOVAL, FLOTAR and all that stuff, read IOPT and so on. Oh, for a decent debugger! The road to recovery: I've deleted all occurrences of NPARR. Code tested OK. Inserted IOPT code, GOPT augmented by DPARA, DPARA1. Code OK. Got rid of FLOPAR, FLOTAR, FLOVAL, FLOOLD, FLONEW. Code OK. Got rid of REYNLD, REYTAR, REYVAL, REYOLD, REYNEW. Code OK. As if I didn't have enough problems, the EDIT program on the SUN seems to suddenly have grown implicit tabs. If you have 17 blank spaces in your text, for instance, they will show up as 12, or 13, or some random but repeatable and totally wrong value. Just UNIX trying to be helpful and customizable, I'm sure. 08 November 1994 Added an input variable, ISHIFT, to control whether or not we try to adjust G after making a prediction. Inserted call to GSHIFT in FLOSOL. I'm stuck with this stupid ISHIFT option. The code is bombing. The ALPHA is horrible with no output, the SUn gives me a singular jacobian. 07 November 1994 I cleaned up the code, so that only one of SENS or GRAD is computed, and whichever is computed is used to computed the derivative GOPT. So of course it's not working now. OK, now it's working. Summary of 636, 637, 638: If I use good direction information (from finite differences) I converge properly in 33 steps. If I use poor information (from sensitivities) then I don't quite converge in 65 steps. Now I have succeeded in using sensitivities or finite difference gradients to make predictions, and it makes a difference. Now I am trying to simplify the use of parameters by having one list, in PARA, for the flow solver, and another, shorter list, in XOPT. A vector IFREE tells me which parameters are actually free to vary, and hence must be passed to the optimizer. This will eliminate all this tedious passing of combinations of scalars and vectors. 06 November 1994 I modified GETDER and DISDER so that GRAD is passed in, and DISDER uses GRAD if IPRED=2. 05 November 1994 Well, the odd thing is that the runs that use the sensitivities don't converge. Now I knew that, and I even count on that other times, but I guess I need to look at perhaps a different problem, so I can get them all to converge, and compare performance. Changing to IBUMP=1 seemed to improve the behavior of 633, so I'm rerunning 634 and 635 now. 635 doesn't converge as well. What's going on? 04 November 1994 OK, strategy. First thing to do is finish the routine that allows you to transfer a solution from one mesh to another, and test it out. Then you can use the sensitivities to compute a new solution, and transfer that solution to the new mesh. Will I be doing what Max wants if I write a routine that predicts the value of the solution, using sensitivities or gradients? Obviously, I would like to start at a point, and then examine the behavior of the flow solution residual along a tangent line. I would want to compare the sensitivities and gradients. I changed SOLCON to use IPRED to determine how to predict the next point. But who ensures that GRAD has a value? Wrote GETGRD2 (don't need the cost gradients. Maybe you could cut that down too?) Now it would be rational to run, say, problem 501 using IPRED=0, IPRED=1 and IPRED=2 and compare them, expecting that first of all GETGRD2 won't work properly, and SOLCON won't use its results properly. Oh, and by the way, you have all that phony code in there for the scalar cases, which will never work. Set up input.633, input.634, input.635, and submitted 633. Finally, I still think you should investigate computing the terms like D2 XSI DALPHA D X, and doing a moving mesh calculation. 03 November 1994 Max asked me to check that, for finer meshes, the approximation between the sensitivities and the finite differences improved for the case where we hold the internal mesh fixed. Results: Case 1: Adjusting the mesh: Run Grid Max(Sen-FD) Max(Sen) Max(FD) 618 11 x 4 .331 .101 .234 627 21 x 7 .481 .275 .267 629 41 x 13 .567 .442 .263 Case 2: Fixing the mesh: Run Grid Max(Sen-FD) Max(Sen) Max(FD) 626 11 x 4 .037 .101 .135 628 21 x 7 .027 .275 .302 630 41 x 13 .017 .442 .460 This suggests that, for finer meshes, the approximation improves, although it's still not great, and the error doesn't drop in a nice way. Modified SETXY so that the "inbetween" columns of interior nodes were given the average Y coordinate of the left and right neighbor nodes. This is required so that when IBUMP=1, we can indeed make most of the elements above the bump nonisoparametric. Once I am convinced this works, I will start running IBUMP=1 runs instead of IBUMP=2. As a minor help to sanity, I am going to eliminate the NNODES variable. 01 November 1994 The results of comparing runs 618 and 626 tell me the following: I thought the sensitivities were very bad, because 618 showed a large discrepancy between the gradient and sensitivity vectors. But 626 showed me that if I held the internal mesh still while computing the sensitivities, that I got very good agreement between them. For some numbers, note that: 618: U-Sensitivity had max value 0.105, U-Finite diff max value 0.234, and maximum difference was 0.3331. 626: U-Sensitivity had max value 0.105, U-Finite diff max value 0.135, and maximum difference was 0.0379. So I have some confidence that the shape sensitivity calculations are correct, at least to within 10%, as long as I compare them to gradients on a fixed internal mesh. Now of course I also plan to use sensitivities to estimate solutions at nearby parameter values, and for that, I now realize that I must compute the sensitivities, and then compute the new mesh, and interpolate values between the two finite element meshes. In order to ease interpolation, I set up a routine GETXSI which, given X and Y and an element number IELEM can compute the XSI, ETA coordinates of that point. Also, I am replacing the variable IT by IELEM. 31 October 1994 In my thesis discussion of the spurious local minimizer, I claimed that the gutters get twice as deep when the mesh is refined. This result was obtained when I was computing with sensitivities, which I know are much less accurate for shape parameters. I recomputed these results using finite difference estimates, and can still claim that the local minimizer is spurious; but I'm not getting nice convergence to a point yet, and the gutters are VERY deep. Max wants me to consider the case where I don't shift the internal nodes as the bump moves, at least when it moves only slightly. This is just to see whether the computation of the shape sensitivities comes closer to the finite difference result. How about IGRID? I'm trying to set this up with INPUT.626, to be compared to INPUT.618. Max also said I should probably run IBUMP=1, that is, use isoparametric elements only just above the bump. 27 October 1994 I've had a fairly miserable time trying to figure out what the appropriate equations are for the shape sensitivities. My first problem is that I am shifting the mesh, so it's natural for me to want to do a lot of calculations in the "master element". These calculations are, in a sense, like the "material derivatives" that you see in fluid calculations. Max's approach is, instead, like Cabuk's, in trying to fix a spatial point in the physical region and derive the appropriate boundary conditions. But I am having trouble believing them. Here is my attempt to work on the continuity equation. Consider a physical point (x,y) which corresponds to the abstract point (xsi,eta), for a given value of alpha. Note that this correspondence depends on alpha, so that we may write x(xsi,eta,alpha)=a1*xsi^2+b1*xsi*eta+c1*eta^2+d1*xsi+e1*eta+f1 y(xsi,eta,alpha)=a2*xsi^2+b2*xsi*eta+c2*eta^2+d2*xsi+e2*eta+f2 where the coefficients a1 through f2 are functions of the x, y coordinates of the images of the nodes of the master element, and as such, depend on alpha. Now consider the continuity equation: dUdx + dVdy = 0 Rewrite this in the master element: dUdxsi dxsidx + dUdeta detadx + dVdxsi dxsidy + dVdeta detady = 0 Now take the partial with respect to alpha, and write u, v for the U and V sensitivities: dudxsi dxsidx + dudeta detadx + dvdxsi dxsidy + dvdeta detady + dUdxsi d2xsidxdalpha + dUdeta d2etadxdalpha + dVdxsi d2xsidydalpha + dVdeta d2etadydalpha = 0 I am now faced with the problem of computing quantities like d2xsidxdalpha. I talked with Max, who at least originally believed that this problem could be handled by computing source terms on the bump (for his fixed spatial point approach), but now he agrees that there may be something more required in his approach, and he encouraged me to continue looking at my approach. More later. 22 October 1994 Now I'm getting ready to finish BMPFX, by first calling BMPSEN to do the right hand side. If that works, then I can proceed to try to nail down the errors. Also, on the trip back from Tullahoma, I got some insight into a better formulation of the boundary conditions for the bump problem, and an explanation of why the boundary condition for u is du/dy dy/deps deps + du/deps = 0 20 October 1994 There must be something wrong with my "BMPFX" routine. If I take the computed sensitivities and multiply them by the matrix, then (at least in an element with zero right hand side) I get 0 for the residual. Well, 1.0E-12, anyway. So toss the old version of BMPFX, and start with FPRIME and build on it. OK, I have BMPFX, which now reports very tiny residuals for the sensitivities and the gradients, at least away from the bump. I haven't inserted the RHS computation into BMPFX yet, but that's next. Also, I should examine the residuals for GRAD and SENS and see whether GRAD behaves as though there's a nonzero RHS away from the bump. I would hope not... Once I have confidence in my sensitivities, I can move on. 18 October 1994 I tracked down the error, which was a very tiny transposition of two arguments to UVAL2. Now I will test out the new quadrature rules. 7 nodes seems to give similar results to 3, but 4 is distinct, and perhaps off, either because of an error I made, or because one of the weights is negative? I have to explain why, on problem 615, I have an enormous discrepancy in the SHAPE sensitivities. Then I will proceed to investigate the use of the sensitivities for new purposes. 17 October 1994 Still trying to track down a discrepancy. I have "AVAR", the factored jacobian lying around. I call FLOSEN to set up the right hand side, which I can represent as ABC*BC, for flow sensitivities, and solve AVAR*SENS=-ABC*BC Then I set up a residual routine that computes AVAR*SENS+ABC*BC and it does not get zero. The discrepancy shows up in the first column of elements, in the horizontal momentum equations. I have to resolve this, so that I know I have accurate sensitivities, so I can go on and do the stuff Max wants. By the way, I am surprised at the low resolving power of my quadrature rule, and want to have an option to increase it. But to do so, I must add a WQUAD array, since I was implicitly assuming that the weights were equal. 13 October 1994 I revised UBUMP to include a nonzero term dPdBMP. Now I have to use it. But I find it will be difficult to do so. BMPSEN sets up the right hand side, but the left hand side is already set, and hard to change, because the matrix A has been evaluated and factored. So we would have to do some delicate adjusting. But no, even that won't work. If the pressure is a boundary condition, then we drop the continuity equation, replacing it by a specification P=VALUE. Could I do something equivalent? 12 October 1994 Set up separate array GRAD to hold the gradient. I think I'd like to write this to the plotfile as well. Noticed that I still can't write the sensitivity values along the bump to the plot file. Verifying that GETGRD can live with calling FLOSOL rather than SOLCON (since we have a solution at an extremely close neighbor!). 09 October 1994 Modified SOLCON to include the cases where REYVAL or FLOVAL were not supposed to be parameters. (Haven't tested this, though). My major concern right now is that I know that the bump sensitivity calculations are wrong, and I don't know how to derive the defining equations so that I can write them in my thesis, and code them in my program. 07 October 1994 Got a new version of SOLCON going, which corrects most of the insufficiencies of the old version. And soon I will add the case where the scalar values of inflow parameter or Reynolds number change, but not right now. Folded March1, 2 and 3 into one marching code called MARCH. I dropped OSOLVE and PSOLVE in favor of QSOLVE. I noticed that QSOLVE didn't do as well on #501 as OSOLVE had, but I don't care right now. 06 October 1994 I wrote up a new version of SOLCON. The old version only continued on the REYNLD parameter, which was fine if we were trying to evaluate a solution at a high Reynolds number, and wanted to use continuation to get out there. But it had no capability for using an old solution plus sensitivities to compute the new solution. Even though I was actually using it as though it had that capability! The current, new version can handle changes in the PARA vector, going from PAROLD to PARNEW and using sensitivities to do so. I want to include the ability to move in the FLOVAL and REYVAL directions as well, so that the new SOLCON can: *) Start from a 0 solution and compute a solution at some set of values of FLOVAL, REYVAL and PARA. *) Move from FLOOLD, REYOLD, PAROLD to FLONEW, REYNEW, PARNEW. *) In particular, move from FLOVAL, REYOLD, PARA to FLOVAL, REYNLD, PARA. (That is, REYNLD continuation, even if REYNLD is not a parameter). 05 October 1994 I got QSOLVE working. It's tempting to consider getting rid of both OSOLVE and PSOLVE, but perhaps instead I should merely retire them. What I want to work on now is improving the prediction of the next point, using IPRED to control this behavior. I'm now passing SENS to SOLCON, and I need to make sure that it is evaluated. (What do MARCH1, 2 and 3 do, for instance?) 04 October 1994 I'm trying to set up a routine QSOLVE which will compute derivatives using finite differences or sensitivities, on my command. I want to do this so that I am in control of the finite difference calculations; I suspect that sometimes X and GRAD(X) do not correspond, because of how 611 tries to economize. But I'm having some problems getting the code to run right now. 30 September 1994 Max gave me four areas to pursue: 1) Calculate u(p+dp)=u(p)+du/dp * dp 2) Fixed jacobian 3) Use coarse grid, switch to fine grid. 4) Only update solution where sensitivities are high. IE, change the domain. 29 September 1994 Split ISHAPE into ISHAPB and ISHAPF, so that bump and flow may be modelled with different basis functions. Also want to add choice ISHAPB=0, allowing smooth polynomial interpolation. Is there a problem computing the basis functions. 24 September 1994 Here's a comparison between #600 and #601. #601 allowed up to 20 steps at grids of 11 by 3, 21 by 7 and 41 by 13: #600 #601 ---- ---- 33 30 opt steps 119 108 calls to Solcon 0 0 continuation steps by Solcon 0 0 calls to Tangnt 119 108 calls to Newton 161 151 Newton iteration steps 280 259 function evaluations. 85 75 finite difference gradient estimations 76 76 LU factors 164 160 calls to the LU solver 164 160 right hand sides for LU solver 372 61 CPU seconds Now a more interesting comparison would be one where the coarse grid approximation was pretty poor, for one of two reasons: high Reynolds number, or non-attainable minimum. Let's go to REYNLD=100 first, for runs #602 and #603. Then consider doing a non-attainable minimum. Drat! At REYNLD=100, I can't get convergence with 11 or 21 grids. Let me think what I want to do, then. Can I run a 61 grid? Let me set that up as #603... 23 September 1994 As a first stab at Max's new questions, I suggest comparing the following two runs: 3 parameter bump problem, started from 0, RE=1. #600: Run 1 uses a 41 by 13 grid all the way. #601: Run 2 uses 11 by 3, 21 by 7 grid, then 41 by 13. A mistake in SOLCON turned up today: if NPARF was zero, then a bug in programming meant that 5 continuation steps were always taken. Now I have to rerun #600 and #601. Just so I can compare, here are the stats for #600 with the SOLCON bug: Bug Fixed ---- ----- 32 33 opt steps 117 119 calls to Solcon 702 0 continuation steps by Solcon 702 0 calls to Tangnt 819 119 calls to Newton 1504 161 Newton iteration steps 2323 280 function evaluations. 84 85 finite difference gradient estimations 427 76 LU factors 2209 164 calls to the LU solver 2209 164 right hand sides for LU solver 2400 372 CPU seconds 22 September 1994 The 578 run came back, and the behavior is a little better. The bump has a peak in the front, and a slight dip before the shoulder at the back. After showing the results to Max, he suggested that he'd seen enough of this line of study, and that we should get together in the next couple of days and talk about a new line of inquiry, which might concentrate on cutting down the amount of work involved in the optimization, including: * using a coarser grid at first, * using higher Newton tolerances, * Taking only one or two Newton steps * using only an Euler approximation to the flow. Question: Do I use an Euler approximation to the Newton starting point now? I don't think so... Also, while reading Max's book, I realized that when I rewrote the FX routine, I used the formulation Integral (dP/dx * dWi) dOmega whereas the standard finite element formulation replaces this by Integral (-P*dWi/dx) dOmega + boundary terms. The boundary terms have the form P*Wi. How would I handle them? 21 September 1994 I found that in the sensitivity calculation I was computing sensitivities for the REYNLDS parameter even if I wasn't using REYNLD as a parameter. Perhaps this is what went wrong on the ALPHA, but I expect there's more... Well, I found it, no thanks to debugging tools on the ALPHA! My Newton iterates occasionally blow up, and so I simply need to check whether the residuals are getting too large, and abort. Now #578 seems to be running along OK. 20 September 1994 I talked with Max about the results of the runs 568 through 576, and he suggested * considering adding the cost of the V discrepancy and seeing if that removed the concavity of the solution bump * trying to enforce bump convexity in some way, by explicit constraints, or rejection. * graph the behavior of the U discrepancy functional as a function of the bump wiggle penalty coefficient. The overview of the runs that I gave Max is as follows: The following runs were all made with a three parameter bump, the target bump running from 3 to 6, the feasible bump from 5 to 6. #568: Starting at 0, minimizing bump has an "S" shape. #572: Same as 568, but add bump wiggle penalty coefficient of 1/1000. Minimizing bump has nice "bluff" shape. #573: Start at solution of #572, but no bump wiggle penalty. Minimizing bump has "M" shape, and is lower in cost than minimizer found by #568. #574: March through S and M bumps, exhibit a "W" shape for the cost function along that line. #575: 2D march in plane of S and M bumps, and solution of 572. Exhibits "hill" between S and M. #576: Sequence of bump wiggle penalties, 1/1000, 1/10,000, 1/100,000, 0. Got 573 solution to 6 digits. While trying to run #578 (adding the V discrepancy) I got a floating point overflow on the ALPHA. I haven't tracked this problem down yet. 19 September 1994 I broke INIT1 into three routines: BMPSPL, FLOSPL and PLTOPN. Before this, I essentially had two copies of the same code, for doing spline set ups for target and feasible solutions. Now that is unified. 15 September 1994 Ray Stell told me how to access the FORTRAN man page on the Sun. I have successfully implemented a scheme for cutting down on the number of jacobian evaluations. If we are using OSOLVE, then during the gradient estimation, which involves a point X-bar near the current point X, I don't re-evaluate the jacobian. This should help a lot. 14 September 1994 To make life easier, I added IMEMRY to FLOW3, which allows me to add new counters, update them, and access their values with ease. I also sent mail to Ray Stell asking where the F77 man page was. 13 September 1994 Very strange behavior on the SUN: No DATE function. No man page for f77. Message of missing routine "cg92_used_" when I tried to link to a compiled library (lib611d.a). Only when I compiled everything together at once could I link. Even then, I couldn't compile with "f77 -c *.f" and then load with "f77 *.o", NO, NO, I had to compile at once: "f77 *.f". What a crock! Well, let's go ahead and try and run something on this heap of horse ordure! The SUN code also crashed, but at least it gave me a warning about what was happening, apparently something involved in CHRCTI, when trying to patch together the name "PARA1(*)". Now I will look more closely at it. Hmmm, more likely problem in CHRDUB. Report: CHRDUB produced its value by using WRITE statements. For the IRIS and the SUN, this was catastrophic, if the value of CHRDUB was actually used immediately in a write statement: write(*,*)chrdub(x) since apparently this invoked the WRITE routine a second time. I have reworded FLOW3 temporarily. The correct conclusions: *) This is another piece of evidence that FUNCTIONS are bad. *) This is another piece of evidence that it's better to do things yourself rather than rely on a WRITE statement. *) Some implementations of FORTRAN are really poor once you try something fancy. 12 September 1994 Added draft code to NEWTON to handle jacobian option with retry. Just submitted 577, with IJAC=10, comparison with 501. It failed, of course, on the SGI, so I have to run it on the ALPHA. I used to be able to run 41 by 13 jobs on the SGI. Has my code gotten that much bigger? 576: Optimize on Cost = U + 0.001*bump, Cost = U + 0.0001*bump, Cost = U + 0.00001*bump, Cost = U 575: 2D march on plane of 568, 573, and 572 solutions. Taking a LONG time to run. 574: 1D march between 568 and 573 solutions. Shows that, at least along this line, there is a rise in the functional in between. 09 September 1994 Discussion of results for 568, 572, 573: 568: Cost depends on U. Solution bump is an S, with negative back. 572: Cost depends on U and bump slope. Solution bump is positive. 573: Cost depends on U, optimization starts from 572 solution. Solution is an M function, with negative part in middle! My conclusion is that we may have very weak dependence on the shape of the bump, though strong dependence on the height. Also, there may be local minima. Max's recommendations: A) Do a 1D march between 568 and 573 solutions, to verify that there seems to be a barrier between the solutions. B) Do a 2D march on the plane containing 568, 573, and 572. C) Do a series of optimizations with WATEB=0.001 on first step, but then decreasing to zero. See if solution goes to 573 solution, or 568, or yet somewhere else. 08 September 1994 The results from 573 came back, convergence to (0.35, -0.20, 0.20), a bump with a slump. The functional was 0.8E-04. I went to compare this with 568 and realized that 568 had not converged since I'd limited it to 30 steps. So I'm rerunning that. Modified FLOW so that IPLOT=1 prints target, point 1 and last point. This way, I can get an idea, in profile plots, of how much the optimization improved the initial functional value. 07 September 1994 The bump penalty (WATEB=0.001) in run 572 smoothed away the negative part of the bump that had been computed in 568. I have to ask Max what to do next for that problem. I inserted piecewise linear code (ISHAPE=1). Max looked at the results, and thought that I had done a run with decreasing values of WATEB, so that 568 would have represented another local minimum. I didn't, and I think 568 is the global minimum (just guessing that). Anyway, he'd like to clarify this point by doing 573: start from the 572 solution, but with EPS=0, and see where you go. Also, perhaps, do a 1 D march between the two solutions. Also, he wanted me to confirm that I have some runs at REYNLD=500, and if not, to make them. 06 September 1994 I inserted code to do piecewise quadratic polynomials rather than splines. I ran the test problem (as INPUT.570), and converged to the correct answer, encouraging. I am now going to rerun a problem (as INPUT.571) which got a concave bump using C1 piecewise cubics (INPUT.568). If the computed bump stays convex, I will be immoderately happy. Max looked at the 568 results, and suggested rerunning the problem with a wiggle penalty. I will do that as soon as 571 completes. Drat! 571 went to a bump described by (.52, -0.5, 0.35). It hadn't actually converged after 30 steps, so I resubmitted with 60, but I'm dubious of the results. I shall add in a piecewise linear option as a final, simplest form. 04 September 1994 After talking with Max, Janet said I should look at the problem where a long shape is modeled by a short one, hoping that the short one will become blunter. I will do so, but I expect I will need to make the bump relatively larger, and with more parameters (certainly more than 1!) and run at a higher Reynolds number. 03 September 1994 I spoke briefly to Janet today, about what problems I'm running. She suggested that the bump problems I'm running may not be the right things to look at; it may not be possible to match when the feasible and target bumps are so far apart, and that maybe instead I should look at comparing a short and a long bump. Also, what if I described the bump shape as a piecewise quadratic? That would be more compatible with the elements, and I could write the code for it myself, instead of having to rely on DeBoor. I still don't understand the not-a-knot condition when you only have three knots, for instance. 31 August 1994 I can't get FLOW3 to run on the IRIS. It's been a long time since I ran it there, but I'm not happy that it's failing now. 30 August 1994 I resurrected FLOW1, which had the old simple iteration code in it, but it couldn't run problem 558 either. It was only when I refined the grid to 41 by 13 that I got convergence, using the FLOW3 code. (Even then, I had some problems). 29 August 1994 I found that the graphics dumps were not being done correctly, except for the target and final points. So I corrected that. I am having troubles with convergence again, and at low Reynolds numbers, like 100. I am sure that this didn't used to happen, and I need to come up with a better strategy for handling it. 27 August 1994 I just spent quite a few days computing test cases for the paper for the Bozeman conference, and I'm done now. I gave copies to Max and Janet to look over. Meanwhile, I'll move on and try to work on the problem where the target is not feasible. 17 August 1994 I modified MARCH1 to update it to the standards of MARCH2. So it projects the derivative GOPT onto the marching direction and returns a line derivative. Also, it computes a finite difference estimate of that same derivative, and writes both to a file. I'm having a problem. I did a MARCH1 from ALPHA=0 to ALPHA=1, printing out COST, and the directional derivatives derived from the sensitivities, and from the finite difference derivatives. Essentially, the directional derivatives differ in sign. I went back and looked at runs 537 (FD) and 538 (SENS) and saw the same problem. 537 converged, and the directional derivatives were correctly negative from ALPHA=0 to ALPHA=0.5. 538 failed to converge to ALPHA=0.5. Instead, it proceeded as long as the directional derivative was negative (til about 0.015), proceeded a little further, as the functional decreased but the directional derivative become positive, and then sputtered out at ALPHA=0.03. So, one possibility, we need to look at the functional with XPROF=3, not XPROF=9. Or does the same problem happen there? No, the problem clears up and the sensitivities are good. 16 August 1994 In order to investigate the scaling problem for the two parameter case, I added an input quantity BMULT which multiplies the bump parameter. 15 August 1994 I added a routine to print out the derivatives of the state variables with respect to the parameters along the profile line, because my printout of the maximum entries of this array was not explaining why the inflow is so much more important than the bump to the functional. Locality is why, and that should show up if we restrict our attention to the profile. 14 August 1994 I added the marching file option to OSOLVE and PSOLVE, only there I write COST, (PARNEW(I),I=1,NPAR) to the file. I'm doing this partly so I can reproduce a graph of the cost function and the iterates in SLAB. 12 August 1994 The SENS vector in OSOLVE is meaningless. It is only carried around for the convenience of GETCST, which computes the COST and GOPT, the derivative of the cost, whether or not SENS makes sens. I propose that we split the cost and cost derivative calculations into GETCST and GETDER. Here I go! I also added a printout to PRSEN which prints out the magnitude of the cost influences. 26 July 1994 Set up my first problem with a nonfeasible target today, INPUT.535. 25 July 1994 In order to allow problems in which the inflow parameter is fixed, I added two input quantities, FLOPAR and FLOTAR. If NPARF is zero, then FLOPAR is the inflow parameter for the feasible solutions, and FLOTAR is the inflow parameter for the target. I still have some discrepancies to clear up, namely, adding MAXPARF to the call lists of SOLCON, UBUMP, UVAL and UVAL2, since I use NPARF+2 to dimension the spline coefficients, but this won't be right if NPARF is zero. Also, I have to rewrite PRSPLN because right now it assumes the declared and used dimensions of a spline array are the same. While running INPUT.534, I saw that the code was having a lot of convergence problems when the inflow was zero. I went into Newton and inserted a patch saying that if the inflow was zero, the first guess should be zero, and the problems vanished. I don't understand how this could be a problem, since the inflow only affects the right hand side of the Newton system... 24 July 1994 For problem 530, I got a strange region of two long thin valleys. I want to plot the flow field for the local minimum, so I'm trying to use MARCH1 to do it. I really should have written MARCH0, but there's so much overhead to maintaining one of those upper level programs! 22 July 1994 I am trying to deal with the bad points by using a signal to the optimizer to tell it I'm rejecting the point. We'll see if this works, and if it helps. If not, I may have to look at a damped Newton method. I was able to help my problems with INPUT.530 by backing away from the outermost points I was computing, which had weird values for the parameters, such as (.53, .64, -0.008, -0.6). This whole issue of wiggly bumps is just waiting to be addressed. I really would prefer to have a problem formulation where the bump stayed convex! While plotting the MARCH2 output of INPUT.531, I noticed that the gradients weren't perpendicular to the contour lines, a sure sign of error. It turned out that, for all these years, I have never taken into account the fact that the two axis can be of different lengths, so that the gradient vectors must be rescaled to show up properly. In other words, the original picture is in a rectangular frame. When I squash it into a square frame, I have to squash the derivative vectors correspondingly. 21 July 1994 I'm having problems, for example, with INPUT.528, because OSolve is getting requests from the optimizer for large steps. When it goes into Newton, the iteration fails. I've tried modifying SolCon to take 5 or 10 continuation steps at least, and I'm still having some problems. This has only occurred because REYNLD=100, but I can see that this can only get worse with larger REYNLD. At this point, I am trying to fix the problem by setting COST=10 when all else fails. 17 July 1994 I'm working on Max #9. I set up a run using sensitivities, and got "false convergence", which I assume to mean that the code is confused: the gradients are nonzero, but going in the gradient direction doesn't result in any further reduction. So I restarted the run, using finite difference gradients, and the code continued for a while, still not getting to the expected target solution. So I'm hoping the cost function contours for this problem are also interesting. I've submitted a run to compute them, but it's having a little trouble. 10 July 1994 Added new input parameter ICOST. If ICOST=1, then the cost weights are renormalized so that the cost of the zero solution is 1, and the component costs that are nonzero are equal. Setting up Max #9, the standard problem at REYNLD=100, I discovered that SOLCON really expects REYNLD to be a parameter. I have to decide whether to go back to always having REYNLD a parameter, and then maybe figure out how to tell the code it's "not in play" sometimes, or figure out how to have SOLCON work in either case. 09 July 1994 The new 521 sensitivities look quite good now. I made the sensitivity movie, although I really need to do some graphs comparing the sensitivities and derivatives directly. Also, for scaling, I just started computing the "cost" of the target bump, and I can then use that to rescale my costs so that each component has equal weight and scale, subject, of course, to user preferences in accord with WATEB and the other factors. 06 July 1994 I think I found the bug that was causing bump sensitivities to be off. The contributions to GOPT should be offset by NPARF entries, but this was not happening. I will recompile, and then resubmit INPUT.521, and see if it does a better job. If so, I need to go back and look at the earlier runs that depended on the penalty parameter. 05 July 1994 Dashed off a very quick 5 minute videotape for Janet to use in the Air Force presentation. 04 July 1994 I went back to computing 21 points on #521, but tried to move the special points in a bit, so that the shape of the banana would show up better. However, now I'm losing definition on the lowest contour levels, as I should have expected. Also, I'm not sure that one of the special points is actually the local minimum, although I've convinced myself that the point drawn in the upper right corner is the global minimum. I went ahead and submitted a job with 21 steps in the Z direction, so that I can go ahead tomorrow and try to make an animation of the sensitivity direction field. 03 July 1994 I was troubled, because the vector plots of the gradient field, or approximations thereof, for the cost, as done by SLAB, weren't crossing the level curves perpendicularly. In fact, they were going all over the place... I finally realized that I had not updated MARCH3 as I had MARCH2, so that instead of writing the projected 2 components of the gradient, and the 2 components of the projected finite difference gradient, I was still writing the NPAR components of the unprojected gradient. Since NPAR=4 for this problem, I wasn't noticing any missing data. I've just fixed that, and resubmitted #521. 01 July 1994 I spent some time writing up, in the files Max1.txt, ..., Max8.txt, my progress so far on the problems of 11 April. I also put the appropriate input, output, data and movie files in directories on the ICAM Mac, and on floppy disks at home. I am working on #8 right now, and rather surprised at how long it's taking to run. It's a 21 by 21 by 21 run, but that should only take 8 times as long as what I'm used to, on ALPHA1. And yet it's taken all day, it seems. It's almost done. I'm puzzled. The gradient vectors I'm getting aren't perpendicular to the contour lines. I can believe that for the sensitivities, but for the finite difference gradients? They even go the wrong way sometimes! I noticed this while staring at the output from INPUT.521. The contours and vector fields by themselves are plausible, but together they're inconsistent. 29 June 1994 I made a few interesting plots with SLAB today, of the sensitivities versus the finite differences. I want to spend the rest of today writing up what I've done so far, and investigating how to go about making the SLAB graphics available on the Macintosh, presumably by making CGM graphics and then having GPLOT turn them into PICT files. 27 June 1994 For MAX #7, I have a plot of the sensitivities on a contour grid. It's interesting, because it makes the problem a little more clear. The long skinny elliptical valley, with the tiny, tiny gradient arrows in the base of the valley. Even the normalized arrows show some interesting features. To proceed, I must first figure out how to train the code to dump the finite difference sensitivities computed by 611. Maybe I'll just compute them myself, for comparison. Then, to do the banana problem, I have to learn how to project the gradients onto that plane. No, just take increments in the plane where you are working. Perfect. OK, but then I still have to project the vectors returned by the sensitivity calculation. I think I have done this, and I will verify it by setting up and running #520, which will examine the banana problem, and write out the cost and the (projected) cost gradient. If this plots out nicely, I'm in good shape. 25 June 1994 For MAX #6, I'm not real happy with my example, since the only reason it's better than an unpenalized problem is that the unpenalized problem gets stuck in a local minimum. Doing almost ANYTHING would knock the iteration out of the local minimum, and set it going to the global minimum. I'd like to find a better example, but for now it will have to do. To compare finite difference sensitivities to the real things, I think I could make some very nice plots. What I want to do is make some 2D runs where I evaluate the function and sensitivities at preselected points. Then I want to plot the vector field of sensitivities, both kinds, and compare them. That means SLAB has to be able to read function values, direction fields, and compute the projection of the direction field onto the plane of reference. Or can FLOW take care of that last part? You need to pull out some vector code. I modified SLAB so that it can accept X and Y components of a vector, and plot the vector field. I modified MARCH2 and MARCH3 to write out the GOPT vector to the plain text output file. I'm running a simple two parameter problem, and we'll see that the vectors point to the minimum, I hope. 24 June 1994 As my first stab at MAX #6, I'm going to use 501 as my base problem, the one with the banana. Then I'm going to compare it to 518, a new problem, where I start out with a small amount of bump penalty, say 0.5, converge weakly, and then restart. Well, I wasted all day on this problem. I'm back to the point where I have to say that the three-parameter bump problem is too ill-conditioned to work with. I needed a problem that was hard to solve, but which a penalty could help with. I tried the banana problem - no luck. The penalized problem got past the local minimum, but then got stuck at a stupid point - something like 0.2, .5, 0.3 instead of 0.375, 0.5, 0.375. If I increase the grid size, both the original and the penalized problem are able to converge. The two parameter problems are too easy, so the penalty has nothing to offer. Well, here's my last idea for today. Let's try and understand why the bump penalty doesn't help the banana problem. After all, if only it didn't get stuck too, it'd be a good example. So let's try it one more time. If it still fails, then let's look at a 3D slice including the local minimum, where the penalized code stops, and the target... OK, on a 21 by 7 grid, using WATEB=0.01 initially, followed by WATEB=0, I was able to converge in about 55 steps total, whereas, as we all know, the regular approach eats the wrong end of the banana. 23 June 1994 The run 515 was made, but the movie is bad, because the ellipses are back to being long and skinny. This is slightly maddening! Let me try a temporary input 516, which drops pressure from the cost. Maybe that will help. Temporary, schmemporary! 516 made a nice movie, with the long narrow loops smoothing out to circles. So now I'm running 517, which will compute the location of the minimum at 11 points, so I can (I hope) mark up the movie with points. Results: WATEB Imin, Bmin 0.50 0.51277, 0.016938 0.45 0.51278, 0.018843 0.40 0.51278, 0.021230 0.35 0.51279, 0.024307 0.30 0.51280, 0.028422 0.25 0.51281, 0.034205 0.20 0.51283, 0.042912 0.15 0.51283, 0.057461 0.10 0.51280, 0.086324 0.05 0.51229, 0.164886 0.00 0.50000, 0.500000 Damn! Does this mean I should cut down the range yet again? No, I think I'll try to work with what I have for now... I'm glad I went with this data. I found out that I COULD do it, but it took some effort. First, I made a PICS file using DICER. Then I read the PICS file into GRAPHIC CONVERTER, and went, frame by frame, marking the location of the minimum with a bright green circle. I saved each frame as a separate PICT file. Then I used CONVERT TO MOVIE to turn the PICT files into a MOOV. It's short and jerky, but it makes the point. 22 June 1994 Today's modest goal: Get the location of the minimum for the 2D problem, MAX#4, as the parameter moves from 0 to 10. I wrote INPUT.514 to do this. Results: WATEB Imin, Bmin 10 0.49130, 0.0052397 9 0.49126, 0.0052793 8 0.49112, 0.0054758 7 0.49085, 0.0059561 6 0.49038, 0.0081341 5 0.49065, 0.011099 4 0.49065, 0.011099 3 0.49112, 0.020401 2 0.49100, 0.020792 1 0.49144, 0.060016 0 0.49593, 0.29901 How do I explain these wretched results? My first question is, what happened at WATEB=0? I'm going to have to rerun this problem with ITYPE=0, and see if I have more confidence in the results. Note that it really seems as though 611 thought it had a minimum for WATEB=0. How can that be? Did the sensitivities crap out again? In that case, time to solve a bigger problem! Secondly, it looks like I should start at WATEB=1 and come in to 0, like I always thought. I recall that didn't seem to make a good movie, because the ellipses weren't nice, but I can't believe that. I'll try that again. If possible, do the same for MAX#1, where the parameter moves from 0 to 1. Second run, WATEB goes from 1 to 0: Results: WATEB Imin, Bmin 1.0 0.49233, 0.074194 0.9 0.49233, 0.074194 0.8 0.49233, 0.074194 0.7 0.49233, 0.074194 0.6 0.49233, 0.074194 0.5 0.49334, 0.10441 0.4 0.49298, 0.10482 0.3 0.49339, 0.13267 0.2 0.49519, 0.22354 0.1 0.49760, 0.37342 0.0 0.49999, 0.49999 Here, at least I have a happy result at the end, even though I have to explain the poor performance of the first few steps! I will rerun this problem, going from WATEB=0.5 to 0. If that turns out OK, then I will try to generate a corresponding movie, so I can label stuff. Results: WATEB Imin, Bmin 0.5 0.49250, 0.093188 0.45 0.49250, 0.093188 0.40 0.49250, 0.11856 0.35 0.49337, 0.13589 0.30 0.49337, 0.13589 0.25 0.49414, 0.16936 0.20 0.49502, 0.19914 0.15 0.49542, 0.22696 0.10 0.49509, 0.30678 0.05 0.49870, 0.41029 0.00 0.49999, 0.49997 Before I completely buy this, let me try two things. A) drop the resetting of the optimization tolerances. (This made almost no difference) B) try OSOLVE. Results: WATEB Imin, Bmin 0.50 0.49319, 0.098633 0.45 0.49340, 0.10878 0.40 0.49364, 0.12110 0.35 0.49395, 0.13633 0.30 0.49433, 0.15549 0.25 0.49481, 0.18005 0.20 0.49544, 0.21216 0.15 0.49626, 0.25487 0.10 0.49730, 0.31242 0.05 0.49858, 0.39087 0.00 0.49999, 0.50000 Clearly, using OSOLVE made a big difference, especially in terms of the "smoothness" of the change in the minima that are calculated. PSOLVE tended to get stuck, and to sit at a point, needlessly recomputing, knowing something was wrong, but unable to move. In particular, OSOLVE used 126 optimization steps, while PSOLVE used 315! OK, now let's look at the movie generated by OSOLVE on this problem, looked at as a 3D problem... I think I should run this as an 11 by 11 by 11 problem! Well, I'll have to wait til tomorrow to pick up the results. 21 June 1994 Worked on SLAB a bit. Fixed the missing lowest color contour region problem. Got the color bar to show up. Modified CHRPLT so it can left, center, or right justify. I then totally wiped out my copy of SLABSUB. Thank you, UNIX, very much. In response to a plea from Jennifer, I resurrected DISPLEE, and am trying to get DRAWCGM running on the SUN. Progress was minimal, given that the FORTRAN compiler on the ICAM Sun can't seem to find the license server, so it won't go anymore, and there seem to be some missing routines from the X libraries on Corona. 20 June 1994 Now, in pursuit of modified MAX#5: redo the movie for MAX#1 but with the minimum (for each slice) drawn in, showing how the minimum moves smoothly. This might be worth a new option on the code? Solve the 3 bump parameter problem with just U in the cost, or solve the 1 bump parameter problem with U, V, P in the cost, with WATEB=1. Use those values of the parameters as a starting point for solving the optimization problem with WATEB=0.75 (or whatever). And so on. I'm dropping the 3 bump parameter problem. It's too singular. I added code to allow the user to specify that REYNLD is a variable (NPARR=1) or is not (NPARR=0). If it is not, then the user may specify a particular value, REYNLD. By the way, what is the situation for inflow? Tomorrow, go back to worrying about making a movie with the minimums marked. 18 June 1994 I'm still having a lot of trouble getting nice smooth ellipses out of the original problem, on MAX#4. Got the ellipses. I did have to use WATEB varying from 0 to 10, which I think is excessive. I probably had a mistake in the earlier, failed, runs, and I now think WATEB going from 0 to 1 would also have ellipses... 17 June 1994 Working on MAX#4. I got an initial 2D slice, and saw the long, thin valleys again. So now I've set up a 3D run (508) and we'll see what we can get. Janet says Max has another list for me! 16 June 1994 I got MAX#2 done. I did slice, stack, and transparency movies of the cost functional involving both horizontal and vertical velocities, which, as expected, did not involve a banana shape. I've submitted a 3D march for pressure now, (MAX#3) and am starting to do the 2 parameter problem (MAX#4). 15 June 1994 Added the ability to read character data to NAMELS, and added FILEG and FILEM to the list. The problem with WATEU=WATEV=1 reached the global minimum. So now I'm running INPUT.504, which will make a 3D picture of the functional around the point which was the local minimum for 502. We'd expect to see that banana shape all cleared up now. (I had a few problems with 504. I had put a few bugs into the code, and the CALVIN server went down, among others. I've submitted the 504 run, but I'll have to come back tomorrow to see it.) 14 June 1994 I'm going to look at MAX #2 now, including the cost of V, and seeing if I get a local minimum (probably) and a banana shape (also probably). It looks to me like I should add a WATER term to allow me to penalize discrepancies in the REYNLD value as well. 13 June 1994 I made some preliminary PICS animations with DICER and Janet suggested I make the pictures bigger, so that's what I did for my final one. It was a fairly involved process. Leaving out the part about getting the binary data to the Mac: I ran DICER, specified an 8 by 12 by 12 dataset of FLOATS, surveyed the range, and then cut it back to 0 to 0.0001, and magnified it 32 times. Then I selected 2 slices (the first and last) and asked DICER to generate a 60 frame animation, saving the result as a PICS file. Then I had to use ConvertToMovie to convert to QuickTime movie, which is playable with SimplePlayer. Unfortunately, the movie was too large, so I had to use Stuffit to squash it down to about 350K. Janet's going to take a copy home to show Max. 11 June 1994 I am running the test with WATEB=0 to WATEB=0.01. 10 June 1994 Some good news: I was able to see the banana shape using DICER. However, it's such a pain to try to read a plain data file! Also, I had to input just a few slices, and be very careful about the index ordering again. Still, the functional values are very small on the first slice, and quite much larger thereafter, so you have to restrict the range. For further work, I propose regenerating the data, and going from WATEB=0 to WATEB=0.01... 09 June 1994 I resurrected FLOW1, and restored the MARCH1 and MARCH2 codes to it, so that I can make comparisons, if necessary. Also, it still has the built in ability to create data files for SLICE, so I may want to rerun problem 20. 08 June 1994 I had a small success. I was able to generate some 3D data on the IRIS in a FORTRAN program, write it to a direct access binary file, transfer that to the Macintosh, and pass it into DICER. Now I am going to try the same operation, but using real data from FLOW3. By the way, the relevant steps included: open(unit=3,status='new',access='direct',file='csc.dat', form='unformatted',recl=1) irec=0 do i=1,nrow do j=1,ncol do k=1,nlay irec=irec+1 write(3,rec=irec)fdat(i,j,k) enddo enddo enddo Use binary mode to transfer data file to Macintosh. Run DICER, choose FILE/OPEN, select FLOAT from data type, then, on next menu, specify data dimensions NROW, NCOL and NLAY. Well, I forgot what a misery DICER is to use. Remember, for instance, that it assumes that the array has a range of -1 to 1, and you have to choose "Survey" AND then "Use Survey" in order to fix that. On top of that, you have to specify "float" every time you first open the file, and you have to specify the dimensions of the array, and the ordering of this is completely unclear. It's either X,Y,Z or Z,Y,X but it never makes any sense whatsoever. Well, entering my dimensions as Z,X,Y, I got a stack of ellipses, which is almost great...except that I was expecting a banana, damn it. So now I have to check my input, check the costs, replot them with a real package, and see if anything else went wrong... Oh, I almost forgot. On top of everything else, the output from the ALPHA was total gibberish. I had to use RECL=4, which makes no sense, and DICER read the data as ranging from -1.0e+38 to +1.0e+38. The RECL=4 is caused by the fact that, if you don't tell it not to, every DEC record has two bytes appended to it. So RECL=1 wouldn't work, and RECL=4 would. (I guess). By inserting the keyword "recordtype=fixed" in the OPEN statement, that error seems to have gone away. I think the other problem is just the dreaded DEC byte swapping that can be cured by the CHRPAK routine BSWAP4. 07 June 1994 I found the original data for the 2D banana contour, and have set up the 3D march to include it. I am running this on the ALPHA. 06 June 1994 The MARCH3 run is working. (I had left NY out of the argument list, so the cost was coming out zero). I made a 4 by 4 by 8 run, and will take the data downstairs into DICER and see what it can do with it. The usual frustrations with the stupid DICER user interface. Right now, it does not read the numbers in correctly, and it only seems to read in one slice, although it never tells me there's a problem. I'll get the documentation from Janet. 05 June 1994 Now I'd like to add a type 4 run (ITYPE=4), corresponding to a 2D march plus regular variation in the penalty function. I will assume the penalty multiplier WATEB goes from 0 to 1, and so all I need to specify further is the number of steps taken in WATEB. I've written the basic MARCH3 code, and compiled it. I added a routine, GETUNT, to keep track of FORTRAN units. This is to make it easy to pick a unit number for MARCH3 to write the cost information data to. I made a run of MARCH3, but it's very fishy... 04 June 1994 I worked on MARCH1 today. I believe I set it up properly, and I ran a test case, which was OK. This is in preparation for MARCH2, which will allow me to look at contour plots of the cost functional. I worked on MARCH2, making the same sort of updates and corrections I had made to MARCH1. That seems to work, too. Now I need to think about how to work on item #1 on Max's list. 03 June 1994 I was despairing, because my simple problem with one inflow and one REYNLD parameter was floundering, getting ridiculous results. I decided that the cheapest test would be to get the "OSOLVE" code running, that is, the stuff that is the same as PSOLVE except that it uses finite difference gradients instead of sensitivities. It turned out that I had a simple sign error in my REYNLD sensitivities. So I've fixed that, and I have OSOLVE running as well. Two small victories. I ran this problem on the ALPHA, to convergence, and this will be run #9. By the way, the ALPHA also accepts the "-fast" switch. Looking now at Max's suggestions, I see that I will need the MARCH codes. So I'm working now on getting MARCH1 running. 02 June 1994 Wu pointed out that you can speed up execution on the SUN by compiling with: f77 -fast -O3 myprog.f I retried the Poiseuille flow problem, but modified the weighting factors so that WATEU=WATEV=WATEP=1, which presumably will help me zero in on the correct solution... 17 May 1994 I may have found the teensy error, namely, that YPROF had not been set. I don't need YPROF, so I'm returning to accessing YC directly when needed, via NPROF. Well, that fixed the problem with the ridiculous costs. Now I will check to see whether the pressure cost has been fixed. No, the pressure cost is still very large. But that's plausible, since the code is converging to a solution with REYNLD=1, instead of 50! What's going on here? Actually, I guess I never tried to optimize over REYNLD and maybe I should look at that first! Poiseuille flow strikes again! I set a target with an inflow parameter of 0.5 and a REYNLD parameter of 50. I set a starting point with inflow 0.5 and REYNLD of 1. The two flows matched, except for the pressure. 16 May 1994 Oh, I'm discouraged. I ran the PSOLVE code, which uses sensitivities to approximate the cost derivatives, rather than finite differences, and didn't get a real great performance. After 20 steps, the code is stuck at a cost of about 0.0013 and moving slowly. More worrisome, though this has no effect on the run, is that the computation of the pressure discrepancy cost seems way out of line with the values FLOW1 gets. So that's what I'm going to concentrate on right now. Well, now I know that the reason the pressure cost is nutty is that the target pressures are nutty! But the velocities are OK. How can this be? I still can't find the problem, but I had a great deal of difficulty justifying the variation between declarations of arrays with maximum and with used dimensions, so I spent two hours setting all the important arrays to use maximum dimensions everywhere. 15 May 1994 I deleted the variable DCDA, which was simply a selection of the entries of SENS. I am still trying to get PSOLVE to work, and I am trying to find out why we aren't computing approximate cost derivatives correctly... 14 May 1994 Working on OSOLVE. I wasn't computing the cost when the optimizer was trying to get the gradients of the cost function. OSOLVE seems to be working now. Run #9 demonstrates the results. Now I want to get PSOLVE working. 13 May 1994 Got rid of NVAR and MAXVAR, which were not needed. With great misgivings, I put together a version of FLOW3D that supposedly includes the capability of optimization. I will rerun problem 1, and then make up a test problem 2 to see if it's still OK. 12 May 1994 I tried to set up MSOLVE as it was. I also printed out PSOLVE and OSOLVE from the old code, so I can compare. Moment of doubt: Does SOLCON do continuation only in the Reynolds parameter? Exactly what does TANGNT compute? 06 May 1994 The results of 26 April suggest that I have a successful starting point strategy. I now propose to isolate that chunk of code in a new subroutine, "GETTAR", and then work on restoring the ability of MSOLVE to do an optimization. Run #8 demonstrates how the new MSOLVE and GETTAR work together. 26 April 1994 Successful run #7, which started with parameters (0, 0, 0, 0, 1) and tried to reach (1/2, 3/8, 1/2, 3/8, 50) with a maximum continuation stepsize of 20. This required 3 continuation steps. The stepsize was 16 1/3. Newton worked fine, and the code reached the target point. So I'm more confident that I have code which can automatically decide to continue if needed, in order to compute the target point, or any point, for that matter. The usual UNIX catastrophe, cp flow3.f flow3sub.f > fred.f Instead of cat flow3.f flow3sub.f > fred.f I will try to smile as people tell me how wonderful UNIX is. I've gotten to the point now where the Newton iteration is the problem. I fixed the fact that the jacobian was singular (because REYNLD had not been set). Now I'm getting the Newton iteration not converging - the residual is doubling each time instead of dropping. Must be adding the wrong things...All is well now! I have printed the run as RUN #6. Next goal, see if continuation works for target points, in the sense that I give a "large" value of REYNLD for the target, requiring SOLCON to take several steps, rather than trying to get there in one jump. 25 April 1994 I'm trying to get SOLCON to compute something now. I think some information is not getting to the right place. Right now, I'm facing a problem where, once I go into the Newton iteration, the jacobian seems to be all zero. 14 April 1994 Run #5 shows that SOLCON works properly, because I was able to get a starting point. However, I need to drop NCON, and switch to something like STPMAX, which measures the maximum difference in parameters allowed before we try continuation. That's my task for today. 13 April 1994 I have written FLOW3, which includes the marching code, and doesn't use it yet. It also needs to add the optimization code. But anyway, it has the MSOLVE routine back in it, which figures out what we're doing. Also, I moved the continuation code into SOLCON, and tried to start setting it up as a continutation solver. Given a solution GSAV and parameters PARSAV, compute solution G at PAR. Just to see if it makes sense, I will now try to run the old problem with FLOW3. This will be run #5. 12 April 1994 Today, I am simply going to try to restore the 1D and 2D marching capability to FLOW3. Question: what happens to local minimum when we refine the mesh? Is it really there? Rerun FLOW run with NX=41, NY=13. This will be run #4 for this new series. Answer: The code "hangs up" at what may be the local minimum, but then escapes from it and gets to the global min. This happens with the functional behaving as follows: .42, .39, .0098, .0095, .0082, .0057, .0109, .0005, .0009, .00029, .00027, .00025, .00017, .00027, .000098, .000053, .000022, .0000053, .0000029, .0000027, .0000027, .0000027, .0000027, .0000027, .0000027, .0000027, .0000026, .0000023, ... and then getting to .237E-11. 11 April 1994 In response to Max's challenge, I first ran input.01 with FLOW (run #1) and input.01 with FLOW2 (run #2). I figured I'd first better make sure I can recompute the boomerang functional picture. As I recall, the parameters were: REYNLD=1, NX=21, NY=7, INFLOW=0.5, BUMP=0.375, 0.5, 0.375 I have submitted this job, and will check it out. The telltale sign is that it does not get to the target value. Yes, the FLOW code indeed fails to reach the target point, stopping instead at a local minimum. This was run #3. I think I should "freeze" NUFLOW, calling it FLOW2, and move on to "FLOW3"... 11 April 1994 Max gave me a list of things to work on: Contour plot of functionals: 1) Functional involving U only, with penalty. Make plots or movie of functional shape versus penalty parameter. Penalty parameter goes to zero, we get boomerang shape. Hope for large enough penalty we get nicer functional. 2) Functional involving just U and V. Do we get the boomerang shape? (no) 3) Functional involving just P. 4) Functional involving U, V, and P, with penalty. Make plots or movie of functional shape versus penalty parameter PP. As PP goes to zero, get thin ellipses, PP increases, fat ellipses. Watch as PP goes to zero, and minimum moves to true minimum. Plot the iterates: 5) Plot the position of the iterates on a contour plot of the functionals, and show how the penalty parameter helps the iteration. Now Max suggests that it would be good enough to simply plot the location of the minimum on each of the slices, and show how the minimum of the "nice" problem smoothly transforms into the minimum of the original problem. I'll have to see if there's an easy way to modify the individual PICT frames of my movie. Then I'd simply have to compute the location of the minimum for a number of values. 6) Show how letting PP go to zero as the iteration proceeds gets the true minimum, and is also efficient. Sensitivities: 7) Compare finite difference sensitivities to the PDE sensitivities. For boomerang functional, solving for PDE sensitivities got us into trouble. Does the PP help? 8) Sensitivity plots (sensitivity versus parameter values) Reynolds Number Effects 9) All of the above at higher Reynolds numbers. Non-Feasible Targets 10) All of the above for nonfeasible targets, eg, compute the target flow using a different length bump (both shorter and longer). Reference Plane 11) Move the reference plane further away from the bump, and show that the sensitivities dramatically decrease in magnitude. The theme of this work is that you have to be careful, and check what you're doing, otherwise your sensitivities won't be sensitive, or you'll be unable to reach your objective, and so on. 07 April 1994 I just spent a day or two trying to make a quick movie, and I want to record all the problems that occurred so that they can be addressed: * First of all, the streamline code still is not correct for isoperimetric elements, and so the most interesting pictures can not be drawn. The code that computes the stream function has to be corrected. * Secondly, if you try to draw the velocity field using NXSKIP and NYSKIP, there is no averaging done. Hence, since the most interesting stuff happens in row 1, it doesn't show up at all. Answer: I have sketched out some replacement code. * The scaling of the velocity vectors does not take into account which vectors will be seen and which not. Hence, if vectors are omitted because the region is shrunk, or because of NXSKIP or NYSKIP, vectors are not scaled up in size properly. * GPLOT was not working on the IRIS. This has been fixed. I contacted Anjana Kar and got a new copy. * GPLOT is not available on the PowerPC Mac, and the MacII + FPU version will not run on it. I have contacted Anjana about this. * My color contour plots are ENORMOUS. For the 61 by 19 region, a single color contour plot was about 7 Megabytes. A simple line contour plot was more like 20 K. (This was PostScript files I believe). One cause of this is the fact that I draw the contours NCON times, overlaying the colors. This is wasteful, it turns out. Secondly (and harder to fix), I treat one element at a time. * I don't know how to turn a PostScript file into a PICT file. There's supposed to be an IRIS utitily, TOPICT, to do it. * The contour levels are set one way for line contours, and another way for color contours. Why? This means that if you naively draw both, they don't line up nicely. * It would be nice to print a color bar. I wrote some code to do this, and stuffed it into DISPLAY. Right now, it does not make sure that there's enough room on the screen to show the color bar! * How come the only simple color table goes from low blue to high yellow? What ever happened to low blue to high red? * The operations after "Setting data window to physical window" are noticeably much longer than simply reading the data. What's going on? Answer: Setting visibility of elements (not needed, and now replaced by setting them all to visible) and computing the location of the boundary. The boundary calculation is very general, because at one time I wanted to also handle Lee's problems. I don't need that anymore, so I can calculate the boundary apriori, just from knowing LONG. Do so, and save some time! * If you change the range of the contours, and the first value you choose is higher than some values in the flow field, they show up as a white region. Ugh. 05 April 1994 DISPLAY can now show contours or color contours of the X or Y component of velocity. 04 April 1994 The anomalies of yesterday are explained: I had a mistake in the UVAL routine, which meant that V and its derivatives were returned as 0. This explains why the code was able to jump to REYNLD=500 in two steps, and why the sensitivities looked so goofy. I am recompiling, and will rerun the two interesting cases of NX=21, REYNLD=10 (quick test) and NX=61, REYNLD=500 (slow test). For REYNLD=10, I get the results I expected. For REYNLD=500, it's not so clear. But one problem is that, with NXSKIP and NYSKIP, I'm just skipping data, so I miss the large values for derivatives, which are just near the bump...So I have to figure out how to average. 03 April 1994 Well this is very strange. On my 61 by 19 run, I forgot to set NCON to 10, leaving it at 2. And yet, I was able to go from a zero solution at REYNLD=0 to a solution at REYNLD=500 in 2 continuation steps! What's going on? 02 April 1994 I have got the sensitivity with respect to the inflow parameters running, I think. I just tacked the REYNLD sensitivity onto the end. DISPLAY plots values for INFLOW and REYNLD sensitivities, although I think the REYNLD sensitivities look wrong. (For a 1 parameter INFLOW, INFLOW and REYNLD sensitivities should be the same. They aren't.) Now I have to figure out the BUMP sensitivities. I've coded them up, but they don't show up as anything but zero. I've found that the routine BMPSEN seems to be computing the right hand side of the linear system as zero, so that's a start. Now I've got the bump guys looking good, and so it's just the REYNLD numbers I don't believe, at least not when I look at them. That's enough for today! Now that I think about it, the REYNLD and INFLOW sensitivities are unrelated, and I believe the picture I saw for REYNLD. So now, on to 500! I will try the following run: NX=61 NY=19 PARTAR(5)=500 NCON=10 IPLOT=1 (only 1 plot at 500, unless I modify the code...) Things to do when you can: If you can, modify the code so that you specify REYSTP, rather than NCON. Also, you need to dump a plot every continuation step, temporarily. 29 March 1994 Max wants some slides from me by next Friday, 8 April. I haven't found out exactly what he wants yet. I propose, though, to try to get the code running on the Iris to REYNLD=500, and to add the sensitivity calculations back in, and to do some plots of bump sensitivity at high REYNLD. Before I go on, I want to get rid of the results files from my attempts to reach REYNLD=1000, so let me make a table first: NX NY REYNLD Time ------------------------------- 31 10 250 174 s = 3 m Iris 41 13 400 2133 s = 35 m Iris 51 16 530 6205 s = 103 m Iris 61 19 680 16000 s = 266 m Iris 71 22 870 736 s = 12 m Cray 209 MFLOPS 81 25 1000 1320 s = 22 m Cray 216 MFLOPS One thing that help me up was the stepsize, so let me try running 61 by 19 with stepsize 25! Result: I was able to run the 21 by 7 problem, from REYNLD=0 to REYNLD=100, using a stepsize of 50. So let me try that stepsize for the 61 problem now! That worked too. I worked out the sensitivity equations, and found that they are not identical to the Navier Stokes equations, because the nonlinear term adds a complication. [LATER NOTE: NO, THEY ARE THE SAME, DON'T WORRY!] I REALLY don't understand where the boundary conditions come from for the bump, though. I want to modify the program so that PAR contains REYNLD as well. This affects PAR, PARTAR, NPAR, and so on. I'll have to be careful. I got this to work. The point of this work is to put all the "parameters" on an equal footing, so that I can easily set up continuation without lots of special cases. The task now is to compute the sensitivities. I can still just continue on REYNLD for now, but I want to be able to compute the sensitivities for various parameters. 23 March 1994 So now I believe that NUFLOW and DISPLAY are working OK together. That is, NUFLOW is back to writing a formatted file, and DISPLAY can read it. The only good news to come out of this is that I now have NXSKIP and NYSKIP working, so that I can draw velocity vectors thinned out as I like. Where does this leave me? * I still wanted to see a plot at REYNLD=1000, which means I need NX=81, NY=?, which means I need to run on the Cray, save a plot file, transfer it to Iris, and view it. * I need to be able to compute the sensitivities. Perhaps I can look at the old NSTOKE code and figure out what's going on. Is this just the tangent vectors? 21 March 1994 I'm having a goofy problem with NXSKIP and NYSKIP, because somehow DISPLAY's value of NX is getting corrupted. Right now, it's 1017, which is impossible. This turns out to be because NVAR and NX got mixed up. So I fixed that, and rather than regenerate the test on the Cray, I just tried to use a DISPLAY.DAT file sitting in the directory. But it failed (it reported the region was 0 by 0, so I didn't wait around for the real bad news). So I recompiled NUFLOW, and ran TEST01, and made a plot file. And guess what, DISPLAY still thinks the region is 0 by 0. So I have to go back in and find out what went wrong... After flailing around trying to find why I can't do graphics anywhere now, I GIVE UP ON BUFPAK! I'm going back to writing stuff out by hand. That means ripping up a lot of stuff in DISPLAY too, but tough. I've had it. 17 March 1994 I've been preoccuppied. So I took the file d_form_cray.dat and ran it through CRAY2IRIS (10 words per record) and made d_uform_iris.dat and ran DISPLAY and saw nice pressure plots. Now I'm going to try to set NPSKIP so I don't get overwhelmed by velocity vectors. Yes, NPSKIP works. But I think I need to be able to control both the row and column, since otherwise I'm just thinning out in one direction. Let me think. OK, I thought. I replaced NPSKIP by NXSKIP and NYSKIP, so I can specify separate row and column skips. I haven't tested this yet. 22 February 1994 So I wrote a little program that reads a Cray formatted "direct access" file (which is really "sequential access!) and converts it to an unformatted direct access file, and I used it on the file I grabbed from the Cray, and it worked. 21 February 1994 It's a waste of time trying to deal with the Cray and graphics. After I spent all that time setting up to write a formatted direct access file, it turns out the Cray pads it all to hell with junk, so why bother? I'm ready to boot the whole of BUFPAK as well, except that it works fine as long as I stay on one machine. So my strategy now should be to write a one time, plug-ugly routine to convert a Cray graphics file to text mode, transfer it, and convert it back to direct access unformatted. 20 February 1994 I modified NUFLOW and NUFLOWSUB by adding the new BUFPAK routines. I also changed NUFLOW so that IWFORM is read from the input file, allowing me to decide at runtime to create formatted or unformatted graphics files. Now I will have to make a simple run on the Iris, and see if DISPLAY can read formatted or unformatted graphics files from NUFLOW, so that I can then repeat the experiment on the Cray. The new DISPLAY can read unformatted files from the new NUFLOW. I got DISPLAY to read a FORMATTED file from NUFLOW, but I thought I would only have to have IRFORM and possibly NWORDR correct. Actually, I also needed WFORMT... This has been a bug for a while. Can't I put NWORDR and NFFAC together and make WFORMT? Now I am trying to generate on the Cray and display on the Iris. I just submitted a Cray job... 19 February 1994 I got my plot file routines rewritten so that it is fairly easy to read or write a formatted or unformatted file on the Cray or Iris. My plan, then, is to rerun the 16 February problem, create a formatted graphics file, transfer that to the Iris, convert that to an unformatted graphics file, and use that with DISPLAY to draw a picture... Naturally, before I run the huge problem, I'll want to test a smaller version... 16 February 1994 I reached REYNLD=1000 with NX=81, NY=25. Now I want to generate a plot file at that point so I can look at stuff. I just realized that, right now, I don't compute the sensitivities, for instance, so I won't be able to plot them. On top of that, I realized that I write the plot file as an unformatted, fixed length file. I can't transfer that from the Cray to the Iris. So I will try to rewrite BUFPAK so that I can choose at run time to write formatted or unformatted files. This should be possible through subroutine MEMORY. 15 February 1994 Larry: for run 71 by 22, got to REYNLD=300 before running out of time. Increased time from 900 to 1000 seconds, cut NCON from 100 to 25. First useful result from Larry. 71 by 22, broke down at REYNLD=880. And I was able to take steps of size 40... 14 February 1994 The usual problems getting going on Larry: they force you to change your password and then log you off; I forgot that QSUB must either be capitalized or start in column 2, or maybe both. I'm having trouble estimating the necessary time for a run. But in any case, I'm trying to get the 71 by 22 problem going. I talked to Max and Janet. Max said that my experiences with the "breakdown" are plausible, since there's something called the cell Reynolds number condition, which says that some spacing h times the Reynolds number has to be less than 1. So my values of NX and NY and breakdown REYNLD should rise roughly linearly together, which they do. He seemed to think that working with a REYNLD value of 500 (which can be done on the Iris) might be reasonable. He suggested that once I get a result at 1,000, I might spend some time doing optimization at REYNLD=500, and do some plots to show what happens to the sensitivities, and, for that matter, the basic solution. Then it would be time to start computing a reduced basis at REYNLD=1 and marching out in REYNLD as far as possible. 09 February 1994 My near-term goal is to be able to compute a flow for REYNLD=1000. I have set up a simple, fixed stepsize continuation process, starting with a zero solution at REYNLD=0, and using the tangent vector for a first approximation at the next value of REYNLD. However, it is clear that I'm going to have to use the Cray to do calculations at high values of REYNLD, since I have to cut down the cell size so much. My Cray account on LARRY is still good. I shall also check PK. I must find my NQS file and prepare a copy of the code to run on the Cray. 08 February 1994 Here are the results of continuing with a step in REYNLD of 10: 21 by 7, breaks down at 260 31 by 10, breaks down at 270 41 by 13, breaks down at 400. 51 by 16, breaks down at 530. 61 by 19, breaks down at 680. 71 by 22, segmentation fault (too big for the Iris!) 81 by 25, segmentation fault (too big for the Iris!) 07 February 1994 Now let's try to compute a target solution at REYNLD=200 by continuation in a fixed stepsize. For simplification, folded INDAT3 into INDAT1, and dropped INDAT2, MAXVEL, SETREY. Sketched a simple continuation method, which seems to be working, for the simple task of starting with a zero solution and computing a solution for a given value of REYNLD (and with the parameters not changing at all). Got solutions at REYNLD=25 (1 step or 2), 200 (4 steps). When I tried to do 1000, though, I failed miserably, even with a stepsize in REYNLD of 10. The Newton iteration kept breaking down around 250-350. I am only using a 21 by 7 grid, so I guess it's time to go up one step to 31 by 10. 06 February 1994 I computed a tangent vector yesterday. Today I'm writing code to check it. TANGNT - The (initial) norm of the tangent vector is 1.0015 TANCHK - Maximum FRPIME*TAN discrepancy= 4E-17 ABS(DNRM2(TAN)-1)) discrepancy= 1E-16 TAN(NVAR)= 0.99846 Now I am happy with the tangent calculation! 04 February 1994 - Removed some routines: GRID2N, GRID4N, UVDUMP, XYDUMP Removed some variables: ISTEP1, ISTEP2, JPLOT, JSTEP1, JSTEP2, MAXSIM, NUMSIM, TOLSIM, Today's goals: get a loop started, set up tangent. 03 February 1994 - I talked with Janet and Max today. First, we agreed that I should try to get the code set up to solve higher Reynolds numbers, using simple continuation with fixed step size, or even extremely simple continuation, using no tangent information. Second, we discussed the reduced basis problem. Max first suggested setting up a bump problem with three parameters, generating the first and second derivatives, and seeing if that was feasible. Because the grid moves, we would have to have some sensible way of interpreting a representation of one flow in terms of another when their grids are different. So we decided that we should look at the inflow problem first. First, choose a high enough Reynolds number so that things are interesting. 400 was suggested, and I was told to run that test where I send in a very sharply varying flow to see how quickly it gets smeared out. Then set up the reduced basis system and solve it. Third, Janet asked that I start taking her through the program, and explaining what's going on. 01 February 1994 The Newton iteration is working now. I can start from a zero point, and get Newton to iterate to a solution. Now I need to look at how this will fit in to the general problem of optimization. Wait, first I'm going to crank up the Reynolds number to 200 and see what happens. The iteration DOES NOT CONVERGE. Now I will need to prepare for continuation. The first steps: Set up the DFDR vector (done). Tack on the Reynolds number as the last variable. (I've started this, by redimensioning some things...) Also, we can assume continuation is always on the Reynolds number, and so the Newton correction step may be written as FP*[-DX] = FX + DFDR * DR, where DR is chosen by me. 31 January 1994 The first run failed, in NEWTON, because DGBTRF returned INFO=1. Perhaps I'm not evaluating the jacobian? No, I'm certainly TRYING to evaluate the jacobian, but the first column comes back all zero. By comparing to FLOW, I'm trying to figure out what's not being set now... I have discovered that I assumed that PHI and PSI had been merged, but I hadn't taken care of that. Working on it... 27 January 1994 Lowered my ambitions. I dropped out the old CSOLVE code, and other extraneous stuff, and vowed that it will be enough if I can just get a single correction. 26 January 1994 Decided to append the Reynolds number to the end of the regular variables in the variable vector. 13 January 1994 I have a new set of routines, NEWTON, FX and FPRIME, which are the nucleus of my new code.