101     type(ex73f90tmodule_type), pointer:: solver
109     PetscCall(SNESGetApplicationContext(mysnes, solver, ierr))
110 … PetscCall(DMCompositeGetAccessArray(solver%da, Xnest, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
112     PetscCall(InitialGuessLocal(solver, Xsub(1), ierr))
118 …PetscCall(DMCompositeRestoreAccessArray(solver%da, Xnest, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ie…
137   subroutine InitialGuessLocal(solver, X1, ierr)  argument
139     type(ex73f90tmodule_type) solver  local
152     hx = one/(solver%mx - 1)
153     hy = one/(solver%my - 1)
154     temp1 = solver%lambda/(solver%lambda + one) + one
159       j = row/solver%mx
160       i = mod(row, solver%mx)
161       temp = min(j, solver%my - j + 1)*hy
162       if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
165         v = temp1*sqrt(min(min(i, solver%mx - i + 1)*hx, temp))
185   subroutine FormJacobian(dummy, X, jac, jac_prec, solver, ierr)  argument
190     type(ex73f90tmodule_type) solver
200     PetscCall(DMCompositeGetAccessArray(solver%da, X, ione, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
203 …PetscCall(MatCreateSubMatrix(jac_prec, solver%isPhi, solver%isPhi, MAT_INITIAL_MATRIX, Amat, ierr))
205     PetscCall(FormJacobianLocal(Xsub(1), Amat, solver, .true., ierr))
207 … PetscCall(DMCompositeRestoreAccessArray(solver%da, X, ione, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
238   subroutine FormJacobianLocal(X1, jac, solver, add_nl_term, ierr)  argument
240     type(ex73f90tmodule_type) solver  local
258     hx = one/(solver%mx - 1)
259     hy = one/(solver%my - 1)
260     sc = solver%lambda
269       j = irow/solver%mx
270       i = mod(irow, solver%mx)
273       if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
287         if (i + 1 == solver%mx - 1) v(4) = 0.0
289         if (j + 1 == solver%my - 1) v(5) = 0.0
290         col(1) = irow - solver%mx
294         col(5) = irow + solver%mx
317   subroutine FormFunction(snesIn, X, F, solver, ierr)  argument
322     type(ex73f90tmodule_type) solver
334     PetscCall(DMCompositeGetAccessArray(solver%da, X, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
335     PetscCall(DMCompositeGetAccessArray(solver%da, F, itwo, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))
337     PetscCall(FormFunctionNLTerm(Xsub(1), Fsub(1), solver, ierr))
338     PetscCall(MatMultAdd(solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr))
341     PetscCall(MatMult(solver%Cmat, Xsub(1), Fsub(2), ierr))
342     PetscCall(MatMultAdd(solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr))
343     PetscCall(MatMultAdd(solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr))
345 … PetscCall(DMCompositeRestoreAccessArray(solver%da, X, itwo, PETSC_NULL_INTEGER_ARRAY, Xsub, ierr))
346 … PetscCall(DMCompositeRestoreAccessArray(solver%da, F, itwo, PETSC_NULL_INTEGER_ARRAY, Fsub, ierr))
364   subroutine FormFunctionNLTerm(X1, F1, solver, ierr)  argument
366     type(ex73f90tmodule_type) solver  local
375     sc = solver%lambda
384       j = irow/solver%mx
385       i = mod(irow, solver%mx)
388       if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
434   type(ex73f90tmodule_type) solver
442   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, solver%rank, ierr))
447   solver%lambda = 6.0
452 …PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', solver%lambda, fl…
453 …PetscCheckA(solver%lambda <= lambda_max .and. solver%lambda >= lambda_min, PETSC_COMM_SELF, PETSC_…
464 …PetscCallA(DMDAGetInfo(daphi, PETSC_NULL_INTEGER, solver%mx, solver%my, PETSC_NULL_INTEGER, PETSC_…
465   N1 = solver%my*solver%mx
466   N2 = solver%my
493   PetscCallA(MatCreate(PETSC_COMM_WORLD, solver%AmatLin, ierr))
494   PetscCallA(MatSetSizes(solver%AmatLin, PETSC_DECIDE, PETSC_DECIDE, N1, N1, ierr))
495   PetscCallA(MatSetUp(solver%AmatLin, ierr))
497   PetscCallA(FormJacobianLocal(x1, solver%AmatLin, solver, .false., ierr))
498   PetscCallA(MatAssemblyBegin(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))
499   PetscCallA(MatAssemblyEnd(solver%AmatLin, MAT_FINAL_ASSEMBLY, ierr))
536     bval(1) = -one/(solver%mx - 2)
540       j = irow/solver%mx   ! row in domain
541       i = mod(irow, solver%mx)
544       if (i == 0 .or. j == 0 .or. i == solver%mx - 1 .or. j == solver%my - 1) then
586   PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD, solver%da, ierr))
587   PetscCallA(DMCompositeAddDM(solver%da, daphi, ierr))
588   PetscCallA(DMCompositeAddDM(solver%da, dalam, ierr))
589   PetscCallA(DMSetFromOptions(solver%da, ierr))
590   PetscCallA(DMSetUp(solver%da, ierr))
591   PetscCallA(DMCompositeGetGlobalISs(solver%da, isglobal, ierr))
592   solver%isPhi = isglobal(1)
593   PetscCallA(PetscObjectReference(solver%isPhi, ierr))
594   solver%isLambda = isglobal(2)
595   PetscCallA(PetscObjectReference(solver%isLambda, ierr))
598   solver%Amat = Amat
599   solver%Bmat = Bmat
600   solver%Cmat = Cmat
601   solver%Dmat = Dmat
609   PetscCallA(DMCompositeRestoreGlobalISs(solver%da, isglobal, ierr))
619   PetscCallA(SNESSetDM(mysnes, solver%da, ierr))
621   PetscCallA(SNESSetApplicationContext(mysnes, solver, ierr))
623   PetscCallA(SNESSetDM(mysnes, solver%da, ierr))
626   PetscCallA(SNESSetFunction(mysnes, r, FormFunction, solver, ierr))
628     PetscCallA(SNESSetObjective(mysnes, MyObjective, solver, ierr))
630   PetscCallA(SNESSetJacobian(mysnes, KKTmat, KKTmat, FormJacobian, solver, ierr))
649   if (solver%rank == 0) then
663   PetscCallA(MatDestroy(solver%AmatLin, ierr))
664   PetscCallA(ISDestroy(solver%isPhi, ierr))
665   PetscCallA(ISDestroy(solver%isLambda, ierr))
673   PetscCallA(DMDestroy(solver%da, ierr))