<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">Barry:</div><div class="gmail_quote"><br></div><div class="gmail_quote">Why the code does not crash in the case of 'nan':</div><div class="gmail_quote"><br></div><div class="gmail_quote"><div class="gmail_quote">Integrating characteristic equations ...</div><div class="gmail_quote">TS: 0.000000</div><div class="gmail_quote"> 0 SNES Function norm 1.490033584870e+02</div><div class="gmail_quote"> 0 KSP Residual norm nan</div><div class="gmail_quote">TS: 0.100000</div><div class="gmail_quote"> 0 SNES Function norm 1.490033584870e+02</div><div class="gmail_quote"> 0 KSP Residual norm nan</div><div class="gmail_quote">...</div><div class="gmail_quote"><br></div><div class="gmail_quote">Previously, when ilu/lu encounters a zero pivot, the code crashes, which asks user to do something about his model or the algorithm; </div><div class="gmail_quote">Now, it gives 'nan', but continues running.</div><div class="gmail_quote">We cannot assume new users know '-ksp_monitor -snes_monitor', and set '-pc_factor_shift'. The frustration we encountered is that the code seems working well, until I checked the solutions at each time steps and found them never change values. Then it took us hours to find out it was caused by shift.</div><div class="gmail_quote"><br></div><div class="gmail_quote">At least for LU, when zero-pivot occurs, we should throw an error.</div><div class="gmail_quote"><br></div></div><div class="gmail_quote">Hong</div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><br>
Hong,<br>
<br>
This is a question for the entire PETSc development team, not just me.<br>
<br>
The reason I made this change was I was sick of<br>
runs with ILU and nearly zero pivots where the preconditioned residual normed ended up initially being like 10^20 and then decreased in a few iterations to 10^10 and so the iteration stopped while, in fact, no progress on the true residual norm at all, meanwhile the user thinks the linear solver is working fine. I figured was better to just "stop" the entire process so the user can do something smarter for preconditioning than using ILU(0). At worse the user can now add the pivot option themselves.<br>
<br>
Was this a good choice, perhaps not? Should we just go back to the old default? Probably we need to at least have better help in terms of telling people to turn on that option?<br>
<span class=""><font color="#888888"><br>
<br>
Barry<br>
</font></span><div class=""><div class="h5"><br>
<br>
<br>
> On Jun 10, 2015, at 4:47 PM, Hong <<a href="mailto:hzhang@mcs.anl.gov">hzhang@mcs.anl.gov</a>> wrote:<br>
><br>
> Daniel,<br>
><br>
> Try to use it with the flag: '-pc_type svd'. Now it works.<br>
> Great, you found the problem!<br>
><br>
> This is very strange<br>
> OK, this must be the change made on the default pc_factor_shift used for ilu/lu.<br>
> In v3.5.4, it is<br>
> -pc_factor_shift_type INBLOCKS<br>
><br>
> now,<br>
> -pc_factor_shift_type NONE<br>
><br>
> Adding the option '-pc_factor_shift_type INBLOCKS'<br>
> or '-pc_factor_shift_type NONZERO'<br>
> the code runs well :-)<br>
><br>
> Barry:<br>
> For ilu, a shift should be used to avoid zero pivot;<br>
> for lu, '-pc_factor_shift_type NONE' should be the default.<br>
> Why do we make such change?<br>
> It seems you made such change:<br>
> <a href="https://bitbucket.org/petsc/petsc/commits/422a814ef4a731b8529cf3a6428db526d183e312" target="_blank">https://bitbucket.org/petsc/petsc/commits/422a814ef4a731b8529cf3a6428db526d183e312</a><br>
><br>
> Hong<br>
><br>
><br>
><br>
> On Wed, Jun 10, 2015 at 3:11 PM, Adrian Maldonado <<a href="mailto:dmaldona@hawk.iit.edu">dmaldona@hawk.iit.edu</a>> wrote:<br>
> I have been able to replicate the error with the new petsc-dev.<br>
><br>
> KSP does not converge. However, the Jacobian matrix is the same!!!<br>
><br>
> This jacobian matrix is not very well conditioned (condition number 1.4814e+03 for a 12x12 matrix) because of the valve... but it didn't give me any problem before.<br>
><br>
> In fact, you can see that I have implemented the same in matlab (matlab folder) and solved the problem with a vanilla newton-rhapson iteration (no preconditioner etc...) and I didn't have any problem at all.<br>
><br>
><br>
><br>
> On Wed, Jun 10, 2015 at 2:37 PM, Adrian Maldonado <<a href="mailto:dmaldona@hawk.iit.edu">dmaldona@hawk.iit.edu</a>> wrote:<br>
> I am upgrading to the new petsc-dev, I will let you know what I find. Maybe is an issue of the finite difference jacobian (how is generated)?<br>
><br>
> I would love to go to to the conference! How can I register?<br>
><br>
> On Wed, Jun 10, 2015 at 12:55 PM, Hong <<a href="mailto:hzhang@mcs.anl.gov">hzhang@mcs.anl.gov</a>> wrote:<br>
> I found which change made difference:<br>
> When I change Makefile<br>
> $ diff old.Makefile new.Makefile<br>
> 1,2c1,2<br>
> < include ${PETSC_DIR}/conf/variables<br>
> < include ${PETSC_DIR}/conf/rules<br>
> ---<br>
> > include ${PETSC_DIR}/lib/petsc/conf/variables<br>
> > include ${PETSC_DIR}/lib/petsc/conf/rules<br>
><br>
> to build tubtrans with petsc-dev, then I get<br>
><br>
> Loaded case file<br>
> 0 KSP Residual norm nan<br>
> SYSTEM SUMMARY<br>
><br>
> Number of variables: 12<br>
> Number of nodes: 2<br>
> Number of pipes: 1<br>
> Integrating characteristic equations ...<br>
> TS: 0.000000<br>
> 0 KSP Residual norm nan<br>
><br>
> Some petsc solvers have changed, but it should not behave like this.<br>
> We have to use a debugger to figure out which step causes this problem :-(<br>
><br>
> Hong<br>
><br>
> On Wed, Jun 10, 2015 at 10:57 AM, Adrian Maldonado <<a href="mailto:dmaldona@hawk.iit.edu">dmaldona@hawk.iit.edu</a>> wrote:<br>
> Sure.<br>
><br>
> This one has the object files, etc..<br>
><br>
> I will check the new petsc-dev this afternoon.<br>
><br>
> On Wed, Jun 10, 2015 at 10:42 AM, Hong <<a href="mailto:hzhang@mcs.anl.gov">hzhang@mcs.anl.gov</a>> wrote:<br>
> How about make a tarball of your current tubtrans and send to me<br>
> (or use google drive/dropbox)?<br>
><br>
> Hong<br>
><br>
><br>
> On Wed, Jun 10, 2015 at 10:28 AM, Adrian Maldonado <<a href="mailto:dmaldona@hawk.iit.edu">dmaldona@hawk.iit.edu</a>> wrote:<br>
> Strange.<br>
><br>
> I am using the master branch from the bitbucket repo (<a href="https://bitbucket.org/damalbel/tubtrans.git" target="_blank">https://bitbucket.org/damalbel/tubtrans.git</a>).<br>
><br>
> The jacobian is the finite difference jacobian (-snes_fd) as I didn't write one (I did write one in the matlab code prototype, though).<br>
><br>
> What I can think now is to check this singular jacobian and see what equation is causing the problem.<br>
><br>
> Let me clone the new petsc release and try to replicate the problem. I have a couple of things to do this morning but I will get back to you this afternoon.<br>
><br>
> On Wed, Jun 10, 2015 at 9:54 AM, Hong <<a href="mailto:hzhang@mcs.anl.gov">hzhang@mcs.anl.gov</a>> wrote:<br>
> Daniel,<br>
><br>
> It seems Jacobian is singular:<br>
> ./tubtrans -snes_monitor -ksp_monitor |more<br>
><br>
> Loaded case file<br>
> 0 SNES Function norm 1.500947612239e+02<br>
> 0 KSP Residual norm -nan<br>
> SYSTEM SUMMARY<br>
><br>
> Number of variables: 24<br>
> Number of nodes: 4<br>
> Number of pipes: 2<br>
> Integrating characteristic equations ...<br>
> TS: 0.000000<br>
> 0 SNES Function norm 1.490033613070e+02<br>
> 0 KSP Residual norm -nan<br>
><br>
> Do you use the same version in<br>
> git clone <a href="https://bitbucket.org/damalbel/tubtrans.git" target="_blank">https://bitbucket.org/damalbel/tubtrans.git</a>?<br>
><br>
> Can you apply my patch to this version and test it with new petsc-3.5 (we are working on its release)?<br>
> The best way to get petsc-dev is using<br>
> • git clone <a href="https://bitbucket.org/petsc/petsc" target="_blank">https://bitbucket.org/petsc/petsc</a><br>
> We will have a PETSc conference next week<br>
> <a href="http://www.mcs.anl.gov/petsc/petsc-20.pdf" target="_blank">http://www.mcs.anl.gov/petsc/petsc-20.pdf</a><br>
><br>
> Would you like join us?<br>
> Hong<br>
><br>
><br>
> On Wed, Jun 10, 2015 at 8:17 AM, Adrian Maldonado <<a href="mailto:dmaldona@hawk.iit.edu">dmaldona@hawk.iit.edu</a>> wrote:<br>
> P.S: I have applied your patch, and using PETSC 3.5.2 the program works all right. I have attached the logs of running test.xml and test2.xml.<br>
><br>
> Again, let me check out the new petsc release.<br>
><br>
> On Wed, Jun 10, 2015 at 8:08 AM, Adrian Maldonado <<a href="mailto:dmaldona@hawk.iit.edu">dmaldona@hawk.iit.edu</a>> wrote:<br>
> I don't see why it should not work. I can see two things that we can check:<br>
><br>
> - Is SNES iterating? Can you check for -snes_monitor and -snes_converged_reason<br>
> - In the program, I initialize the the solutions vector 'x'to ones as an initial guess (line 102), maybe somehow petsc is unable to write on it?<br>
><br>
> Let me clone the new petsc-dev version and I will try to replicate the error.<br>
><br>
> On Tue, Jun 9, 2015 at 4:56 PM, Hong <<a href="mailto:hzhang@mcs.anl.gov">hzhang@mcs.anl.gov</a>> wrote:<br>
> Danniel,<br>
> Thanks for detailed instruction. I'm using petsc-dev, which will be released today or tomorrow. I'm building a release version to see if that causes the problem.<br>
><br>
> However, your code should generate same solution with my patch (attached) because it only replace<br>
> VecGetArray() with VecGetArrayRead(), and fix a memory leak problem.<br>
><br>
> Here is what I did:<br>
> 1. install petsc-dev (see <a href="http://www.mcs.anl.gov/petsc/developers/index.html" target="_blank">http://www.mcs.anl.gov/petsc/developers/index.html</a>)<br>
> git clone <a href="https://bitbucket.org/petsc/petsc" target="_blank">https://bitbucket.org/petsc/petsc</a><br>
><br>
> 2. build petsc-dev<br>
><br>
> 3. git clone <a href="https://bitbucket.org/damalbel/tubtrans.git" target="_blank">https://bitbucket.org/damalbel/tubtrans.git</a><br>
> patch it:<br>
> cd tubtrans<br>
> patch -Np1 < patch_tubtrans<br>
> make tubtrans<br>
><br>
> 4. run tubtrans<br>
> ./tubtrans<br>
><br>
> Loaded case file<br>
> SYSTEM SUMMARY<br>
><br>
> Number of variables: 24<br>
> Number of nodes: 4<br>
> Number of pipes: 2<br>
> Integrating characteristic equations ...<br>
> TS: 0.000000<br>
> Printing solution vector...<br>
> Vec Object: 1 MPI processes<br>
> type: seq<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> ...<br>
><br>
> Changing to different input file:<br>
> ./tubtrans -f input/test.xml |more<br>
> Loaded case file<br>
> SYSTEM SUMMARY<br>
><br>
> Number of variables: 12<br>
> Number of nodes: 2<br>
> Number of pipes: 1<br>
> Integrating characteristic equations ...<br>
> TS: 0.000000<br>
> Printing solution vector...<br>
> Vec Object: 1 MPI processes<br>
> type: seq<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> 1<br>
> TS: 0.100000<br>
> Printing solution vector...<br>
> Vec Object: 1 MPI processes<br>
> type: seq<br>
> 1<br>
> 1<br>
> 1<br>
><br>
> Do I do something wrong?<br>
><br>
> Hong<br>
><br>
><br>
> On Tue, Jun 9, 2015 at 12:15 AM, Adrian Maldonado <<a href="mailto:dmaldona@hawk.iit.edu">dmaldona@hawk.iit.edu</a>> wrote:<br>
> Professor Hong,<br>
><br>
> I am glad to hear from you! I am happy that you are finding the code useful. As for the solution vector being constant ones - that doesn't seem right.<br>
><br>
> The master branch of the bitbucket repository with the Petsc Release Version 3.5.2 seems to work fine in my laptop.<br>
><br>
> I have written a small walk through below.<br>
><br>
> Please, let me know if you want me to take a look at the code or I can help you in any other way.<br>
><br>
> Walk through:<br>
><br>
> - The program uses XML (libxml2) to parse data for different case scenarios. The cases are stored in the directory 'input/' and you can select the case to run in the petscopt file (e.g -f input/test.xml)<br>
><br>
> - The file test.xml, for example, includes a simple one-pipe test scenario from the textbook 'Fluid Transients' by Wyle. The content is shown below:<br>
><br>
> <?xml version="1.0"?><br>
> <case><br>
> <nnode>2</nnode><br>
> <pipe id="1"><br>
> <fr>1</fr><br>
> <to>2</to><br>
> <fric>0.018</fric><br>
> <diam>0.5</diam><br>
> <a>1200</a><br>
> <len>600</len><br>
> </pipe><br>
> <reservoir id="1"><br>
> <node>1</node><br>
> <head>150</head><br>
> </reservoir><br>
><br>
> <valve id="1"><br>
> <node>2</node><br>
> <lossc>0.009</lossc><br>
> </valve><br>
> </case><br>
><br>
> In this case we have a conduct of two nodes (boundary equations). The only conduct (or pipe) goes from node 1 to node 2. The entries 'fric', 'diam', 'a' and 'len' are the parameters of the pipeline. Then, the file includes a reservoir and a valve as boundary equations.<br>
><br>
> The reservoir defines a constant pressure (water head) in the node 1 whereas the valve defines a relationship between the pressure and the flow with friction losses.<br>
><br>
> The program, after creating the necessary data structures and discretization (according to Courant number) will proceed to the initialization. The initialization is solving the system state for the steady-state conditions with SNES.<br>
><br>
> Running the program in gdb for 'test.xml' and printing the 'x' vector will produce this:<br>
><br>
> Breakpoint 1, main (argc=1, argv=0x7fffffffde68) at main.c:115<br>
> 115 ierr = SNESSolve(snes, NULL, x);CHKERRQ(ierr);<br>
> (gdb) n<br>
> 120 ierr = MatDenseGetArray(SOL, &ss);CHKERRQ(ierr);<br>
> (gdb) p VecView(x, 0)<br>
> Vec Object: 1 MPI processes<br>
> type: seq<br>
> 0.477432<br>
> 150<br>
> 0.477432<br>
> 148.698<br>
> 0.477432<br>
> 147.395<br>
> 0.477432<br>
> 146.093<br>
> 0.477432<br>
> 144.791<br>
> 0.477432<br>
> 143.488<br>
><br>
> The ordering here is, according to the discretization x = [Q0 H0 Q1 H1 ... Qn Hn]' . Where Q is the water flow and H is the water head or pressure. You can see how the water flow is constant in steady state along the pipeline and the pressure in the first node is 150, as the reservoir sets this as boundary equation.<br>
><br>
> For the transient we simulate the sudden closure of the valve (dirac delta) although we can implement different closure dynamics.<br>
><br>
> I have compiled the debug mode of the program (make debug) which prints the 'x' vector and each time step. I am attaching the results in the file log.txt so you can compare it. Each time step solution is stored in the matrix SOL and finally written to a binary file 'output.out'.<br>
><br>
> The first integration time-step looks like:<br>
><br>
> 0.477432<br>
> 150<br>
> 0.477432<br>
> 148.698<br>
> 0.477432<br>
> 147.395<br>
> 0.477432<br>
> 146.093<br>
> 0.477432<br>
> 144.791<br>
> -1.11838e-12<br>
> 441.046<br>
><br>
> The pen-ultimate entry of the x vector is the flow at the valve, which is 0 after the valve closes. At the end of the valve an steep rise in the pressure happens, as expected. This pressure wave will oscillate along the pipe according to the method of characteristics.<br>
><br>
> More complex scenarios are included in '/input' and a script with some python macros 'intPlot.py' is included to produce plots of the pressure and flow profiles within the system.<br>
><br>
> I hope this helped you and please, let me know in what can I help.<br>
><br>
> Yours,<br>
><br>
><br>
> On Mon, Jun 8, 2015 at 10:06 PM, Hong <<a href="mailto:hzhang@mcs.anl.gov">hzhang@mcs.anl.gov</a>> wrote:<br>
> Danniel,<br>
><br>
> Hello, this is Hong, your professor in CS595 :-)<br>
><br>
> We are working on a project for simulating water cycle. The math model to be started is a shallow water river, very similar to your project. I start looking at your codes. After small update with latest petsc-dev, I'm able to run it.<br>
><br>
> Look at the solutions at each time step, I found that solution vector x generated in main.c are alway a constant vector, x=1.<br>
> Am I doing the right thing?<br>
><br>
> Hong<br>
><br>
><br>
><br>
> --<br>
> D. Adrian Maldonado, PhD Candidate<br>
> Electrical & Computer Engineering Dept.<br>
> Illinois Institute of Technology<br>
> 3301 S. Dearborn Street, Chicago, IL 60616<br>
><br>
><br>
><br>
><br>
> --<br>
> D. Adrian Maldonado, PhD Candidate<br>
> Electrical & Computer Engineering Dept.<br>
> Illinois Institute of Technology<br>
> 3301 S. Dearborn Street, Chicago, IL 60616<br>
><br>
><br>
><br>
> --<br>
> D. Adrian Maldonado, PhD Candidate<br>
> Electrical & Computer Engineering Dept.<br>
> Illinois Institute of Technology<br>
> 3301 S. Dearborn Street, Chicago, IL 60616<br>
><br>
><br>
><br>
><br>
> --<br>
> D. Adrian Maldonado, PhD Candidate<br>
> Electrical & Computer Engineering Dept.<br>
> Illinois Institute of Technology<br>
> 3301 S. Dearborn Street, Chicago, IL 60616<br>
><br>
><br>
><br>
><br>
> --<br>
> D. Adrian Maldonado, PhD Candidate<br>
> Electrical & Computer Engineering Dept.<br>
> Illinois Institute of Technology<br>
> 3301 S. Dearborn Street, Chicago, IL 60616<br>
><br>
><br>
><br>
><br>
> --<br>
> D. Adrian Maldonado, PhD Candidate<br>
> Electrical & Computer Engineering Dept.<br>
> Illinois Institute of Technology<br>
> 3301 S. Dearborn Street, Chicago, IL 60616<br>
><br>
><br>
><br>
> --<br>
> D. Adrian Maldonado, PhD Candidate<br>
> Electrical & Computer Engineering Dept.<br>
> Illinois Institute of Technology<br>
> 3301 S. Dearborn Street, Chicago, IL 60616<br>
><br>
><br>
><br>
> --<br>
> D. Adrian Maldonado, PhD Candidate<br>
> Electrical & Computer Engineering Dept.<br>
> Illinois Institute of Technology<br>
> 3301 S. Dearborn Street, Chicago, IL 60616<br>
><br>
<br>
</div></div></blockquote></div><br></div></div>