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