1c4762a1bSJed Brown! 2c4762a1bSJed Brown! Description: Solves a nonlinear system in parallel with SNES. 3c4762a1bSJed Brown! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular 4c4762a1bSJed Brown! domain, using distributed arrays (DMDAs) to partition the parallel grid. 5c4762a1bSJed Brown! The command line options include: 6c4762a1bSJed Brown! -par <parameter>, where <parameter> indicates the nonlinearity of the problem 7c4762a1bSJed Brown! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81) 8c4762a1bSJed Brown! 9c4762a1bSJed Brown! This system (A) is augmented with constraints: 10c4762a1bSJed Brown! 11c4762a1bSJed Brown! A -B * phi = rho 12c4762a1bSJed Brown! -C I lam = 0 13c4762a1bSJed Brown! 14c4762a1bSJed Brown! where I is the identity, A is the "normal" Poisson equation, B is the "distributor" of the 15c4762a1bSJed Brown! total flux (the first block equation is the flux surface averaging equation). The second 16c4762a1bSJed Brown! equation lambda = C * x enforces the surface flux auxiliary equation. B and C have all 17c4762a1bSJed Brown! positive entries, areas in C and fraction of area in B. 18c4762a1bSJed Brown! 19c4762a1bSJed Brown!!/*T 20c4762a1bSJed Brown! Concepts: SNES^parallel Bratu example 21c4762a1bSJed Brown! Concepts: MatNest 22c4762a1bSJed Brown! Processors: n 23c4762a1bSJed Brown!T*/ 24c4762a1bSJed Brown 25c4762a1bSJed Brown! 26c4762a1bSJed Brown! -------------------------------------------------------------------------- 27c4762a1bSJed Brown! 28c4762a1bSJed Brown! Solid Fuel Ignition (SFI) problem. This problem is modeled by 29c4762a1bSJed Brown! the partial differential equation 30c4762a1bSJed Brown! 31c4762a1bSJed Brown! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 32c4762a1bSJed Brown! 33c4762a1bSJed Brown! with boundary conditions 34c4762a1bSJed Brown! 35c4762a1bSJed Brown! u = 0 for x = 0, x = 1, y = 0, y = 1. 36c4762a1bSJed Brown! 37c4762a1bSJed Brown! A finite difference approximation with the usual 5-point stencil 38c4762a1bSJed Brown! is used to discretize the boundary value problem to obtain a nonlinear 39c4762a1bSJed Brown! system of equations. 40c4762a1bSJed Brown! 41c4762a1bSJed Brown! -------------------------------------------------------------------------- 42c4762a1bSJed Brown! The following define must be used before including any PETSc include files 43c4762a1bSJed Brown! into a module or interface. This is because they can't handle declarations 44c4762a1bSJed Brown! in them 45c4762a1bSJed Brown! 46c4762a1bSJed Brown module petsc_kkt_solver 47c4762a1bSJed Brown#include <petsc/finclude/petscdm.h> 48c4762a1bSJed Brown#include <petsc/finclude/petscmat.h> 49c4762a1bSJed Brown use petscdm 50c4762a1bSJed Brown use petscmat 51c4762a1bSJed Brown type petsc_kkt_solver_type 52c4762a1bSJed Brown DM::da 53c4762a1bSJed Brown! temp A block stuff 54c4762a1bSJed Brown PetscInt mx,my 55c4762a1bSJed Brown PetscMPIInt rank 56c4762a1bSJed Brown PetscReal lambda 57c4762a1bSJed Brown! Mats 58c4762a1bSJed Brown Mat::Amat,AmatLin,Bmat,CMat,Dmat 59c4762a1bSJed Brown IS::isPhi,isLambda 60c4762a1bSJed Brown end type petsc_kkt_solver_type 61c4762a1bSJed Brown 62c4762a1bSJed Brown end module petsc_kkt_solver 63c4762a1bSJed Brown 64c4762a1bSJed Brown module petsc_kkt_solver_interfaces 65c4762a1bSJed Brown use petsc_kkt_solver 66c4762a1bSJed Brown 67c4762a1bSJed Brown Interface SNESSetApplicationContext 68c4762a1bSJed Brown Subroutine SNESSetApplicationContext(snesIn,ctx,ierr) 69c4762a1bSJed Brown#include <petsc/finclude/petscsnes.h> 70c4762a1bSJed Brown use petscsnes 71c4762a1bSJed Brown use petsc_kkt_solver 72c4762a1bSJed Brown SNES:: snesIn 73c4762a1bSJed Brown type(petsc_kkt_solver_type) ctx 74c4762a1bSJed Brown PetscErrorCode ierr 75c4762a1bSJed Brown End Subroutine 76c4762a1bSJed Brown End Interface SNESSetApplicationContext 77c4762a1bSJed Brown 78c4762a1bSJed Brown Interface SNESGetApplicationContext 79c4762a1bSJed Brown Subroutine SNESGetApplicationContext(snesIn,ctx,ierr) 80c4762a1bSJed Brown#include <petsc/finclude/petscsnes.h> 81c4762a1bSJed Brown use petscsnes 82c4762a1bSJed Brown use petsc_kkt_solver 83c4762a1bSJed Brown SNES:: snesIn 84c4762a1bSJed Brown type(petsc_kkt_solver_type), pointer :: ctx 85c4762a1bSJed Brown PetscErrorCode ierr 86c4762a1bSJed Brown End Subroutine 87c4762a1bSJed Brown End Interface SNESGetApplicationContext 88c4762a1bSJed Brown end module petsc_kkt_solver_interfaces 89c4762a1bSJed Brown 90c4762a1bSJed Brown program main 91c4762a1bSJed Brown#include <petsc/finclude/petscdm.h> 92c4762a1bSJed Brown#include <petsc/finclude/petscsnes.h> 93c4762a1bSJed Brown use petscdm 94c4762a1bSJed Brown use petscdmda 95c4762a1bSJed Brown use petscsnes 96c4762a1bSJed Brown use petsc_kkt_solver 97c4762a1bSJed Brown use petsc_kkt_solver_interfaces 98c4762a1bSJed Brown implicit none 99c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 100c4762a1bSJed Brown! Variable declarations 101c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102c4762a1bSJed Brown! 103c4762a1bSJed Brown! Variables: 104c4762a1bSJed Brown! mysnes - nonlinear solver 105c4762a1bSJed Brown! x, r - solution, residual vectors 106c4762a1bSJed Brown! J - Jacobian matrix 107c4762a1bSJed Brown! its - iterations for convergence 108c4762a1bSJed Brown! Nx, Ny - number of preocessors in x- and y- directions 109c4762a1bSJed Brown! 110c4762a1bSJed Brown SNES:: mysnes 111c4762a1bSJed Brown Vec:: x,r,x2,x1,x1loc,x2loc 112c4762a1bSJed Brown Mat:: Amat,Bmat,Cmat,Dmat,KKTMat,matArray(4) 113c4762a1bSJed Brown! Mat:: tmat 114c4762a1bSJed Brown DM:: daphi,dalam 115c4762a1bSJed Brown IS:: isglobal(2) 116c4762a1bSJed Brown PetscErrorCode ierr 117c4762a1bSJed Brown PetscInt its,N1,N2,i,j,irow,row(1) 118c4762a1bSJed Brown PetscInt col(1),low,high,lamlow,lamhigh 119c4762a1bSJed Brown PetscBool flg 120c4762a1bSJed Brown PetscInt ione,nfour,itwo,nloc,nloclam 121c4762a1bSJed Brown PetscReal lambda_max,lambda_min 122c4762a1bSJed Brown type(petsc_kkt_solver_type) solver 123c4762a1bSJed Brown PetscScalar bval(1),cval(1),one 124c4762a1bSJed Brown 125c4762a1bSJed Brown! Note: Any user-defined Fortran routines (such as FormJacobian) 126c4762a1bSJed Brown! MUST be declared as external. 127c4762a1bSJed Brown external FormInitialGuess,FormJacobian,FormFunction 128c4762a1bSJed Brown 129c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130c4762a1bSJed Brown! Initialize program 131c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 132c4762a1bSJed Brown call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 133c4762a1bSJed Brown if (ierr .ne. 0) then 134c4762a1bSJed Brown print*,'PetscInitialize failed' 135c4762a1bSJed Brown stop 136c4762a1bSJed Brown endif 137c4762a1bSJed Brown call MPI_Comm_rank(PETSC_COMM_WORLD,solver%rank,ierr);CHKERRA(ierr) 138c4762a1bSJed Brown 139c4762a1bSJed Brown! Initialize problem parameters 140c4762a1bSJed Brown lambda_max = 6.81_PETSC_REAL_KIND 141c4762a1bSJed Brown lambda_min = 0.0 142c4762a1bSJed Brown solver%lambda = 6.0 143c4762a1bSJed Brown ione = 1 144c4762a1bSJed Brown nfour = 4 145c4762a1bSJed Brown itwo = 2 146c4762a1bSJed Brown call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par', solver%lambda,flg,ierr);CHKERRA(ierr) 147c4762a1bSJed Brown if (solver%lambda .ge. lambda_max .or. solver%lambda .lt. lambda_min) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_USER,'Lambda provided with -par is out of range'); endif 148c4762a1bSJed Brown 149c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150c4762a1bSJed Brown! Create vector data structures; set function evaluation routine 151c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 152c4762a1bSJed Brown 153c4762a1bSJed Brown! just get size 154c4762a1bSJed Brown call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nfour,nfour,PETSC_DECIDE,PETSC_DECIDE,ione,ione, & 155c4762a1bSJed Brown & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,daphi,ierr);CHKERRA(ierr) 156c4762a1bSJed Brown call DMSetFromOptions(daphi,ierr);CHKERRA(ierr) 157c4762a1bSJed Brown call DMSetUp(daphi,ierr);CHKERRA(ierr) 158c4762a1bSJed Brown call DMDAGetInfo(daphi,PETSC_NULL_INTEGER,solver%mx,solver%my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 159c4762a1bSJed Brown & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr) 160c4762a1bSJed Brown N1 = solver%my*solver%mx 161c4762a1bSJed Brown N2 = solver%my 162c4762a1bSJed Brown flg = .false. 163c4762a1bSJed Brown call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_constraints',flg,flg,ierr);CHKERRA(ierr) 164c4762a1bSJed Brown if (flg) then 165c4762a1bSJed Brown N2 = 0 166c4762a1bSJed Brown endif 167c4762a1bSJed Brown 168c4762a1bSJed Brown call DMDestroy(daphi,ierr);CHKERRA(ierr) 169c4762a1bSJed Brown 170c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 171c4762a1bSJed Brown! Create matrix data structure; set Jacobian evaluation routine 172c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 173c4762a1bSJed Brown call DMShellCreate(PETSC_COMM_WORLD,daphi,ierr);CHKERRA(ierr) 174c4762a1bSJed Brown call DMSetOptionsPrefix(daphi,'phi_',ierr);CHKERRA(ierr) 175c4762a1bSJed Brown call DMSetFromOptions(daphi,ierr);CHKERRA(ierr) 176c4762a1bSJed Brown 177c4762a1bSJed Brown call VecCreate(PETSC_COMM_WORLD,x1,ierr);CHKERRA(ierr) 178c4762a1bSJed Brown call VecSetSizes(x1,PETSC_DECIDE,N1,ierr);CHKERRA(ierr) 179c4762a1bSJed Brown call VecSetFromOptions(x1,ierr);CHKERRA(ierr) 180c4762a1bSJed Brown 181c4762a1bSJed Brown call VecGetOwnershipRange(x1,low,high,ierr);CHKERRA(ierr) 182c4762a1bSJed Brown nloc = high - low 183c4762a1bSJed Brown 184c4762a1bSJed Brown call MatCreate(PETSC_COMM_WORLD,Amat,ierr);CHKERRA(ierr) 185c4762a1bSJed Brown call MatSetSizes(Amat,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr);CHKERRA(ierr) 186c4762a1bSJed Brown call MatSetUp(Amat,ierr);CHKERRA(ierr) 187c4762a1bSJed Brown 188c4762a1bSJed Brown call MatCreate(PETSC_COMM_WORLD,solver%AmatLin,ierr);CHKERRA(ierr) 189c4762a1bSJed Brown call MatSetSizes(solver%AmatLin,PETSC_DECIDE,PETSC_DECIDE,N1,N1,ierr);CHKERRA(ierr) 190c4762a1bSJed Brown call MatSetUp(solver%AmatLin,ierr);CHKERRA(ierr) 191c4762a1bSJed Brown 192c4762a1bSJed Brown call FormJacobianLocal(x1,solver%AmatLin,solver,.false.,ierr);CHKERRA(ierr) 193c4762a1bSJed Brown call MatAssemblyBegin(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr) 194c4762a1bSJed Brown call MatAssemblyEnd(solver%AmatLin,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr) 195c4762a1bSJed Brown 196c4762a1bSJed Brown call DMShellSetGlobalVector(daphi,x1,ierr);CHKERRA(ierr) 197c4762a1bSJed Brown call DMShellSetMatrix(daphi,Amat,ierr);CHKERRA(ierr) 198c4762a1bSJed Brown 199c4762a1bSJed Brown call VecCreate(PETSC_COMM_SELF,x1loc,ierr);CHKERRA(ierr) 200c4762a1bSJed Brown call VecSetSizes(x1loc,nloc,nloc,ierr);CHKERRA(ierr) 201c4762a1bSJed Brown call VecSetFromOptions(x1loc,ierr);CHKERRA(ierr) 202c4762a1bSJed Brown call DMShellSetLocalVector(daphi,x1loc,ierr);CHKERRA(ierr) 203c4762a1bSJed Brown 204c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 205c4762a1bSJed Brown! Create B, C, & D matrices 206c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 207c4762a1bSJed Brown call MatCreate(PETSC_COMM_WORLD,Cmat,ierr);CHKERRA(ierr) 208c4762a1bSJed Brown call MatSetSizes(Cmat,PETSC_DECIDE,PETSC_DECIDE,N2,N1,ierr);CHKERRA(ierr) 209c4762a1bSJed Brown call MatSetUp(Cmat,ierr);CHKERRA(ierr) 210c4762a1bSJed Brown! create data for C and B 211c4762a1bSJed Brown call MatCreate(PETSC_COMM_WORLD,Bmat,ierr);CHKERRA(ierr) 212c4762a1bSJed Brown call MatSetSizes(Bmat,PETSC_DECIDE,PETSC_DECIDE,N1,N2,ierr);CHKERRA(ierr) 213c4762a1bSJed Brown call MatSetUp(Bmat,ierr);CHKERRA(ierr) 214c4762a1bSJed Brown! create data for D 215c4762a1bSJed Brown call MatCreate(PETSC_COMM_WORLD,Dmat,ierr);CHKERRA(ierr) 216c4762a1bSJed Brown call MatSetSizes(Dmat,PETSC_DECIDE,PETSC_DECIDE,N2,N2,ierr);CHKERRA(ierr) 217c4762a1bSJed Brown call MatSetUp(Dmat,ierr);CHKERRA(ierr) 218c4762a1bSJed Brown 219c4762a1bSJed Brown call VecCreate(PETSC_COMM_WORLD,x2,ierr);CHKERRA(ierr) 220c4762a1bSJed Brown call VecSetSizes(x2,PETSC_DECIDE,N2,ierr);CHKERRA(ierr) 221c4762a1bSJed Brown call VecSetFromOptions(x2,ierr);CHKERRA(ierr) 222c4762a1bSJed Brown 223c4762a1bSJed Brown call VecGetOwnershipRange(x2,lamlow,lamhigh,ierr);CHKERRA(ierr) 224c4762a1bSJed Brown nloclam = lamhigh-lamlow 225c4762a1bSJed Brown 226c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 227c4762a1bSJed Brown! Set fake B and C 228c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 229c4762a1bSJed Brown one = 1.0 230c4762a1bSJed Brown if (N2 .gt. 0) then 231c4762a1bSJed Brown bval(1) = -one/(solver%mx-2) 232c4762a1bSJed Brown! cval = -one/(solver%my*solver%mx) 233c4762a1bSJed Brown cval(1) = -one 234c4762a1bSJed Brown do 20 irow=low,high-1 235c4762a1bSJed Brown j = irow/solver%mx ! row in domain 236c4762a1bSJed Brown i = mod(irow,solver%mx) 237c4762a1bSJed Brown row(1) = irow 238c4762a1bSJed Brown col(1) = j 239c4762a1bSJed Brown if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 240c4762a1bSJed Brown ! no op 241c4762a1bSJed Brown else 242c4762a1bSJed Brown call MatSetValues(Bmat,ione,row,ione,col,bval,INSERT_VALUES,ierr);CHKERRA(ierr) 243c4762a1bSJed Brown endif 244c4762a1bSJed Brown row(1) = j 245c4762a1bSJed Brown call MatSetValues(Cmat,ione,row,ione,row,cval,INSERT_VALUES,ierr);CHKERRA(ierr) 246c4762a1bSJed Brown 20 continue 247c4762a1bSJed Brown endif 248c4762a1bSJed Brown call MatAssemblyBegin(Bmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr) 249c4762a1bSJed Brown call MatAssemblyEnd(Bmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr) 250c4762a1bSJed Brown call MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr) 251c4762a1bSJed Brown call MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr) 252c4762a1bSJed Brown 253c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 254c4762a1bSJed Brown! Set D (indentity) 255c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 256c4762a1bSJed Brown do 30 j=lamlow,lamhigh-1 257c4762a1bSJed Brown row(1) = j 258c4762a1bSJed Brown cval(1) = one 259c4762a1bSJed Brown call MatSetValues(Dmat,ione,row,ione,row,cval,INSERT_VALUES,ierr);CHKERRA(ierr) 260c4762a1bSJed Brown 30 continue 261c4762a1bSJed Brown call MatAssemblyBegin(Dmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr) 262c4762a1bSJed Brown call MatAssemblyEnd(Dmat,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr) 263c4762a1bSJed Brown 264c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 265c4762a1bSJed Brown! DM for lambda (dalam) : temp driver for A block, setup A block solver data 266c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 267c4762a1bSJed Brown call DMShellCreate(PETSC_COMM_WORLD,dalam,ierr);CHKERRA(ierr) 268c4762a1bSJed Brown call DMShellSetGlobalVector(dalam,x2,ierr);CHKERRA(ierr) 269c4762a1bSJed Brown call DMShellSetMatrix(dalam,Dmat,ierr);CHKERRA(ierr) 270c4762a1bSJed Brown 271c4762a1bSJed Brown call VecCreate(PETSC_COMM_SELF,x2loc,ierr);CHKERRA(ierr) 272c4762a1bSJed Brown call VecSetSizes(x2loc,nloclam,nloclam,ierr);CHKERRA(ierr) 273c4762a1bSJed Brown call VecSetFromOptions(x2loc,ierr);CHKERRA(ierr) 274c4762a1bSJed Brown call DMShellSetLocalVector(dalam,x2loc,ierr);CHKERRA(ierr) 275c4762a1bSJed Brown 276c4762a1bSJed Brown call DMSetOptionsPrefix(dalam,'lambda_',ierr);CHKERRA(ierr) 277c4762a1bSJed Brown call DMSetFromOptions(dalam,ierr);CHKERRA(ierr) 278c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 279c4762a1bSJed Brown! Create field split DA 280c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 281c4762a1bSJed Brown call DMCompositeCreate(PETSC_COMM_WORLD,solver%da,ierr);CHKERRA(ierr) 282c4762a1bSJed Brown call DMCompositeAddDM(solver%da,daphi,ierr);CHKERRA(ierr) 283c4762a1bSJed Brown call DMCompositeAddDM(solver%da,dalam,ierr);CHKERRA(ierr) 284c4762a1bSJed Brown call DMSetFromOptions(solver%da,ierr);CHKERRA(ierr) 285c4762a1bSJed Brown call DMSetUp(solver%da,ierr);CHKERRA(ierr) 286c4762a1bSJed Brown call DMCompositeGetGlobalISs(solver%da,isglobal,ierr);CHKERRA(ierr) 287c4762a1bSJed Brown solver%isPhi = isglobal(1) 288c4762a1bSJed Brown solver%isLambda = isglobal(2) 289c4762a1bSJed Brown 290c4762a1bSJed Brown! cache matrices 291c4762a1bSJed Brown solver%Amat = Amat 292c4762a1bSJed Brown solver%Bmat = Bmat 293c4762a1bSJed Brown solver%Cmat = Cmat 294c4762a1bSJed Brown solver%Dmat = Dmat 295c4762a1bSJed Brown 296c4762a1bSJed Brown matArray(1) = Amat 297c4762a1bSJed Brown matArray(2) = Bmat 298c4762a1bSJed Brown matArray(3) = Cmat 299c4762a1bSJed Brown matArray(4) = Dmat 300c4762a1bSJed Brown 301c4762a1bSJed Brown call MatCreateNest(PETSC_COMM_WORLD,itwo,isglobal,itwo,isglobal,matArray,KKTmat,ierr);CHKERRA(ierr) 302c4762a1bSJed Brown call MatSetFromOptions(KKTmat,ierr);CHKERRA(ierr) 303c4762a1bSJed Brown 304c4762a1bSJed Brown! Extract global and local vectors from DMDA; then duplicate for remaining 305c4762a1bSJed Brown! vectors that are the same types 306c4762a1bSJed Brown call MatCreateVecs(KKTmat,x,PETSC_NULL_VEC,ierr);CHKERRA(ierr) 307c4762a1bSJed Brown call VecDuplicate(x,r,ierr);CHKERRA(ierr) 308c4762a1bSJed Brown 309c4762a1bSJed Brown call SNESCreate(PETSC_COMM_WORLD,mysnes,ierr);CHKERRA(ierr) 310c4762a1bSJed Brown 311c4762a1bSJed Brown call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr) 312c4762a1bSJed Brown 313c4762a1bSJed Brown call SNESSetApplicationContext(mysnes,solver,ierr);CHKERRA(ierr) 314c4762a1bSJed Brown 315c4762a1bSJed Brown call SNESSetDM(mysnes,solver%da,ierr);CHKERRA(ierr) 316c4762a1bSJed Brown 317c4762a1bSJed Brown! Set function evaluation routine and vector 318c4762a1bSJed Brown call SNESSetFunction(mysnes,r,FormFunction,solver,ierr);CHKERRA(ierr) 319c4762a1bSJed Brown 320c4762a1bSJed Brown call SNESSetJacobian(mysnes,KKTmat,KKTmat,FormJacobian,solver,ierr);CHKERRA(ierr) 321c4762a1bSJed Brown 322c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 323c4762a1bSJed Brown! Customize nonlinear solver; set runtime options 324c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 325c4762a1bSJed Brown! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 326c4762a1bSJed Brown call SNESSetFromOptions(mysnes,ierr);CHKERRA(ierr) 327c4762a1bSJed Brown 328c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 329c4762a1bSJed Brown! Evaluate initial guess; then solve nonlinear system. 330c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 331c4762a1bSJed Brown! Note: The user should initialize the vector, x, with the initial guess 332c4762a1bSJed Brown! for the nonlinear solver prior to calling SNESSolve(). In particular, 333c4762a1bSJed Brown! to employ an initial guess of zero, the user should explicitly set 334c4762a1bSJed Brown! this vector to zero by calling VecSet(). 335c4762a1bSJed Brown 336c4762a1bSJed Brown call FormInitialGuess(mysnes,x,ierr);CHKERRA(ierr) 337c4762a1bSJed Brown call SNESSolve(mysnes,PETSC_NULL_VEC,x,ierr);CHKERRA(ierr) 338c4762a1bSJed Brown call SNESGetIterationNumber(mysnes,its,ierr);CHKERRA(ierr) 339c4762a1bSJed Brown if (solver%rank .eq. 0) then 340c4762a1bSJed Brown write(6,100) its 341c4762a1bSJed Brown endif 342c4762a1bSJed Brown 100 format('Number of SNES iterations = ',i5) 343c4762a1bSJed Brown 344c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 345c4762a1bSJed Brown! Free work space. All PETSc objects should be destroyed when they 346c4762a1bSJed Brown! are no longer needed. 347c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 348c4762a1bSJed Brown call MatDestroy(KKTmat,ierr);CHKERRA(ierr) 349c4762a1bSJed Brown call MatDestroy(Amat,ierr);CHKERRA(ierr) 350c4762a1bSJed Brown call MatDestroy(Dmat,ierr);CHKERRA(ierr) 351c4762a1bSJed Brown call MatDestroy(Bmat,ierr);CHKERRA(ierr) 352c4762a1bSJed Brown call MatDestroy(Cmat,ierr);CHKERRA(ierr) 353c4762a1bSJed Brown call MatDestroy(solver%AmatLin,ierr);CHKERRA(ierr) 354c4762a1bSJed Brown call ISDestroy(solver%isPhi,ierr);CHKERRA(ierr) 355c4762a1bSJed Brown call ISDestroy(solver%isLambda,ierr);CHKERRA(ierr) 356c4762a1bSJed Brown call VecDestroy(x,ierr);CHKERRA(ierr) 357c4762a1bSJed Brown call VecDestroy(x2,ierr);CHKERRA(ierr) 358c4762a1bSJed Brown call VecDestroy(x1,ierr);CHKERRA(ierr) 359c4762a1bSJed Brown call VecDestroy(x1loc,ierr);CHKERRA(ierr) 360c4762a1bSJed Brown call VecDestroy(x2loc,ierr);CHKERRA(ierr) 361c4762a1bSJed Brown call VecDestroy(r,ierr);CHKERRA(ierr) 362c4762a1bSJed Brown call SNESDestroy(mysnes,ierr);CHKERRA(ierr) 363c4762a1bSJed Brown call DMDestroy(solver%da,ierr);CHKERRA(ierr) 364c4762a1bSJed Brown call DMDestroy(daphi,ierr);CHKERRA(ierr) 365c4762a1bSJed Brown call DMDestroy(dalam,ierr);CHKERRA(ierr) 366c4762a1bSJed Brown 367c4762a1bSJed Brown call PetscFinalize(ierr) 368c4762a1bSJed Brown end 369c4762a1bSJed Brown 370c4762a1bSJed Brown! --------------------------------------------------------------------- 371c4762a1bSJed Brown! 372c4762a1bSJed Brown! FormInitialGuess - Forms initial approximation. 373c4762a1bSJed Brown! 374c4762a1bSJed Brown! Input Parameters: 375c4762a1bSJed Brown! X - vector 376c4762a1bSJed Brown! 377c4762a1bSJed Brown! Output Parameter: 378c4762a1bSJed Brown! X - vector 379c4762a1bSJed Brown! 380c4762a1bSJed Brown! Notes: 381c4762a1bSJed Brown! This routine serves as a wrapper for the lower-level routine 382c4762a1bSJed Brown! "InitialGuessLocal", where the actual computations are 383c4762a1bSJed Brown! done using the standard Fortran style of treating the local 384c4762a1bSJed Brown! vector data as a multidimensional array over the local mesh. 385c4762a1bSJed Brown! This routine merely handles ghost point scatters and accesses 386c4762a1bSJed Brown! the local vector data via VecGetArrayF90() and VecRestoreArrayF90(). 387c4762a1bSJed Brown! 388c4762a1bSJed Brown subroutine FormInitialGuess(mysnes,Xnest,ierr) 389c4762a1bSJed Brown#include <petsc/finclude/petscsnes.h> 390c4762a1bSJed Brown use petscsnes 391c4762a1bSJed Brown use petsc_kkt_solver 392c4762a1bSJed Brown use petsc_kkt_solver_interfaces 393c4762a1bSJed Brown implicit none 394c4762a1bSJed Brown! Input/output variables: 395c4762a1bSJed Brown SNES:: mysnes 396c4762a1bSJed Brown Vec:: Xnest 397c4762a1bSJed Brown PetscErrorCode ierr 398c4762a1bSJed Brown 399c4762a1bSJed Brown! Declarations for use with local arrays: 400c4762a1bSJed Brown type(petsc_kkt_solver_type), pointer:: solver 401c4762a1bSJed Brown Vec:: Xsub(2) 402c4762a1bSJed Brown PetscInt:: izero,ione,itwo 403c4762a1bSJed Brown 404c4762a1bSJed Brown izero = 0 405c4762a1bSJed Brown ione = 1 406c4762a1bSJed Brown itwo = 2 407c4762a1bSJed Brown ierr = 0 408c4762a1bSJed Brown call SNESGetApplicationContext(mysnes,solver,ierr);CHKERRQ(ierr) 409c4762a1bSJed Brown call DMCompositeGetAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr) 410c4762a1bSJed Brown 411c4762a1bSJed Brown call InitialGuessLocal(solver,Xsub(1),ierr);CHKERRQ(ierr) 412c4762a1bSJed Brown call VecAssemblyBegin(Xsub(1),ierr);CHKERRQ(ierr) 413c4762a1bSJed Brown call VecAssemblyEnd(Xsub(1),ierr);CHKERRQ(ierr) 414c4762a1bSJed Brown 415c4762a1bSJed Brown! zero out lambda 416c4762a1bSJed Brown call VecZeroEntries(Xsub(2),ierr);CHKERRQ(ierr) 417c4762a1bSJed Brown call DMCompositeRestoreAccessArray(solver%da,Xnest,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr) 418c4762a1bSJed Brown 419c4762a1bSJed Brown return 420c4762a1bSJed Brown end subroutine FormInitialGuess 421c4762a1bSJed Brown 422c4762a1bSJed Brown! --------------------------------------------------------------------- 423c4762a1bSJed Brown! 424c4762a1bSJed Brown! InitialGuessLocal - Computes initial approximation, called by 425c4762a1bSJed Brown! the higher level routine FormInitialGuess(). 426c4762a1bSJed Brown! 427c4762a1bSJed Brown! Input Parameter: 428c4762a1bSJed Brown! X1 - local vector data 429c4762a1bSJed Brown! 430c4762a1bSJed Brown! Output Parameters: 431c4762a1bSJed Brown! x - local vector data 432c4762a1bSJed Brown! ierr - error code 433c4762a1bSJed Brown! 434c4762a1bSJed Brown! Notes: 435c4762a1bSJed Brown! This routine uses standard Fortran-style computations over a 2-dim array. 436c4762a1bSJed Brown! 437c4762a1bSJed Brown subroutine InitialGuessLocal(solver,X1,ierr) 438c4762a1bSJed Brown#include <petsc/finclude/petscsys.h> 439c4762a1bSJed Brown use petscsys 440c4762a1bSJed Brown use petsc_kkt_solver 441c4762a1bSJed Brown implicit none 442c4762a1bSJed Brown! Input/output variables: 443c4762a1bSJed Brown type (petsc_kkt_solver_type) solver 444c4762a1bSJed Brown Vec:: X1 445c4762a1bSJed Brown PetscErrorCode ierr 446c4762a1bSJed Brown 447c4762a1bSJed Brown! Local variables: 448c4762a1bSJed Brown PetscInt row,i,j,ione,low,high 449c4762a1bSJed Brown PetscReal temp1,temp,hx,hy,v 450c4762a1bSJed Brown PetscReal one 451c4762a1bSJed Brown 452c4762a1bSJed Brown! Set parameters 453c4762a1bSJed Brown ione = 1 454c4762a1bSJed Brown ierr = 0 455c4762a1bSJed Brown one = 1.0 456c4762a1bSJed Brown hx = one/(solver%mx-1) 457c4762a1bSJed Brown hy = one/(solver%my-1) 458c4762a1bSJed Brown temp1 = solver%lambda/(solver%lambda + one) + one 459c4762a1bSJed Brown 460c4762a1bSJed Brown call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr) 461c4762a1bSJed Brown 462c4762a1bSJed Brown do 20 row=low,high-1 463c4762a1bSJed Brown j = row/solver%mx 464c4762a1bSJed Brown i = mod(row,solver%mx) 465c4762a1bSJed Brown temp = min(j,solver%my-j+1)*hy 466c4762a1bSJed Brown if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 467c4762a1bSJed Brown v = 0.0 468c4762a1bSJed Brown else 469c4762a1bSJed Brown v = temp1 * sqrt(min(min(i,solver%mx-i+1)*hx,temp)) 470c4762a1bSJed Brown endif 471c4762a1bSJed Brown call VecSetValues(X1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr) 472c4762a1bSJed Brown 20 continue 473c4762a1bSJed Brown 474c4762a1bSJed Brown return 475c4762a1bSJed Brown end subroutine InitialGuessLocal 476c4762a1bSJed Brown 477c4762a1bSJed Brown! --------------------------------------------------------------------- 478c4762a1bSJed Brown! 479c4762a1bSJed Brown! FormJacobian - Evaluates Jacobian matrix. 480c4762a1bSJed Brown! 481c4762a1bSJed Brown! Input Parameters: 482c4762a1bSJed Brown! dummy - the SNES context 483c4762a1bSJed Brown! x - input vector 484c4762a1bSJed Brown! solver - solver data 485c4762a1bSJed Brown! 486c4762a1bSJed Brown! Output Parameters: 487c4762a1bSJed Brown! jac - Jacobian matrix 488c4762a1bSJed Brown! jac_prec - optionally different preconditioning matrix (not used here) 489c4762a1bSJed Brown! flag - flag indicating matrix structure 490c4762a1bSJed Brown! 491c4762a1bSJed Brown subroutine FormJacobian(dummy,X,jac,jac_prec,solver,ierr) 492c4762a1bSJed Brown#include <petsc/finclude/petscsnes.h> 493c4762a1bSJed Brown use petscsnes 494c4762a1bSJed Brown use petsc_kkt_solver 495c4762a1bSJed Brown implicit none 496c4762a1bSJed Brown! Input/output variables: 497c4762a1bSJed Brown SNES:: dummy 498c4762a1bSJed Brown Vec:: X 499c4762a1bSJed Brown Mat:: jac,jac_prec 500c4762a1bSJed Brown type(petsc_kkt_solver_type) solver 501c4762a1bSJed Brown PetscErrorCode ierr 502c4762a1bSJed Brown 503c4762a1bSJed Brown! Declarations for use with local arrays: 504c4762a1bSJed Brown Vec:: Xsub(1) 505c4762a1bSJed Brown Mat:: Amat 506*2a27bf02SStefano Zampini PetscInt ione 507c4762a1bSJed Brown 508c4762a1bSJed Brown ione = 1 509c4762a1bSJed Brown 510c4762a1bSJed Brown call DMCompositeGetAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr) 511c4762a1bSJed Brown 512c4762a1bSJed Brown! Compute entries for the locally owned part of the Jacobian preconditioner. 513c4762a1bSJed Brown call MatCreateSubMatrix(jac_prec,solver%isPhi,solver%isPhi,MAT_INITIAL_MATRIX,Amat,ierr);CHKERRQ(ierr) 514c4762a1bSJed Brown 515c4762a1bSJed Brown call FormJacobianLocal(Xsub(1),Amat,solver,.true.,ierr);CHKERRQ(ierr) 516c4762a1bSJed Brown call MatDestroy(Amat,ierr);CHKERRQ(ierr) ! discard our reference 517c4762a1bSJed Brown call DMCompositeRestoreAccessArray(solver%da,X,ione,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr) 518c4762a1bSJed Brown 519c4762a1bSJed Brown ! the rest of the matrix is not touched 520c4762a1bSJed Brown call MatAssemblyBegin(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr) 521c4762a1bSJed Brown call MatAssemblyEnd(jac_prec,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr) 522c4762a1bSJed Brown if (jac .ne. jac_prec) then 523c4762a1bSJed Brown call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr) 524c4762a1bSJed Brown call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr) 525c4762a1bSJed Brown end if 526c4762a1bSJed Brown 527c4762a1bSJed Brown! Tell the matrix we will never add a new nonzero location to the 528c4762a1bSJed Brown! matrix. If we do it will generate an error. 529c4762a1bSJed Brown call MatSetOption(jac_prec,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr);CHKERRQ(ierr) 530c4762a1bSJed Brown 531c4762a1bSJed Brown return 532c4762a1bSJed Brown end subroutine FormJacobian 533c4762a1bSJed Brown 534c4762a1bSJed Brown! --------------------------------------------------------------------- 535c4762a1bSJed Brown! 536c4762a1bSJed Brown! FormJacobianLocal - Computes Jacobian preconditioner matrix, 537c4762a1bSJed Brown! called by the higher level routine FormJacobian(). 538c4762a1bSJed Brown! 539c4762a1bSJed Brown! Input Parameters: 540c4762a1bSJed Brown! x - local vector data 541c4762a1bSJed Brown! 542c4762a1bSJed Brown! Output Parameters: 543c4762a1bSJed Brown! jac - Jacobian preconditioner matrix 544c4762a1bSJed Brown! ierr - error code 545c4762a1bSJed Brown! 546c4762a1bSJed Brown! Notes: 547c4762a1bSJed Brown! This routine uses standard Fortran-style computations over a 2-dim array. 548c4762a1bSJed Brown! 549c4762a1bSJed Brown subroutine FormJacobianLocal(X1,jac,solver,add_nl_term,ierr) 550c4762a1bSJed Brown#include <petsc/finclude/petscmat.h> 551c4762a1bSJed Brown use petscmat 552c4762a1bSJed Brown use petsc_kkt_solver 553c4762a1bSJed Brown implicit none 554c4762a1bSJed Brown! Input/output variables: 555c4762a1bSJed Brown type (petsc_kkt_solver_type) solver 556c4762a1bSJed Brown Vec:: X1 557c4762a1bSJed Brown Mat:: jac 558c4762a1bSJed Brown logical add_nl_term 559c4762a1bSJed Brown PetscErrorCode ierr 560c4762a1bSJed Brown 561c4762a1bSJed Brown! Local variables: 562c4762a1bSJed Brown PetscInt irow,row(1),col(5),i,j 563c4762a1bSJed Brown PetscInt ione,ifive,low,high,ii 564c4762a1bSJed Brown PetscScalar two,one,hx,hy,hy2inv 565c4762a1bSJed Brown PetscScalar hx2inv,sc,v(5) 566c4762a1bSJed Brown PetscScalar,pointer :: lx_v(:) 567c4762a1bSJed Brown 568c4762a1bSJed Brown! Set parameters 569c4762a1bSJed Brown ione = 1 570c4762a1bSJed Brown ifive = 5 571c4762a1bSJed Brown one = 1.0 572c4762a1bSJed Brown two = 2.0 573c4762a1bSJed Brown hx = one/(solver%mx-1) 574c4762a1bSJed Brown hy = one/(solver%my-1) 575c4762a1bSJed Brown sc = solver%lambda 576c4762a1bSJed Brown hx2inv = one/(hx*hx) 577c4762a1bSJed Brown hy2inv = one/(hy*hy) 578c4762a1bSJed Brown 579c4762a1bSJed Brown call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr) 580c4762a1bSJed Brown call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr) 581c4762a1bSJed Brown 582c4762a1bSJed Brown ii = 0 583c4762a1bSJed Brown do 20 irow=low,high-1 584c4762a1bSJed Brown j = irow/solver%mx 585c4762a1bSJed Brown i = mod(irow,solver%mx) 586c4762a1bSJed Brown ii = ii + 1 ! one based local index 587c4762a1bSJed Brown! boundary points 588c4762a1bSJed Brown if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 589c4762a1bSJed Brown col(1) = irow 590c4762a1bSJed Brown row(1) = irow 591c4762a1bSJed Brown v(1) = one 592c4762a1bSJed Brown call MatSetValues(jac,ione,row,ione,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr) 593c4762a1bSJed Brown! interior grid points 594c4762a1bSJed Brown else 595c4762a1bSJed Brown v(1) = -hy2inv 596c4762a1bSJed Brown if (j-1==0) v(1) = 0.0 597c4762a1bSJed Brown v(2) = -hx2inv 598c4762a1bSJed Brown if (i-1==0) v(2) = 0.0 599c4762a1bSJed Brown v(3) = two*(hx2inv + hy2inv) 600c4762a1bSJed Brown if (add_nl_term) v(3) = v(3) - sc*exp(lx_v(ii)) 601c4762a1bSJed Brown v(4) = -hx2inv 602c4762a1bSJed Brown if (i+1==solver%mx-1) v(4) = 0.0 603c4762a1bSJed Brown v(5) = -hy2inv 604c4762a1bSJed Brown if (j+1==solver%my-1) v(5) = 0.0 605c4762a1bSJed Brown col(1) = irow - solver%mx 606c4762a1bSJed Brown col(2) = irow - 1 607c4762a1bSJed Brown col(3) = irow 608c4762a1bSJed Brown col(4) = irow + 1 609c4762a1bSJed Brown col(5) = irow + solver%mx 610c4762a1bSJed Brown row(1) = irow 611c4762a1bSJed Brown call MatSetValues(jac,ione,row,ifive,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr) 612c4762a1bSJed Brown endif 613c4762a1bSJed Brown 20 continue 614c4762a1bSJed Brown 615c4762a1bSJed Brown call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr) 616c4762a1bSJed Brown 617c4762a1bSJed Brown return 618c4762a1bSJed Brown end subroutine FormJacobianLocal 619c4762a1bSJed Brown 620c4762a1bSJed Brown! --------------------------------------------------------------------- 621c4762a1bSJed Brown! 622c4762a1bSJed Brown! FormFunction - Evaluates nonlinear function, F(x). 623c4762a1bSJed Brown! 624c4762a1bSJed Brown! Input Parameters: 625c4762a1bSJed Brown! snes - the SNES context 626c4762a1bSJed Brown! X - input vector 627c4762a1bSJed Brown! dummy - optional user-defined context, as set by SNESSetFunction() 628c4762a1bSJed Brown! (not used here) 629c4762a1bSJed Brown! 630c4762a1bSJed Brown! Output Parameter: 631c4762a1bSJed Brown! F - function vector 632c4762a1bSJed Brown! 633c4762a1bSJed Brown subroutine FormFunction(snesIn,X,F,solver,ierr) 634c4762a1bSJed Brown#include <petsc/finclude/petscsnes.h> 635c4762a1bSJed Brown use petscsnes 636c4762a1bSJed Brown use petsc_kkt_solver 637c4762a1bSJed Brown implicit none 638c4762a1bSJed Brown! Input/output variables: 639c4762a1bSJed Brown SNES:: snesIn 640c4762a1bSJed Brown Vec:: X,F 641c4762a1bSJed Brown PetscErrorCode ierr 642c4762a1bSJed Brown type (petsc_kkt_solver_type) solver 643c4762a1bSJed Brown 644c4762a1bSJed Brown! Declarations for use with local arrays: 645c4762a1bSJed Brown Vec:: Xsub(2),Fsub(2) 646*2a27bf02SStefano Zampini PetscInt itwo 647c4762a1bSJed Brown 648c4762a1bSJed Brown! Scatter ghost points to local vector, using the 2-step process 649c4762a1bSJed Brown! DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 650c4762a1bSJed Brown! By placing code between these two statements, computations can 651c4762a1bSJed Brown! be done while messages are in transition. 652c4762a1bSJed Brown 653c4762a1bSJed Brown itwo = 2 654c4762a1bSJed Brown call DMCompositeGetAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr) 655c4762a1bSJed Brown call DMCompositeGetAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr) 656c4762a1bSJed Brown 657c4762a1bSJed Brown call FormFunctionNLTerm( Xsub(1), Fsub(1), solver, ierr);CHKERRQ(ierr) 658c4762a1bSJed Brown call MatMultAdd( solver%AmatLin, Xsub(1), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr) 659c4762a1bSJed Brown 660c4762a1bSJed Brown! do rest of operator (linear) 661c4762a1bSJed Brown call MatMult( solver%Cmat, Xsub(1), Fsub(2), ierr);CHKERRQ(ierr) 662c4762a1bSJed Brown call MatMultAdd( solver%Bmat, Xsub(2), Fsub(1), Fsub(1), ierr);CHKERRQ(ierr) 663c4762a1bSJed Brown call MatMultAdd( solver%Dmat, Xsub(2), Fsub(2), Fsub(2), ierr);CHKERRQ(ierr) 664c4762a1bSJed Brown 665c4762a1bSJed Brown call DMCompositeRestoreAccessArray(solver%da,X,itwo,PETSC_NULL_INTEGER,Xsub,ierr);CHKERRQ(ierr) 666c4762a1bSJed Brown call DMCompositeRestoreAccessArray(solver%da,F,itwo,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr) 667c4762a1bSJed Brown return 668c4762a1bSJed Brown end subroutine formfunction 669c4762a1bSJed Brown 670c4762a1bSJed Brown! --------------------------------------------------------------------- 671c4762a1bSJed Brown! 672c4762a1bSJed Brown! FormFunctionNLTerm - Computes nonlinear function, called by 673c4762a1bSJed Brown! the higher level routine FormFunction(). 674c4762a1bSJed Brown! 675c4762a1bSJed Brown! Input Parameter: 676c4762a1bSJed Brown! x - local vector data 677c4762a1bSJed Brown! 678c4762a1bSJed Brown! Output Parameters: 679c4762a1bSJed Brown! f - local vector data, f(x) 680c4762a1bSJed Brown! ierr - error code 681c4762a1bSJed Brown! 682c4762a1bSJed Brown! Notes: 683c4762a1bSJed Brown! This routine uses standard Fortran-style computations over a 2-dim array. 684c4762a1bSJed Brown! 685c4762a1bSJed Brown subroutine FormFunctionNLTerm(X1,F1,solver,ierr) 686c4762a1bSJed Brown#include <petsc/finclude/petscvec.h> 687c4762a1bSJed Brown use petscvec 688c4762a1bSJed Brown use petsc_kkt_solver 689c4762a1bSJed Brown implicit none 690c4762a1bSJed Brown! Input/output variables: 691c4762a1bSJed Brown type (petsc_kkt_solver_type) solver 692c4762a1bSJed Brown Vec:: X1,F1 693c4762a1bSJed Brown PetscErrorCode ierr 694c4762a1bSJed Brown! Local variables: 695*2a27bf02SStefano Zampini PetscScalar sc 696c4762a1bSJed Brown PetscScalar u,v(1) 697c4762a1bSJed Brown PetscInt i,j,low,high,ii,ione,irow,row(1) 698c4762a1bSJed Brown PetscScalar,pointer :: lx_v(:) 699c4762a1bSJed Brown 700c4762a1bSJed Brown sc = solver%lambda 701c4762a1bSJed Brown ione = 1 702c4762a1bSJed Brown 703c4762a1bSJed Brown call VecGetArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr) 704c4762a1bSJed Brown call VecGetOwnershipRange(X1,low,high,ierr);CHKERRQ(ierr) 705c4762a1bSJed Brown 706c4762a1bSJed Brown! Compute function over the locally owned part of the grid 707c4762a1bSJed Brown ii = 0 708c4762a1bSJed Brown do 20 irow=low,high-1 709c4762a1bSJed Brown j = irow/solver%mx 710c4762a1bSJed Brown i = mod(irow,solver%mx) 711c4762a1bSJed Brown ii = ii + 1 ! one based local index 712c4762a1bSJed Brown row(1) = irow 713c4762a1bSJed Brown if (i .eq. 0 .or. j .eq. 0 .or. i .eq. solver%mx-1 .or. j .eq. solver%my-1) then 714c4762a1bSJed Brown v(1) = 0.0 715c4762a1bSJed Brown else 716c4762a1bSJed Brown u = lx_v(ii) 717c4762a1bSJed Brown v(1) = -sc*exp(u) 718c4762a1bSJed Brown endif 719c4762a1bSJed Brown call VecSetValues(F1,ione,row,v,INSERT_VALUES,ierr);CHKERRQ(ierr) 720c4762a1bSJed Brown 20 continue 721c4762a1bSJed Brown 722c4762a1bSJed Brown call VecRestoreArrayReadF90(X1,lx_v,ierr);CHKERRQ(ierr) 723c4762a1bSJed Brown 724c4762a1bSJed Brown call VecAssemblyBegin(F1,ierr);CHKERRQ(ierr) 725c4762a1bSJed Brown call VecAssemblyEnd(F1,ierr);CHKERRQ(ierr) 726c4762a1bSJed Brown 727c4762a1bSJed Brown ierr = 0 728c4762a1bSJed Brown return 729c4762a1bSJed Brown end subroutine FormFunctionNLTerm 730c4762a1bSJed Brown 731c4762a1bSJed Brown!/*TEST 732c4762a1bSJed Brown! 733c4762a1bSJed Brown! build: 7349e489221SSatish Balay! requires: !single !complex 735c4762a1bSJed Brown! 736c4762a1bSJed Brown! test: 737c4762a1bSJed Brown! nsize: 4 738c4762a1bSJed Brown! args: -par 5.0 -da_grid_x 10 -da_grid_y 10 -snes_monitor_short -snes_linesearch_type basic -snes_converged_reason -ksp_type fgmres -ksp_norm_type unpreconditioned -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type upper -ksp_monitor_short -fieldsplit_lambda_ksp_type preonly -fieldsplit_lambda_pc_type jacobi -fieldsplit_phi_pc_type gamg -fieldsplit_phi_pc_gamg_agg_nsmooths 1 -fieldsplit_phi_pc_gamg_threshold 0. -fieldsplit_phi_pc_gamg_esteig_ksp_type cg 739c4762a1bSJed Brown! 740c4762a1bSJed Brown!TEST*/ 741