[petsc-users] What is the difference between these two approaches?
(Rebecca) Xuefei YUAN
xy2102 at columbia.edu
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
More information about the petsc-users
mailing list