Matrix-vector multiplication in weak formulation

[andryr (Andry Rafaralahy)]

Hello everyone

I am having issues trying to solve a problem that has a matrix-vector product in its weak formulation.
I defined a macro epsilon as follows :

 macro epsilon(u1, u2) [dx(u1), dy(u2), dy(u1)+dx(u2)] //

and there is a term in the weak formulation that looks like this :

int2d(Th)(
        epsilon(uu2, vv2)' * C * epsilon(w, s) 
    )

where C is a 3x3 real matrix. The issue is that this code won’t compile and this is the error message I get

   69 :         (C * epsilon(uu2, vv2)     [dx(uu2), dy( vv2), dy(uu2)+dx( vv2)] )
 error operator *  <14Matrice_CreuseIdE>, <10LinearCombI7MGauche4C_F0E> 

However, when I use remove the matrix C from the weak formulation it compiles and run fine.

int2d(Th)(
        epsilon(uu2, vv2)' * epsilon(w, s) 
    )

Any help would be appreciated

Thanks

[prj]

Sadly, I don’t think it’s quite possible. However, you can trick the varf with a hand-defined macro.

mesh Th;
real[int, int] C(2,2);
macro Cdense [[C(0,0), C(0,1)],[C(1,0), C(1,1)]]//
varf vPb(u, v)= int2d(Th)(
        [2*v, v]'*Cdense*[3*u, u]
    );

[andryr (Andry Rafaralahy)]

It works ! Thank you very much.

[Lee]

Hi Andry

I have faced a similar problem like you. But the elements of C matrix in my code are all functions of other field variables, like C(i,j) = Cij(m,n), here m and n are known field variables. And I use
func C=[ [C11,C12],[C21,C22] ];
to define C, but the same error of yours appears.

Could you please lend me a hand with this? Thank you very much.

Lee

[Lee]

Dear prj

would you help me out with the similar problem? the short description is presented in this discussion part.

And by the way, the varf equation I use is written in a very long formate and it takes a very long time to finish the computation of matrix A=equ(Th,Th).
I want to retrive its matrix calculation form to speed up it. Could you also give me some suggestions on how to speed up the process?

Thanks a lot.

[prj]

You cannot use a func, use a macro instead.
For speeding up your program, just use PETSc, see one of the many examples + this tutorial.
Basically, replace

mesh Th;
func Pk = P1;
varf vPb(u, v) = int2d(Th)(...);
fespace Vh(Th, Pk);
matrix A = equ(Vh, Vh);
real[int] rhs = equ(0, Vh);
real[int] prod = A * rhs;
set(A, solver = sparsesolver);
real[int] sol = A^-1 * rhs;

By

mesh Th;
func Pk = P1;
load "PETSc"
Mat A;
createMat(Th, A, Pk)
varf vPb(u, v) = int2d(Th)(...);
fespace Vh(Th, Pk);
Mat A = equ(Vh, Vh);
real[int] rhs = equ(0, Vh);
real[int] prod = A * rhs;
set(A, sparams = "-pc_type lu");
real[int] sol = A^-1 * rhs;

You should get near-linear speedup for the assembly phase, and super-linear speedup for the solution phase (at least up until ~4 to 8 processes).

[Lee]

Dear prj

Thanks for your reply, I will try PETSc.

[Lee]

Dear prj

My problem is a constrained non-linear problem, It works well using IPOPT provided in FreeFem++.
But I have not found the parallel version of IPOPT in FreeFem++ so it take a very long time in sequence form, especially the formation of matrix HJ .
Now I have exactly the form of J, DJ and HJ, can my problem be able to transformed using PETSc? Can you give me a general guideline for this transform process?

Thanks a lot.

[prj]

Sure! You can use Tao, the Toolkit for Advanced Optimization, which is part of PETSc. The interface ressembles what you’d do with Ipopt, but of course, everything is parallel
You can have a look at one of these examples [1,2,3]. Also, example12.edp from this tutorial.

[Lee]

Thank you very much.

[Lee]

Dear prj

I have tried to read the example https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp using TAO and run it on my PC. I am confused by the problems below:

1. The output massege always prints “Solver terminated: -6 Line Search Failure” after several iterations. Does it mean the solving process is unsuccessful? If so, what should I do to make the line search successful?

2. Where can I found the detail explaination of the meaning of some functions like “changeNumbering”? I am not sure my guess is ture or not. This makes me confused on retrieving the desired global field variables.

3. My own problem is a little bit different from normal constrained nonlinear optimization, DJ is not the direct partial derivative of J and can not form a exact corresponding J. In my original IPOPT code, I must turn off the linesearch function (with the code “linesearch=false” in IPOPT) to get the results or it will trapped in an endless line search process. How should I deal with this problem with TAO?

Thanks for your reply.

[prj]

1. sorry, there was a bug in minimal-surface-Tao-2d-PETSc.edp, I fixed it in this commit. With quasi-Newton, you should now get Solver terminated: -2 Maximum Iterations, and with Newton Solution converged: ||g(X)|| <= gatol.
2. please have a look a section 8 of this tutorial. Also, you can start with a sequential example, and we’ll parallelize it afterwards. There is no need to bother with more than one process in the beginning.
3. you can try one of these methods, surely all of them don’t do linesearch. If you formulate exactly your problem, I may gave you a better answer.

[Lee]

Dear prj

Thank you for your timely reply. I have managed to greatly simplify my problem and I tried to rewrite it in the same form of [https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/minimal-surface-Tao-2d-PETSc.edp].
The original problem using IPOPT and the rewrite one are uploaded as follows.
Original_Ipopt.edp (3.1 KB) TAO_test.edp (4.3 KB)

The problem I want to solve is provided in Original_Ipopt.edp. It is a simplified bound constrained-nonlinear-time dependent-multiple variable problem. It works fine using IPOPT. But I have not found a TAO example that meet all the requirement I need. So I tried to rewrite it using TAO as in TAO_test.edp.

In the TAO_test.edp I uploaded, an error is encountered in the line 28 as createMat(Th, H, Pk);
The message says that “meshN The Identifier meshN does not exist”. How to correct it?

Could you please point out the incorrect parts in the TAO_test.edp and give me your valuable suggestions on how to correct them?

Expecting your reply.

[prj]

Here, I’ve fixed part of the script:

1. macro def and init have to be scoped
2. macro dimension must be defined before including macro_ddm.idp
3. can’t use both plotMPI and fespace Xh (I’ll fix that in FreeFEM)
4. some functions were not defined at the right place

The code runs, but I don’t know if it’s working properly. TAO_test.edp (4.3 KB)

[Lee]

Thank you very much

[Lee]

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.

[prj]

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.

[Lee]

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.

[prj]

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.

[Lee]

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.