From: Kameeko K. <ka...@gm...> - 2015-01-07 03:56:30
|
Hello, I've attached a simplified version of my code that reproduces the problem I'm having, namely that the adjoint calculated has a very nonzero residual when the residual is calculated manually from the transpose of system.matrix and adjoint_rhs. The simplified system has two variables; I've not yet been able to replicate the phenomenon with a single-variable system, though I am not sure whether that is a coincidence or related to the problem. Thanks, Harriet On Sun, Jan 4, 2015 at 2:54 PM, Kameeko Kiwi <ka...@gm...> wrote: > 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 >> > > |