[petsc-users] What is the difference between these two approaches?

Tue Feb 9 12:56:52 CST 2010

Dear Jed,

'Perhaps a wrong or low-accuracy Jacobian?  Are you using an analytic
Jacobian?  Does the history change if you tighten the linear solve
tolerance?  If the system is very ill-conditioned, you might only be
able to converge this far.'

There is no analytic Jacobian for this code(approach 1), so I ran with  
-snes_mf or -dmmg_jacobian_mf_fd options. Both give me similar  
behaviors of residual history. The problem is solved via approach 2,  
and the residual is able to go down as small as 1e-10.

>> Dear all,
>>
>> I do have an PDE-constrained optimization problem:
>>
>> F1(phi(x,y),c) = L(phi(x,y)) + rho(s1(x,y),s2(x,y)) = 0 (1)
>> F2(phi(x,y),c) = int_{domain}(c*rho(s1(x,y),s2(x,y)) - 1.0)dxdy = 0 (2)
>>
>> where phi(x,y) is defined at each grid point and c is a scalar
>> parameter that satisfying the equation (2).
>
> Is c the only parameter you are optimizing?  What does your system look
> like in "minimize functional subject to PDE" form?

Yes, c is the only parameter. The "optimizing" shows up here because  
there is a parameter c in F1, and c satisfies F2. Actually I am  
solving F1 and need F2 to be satisfied, thus I can solve F1&F2 as in  
approach 1(in this case, F2 is not strictly satisfied) together or  
solve F1 first and then get c from F2 as in approach 2(in this case,  
F2 is strictly satisfied).

>> I have two pieces of codes to solve them by different approaches.
>>
>> The first approach is to use DMComposite() to manage unknowns and
>> solve F1&F2 at the same time, and it calls DMMGSolve() once.
>>
>> The second approach is to solve F1 with a guess c for some number(i.e.
>> 1) times of nonlinear and linear iteration, and then update c from F2.
>> Solve F1 again with updated c till some conditions satisfied, for
>> example,
>
Perhaps the iteration has stagnated or found a different local minimum?
Have you compared the objective functions?
I thought about Null space for the problem. Because in approach 2, a  
constant vector is a null space of the solution, so I add

ierr = DMMGSetNullSpace(dmmg, PETSC_TRUE, 0, PETSC_NULL);CHKERRQ(ierr);

but for approach 1, I wrote a CreateNullSpace() as in

ierr = DMMGSetNullSpace(dmmg, PETSC_FALSE, 1, CreateNullSpace);CHKERRQ(ierr);

Will this be the source of difference?

>> d) the shortcoming of approach one is PETSc has not ready for assemble
>> a Jacobian for object DMComposite, so the second approach is able to
>> write an analytic Jacobian.
>
> Does DMCompositeGetMatrix() not work for you?

I tried DMCompositeGetMatrix() last July, the conversations are:

Hi,

I am working on an optimization problem, in which I would like to  
assemble a Jacobian matrix. Thus  
DMMGSetSNES(dmmg,FormFunction,FormJacobian) is called.

In damgsnes.c:637, in calling DMGetMatrix(), it calls  
DMCompositeGetMatrix() where the temp matrix Atmp has been freed  
before it passes any information to J at pack.c:1722 and 1774.

So after calling DMGetMatrix() in DMMGSetSNES, the stencil of the  
dmmg[i]->B has unchanged, i.e.,

(gdb) p dmmg[0]->B->stencil
$107 = {dim = 0, dims = {0, 0, 0, 0}, starts = {0, 0, 0, 0}, noc =  
PETSC_FALSE}
(gdb) where
#0  DMMGSetSNES (dmmg=0x8856208, function=0x804c84f <FormFunction>,
    jacobian=0x8052932 <FormJacobian>) at damgsnes.c:641
#1  0x0804c246 in main (argc=Cannot access memory at address 0x0
) at tworeggt.c:126

I compare this with
http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/snes/examples/tutorials/ex18.c.html

and it shows that the stencil has been carried out and passed to  
dmmg[0]->B as follows:

(gdb) p dmmg[i]->B->stencil
$80 = {dim = 2, dims = {5, 5, 1, 0}, starts = {0, 0, 0, 0}, noc = PETSC_TRUE}
(gdb) where
#0  DMMGSetSNES (dmmg=0x884b530, function=0x804c364 <FormFunction>,
    jacobian=0x804d34d <FormJacobian>) at damgsnes.c:642
#1  0x0804b969 in main (argc=Cannot access memory at address 0x2
) at ex18.c:100

Because of this missing stencil of Jacobian matrix, I get the error  
code as follows:
Program received signal SIGSEGV, Segmentation fault.
0x082447c2 in ISLocalToGlobalMappingApply (mapping=0x0, N=1, in=0xbff8f250,
    out=0xbff8ce14) at /home/rebecca/soft/petsc-3.0.0-p1/include/petscis.h:129
129       PetscInt i,*idx = mapping->indices,Nmax = mapping->n;
(gdb) where
#0  0x082447c2 in ISLocalToGlobalMappingApply (mapping=0x0, N=1,
    in=0xbff8f250, out=0xbff8ce14)
    at /home/rebecca/soft/petsc-3.0.0-p1/include/petscis.h:129
#1  0x0824440c in MatSetValuesLocal (mat=0x88825e8, nrow=1, irow=0xbff8f250,
    ncol=4, icol=0xbff8ee50, y=0xbff8f628, addv=INSERT_VALUES) at matrix.c:1583
#2  0x08240aae in MatSetValuesStencil (mat=0x88825e8, m=1, idxm=0xbff8f6b8,
    n=4, idxn=0xbff8f4b4, v=0xbff8f628, addv=INSERT_VALUES) at matrix.c:1099
#3  0x08053835 in FormJacobian (snes=0x8874700, X=0x8856778, J=0x88747d0,
    B=0x88747d4, flg=0xbff8f8d4, ptr=0x8856338) at tworeggt.c:937
#4  0x0805a5cf in DMMGComputeJacobian_Multigrid (snes=0x8874700, X=0x8856778,
    J=0x88747d0, B=0x88747d4, flag=0xbff8f8d4, ptr=0x8856208) at damgsnes.c:60
#5  0x0806b18a in SNESComputeJacobian (snes=0x8874700, X=0x8856778,
    A=0x88747d0, B=0x88747d4, flg=0xbff8f8d4) at snes.c:1111
#6  0x08084945 in SNESSolve_LS (snes=0x8874700) at ls.c:189
#7  0x08073198 in SNESSolve (snes=0x8874700, b=0x0, x=0x8856778) at  
snes.c:2221
#8  0x0805d5f9 in DMMGSolveSNES (dmmg=0x8856208, level=0) at damgsnes.c:510
#9  0x08056e38 in DMMGSolve (dmmg=0x8856208) at damg.c:372
#10 0x0804c3fe in main (argc=128, argv=0xbff90c04) at tworeggt.c:131

I think there might be a bug in DMCompositeGetMatrix().

And Barry said something about it as

'You cannot use the stencil operations to put values into a "composite  
matrix". The numbering of rows and columns of the composite matrix  
reflect all the different variables (unknowns) sp do not match what  
they are for a single component.'

I did not quite get it at that time...

Thanks so much!

Rebecca
> Jed
>
>

-- 
(Rebecca) Xuefei YUAN
Department of Applied Physics and Applied Mathematics
Columbia University
Tel:917-399-8032
www.columbia.edu/~xy2102