Thank you very much 
Dear prj
I have run the TAO_test.edp, but the it stucked at line 36 and 37 as
Line 36 plot(part, fill=1, value=1,wait=1,cmm=“part befor creatPartition”);
Line 37 createPartition(Th, part[], P0)
These two are used in the build up block of ThNo. The output is as below:
— partition of unity built (in 2.956960e-02)
— global numbering created (in 3.081115e-03)
— global CSR created (in 4.659305e-03)
Warning: May be a bug in your script,
a part of the plot is wrong t (mesh or FE function, curve) => skip the item 1 in plot command
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
Error plot item empty 0
The message “Error plot item empty 0” keep showing up and the next plot function in line 37did not happen.
Could you please give me hint on how to fix it?
Thank you very much.
Sorry, you need to add buildDmesh(Th) after the plot line 21. Then, Th will be distributed and both createMat and createPartition will work as expected.
It works! And the results seem to be the same as the Ipopt one I uploaded and the speed is extremely fast!
Thank you very much!
Do you have a full tutorial for the TAO in FreeFem++ to explain the functions like “createPartition”? The “Miscellaneous” part is FreeFem++ Doc is not clear enough for me when compared with explaination of IPOPT.
Thank you again.
Tutorial about FreeFEM + parallelism is available there, section 8, read the PDF about PETSc, there are some explanation. If something is not clear, please let me know.
List of Tao solvers are available there, if you see something you like that is not interfaced yet in FreeFEM, please let me know.
Dear prj
I have implemented my complete problem using TAO, but the due to I can not find a correct optimizing function J for my DJ and HJ, almost each Newton iteration needs very long CG iteration computing J and DJ and can not reach convergence.
In my original IPOPT program, I turned off the linesearch function and the J and DJ would be computed only once in each Newton iteration. In FreeFem++4.6 doc, it says that when linesearch is turned off, the method becomes a standard Newton algorithm instead of a primal-dual system(as in page 234). Would you please give me some suggestions on how to achieve this using TAO? I have tried all the possible TAO solvers for bound-constrained nonlinear problem and the problem continue to show up.
I also tried the SNESSolve in PETSc, since it do not need a J. But I have questions on the meaning of parameters bPETSc and xPETSc in the SNESSolve implement line below:
SNESSolve(A, funcJ, funcRes, bPETSc, xPETSc,xl = xlPETSc, xu = xuPETSc, sparams = “-snes_monitor -ksp_converged_reason -snes_view -snes_vi_monitor -snes_type vinewtonrsls -snes_rtol 1.0e-6 -pc_type lu”);
It seems bPETSc is the initial value and xPETSc is the output, but I found that before entering SNESSolve, some use xPETSc=bPETSc, some use xPETSc=0 and some use xPETSc=0.1 in the provided examples using SNESSolve. Even for the same example, xPETSc=0.1 and xPETSc=0.2 would leads to different solutions. So would you please also tell me the exact meaning of this two parameters and their values? And by the way, the boundary condition is applied using the same way as TAO in a bound-constrained form, is this correct in SNESSolve or must be applied in the varf form?
Thank you very much for your reply.
I’m not a TAO expert, you’ll have to try different methods. Also, I still don’t know what you are trying to solve exactly, so it’s difficult for me to give you precise pointers. SNESSolve either solves F(x) = 0, or F(x) = b. You are free to put what you want in the output vector, but this acts as an initial guess for nonlinear solves, so it’s best to use something not too far from the solution otherwise your Newton may never converge.
Thank you for your reply, prj.
Here I uploaded my problem using SNESSolve, all the macro and constants are deleted to reduce the program from 2300 lines to 270+ lines. This is a nonlinear problem F(x)=0, the F is given in funcRes part and the correct Jacobian matrix is obtained in funcJ part. Their forms have been proved in my original IPOPT program so don’t worry about that. The solution vector is [u1,u2,u3,eta,eta0], boundary condition is given for u1, u2 and u3, bound Inequality constrain [0,1] is given for eta and eta0. Also, the eta and eta0 are time dependent so their values of previous time step is defined as etaOLD and eta0OLD.PETSc_Example.edp (18.0 KB)
Would you please help me point out the incorrect part in this program? Thanke you very much. I’m nont sure whether I transform the code in a correct way.
Your script doesn’t even run.
Sorry, I don’t have a simplified version that can run, sorry.
No, I mean that your script is wrong, I can’t even launch it on my machine.
Is the way of adding boundary condition in SNESSolve the same with TAO? If the initial value of boundary nodes are set to meet the requirements of boundary condition, can it be maintained using the constrain of upper and lower bound or using the form as on(1,2,3,4, du1=0,du2=0,du3=0) in varf function?
Thank you for your reply.
The boundary conditions have to be set in your initial guess.
Thans for your timely reply.
I have found that in the definition of funcRes and funcJ, “tgv” is used in the vector and matrix forming such as “Loc = vPb(Wh, Wh, tgv = -2)” in SNESSolve. The tgv in FreeFem++ is always represent 1e30 to implement Dirichlet boundary conditions and in most situation can be neglected.
So is it necessary to define funcRes and funcJ with the form “Loc = vPb(Wh, Wh, tgv = -2)” or with the form in which tgv can be neglected as “Loc = vPb(Wh, Wh)”?
If correct form is Loc = vPb(Wh, Wh, tgv = -2), which value of tgv should be used, -1, -2 or 1e30 and why?
Thank you.
I always advise against using the default value of 10e+30. -1 (or -2 if your problem has homogeneous boundary conditions and is symmetric) is always better when you do iterative algorithms, such as a Newton loop in SNESSolve/TaoSolve.
Thanke you. I’m now using tgv=-2 in a SNESSolve iteration.
Dear prj
Here I have uploaded a runable simplified version of my problem.
This is a nonlinear problem F(x)=0, the F is given in funcRes part and the correct Jacobian matrix is obtained in funcJ part. Their forms have been proved in my original IPOPT program so don’t worry about that. The solution vector is [u1,u2,u3,eta,eta0], boundary condition is given for u1, u2 and u3, bound Inequality constrain [0,1] is given for eta and eta0. Also, the eta and eta0 are time dependent so their values of previous time step is defined as etaOLD and eta0OLD.
Could you please point out the errors it contains and give me some suggestions on how to correct it? My own version does not run well and can not get the correct results as IPOPT does,the boundary constrain are some times violated in the line search iterations in Newton computation with the default solver option as " -snes_monitor -ksp_converged_reason -snes_view -snes_vi_monitor -snes_type vinewtonrsls -snes_max_it 1000 -snes_rtol 1.0e-6 -pc_type lu ". And moreover, if I turn off the line search using -snes_linesearch_type basec or change the preconditioner, the results converge would fail.
PETSc_exa.edp (103.6 KB)
Thank you very much.
Your script has 1.5k source lines of code. Start with something small, if it doesn’t give you what you want, I’ll help you debug. Then, incrementally add complexity on top of functioning code.
Thank you. I will try that.
Do you know how to use the PDIPM solver (-tao type pdipm) in TAO? I found it may be the same solver used in IPOPT, but there is not much information about it. I have try several time with the command
“TaoSolve(H, J, DJ, uPETSc, xl=lbPETSc, xu=ubPETSc, sparams="-tao_monitor -tao_view -tao_type pdipm -tao_max_it 1000 -tao_gatol 1e-4 -pc_type lu", HessianRoutine=HJ)”
but failed, the iteration can not even start.
Does it fail on a complicated problem, or on a very basic one? If it’s again on your 1.5k lines of code solver, then it’s not surprise, I guess you have other issues that you need to figure out first. And that’s why I recommend that you start with a much simpler version of your solver and gradually transition into the complete one.