1! --------------------------------------------------------------------- 2! 3! Solid Fuel Ignition (SFI) problem. This problem is modeled by 4! the partial differential equation 5! 6! -Laplacian(u) - lambda * exp(u) = 0, 0 < x,y < 1, 7! 8! with boundary conditions 9! 10! u = 0 for x = 0, x = 1, y = 0, y = 1, 11! 12! A finite difference approximation with the usual 5-point stencil 13! is used to discretize the boundary value problem to obtain a 14! nonlinear system of equations. The problem is solved in a 2D 15! rectangular domain, using distributed arrays (DAs) to partition 16! the parallel grid. 17! 18! -------------------------------------------------------------------- 19 20#include <petsc/finclude/petscdm.h> 21 22module Bratu2D 23 24 use petsc 25 implicit none 26 27 type gridinfo 28 PetscInt mx,xs,xe,xm,gxs,gxe,gxm 29 PetscInt my,ys,ye,ym,gys,gye,gym 30 end type gridinfo 31 32contains 33 34 subroutine GetGridInfo(da, grd, ierr) 35 implicit none 36 DM da 37 type(gridinfo) grd 38 PetscErrorCode ierr 39 ! 40 call DMDAGetInfo(da, PETSC_NULL_INTEGER, & 41 & grd%mx, grd%my, PETSC_NULL_INTEGER, & 42 & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 43 & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 44 & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 45 & PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr) 46 call DMDAGetCorners(da, & 47 & grd%xs,grd%ys,PETSC_NULL_INTEGER, & 48 & grd%xm,grd%ym,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr) 49 call DMDAGetGhostCorners(da, & 50 & grd%gxs,grd%gys,PETSC_NULL_INTEGER, & 51 & grd%gxm,grd%gym,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr) 52 53 grd%xs = grd%xs+1 54 grd%ys = grd%ys+1 55 grd%gxs = grd%gxs+1 56 grd%gys = grd%gys+1 57 58 grd%ye = grd%ys+grd%ym-1 59 grd%xe = grd%xs+grd%xm-1 60 grd%gye = grd%gys+grd%gym-1 61 grd%gxe = grd%gxs+grd%gxm-1 62 63 end subroutine GetGridInfo 64 65 subroutine InitGuessLocal(grd, x, lambda, ierr) 66 implicit none 67 type(gridinfo) grd 68 PetscScalar x(grd%xs:grd%xe,grd%ys:grd%ye) 69 PetscReal lambda 70 PetscErrorCode ierr 71 ! 72 PetscInt i, j 73 PetscReal hx,hy,temp,temp1,one 74 75 one = 1.0 76 hx = one/(dble(grd%mx-1)) 77 hy = one/(dble(grd%my-1)) 78 temp1 = lambda/(lambda+one) 79 80 do j=grd%ys,grd%ye 81 temp = dble(min(j-1,grd%my-j))*hy 82 do i=grd%xs,grd%xe 83 if (i==1 .or. j==1 .or. i==grd%mx .or. j==grd%my) then 84 ! boundary points 85 x(i,j) = 0.0 86 else 87 ! interior grid points 88 x(i,j) = temp1*sqrt(min(dble(min(i-1,grd%mx-i)*hx),dble(temp))) 89 end if 90 end do 91 end do 92 ierr = 0 93 94 end subroutine InitGuessLocal 95 96 subroutine FunctionLocal(grd, x, f, lambda, ierr) 97 implicit none 98 type(gridinfo) grd 99 PetscScalar x(grd%gxs:grd%gxe,grd%gys:grd%gye) 100 PetscScalar f(grd%xs:grd%xe,grd%ys:grd%ye) 101 PetscReal lambda 102 PetscErrorCode ierr 103 ! 104 PetscInt i,j 105 PetscReal hx,hy,hxdhy,hydhx,sc,one,two 106 PetscScalar u,uxx,uyy 107 108 one = 1.0 109 two = 2.0 110 hx = one/dble(grd%mx-1) 111 hy = one/dble(grd%my-1) 112 sc = hx*hy 113 hxdhy = hx/hy 114 hydhx = hy/hx 115 116 do j=grd%ys,grd%ye 117 do i=grd%xs,grd%xe 118 if (i==1 .or. j==1 .or. i==grd%mx .or. j==grd%my) then 119 ! boundary points 120 f(i,j) = x(i,j) - 0.0 121 else 122 ! interior grid points 123 u = x(i,j) 124 uxx = (two*u - x(i-1,j) - x(i+1,j)) * hydhx 125 uyy = (two*u - x(i,j-1) - x(i,j+1)) * hxdhy 126 f(i,j) = uxx + uyy - lambda*exp(u)*sc 127 end if 128 end do 129 end do 130 ierr = 0 131 132 end subroutine FunctionLocal 133 134 subroutine JacobianLocal(grd, x, Jac, lambda, ierr) 135 implicit none 136 type(gridinfo) grd 137 PetscScalar x(grd%gxs:grd%gxe,grd%gys:grd%gye) 138 Mat Jac 139 PetscReal lambda 140 PetscErrorCode ierr 141 ! 142 PetscInt i,j,row(1),col(5) 143 PetscInt ione,ifive 144 PetscReal hx,hy,hxdhy,hydhx,sc,v(5),one,two 145 146 ione = 1 147 ifive = 5 148 one = 1.0 149 two = 2.0 150 hx = one/dble(grd%mx-1) 151 hy = one/dble(grd%my-1) 152 sc = hx*hy 153 hxdhy = hx/hy 154 hydhx = hy/hx 155 156 do j=grd%ys,grd%ye 157 row = (j - grd%gys)*grd%gxm + grd%xs - grd%gxs - 1 158 do i=grd%xs,grd%xe 159 row = row + 1 160 if (i==1 .or. j==1 .or. i==grd%mx .or. j==grd%my) then 161 ! boundary points 162 col(1) = row(1) 163 v(1) = one 164 call MatSetValuesLocal(Jac,ione,row,ione,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr) 165 else 166 ! interior grid points 167 v(1) = -hxdhy 168 v(2) = -hydhx 169 v(3) = two*(hydhx + hxdhy) - lambda*exp(x(i,j))*sc 170 v(4) = -hydhx 171 v(5) = -hxdhy 172 col(1) = row(1) - grd%gxm 173 col(2) = row(1) - 1 174 col(3) = row(1) 175 col(4) = row(1) + 1 176 col(5) = row(1) + grd%gxm 177 call MatSetValuesLocal(Jac,ione,row,ifive,col,v,INSERT_VALUES,ierr);CHKERRQ(ierr) 178 end if 179 end do 180 end do 181 182 end subroutine JacobianLocal 183 184end module Bratu2D 185 186! -------------------------------------------------------------------- 187 188subroutine FormInitGuess(da, X, lambda, ierr) 189 use Bratu2D 190 implicit none 191 DM da 192 Vec X 193 PetscReal lambda 194 PetscErrorCode ierr 195 ! 196 type(gridinfo) :: grd 197 PetscScalar,pointer :: xx(:) 198 199 call VecGetArrayF90(X,xx,ierr);CHKERRQ(ierr) 200 201 call GetGridInfo(da,grd,ierr);CHKERRQ(ierr) 202 call InitGuessLocal(grd,xx,lambda,ierr);CHKERRQ(ierr) 203 204 call VecRestoreArrayF90(X,xx,ierr);CHKERRQ(ierr) 205 206end subroutine FormInitGuess 207 208 209subroutine FormFunction(da, X, F, lambda, ierr) 210 use Bratu2D 211 implicit none 212 DM da 213 Vec X 214 Vec F 215 PetscReal lambda 216 PetscErrorCode ierr 217 ! 218 type(gridinfo) :: grd 219 Vec :: localX 220 PetscScalar,pointer :: xx(:) 221 PetscScalar,pointer :: ff(:) 222 223 call DMGetLocalVector(da,localX,ierr);CHKERRQ(ierr) 224 call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr) 225 call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr) 226 227 call VecGetArrayF90(localX,xx,ierr);CHKERRQ(ierr) 228 call VecGetArrayF90(F,ff,ierr);CHKERRQ(ierr) 229 230 call GetGridInfo(da,grd,ierr);CHKERRQ(ierr) 231 call FunctionLocal(grd,xx,ff,lambda,ierr);CHKERRQ(ierr) 232 233 call VecRestoreArrayF90(F,ff,ierr);CHKERRQ(ierr) 234 call VecRestoreArrayF90(localX,xx,ierr);CHKERRQ(ierr) 235 236 call DMRestoreLocalVector(da,localX,ierr);CHKERRQ(ierr) 237 238end subroutine FormFunction 239 240 241subroutine FormJacobian(da, X, J, lambda, ierr) 242 use Bratu2D 243 implicit none 244 DM da 245 Vec X 246 Mat J 247 PetscReal lambda 248 PetscErrorCode ierr 249 ! 250 type(gridinfo) :: grd 251 Vec :: localX 252 PetscScalar,pointer :: xx(:) 253 254 call DMGetLocalVector(da,localX,ierr);CHKERRQ(ierr) 255 call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr) 256 call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr) 257 call VecGetArrayF90(localX,xx,ierr);CHKERRQ(ierr) 258 259 call GetGridInfo(da,grd,ierr);CHKERRQ(ierr) 260 call JacobianLocal(grd,xx,J,lambda,ierr);CHKERRQ(ierr) 261 262 call VecRestoreArrayF90(localX,xx,ierr);CHKERRQ(ierr) 263 call DMRestoreLocalVector(da,localX,ierr);CHKERRQ(ierr) 264 265 call MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr) 266 call MatAssemblyEnd (J,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr) 267 268end subroutine FormJacobian 269 270! -------------------------------------------------------------------- 271 272 273! Local Variables: 274! mode: f90 275! End: 276