Hello everyone, I am running into an issue with SNES solver where the residual in the first iteration gets calculated fine, but on the 2nd iteration, even though it is getting calculated locally, it doesn't seem to be "assigned" to the residual vector correctly. The residual vector here is referred to as petsc_PVector which is defined in a type (struct). It seems like the first iteration, it behaves as expected and that petsc_PVector gets updated, 2nd iteration however seems like there are 2 versions of it; one that is being passed into FormFunctionU, and another that exists in the type (struct). This result in the residual on the 2nd iteration to be exactly = 0 and convergence. Any insights, suggestions, or experiences you can share regarding this issue would be immensely appreciated. I'm particularly interested in understanding why the residual vector behaves differently in the second iteration. Best, Khaled The following is a snapshot of how I initialize and use SNES: subroutine set_function_jacobian_to_solver_u(nlU) type(variable_node_u),target,intent(inout)::nlU PetscErrorCode::ierr type(variable_node_part),pointer::var_nd_part integer::ipart ! external FormFunctionU,FormJacobianU do ipart=1,nlU%variable_node%npart var_nd_part=>nlU%variable_node%parts(ipart) PetscCallA(SNESSetFunction(var_nd_part%petsc_snes,var_nd_part%petsc_PVector,FormFunctionU,nlU,ierr)) PetscCallA(SNESSetJacobian(var_nd_part%petsc_snes,var_nd_part%petsc_KMatrix,var_nd_part%petsc_KMatrix,FormJacobianU,nlU,ierr)) PetscCallA(SNESSetFromOptions(var_nd_part%petsc_snes,ierr)) end do var_nd_part=>null() end subroutine set_function_jacobian_to_solver_u subroutine FormFunctionU(petsc_snes,petsc_vector,petsc_PVector,nlU,ierr) type(tSNES)::petsc_snes type(tVec)::petsc_vector type(tVec)::petsc_PVector type(PetscScalar),pointer::vec_ptr_a(:),vec_ptr_b(:) PetscErrorCode::ierr type(variable_node_u),target::nlU type(PetscScalar),pointer::val_petsc_PVector(:) type(parameter),pointer::para_data type(mesh),pointer::mesh_data type(element),pointer::eles(:) integer::k,i,flag,pid PetscInt::its para_data=>nlU%para_data mesh_data=>nlU%mesh_data eles=>nlU%eles call MPI_Comm_rank(MPI_COMM_WORLD, pid, ierr) call update_data_by_petsc_vector(nlU%variable_node,petsc_vector) call set_petsc_PVector_to_zero(nlU%variable_node) call update_petsc_PVector(nlU%variable_node) flag = 0 do k=1,mesh_data%nel call get_element_value(nlU%data_e, eles(k)%nnd_e, mesh_data%elements(:,k), nlU%nvar_total, nlU%data) call cal_element_ResidualU(nlU%PVector_e, eles(k), nlU%data_e, para_data) call set_petsc_PVector_values_by_element(nlU%variable_node, nlU%PVector_e, eles(k)%nnd_e, mesh_data%elements(:,k)) end do call update_petsc_PVector(nlU%variable_node) para_data=>null() mesh_data=>null() eles=>null() return end subroutine FormFunctionU call SNESSolve(nlU%variable_node%petsc_snes,PETSC_NULL_VEC, nlU%variable_node%petsc_vector,error_petsc)