Matrix-vector multiplication in weak formulation

General Discussion
41

Thank you for timely reply.

When BNLS is used in my original problem, it’s very hard to converge in the second time step iteration. The residual output information indicate the iterations seems stuck in a small interval and repeat again and again. This phenomenon also appears in IPOPT, but once I turn off line search with command linesearch=false, the IPOPT can give correct results.

The results for the first time step is almost the same with IPOPT. Since IPOPT uses the primal-dual interior point method, I think maybe ipm or PDIPM may help to retrieve the correct IPOPT results. If PDIPM is not readily available, how about ipm?

Meanwhile, I’m trying other Tao solvers, but till now, the problem is not solved.

42

You can try both IPM and PDIPM, I was just implying that you may run into a bug. But there are indeed many available options, I don’t really know what you are solving and I’m not really a Tao expert, so I’m not very useful.

43

I tried to use IPM in a much smaller bound constrained example, but it can not run. So I uploaded the example (as the TAOexample_use_ipm.edp above) if you can help me point out the error.

Thank you very much.

44

Dear prj

Here I have build a small example using tao_ipm solver. But the error continue to show in the cmd window as:

“[0]IPETSC ERROR: Caught signal number 11 SEGV: Segementation Violation, probably memory access out of range”

I totally have no idea how to fix it. Could you please help me correct it? The edp is uploaded as follows.
ipmExa.edp (5.9 KB)

Thank you very much.

45

The equality constraint also appears in the Hessian, so you need to change your HJ to something like:

func int HJ(real[int]& X, real[int]& E) {
[...]
}
46

The problem only has simple inequality constraint as xlPETSc and xuPETSc. But if I delete the funcE and funcJE, how should I write the TaoSolve information?

If I write it as follows:
TaoSolve(H, J, DJ, uPETSc, xl = lbPETSc, xu = ubPETSc, sparams = “-tao_monitor -tao_view -tao_type ipm -pc_type lu -pc_factor_mat_solver_type mumps -tao_max_it 1000 -tao_gatol 1e-4”, HessianRoutine = HJ);
The error I mentioned continues to appear.

47

I can’t reproduce the previous error without equality constraint.

48

The error message is:
[0] PETSC EROOR: Out of memory. This could be due to allocating too large an object or bleeding by not properly destroyubg unneeded objects.
It still can not run properly for multiple processors? How should I correct it?

49

PDIPM and IPM are bugged when using no inequality or equality constraint. This will be fixed in the next PETSc release. I suggest you use other Tao solvers in the meantime.

50

Thank you for your reply.

Could you recommend some Tao solver which has the same function as interior point method?

51

No, I’d suggest you try different alternative. You can tune the line search parameters by looking at the PETSc/Tao documentation that I’ve already pasted above.

52

I try to reformulate the example in Transmission Problem (Transmission problem) into the type where the diffusion coefficient is a matrix, which depends on the region. When the diffusion coefficient is constant matrix, the results are satisfying. When I add region-dependent formulation the problem diverges, even if the matrix has the same values in both regions. Here is the code I tried (works with line 15, doesn’t work with line 17):

border a(t=0, 1){x=t; y=0;};
border b(t=0, 1){x=1; y=t;};
border c(t=0, 1){x=1-t; y=1;};
border d(t=0, 1){x=0; y=1-t;};
border i1(t=0, 1){x=t; y=0.5;};
mesh th = buildmesh(a(10) + b(10) + c(10) +d(10)+i1(10));
fespace Ph(th, P2);

Ph reg=region;
int nupper = reg(0.25, 0.75);
int nlower = reg(0.25, 0.25);

matrix A = [[1,0],[0,0.01]];
//this part works
// matrix nu = A;
//this part doesnt
matrix nu = A*(region==nlower) + A*(region==nupper);

Ph u,v;
solve lap (u, v)
    = int2d(th)(
         [dx(v), dy(v)]'*([[nu(0,0), nu(0,1)],[nu(1,0), nu(1,1)]]*[dx(u), dy(u)])
    )
    + int2d(th)(
        -1*v
    )
    + on(a, b, c, d, u=0)
    ;
plot(u);

I was able to make it converge by changing line 17 but the results are the same as with constant matrix in both regions :

matrix B = [[0.01,0],[0,0.01]];
matrix nu = B*(region!=nlower&&region!=nupper)+ A*(region==nlower);