Re: [Libmesh-users] Adjoint Solve

[Kameeko K. <ka...@gm...>]

Hello,

In case there were initialization issues arising from my not having
performed a forward solve first, I added print lines after line 824 in
petsc_linear_solver.C, after the matrices and vectors are closed, to see if
PetscLinearSolver::adjoint_solve was trying to solve the correct linear
system. The matrix_in and rhs_in are what I expect them to be and match the
ones used to manually calculate the residual that turned out to be grossly
nonzero.

Thanks,
Harriet

On Tue, Dec 30, 2014 at 9:47 PM, Kameeko Kiwi <ka...@gm...> wrote:

> Hello,
>
> I'm trying to solve an adjoint problem for a linear system without first
> solving the forward system. libMesh is giving me an adjoint solution which
> results in a grossly non-zero residual (calculated using system.matrix and
> system.get_adjoint_rhs, since the Jacobian is symmetric for now) , O(e-4)
> in L2. However, the final residual output given by adjoint_solve() is
> O(e-11). This occurred on two separate computers, so I don't believe it is
> a product of my particular petsc installation.
>
> If I use the debugging options to print out the jacobian and adjoint rhs,
> copy them into Matlab and calculate the adjoint solution there, I get a
> solution for which the manually calculated residual (calculated in libMesh
> by reading in the Matlab solution from a text file) is much smaller
> (O(e-15)).
>
> If I implement the adjoint system as the forward problem and solve it
> using system.solve(), I get a solution that agrees with that produced by
> Matlab.
>
> It seems that the adjoint is not being solved to the user specified
> accuracy. Shrinking solver->initial_linear_tolerance (currently at e-10)
> did not help. I can send my code if needed.
>
> Thanks,
> Harriet
>