xref: /petsc/src/binding/petsc4py/demo/legacy/wrap-f2py/Bratu2D.F90 (revision edb0e59d3c097acd4a4005a4e51d4daa5c739255) 
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    PetscCall(DMDAGetInfo(da, PETSC_NULL_INTEGER, grd%mx, grd%my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
41    PetscCall(DMDAGetCorners(da, grd%xs, grd%ys, PETSC_NULL_INTEGER, grd%xm, grd%ym, PETSC_NULL_INTEGER, ierr))
42    PetscCall(DMDAGetGhostCorners(da, grd%gxs, grd%gys, PETSC_NULL_INTEGER, grd%gxm, grd%gym, PETSC_NULL_INTEGER, ierr))
43
44    grd%xs = grd%xs + 1
45    grd%ys = grd%ys + 1
46    grd%gxs = grd%gxs + 1
47    grd%gys = grd%gys + 1
48
49    grd%ye = grd%ys + grd%ym - 1
50    grd%xe = grd%xs + grd%xm - 1
51    grd%gye = grd%gys + grd%gym - 1
52    grd%gxe = grd%gxs + grd%gxm - 1
53
54  end subroutine GetGridInfo
55
56  subroutine InitGuessLocal(grd, x, lambda, ierr)
57    implicit none
58    type(gridinfo) grd
59    PetscScalar x(grd%xs:grd%xe, grd%ys:grd%ye)
60    PetscReal lambda
61    PetscErrorCode ierr
62    !
63    PetscInt i, j
64    PetscReal hx, hy, temp, temp1, one
65
66    one = 1.0
67    hx = one/(dble(grd%mx - 1))
68    hy = one/(dble(grd%my - 1))
69    temp1 = lambda/(lambda + one)
70
71    do j = grd%ys, grd%ye
72      temp = dble(min(j - 1, grd%my - j))*hy
73      do i = grd%xs, grd%xe
74        if (i == 1 .or. j == 1 .or. i == grd%mx .or. j == grd%my) then
75          ! boundary points
76          x(i, j) = 0.0
77        else
78          ! interior grid points
79          x(i, j) = temp1*sqrt(min(dble(min(i - 1, grd%mx - i)*hx), dble(temp)))
80        end if
81      end do
82    end do
83    ierr = 0
84
85  end subroutine InitGuessLocal
86
87  subroutine FunctionLocal(grd, x, f, lambda, ierr)
88    implicit none
89    type(gridinfo) grd
90    PetscScalar x(grd%gxs:grd%gxe, grd%gys:grd%gye)
91    PetscScalar f(grd%xs:grd%xe, grd%ys:grd%ye)
92    PetscReal lambda
93    PetscErrorCode ierr
94    !
95    PetscInt i, j
96    PetscReal hx, hy, hxdhy, hydhx, sc, one, two
97    PetscScalar u, uxx, uyy
98
99    one = 1.0
100    two = 2.0
101    hx = one/dble(grd%mx - 1)
102    hy = one/dble(grd%my - 1)
103    sc = hx*hy
104    hxdhy = hx/hy
105    hydhx = hy/hx
106
107    do j = grd%ys, grd%ye
108      do i = grd%xs, grd%xe
109        if (i == 1 .or. j == 1 .or. i == grd%mx .or. j == grd%my) then
110          ! boundary points
111          f(i, j) = x(i, j) - 0.0
112        else
113          ! interior grid points
114          u = x(i, j)
115          uxx = (two*u - x(i - 1, j) - x(i + 1, j))*hydhx
116          uyy = (two*u - x(i, j - 1) - x(i, j + 1))*hxdhy
117          f(i, j) = uxx + uyy - lambda*exp(u)*sc
118        end if
119      end do
120    end do
121    ierr = 0
122
123  end subroutine FunctionLocal
124
125  subroutine JacobianLocal(grd, x, Jac, lambda, ierr)
126    implicit none
127    type(gridinfo) grd
128    PetscScalar x(grd%gxs:grd%gxe, grd%gys:grd%gye)
129    Mat Jac
130    PetscReal lambda
131    PetscErrorCode ierr
132    !
133    PetscInt i, j, row(1), col(5)
134    PetscInt ione, ifive
135    PetscReal hx, hy, hxdhy, hydhx, sc, v(5), one, two
136
137    ione = 1
138    ifive = 5
139    one = 1.0
140    two = 2.0
141    hx = one/dble(grd%mx - 1)
142    hy = one/dble(grd%my - 1)
143    sc = hx*hy
144    hxdhy = hx/hy
145    hydhx = hy/hx
146
147    do j = grd%ys, grd%ye
148      row = (j - grd%gys)*grd%gxm + grd%xs - grd%gxs - 1
149      do i = grd%xs, grd%xe
150        row = row + 1
151        if (i == 1 .or. j == 1 .or. i == grd%mx .or. j == grd%my) then
152          ! boundary points
153          col(1) = row(1)
154          v(1) = one
155          PetscCall(MatSetValuesLocal(Jac, ione, row, ione, col, v, INSERT_VALUES, ierr))
156        else
157          ! interior grid points
158          v(1) = -hxdhy
159          v(2) = -hydhx
160          v(3) = two*(hydhx + hxdhy) - lambda*exp(x(i, j))*sc
161          v(4) = -hydhx
162          v(5) = -hxdhy
163          col(1) = row(1) - grd%gxm
164          col(2) = row(1) - 1
165          col(3) = row(1)
166          col(4) = row(1) + 1
167          col(5) = row(1) + grd%gxm
168          PetscCall(MatSetValuesLocal(Jac, ione, row, ifive, col, v, INSERT_VALUES, ierr))
169        end if
170      end do
171    end do
172
173  end subroutine JacobianLocal
174
175end module Bratu2D
176
177! --------------------------------------------------------------------
178
179subroutine FormInitGuess(da, X, lambda, ierr)
180  use Bratu2D
181  implicit none
182  DM da
183  Vec X
184  PetscReal lambda
185  PetscErrorCode ierr
186  !
187  type(gridinfo)      :: grd
188  PetscScalar, pointer :: xx(:)
189
190  PetscCall(VecGetArray(X, xx, ierr))
191  PetscCall(GetGridInfo(da, grd, ierr))
192  PetscCall(InitGuessLocal(grd, xx, lambda, ierr))
193  PetscCall(VecRestoreArray(X, xx, ierr))
194
195end subroutine FormInitGuess
196
197subroutine FormFunction(da, X, F, lambda, ierr)
198  use Bratu2D
199  implicit none
200  DM da
201  Vec X
202  Vec F
203  PetscReal lambda
204  PetscErrorCode ierr
205  !
206  type(gridinfo)      :: grd
207  Vec                 :: localX
208  PetscScalar, pointer :: xx(:)
209  PetscScalar, pointer :: ff(:)
210
211  PetscCall(DMGetLocalVector(da, localX, ierr))
212  PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr))
213  PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr))
214
215  PetscCall(VecGetArray(localX, xx, ierr))
216  PetscCall(VecGetArray(F, ff, ierr))
217
218  PetscCall(GetGridInfo(da, grd, ierr))
219  PetscCall(FunctionLocal(grd, xx, ff, lambda, ierr))
220
221  PetscCall(VecRestoreArray(F, ff, ierr))
222  PetscCall(VecRestoreArray(localX, xx, ierr))
223  PetscCall(DMRestoreLocalVector(da, localX, ierr))
224
225end subroutine FormFunction
226
227subroutine FormJacobian(da, X, J, lambda, ierr)
228  use Bratu2D
229  implicit none
230  DM da
231  Vec X
232  Mat J
233  PetscReal lambda
234  PetscErrorCode ierr
235  !
236  type(gridinfo)      :: grd
237  Vec                 :: localX
238  PetscScalar, pointer :: xx(:)
239
240  PetscCall(DMGetLocalVector(da, localX, ierr))
241  PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr))
242  PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr))
243  PetscCall(VecGetArray(localX, xx, ierr))
244
245  PetscCall(GetGridInfo(da, grd, ierr))
246  PetscCall(JacobianLocal(grd, xx, J, lambda, ierr))
247
248  PetscCall(VecRestoreArray(localX, xx, ierr))
249  PetscCall(DMRestoreLocalVector(da, localX, ierr))
250
251  PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY, ierr))
252  PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY, ierr))
253
254end subroutine FormJacobian
255
256! --------------------------------------------------------------------
257
258! Local Variables:
259! mode: f90
260! End:
261