xref: /petsc/src/snes/tutorials/ex5f90.F90 (revision 4e8208cbcbc709572b8abe32f33c78b69c819375)
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
10c4762a1bSJed Brown!
11c4762a1bSJed Brown!  --------------------------------------------------------------------------
12c4762a1bSJed Brown!
13c4762a1bSJed Brown!  Solid Fuel Ignition (SFI) problem.  This problem is modeled by
14c4762a1bSJed Brown!  the partial differential equation
15c4762a1bSJed Brown!
16c4762a1bSJed Brown!          -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
17c4762a1bSJed Brown!
18c4762a1bSJed Brown!  with boundary conditions
19c4762a1bSJed Brown!
20c4762a1bSJed Brown!           u = 0  for  x = 0, x = 1, y = 0, y = 1.
21c4762a1bSJed Brown!
22c4762a1bSJed Brown!  A finite difference approximation with the usual 5-point stencil
23c4762a1bSJed Brown!  is used to discretize the boundary value problem to obtain a nonlinear
24c4762a1bSJed Brown!  system of equations.
25c4762a1bSJed Brown!
26c4762a1bSJed Brown!  The uniprocessor version of this code is snes/tutorials/ex4f.F
27c4762a1bSJed Brown!
28c4762a1bSJed Brown!  --------------------------------------------------------------------------
29c4762a1bSJed Brown!  The following define must be used before including any PETSc include files
30c4762a1bSJed Brown!  into a module or interface. This is because they can't handle declarations
31c4762a1bSJed Brown!  in them
32c4762a1bSJed Brown!
33ce78bad3SBarry Smith#include <petsc/finclude/petscsnes.h>
34ce78bad3SBarry Smith#include <petsc/finclude/petscdmda.h>
35b06eb4cdSBarry Smithmodule ex5module
36c4762a1bSJed Brown  use petscsnes
37dfbbaf82SBarry Smith  use petscdmda
38e7a95102SMartin Diehl  implicit none
39*2a8381b2SBarry Smith  type AppCtx
40c4762a1bSJed Brown    PetscInt xs, xe, xm, gxs, gxe, gxm
41c4762a1bSJed Brown    PetscInt ys, ye, ym, gys, gye, gym
42c4762a1bSJed Brown    PetscInt mx, my
43c4762a1bSJed Brown    PetscMPIInt rank
44c4762a1bSJed Brown    PetscReal lambda
45*2a8381b2SBarry Smith  end type AppCtx
46e7a95102SMartin Diehl
47c4762a1bSJed Browncontains
48c4762a1bSJed Brown! ---------------------------------------------------------------------
49c4762a1bSJed Brown!
50c4762a1bSJed Brown!  FormFunction - Evaluates nonlinear function, F(x).
51c4762a1bSJed Brown!
52c4762a1bSJed Brown!  Input Parameters:
53c4762a1bSJed Brown!  snes - the SNES context
54c4762a1bSJed Brown!  X - input vector
55c4762a1bSJed Brown!  dummy - optional user-defined context, as set by SNESSetFunction()
56c4762a1bSJed Brown!          (not used here)
57c4762a1bSJed Brown!
58c4762a1bSJed Brown!  Output Parameter:
59c4762a1bSJed Brown!  F - function vector
60c4762a1bSJed Brown!
61c4762a1bSJed Brown!  Notes:
62c4762a1bSJed Brown!  This routine serves as a wrapper for the lower-level routine
63c4762a1bSJed Brown!  "FormFunctionLocal", where the actual computations are
64c4762a1bSJed Brown!  done using the standard Fortran style of treating the local
65c4762a1bSJed Brown!  vector data as a multidimensional array over the local mesh.
66c4762a1bSJed Brown!  This routine merely handles ghost point scatters and accesses
67ce78bad3SBarry Smith!  the local vector data via VecGetArray() and VecRestoreArray().
68c4762a1bSJed Brown!
69*2a8381b2SBarry Smith  subroutine FormFunction(snes, X, F, ctx, ierr)
70c4762a1bSJed Brown    implicit none
71c4762a1bSJed Brown
72c4762a1bSJed Brown!  Input/output variables:
73c4762a1bSJed Brown    SNES snes
74c4762a1bSJed Brown    Vec X, F
75c4762a1bSJed Brown    PetscErrorCode ierr
76*2a8381b2SBarry Smith    type(AppCtx) ctx
77c4762a1bSJed Brown    DM da
78c4762a1bSJed Brown
79c4762a1bSJed Brown!  Declarations for use with local arrays:
80c4762a1bSJed Brown    PetscScalar, pointer :: lx_v(:), lf_v(:)
81c4762a1bSJed Brown    Vec localX
82c4762a1bSJed Brown
83c4762a1bSJed Brown!  Scatter ghost points to local vector, using the 2-step process
84c4762a1bSJed Brown!     DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
85c4762a1bSJed Brown!  By placing code between these two statements, computations can
86c4762a1bSJed Brown!  be done while messages are in transition.
87d8606c27SBarry Smith    PetscCall(SNESGetDM(snes, da, ierr))
88d8606c27SBarry Smith    PetscCall(DMGetLocalVector(da, localX, ierr))
89d8606c27SBarry Smith    PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr))
90d8606c27SBarry Smith    PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr))
91c4762a1bSJed Brown
92c4762a1bSJed Brown!  Get a pointer to vector data.
93ce78bad3SBarry Smith!    - For default PETSc vectors, VecGetArray() returns a pointer to
94c4762a1bSJed Brown!      the data array. Otherwise, the routine is implementation dependent.
95ce78bad3SBarry Smith!    - You MUST call VecRestoreArray() when you no longer need access to
96c4762a1bSJed Brown!      the array.
97ce78bad3SBarry Smith!    - Note that the interface to VecGetArray() differs from VecGetArray().
98c4762a1bSJed Brown
99ce78bad3SBarry Smith    PetscCall(VecGetArray(localX, lx_v, ierr))
100ce78bad3SBarry Smith    PetscCall(VecGetArray(F, lf_v, ierr))
101c4762a1bSJed Brown
102c4762a1bSJed Brown!  Compute function over the locally owned part of the grid
103*2a8381b2SBarry Smith    PetscCall(FormFunctionLocal(lx_v, lf_v, ctx, ierr))
104c4762a1bSJed Brown
105c4762a1bSJed Brown!  Restore vectors
106ce78bad3SBarry Smith    PetscCall(VecRestoreArray(localX, lx_v, ierr))
107ce78bad3SBarry Smith    PetscCall(VecRestoreArray(F, lf_v, ierr))
108c4762a1bSJed Brown
109c4762a1bSJed Brown!  Insert values into global vector
110c4762a1bSJed Brown
111d8606c27SBarry Smith    PetscCall(DMRestoreLocalVector(da, localX, ierr))
112*2a8381b2SBarry Smith    PetscCall(PetscLogFlops(11.0d0*ctx%ym*ctx%xm, ierr))
113c4762a1bSJed Brown
114d8606c27SBarry Smith!      PetscCallA(VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr))
115d8606c27SBarry Smith!      PetscCallA(VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr))
116c4762a1bSJed Brown  end subroutine formfunction
117e7a95102SMartin Diehl
118e7a95102SMartin Diehl! ---------------------------------------------------------------------
119e7a95102SMartin Diehl!
120e7a95102SMartin Diehl!  FormInitialGuess - Forms initial approximation.
121e7a95102SMartin Diehl!
122e7a95102SMartin Diehl!  Input Parameters:
123e7a95102SMartin Diehl!  X - vector
124e7a95102SMartin Diehl!
125e7a95102SMartin Diehl!  Output Parameter:
126e7a95102SMartin Diehl!  X - vector
127e7a95102SMartin Diehl!
128e7a95102SMartin Diehl!  Notes:
129e7a95102SMartin Diehl!  This routine serves as a wrapper for the lower-level routine
130e7a95102SMartin Diehl!  "InitialGuessLocal", where the actual computations are
131e7a95102SMartin Diehl!  done using the standard Fortran style of treating the local
132e7a95102SMartin Diehl!  vector data as a multidimensional array over the local mesh.
133e7a95102SMartin Diehl!  This routine merely handles ghost point scatters and accesses
134e7a95102SMartin Diehl!  the local vector data via VecGetArray() and VecRestoreArray().
135e7a95102SMartin Diehl!
136e7a95102SMartin Diehl  subroutine FormInitialGuess(snes, X, ierr)
137e7a95102SMartin Diehl!  Input/output variables:
138e7a95102SMartin Diehl    SNES snes
139*2a8381b2SBarry Smith    type(AppCtx), pointer:: ctx
140e7a95102SMartin Diehl    Vec X
141e7a95102SMartin Diehl    PetscErrorCode ierr
142e7a95102SMartin Diehl    DM da
143e7a95102SMartin Diehl
144e7a95102SMartin Diehl!  Declarations for use with local arrays:
145e7a95102SMartin Diehl    PetscScalar, pointer :: lx_v(:)
146e7a95102SMartin Diehl
147e7a95102SMartin Diehl    ierr = 0
148e7a95102SMartin Diehl    PetscCallA(SNESGetDM(snes, da, ierr))
149*2a8381b2SBarry Smith    PetscCallA(SNESGetApplicationContext(snes, ctx, ierr))
150e7a95102SMartin Diehl!  Get a pointer to vector data.
151e7a95102SMartin Diehl!    - For default PETSc vectors, VecGetArray() returns a pointer to
152e7a95102SMartin Diehl!      the data array. Otherwise, the routine is implementation dependent.
153e7a95102SMartin Diehl!    - You MUST call VecRestoreArray() when you no longer need access to
154e7a95102SMartin Diehl!      the array.
155e7a95102SMartin Diehl!    - Note that the interface to VecGetArray() differs from VecGetArray().
156e7a95102SMartin Diehl
157e7a95102SMartin Diehl    PetscCallA(VecGetArray(X, lx_v, ierr))
158e7a95102SMartin Diehl
159e7a95102SMartin Diehl!  Compute initial guess over the locally owned part of the grid
160*2a8381b2SBarry Smith    PetscCallA(InitialGuessLocal(ctx, lx_v, ierr))
161e7a95102SMartin Diehl
162e7a95102SMartin Diehl!  Restore vector
163e7a95102SMartin Diehl    PetscCallA(VecRestoreArray(X, lx_v, ierr))
164e7a95102SMartin Diehl
165e7a95102SMartin Diehl!  Insert values into global vector
166e7a95102SMartin Diehl
167e7a95102SMartin Diehl  end
168e7a95102SMartin Diehl
169e7a95102SMartin Diehl! ---------------------------------------------------------------------
170e7a95102SMartin Diehl!
171e7a95102SMartin Diehl!  InitialGuessLocal - Computes initial approximation, called by
172e7a95102SMartin Diehl!  the higher level routine FormInitialGuess().
173e7a95102SMartin Diehl!
174e7a95102SMartin Diehl!  Input Parameter:
175e7a95102SMartin Diehl!  x - local vector data
176e7a95102SMartin Diehl!
177e7a95102SMartin Diehl!  Output Parameters:
178e7a95102SMartin Diehl!  x - local vector data
179e7a95102SMartin Diehl!  ierr - error code
180e7a95102SMartin Diehl!
181e7a95102SMartin Diehl!  Notes:
182e7a95102SMartin Diehl!  This routine uses standard Fortran-style computations over a 2-dim array.
183e7a95102SMartin Diehl!
184*2a8381b2SBarry Smith  subroutine InitialGuessLocal(ctx, x, ierr)
185e7a95102SMartin Diehl!  Input/output variables:
186*2a8381b2SBarry Smith    type(AppCtx) ctx
187*2a8381b2SBarry Smith    PetscScalar x(ctx%xs:ctx%xe, ctx%ys:ctx%ye)
188e7a95102SMartin Diehl    PetscErrorCode ierr
189e7a95102SMartin Diehl
190e7a95102SMartin Diehl!  Local variables:
191e7a95102SMartin Diehl    PetscInt i, j
192e7a95102SMartin Diehl    PetscReal temp1, temp, hx, hy
193e7a95102SMartin Diehl    PetscReal one
194e7a95102SMartin Diehl
195e7a95102SMartin Diehl!  Set parameters
196e7a95102SMartin Diehl
197e7a95102SMartin Diehl    ierr = 0
198e7a95102SMartin Diehl    one = 1.0
199*2a8381b2SBarry Smith    hx = one/(ctx%mx - 1)
200*2a8381b2SBarry Smith    hy = one/(ctx%my - 1)
201*2a8381b2SBarry Smith    temp1 = ctx%lambda/(ctx%lambda + one)
202e7a95102SMartin Diehl
203*2a8381b2SBarry Smith    do j = ctx%ys, ctx%ye
204*2a8381b2SBarry Smith      temp = min(j - 1, ctx%my - j)*hy
205*2a8381b2SBarry Smith      do i = ctx%xs, ctx%xe
206*2a8381b2SBarry Smith        if (i == 1 .or. j == 1 .or. i == ctx%mx .or. j == ctx%my) then
207e7a95102SMartin Diehl          x(i, j) = 0.0
208e7a95102SMartin Diehl        else
209*2a8381b2SBarry Smith          x(i, j) = temp1*sqrt(min(hx*min(i - 1, ctx%mx - i), temp))
210e7a95102SMartin Diehl        end if
211e7a95102SMartin Diehl      end do
212e7a95102SMartin Diehl    end do
213e7a95102SMartin Diehl
214e7a95102SMartin Diehl  end
215e7a95102SMartin Diehl
216e7a95102SMartin Diehl! ---------------------------------------------------------------------
217e7a95102SMartin Diehl!
218e7a95102SMartin Diehl!  FormFunctionLocal - Computes nonlinear function, called by
219e7a95102SMartin Diehl!  the higher level routine FormFunction().
220e7a95102SMartin Diehl!
221e7a95102SMartin Diehl!  Input Parameter:
222e7a95102SMartin Diehl!  x - local vector data
223e7a95102SMartin Diehl!
224e7a95102SMartin Diehl!  Output Parameters:
225e7a95102SMartin Diehl!  f - local vector data, f(x)
226e7a95102SMartin Diehl!  ierr - error code
227e7a95102SMartin Diehl!
228e7a95102SMartin Diehl!  Notes:
229e7a95102SMartin Diehl!  This routine uses standard Fortran-style computations over a 2-dim array.
230e7a95102SMartin Diehl!
231*2a8381b2SBarry Smith  subroutine FormFunctionLocal(x, f, ctx, ierr)
232e7a95102SMartin Diehl!  Input/output variables:
233*2a8381b2SBarry Smith    type(AppCtx) ctx
234*2a8381b2SBarry Smith    PetscScalar x(ctx%gxs:ctx%gxe, ctx%gys:ctx%gye)
235*2a8381b2SBarry Smith    PetscScalar f(ctx%xs:ctx%xe, ctx%ys:ctx%ye)
236e7a95102SMartin Diehl    PetscErrorCode ierr
237e7a95102SMartin Diehl
238e7a95102SMartin Diehl!  Local variables:
239e7a95102SMartin Diehl    PetscScalar two, one, hx, hy, hxdhy, hydhx, sc
240e7a95102SMartin Diehl    PetscScalar u, uxx, uyy
241e7a95102SMartin Diehl    PetscInt i, j
242e7a95102SMartin Diehl
243e7a95102SMartin Diehl    one = 1.0
244e7a95102SMartin Diehl    two = 2.0
245*2a8381b2SBarry Smith    hx = one/(ctx%mx - 1)
246*2a8381b2SBarry Smith    hy = one/(ctx%my - 1)
247*2a8381b2SBarry Smith    sc = hx*hy*ctx%lambda
248e7a95102SMartin Diehl    hxdhy = hx/hy
249e7a95102SMartin Diehl    hydhx = hy/hx
250e7a95102SMartin Diehl
251e7a95102SMartin Diehl!  Compute function over the locally owned part of the grid
252e7a95102SMartin Diehl
253*2a8381b2SBarry Smith    do j = ctx%ys, ctx%ye
254*2a8381b2SBarry Smith      do i = ctx%xs, ctx%xe
255*2a8381b2SBarry Smith        if (i == 1 .or. j == 1 .or. i == ctx%mx .or. j == ctx%my) then
256e7a95102SMartin Diehl          f(i, j) = x(i, j)
257e7a95102SMartin Diehl        else
258e7a95102SMartin Diehl          u = x(i, j)
259e7a95102SMartin Diehl          uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j))
260e7a95102SMartin Diehl          uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1))
261e7a95102SMartin Diehl          f(i, j) = uxx + uyy - sc*exp(u)
262e7a95102SMartin Diehl        end if
263e7a95102SMartin Diehl      end do
264e7a95102SMartin Diehl    end do
265e7a95102SMartin Diehl
266e7a95102SMartin Diehl  end
267e7a95102SMartin Diehl
268e7a95102SMartin Diehl! ---------------------------------------------------------------------
269e7a95102SMartin Diehl!
270e7a95102SMartin Diehl!  FormJacobian - Evaluates Jacobian matrix.
271e7a95102SMartin Diehl!
272e7a95102SMartin Diehl!  Input Parameters:
273e7a95102SMartin Diehl!  snes     - the SNES context
274e7a95102SMartin Diehl!  x        - input vector
275e7a95102SMartin Diehl!  dummy    - optional user-defined context, as set by SNESSetJacobian()
276e7a95102SMartin Diehl!             (not used here)
277e7a95102SMartin Diehl!
278e7a95102SMartin Diehl!  Output Parameters:
279e7a95102SMartin Diehl!  jac      - Jacobian matrix
280e7a95102SMartin Diehl!  jac_prec - optionally different matrix used to construct the preconditioner (not used here)
281e7a95102SMartin Diehl!
282e7a95102SMartin Diehl!  Notes:
283e7a95102SMartin Diehl!  This routine serves as a wrapper for the lower-level routine
284e7a95102SMartin Diehl!  "FormJacobianLocal", where the actual computations are
285e7a95102SMartin Diehl!  done using the standard Fortran style of treating the local
286e7a95102SMartin Diehl!  vector data as a multidimensional array over the local mesh.
287e7a95102SMartin Diehl!  This routine merely accesses the local vector data via
288e7a95102SMartin Diehl!  VecGetArray() and VecRestoreArray().
289e7a95102SMartin Diehl!
290e7a95102SMartin Diehl!  Notes:
291e7a95102SMartin Diehl!  Due to grid point reordering with DMDAs, we must always work
292e7a95102SMartin Diehl!  with the local grid points, and then transform them to the new
293e7a95102SMartin Diehl!  global numbering with the "ltog" mapping
294e7a95102SMartin Diehl!  We cannot work directly with the global numbers for the original
295e7a95102SMartin Diehl!  uniprocessor grid!
296e7a95102SMartin Diehl!
297e7a95102SMartin Diehl!  Two methods are available for imposing this transformation
298e7a95102SMartin Diehl!  when setting matrix entries:
299e7a95102SMartin Diehl!    (A) MatSetValuesLocal(), using the local ordering (including
300e7a95102SMartin Diehl!        ghost points!)
301e7a95102SMartin Diehl!        - Set matrix entries using the local ordering
302e7a95102SMartin Diehl!          by calling MatSetValuesLocal()
303e7a95102SMartin Diehl!    (B) MatSetValues(), using the global ordering
304e7a95102SMartin Diehl
305e7a95102SMartin Diehl!        - Set matrix entries using the global ordering by calling
306e7a95102SMartin Diehl!          MatSetValues()
307e7a95102SMartin Diehl!  Option (A) seems cleaner/easier in many cases, and is the procedure
308e7a95102SMartin Diehl!  used in this example.
309e7a95102SMartin Diehl!
310*2a8381b2SBarry Smith  subroutine FormJacobian(snes, X, jac, jac_prec, ctx, ierr)
311e7a95102SMartin Diehl!  Input/output variables:
312e7a95102SMartin Diehl    SNES snes
313e7a95102SMartin Diehl    Vec X
314e7a95102SMartin Diehl    Mat jac, jac_prec
315*2a8381b2SBarry Smith    type(AppCtx) ctx
316e7a95102SMartin Diehl    PetscErrorCode ierr
317e7a95102SMartin Diehl    DM da
318e7a95102SMartin Diehl
319e7a95102SMartin Diehl!  Declarations for use with local arrays:
320e7a95102SMartin Diehl    PetscScalar, pointer :: lx_v(:)
321e7a95102SMartin Diehl    Vec localX
322e7a95102SMartin Diehl
323e7a95102SMartin Diehl!  Scatter ghost points to local vector, using the 2-step process
324e7a95102SMartin Diehl!     DMGlobalToLocalBegin(), DMGlobalToLocalEnd()
325e7a95102SMartin Diehl!  Computations can be done while messages are in transition,
326e7a95102SMartin Diehl!  by placing code between these two statements.
327e7a95102SMartin Diehl
328e7a95102SMartin Diehl    PetscCallA(SNESGetDM(snes, da, ierr))
329e7a95102SMartin Diehl    PetscCallA(DMGetLocalVector(da, localX, ierr))
330e7a95102SMartin Diehl    PetscCallA(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr))
331e7a95102SMartin Diehl    PetscCallA(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr))
332e7a95102SMartin Diehl
333e7a95102SMartin Diehl!  Get a pointer to vector data
334e7a95102SMartin Diehl    PetscCallA(VecGetArray(localX, lx_v, ierr))
335e7a95102SMartin Diehl
336e7a95102SMartin Diehl!  Compute entries for the locally owned part of the Jacobian preconditioner.
337*2a8381b2SBarry Smith    PetscCallA(FormJacobianLocal(lx_v, jac_prec, ctx, ierr))
338e7a95102SMartin Diehl
339e7a95102SMartin Diehl!  Assemble matrix, using the 2-step process:
340e7a95102SMartin Diehl!     MatAssemblyBegin(), MatAssemblyEnd()
341e7a95102SMartin Diehl!  Computations can be done while messages are in transition,
342e7a95102SMartin Diehl!  by placing code between these two statements.
343e7a95102SMartin Diehl
344e7a95102SMartin Diehl    PetscCallA(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
345e7a95102SMartin Diehl    if (jac /= jac_prec) then
346e7a95102SMartin Diehl      PetscCallA(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
347e7a95102SMartin Diehl    end if
348e7a95102SMartin Diehl    PetscCallA(VecRestoreArray(localX, lx_v, ierr))
349e7a95102SMartin Diehl    PetscCallA(DMRestoreLocalVector(da, localX, ierr))
350e7a95102SMartin Diehl    PetscCallA(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
351e7a95102SMartin Diehl    if (jac /= jac_prec) then
352e7a95102SMartin Diehl      PetscCallA(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
353e7a95102SMartin Diehl    end if
354e7a95102SMartin Diehl
355e7a95102SMartin Diehl!  Tell the matrix we will never add a new nonzero location to the
356e7a95102SMartin Diehl!  matrix. If we do it will generate an error.
357e7a95102SMartin Diehl
358e7a95102SMartin Diehl    PetscCallA(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr))
359e7a95102SMartin Diehl
360e7a95102SMartin Diehl  end
361e7a95102SMartin Diehl
362e7a95102SMartin Diehl! ---------------------------------------------------------------------
363e7a95102SMartin Diehl!
364e7a95102SMartin Diehl!  FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner,
365e7a95102SMartin Diehl!  called by the higher level routine FormJacobian().
366e7a95102SMartin Diehl!
367e7a95102SMartin Diehl!  Input Parameters:
368e7a95102SMartin Diehl!  x        - local vector data
369e7a95102SMartin Diehl!
370e7a95102SMartin Diehl!  Output Parameters:
371e7a95102SMartin Diehl!  jac_prec - Jacobian matrix used to compute the preconditioner
372e7a95102SMartin Diehl!  ierr     - error code
373e7a95102SMartin Diehl!
374e7a95102SMartin Diehl!  Notes:
375e7a95102SMartin Diehl!  This routine uses standard Fortran-style computations over a 2-dim array.
376e7a95102SMartin Diehl!
377e7a95102SMartin Diehl!  Notes:
378e7a95102SMartin Diehl!  Due to grid point reordering with DMDAs, we must always work
379e7a95102SMartin Diehl!  with the local grid points, and then transform them to the new
380e7a95102SMartin Diehl!  global numbering with the "ltog" mapping
381e7a95102SMartin Diehl!  We cannot work directly with the global numbers for the original
382e7a95102SMartin Diehl!  uniprocessor grid!
383e7a95102SMartin Diehl!
384e7a95102SMartin Diehl!  Two methods are available for imposing this transformation
385e7a95102SMartin Diehl!  when setting matrix entries:
386e7a95102SMartin Diehl!    (A) MatSetValuesLocal(), using the local ordering (including
387e7a95102SMartin Diehl!        ghost points!)
388e7a95102SMartin Diehl!        - Set matrix entries using the local ordering
389e7a95102SMartin Diehl!          by calling MatSetValuesLocal()
390e7a95102SMartin Diehl!    (B) MatSetValues(), using the global ordering
391e7a95102SMartin Diehl!        - Then apply this map explicitly yourself
392e7a95102SMartin Diehl!        - Set matrix entries using the global ordering by calling
393e7a95102SMartin Diehl!          MatSetValues()
394e7a95102SMartin Diehl!  Option (A) seems cleaner/easier in many cases, and is the procedure
395e7a95102SMartin Diehl!  used in this example.
396e7a95102SMartin Diehl!
397*2a8381b2SBarry Smith  subroutine FormJacobianLocal(x, jac_prec, ctx, ierr)
398e7a95102SMartin Diehl!  Input/output variables:
399*2a8381b2SBarry Smith    type(AppCtx) ctx
400*2a8381b2SBarry Smith    PetscScalar x(ctx%gxs:ctx%gxe, ctx%gys:ctx%gye)
401e7a95102SMartin Diehl    Mat jac_prec
402e7a95102SMartin Diehl    PetscErrorCode ierr
403e7a95102SMartin Diehl
404e7a95102SMartin Diehl!  Local variables:
405e7a95102SMartin Diehl    PetscInt row, col(5), i, j
406e7a95102SMartin Diehl    PetscInt ione, ifive
407e7a95102SMartin Diehl    PetscScalar two, one, hx, hy, hxdhy
408e7a95102SMartin Diehl    PetscScalar hydhx, sc, v(5)
409e7a95102SMartin Diehl
410e7a95102SMartin Diehl!  Set parameters
411e7a95102SMartin Diehl    ione = 1
412e7a95102SMartin Diehl    ifive = 5
413e7a95102SMartin Diehl    one = 1.0
414e7a95102SMartin Diehl    two = 2.0
415*2a8381b2SBarry Smith    hx = one/(ctx%mx - 1)
416*2a8381b2SBarry Smith    hy = one/(ctx%my - 1)
417e7a95102SMartin Diehl    sc = hx*hy
418e7a95102SMartin Diehl    hxdhy = hx/hy
419e7a95102SMartin Diehl    hydhx = hy/hx
420e7a95102SMartin Diehl
421e7a95102SMartin Diehl!  Compute entries for the locally owned part of the Jacobian.
422e7a95102SMartin Diehl!   - Currently, all PETSc parallel matrix formats are partitioned by
423e7a95102SMartin Diehl!     contiguous chunks of rows across the processors.
424e7a95102SMartin Diehl!   - Each processor needs to insert only elements that it owns
425e7a95102SMartin Diehl!     locally (but any non-local elements will be sent to the
426e7a95102SMartin Diehl!     appropriate processor during matrix assembly).
427e7a95102SMartin Diehl!   - Here, we set all entries for a particular row at once.
428e7a95102SMartin Diehl!   - We can set matrix entries either using either
429e7a95102SMartin Diehl!     MatSetValuesLocal() or MatSetValues(), as discussed above.
430e7a95102SMartin Diehl!   - Note that MatSetValues() uses 0-based row and column numbers
431e7a95102SMartin Diehl!     in Fortran as well as in C.
432e7a95102SMartin Diehl
433*2a8381b2SBarry Smith    do j = ctx%ys, ctx%ye
434*2a8381b2SBarry Smith      row = (j - ctx%gys)*ctx%gxm + ctx%xs - ctx%gxs - 1
435*2a8381b2SBarry Smith      do i = ctx%xs, ctx%xe
436e7a95102SMartin Diehl        row = row + 1
437e7a95102SMartin Diehl!           boundary points
438*2a8381b2SBarry Smith        if (i == 1 .or. j == 1 .or. i == ctx%mx .or. j == ctx%my) then
439e7a95102SMartin Diehl          col(1) = row
440e7a95102SMartin Diehl          v(1) = one
441e7a95102SMartin Diehl          PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ione, col, v, INSERT_VALUES, ierr))
442e7a95102SMartin Diehl!           interior grid points
443e7a95102SMartin Diehl        else
444e7a95102SMartin Diehl          v(1) = -hxdhy
445e7a95102SMartin Diehl          v(2) = -hydhx
446*2a8381b2SBarry Smith          v(3) = two*(hydhx + hxdhy) - sc*ctx%lambda*exp(x(i, j))
447e7a95102SMartin Diehl          v(4) = -hydhx
448e7a95102SMartin Diehl          v(5) = -hxdhy
449*2a8381b2SBarry Smith          col(1) = row - ctx%gxm
450e7a95102SMartin Diehl          col(2) = row - 1
451e7a95102SMartin Diehl          col(3) = row
452e7a95102SMartin Diehl          col(4) = row + 1
453*2a8381b2SBarry Smith          col(5) = row + ctx%gxm
454e7a95102SMartin Diehl          PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ifive, col, v, INSERT_VALUES, ierr))
455e7a95102SMartin Diehl        end if
456e7a95102SMartin Diehl      end do
457e7a95102SMartin Diehl    end do
458e7a95102SMartin Diehl
459e7a95102SMartin Diehl  end
460e7a95102SMartin Diehl
461b06eb4cdSBarry Smithend module ex5module
462c4762a1bSJed Brown
463c4762a1bSJed Brownprogram main
464b06eb4cdSBarry Smith  use ex5module
465c4762a1bSJed Brown  implicit none
466c4762a1bSJed Brown!
467c4762a1bSJed Brown
468c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
469c4762a1bSJed Brown!                   Variable declarations
470c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
471c4762a1bSJed Brown!
472c4762a1bSJed Brown!  Variables:
473c4762a1bSJed Brown!     snes        - nonlinear solver
474c4762a1bSJed Brown!     x, r        - solution, residual vectors
475c4762a1bSJed Brown!     J           - Jacobian matrix
476c4762a1bSJed Brown!     its         - iterations for convergence
477c4762a1bSJed Brown!     Nx, Ny      - number of preocessors in x- and y- directions
478c4762a1bSJed Brown!     matrix_free - flag - 1 indicates matrix-free version
479c4762a1bSJed Brown!
480c4762a1bSJed Brown  SNES snes
481c4762a1bSJed Brown  Vec x, r
482c4762a1bSJed Brown  Mat J
483c4762a1bSJed Brown  PetscErrorCode ierr
484c4762a1bSJed Brown  PetscInt its
485c4762a1bSJed Brown  PetscBool flg, matrix_free
486c4762a1bSJed Brown  PetscInt ione, nfour
487c4762a1bSJed Brown  PetscReal lambda_max, lambda_min
488*2a8381b2SBarry Smith  type(AppCtx) ctx
489c4762a1bSJed Brown  DM da
490c4762a1bSJed Brown
491c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
492c4762a1bSJed Brown!  Initialize program
493c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
494d8606c27SBarry Smith  PetscCallA(PetscInitialize(ierr))
495*2a8381b2SBarry Smith  PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, ctx%rank, ierr))
496c4762a1bSJed Brown
497c4762a1bSJed Brown!  Initialize problem parameters
498c4762a1bSJed Brown  lambda_max = 6.81
499c4762a1bSJed Brown  lambda_min = 0.0
500*2a8381b2SBarry Smith  ctx%lambda = 6.0
501c4762a1bSJed Brown  ione = 1
502c4762a1bSJed Brown  nfour = 4
503*2a8381b2SBarry Smith  PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', ctx%lambda, flg, ierr))
504*2a8381b2SBarry Smith  PetscCheckA(ctx%lambda < lambda_max .and. ctx%lambda > lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda provided with -par is out of range')
505c4762a1bSJed Brown
506c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
507c4762a1bSJed Brown!  Create nonlinear solver context
508c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
509d8606c27SBarry Smith  PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))
510c4762a1bSJed Brown
511c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
512c4762a1bSJed Brown!  Create vector data structures; set function evaluation routine
513c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
514c4762a1bSJed Brown
515c4762a1bSJed Brown!  Create distributed array (DMDA) to manage parallel grid and vectors
516c4762a1bSJed Brown
517c4762a1bSJed Brown! This really needs only the star-type stencil, but we use the box
518c4762a1bSJed Brown! stencil temporarily.
5195d83a8b1SBarry Smith  PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, nfour, nfour, PETSC_DECIDE, PETSC_DECIDE, ione, ione, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, da, ierr))
520d8606c27SBarry Smith  PetscCallA(DMSetFromOptions(da, ierr))
521d8606c27SBarry Smith  PetscCallA(DMSetUp(da, ierr))
522c4762a1bSJed Brown
523*2a8381b2SBarry Smith  PetscCallA(DMDAGetInfo(da, PETSC_NULL_INTEGER, ctx%mx, ctx%my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))
524c4762a1bSJed Brown
525c4762a1bSJed Brown!
526c4762a1bSJed Brown!   Visualize the distribution of the array across the processors
527c4762a1bSJed Brown!
528d8606c27SBarry Smith!     PetscCallA(DMView(da,PETSC_VIEWER_DRAW_WORLD,ierr))
529c4762a1bSJed Brown
530c4762a1bSJed Brown!  Extract global and local vectors from DMDA; then duplicate for remaining
531c4762a1bSJed Brown!  vectors that are the same types
532d8606c27SBarry Smith  PetscCallA(DMCreateGlobalVector(da, x, ierr))
533d8606c27SBarry Smith  PetscCallA(VecDuplicate(x, r, ierr))
534c4762a1bSJed Brown
535c4762a1bSJed Brown!  Get local grid boundaries (for 2-dimensional DMDA)
536*2a8381b2SBarry Smith  PetscCallA(DMDAGetCorners(da, ctx%xs, ctx%ys, PETSC_NULL_INTEGER, ctx%xm, ctx%ym, PETSC_NULL_INTEGER, ierr))
537*2a8381b2SBarry Smith  PetscCallA(DMDAGetGhostCorners(da, ctx%gxs, ctx%gys, PETSC_NULL_INTEGER, ctx%gxm, ctx%gym, PETSC_NULL_INTEGER, ierr))
538c4762a1bSJed Brown
539c4762a1bSJed Brown!  Here we shift the starting indices up by one so that we can easily
540c4762a1bSJed Brown!  use the Fortran convention of 1-based indices (rather 0-based indices).
541*2a8381b2SBarry Smith  ctx%xs = ctx%xs + 1
542*2a8381b2SBarry Smith  ctx%ys = ctx%ys + 1
543*2a8381b2SBarry Smith  ctx%gxs = ctx%gxs + 1
544*2a8381b2SBarry Smith  ctx%gys = ctx%gys + 1
545c4762a1bSJed Brown
546*2a8381b2SBarry Smith  ctx%ye = ctx%ys + ctx%ym - 1
547*2a8381b2SBarry Smith  ctx%xe = ctx%xs + ctx%xm - 1
548*2a8381b2SBarry Smith  ctx%gye = ctx%gys + ctx%gym - 1
549*2a8381b2SBarry Smith  ctx%gxe = ctx%gxs + ctx%gxm - 1
550c4762a1bSJed Brown
551*2a8381b2SBarry Smith  PetscCallA(SNESSetApplicationContext(snes, ctx, ierr))
552c4762a1bSJed Brown
553c4762a1bSJed Brown!  Set function evaluation routine and vector
554*2a8381b2SBarry Smith  PetscCallA(SNESSetFunction(snes, r, FormFunction, ctx, ierr))
555c4762a1bSJed Brown
556c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
557c4762a1bSJed Brown!  Create matrix data structure; set Jacobian evaluation routine
558c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
559c4762a1bSJed Brown
560c4762a1bSJed Brown!  Set Jacobian matrix data structure and default Jacobian evaluation
561c4762a1bSJed Brown!  routine. User can override with:
562c4762a1bSJed Brown!     -snes_fd : default finite differencing approximation of Jacobian
563c4762a1bSJed Brown!     -snes_mf : matrix-free Newton-Krylov method with no preconditioning
564c4762a1bSJed Brown!                (unless user explicitly sets preconditioner)
5657addb90fSBarry Smith!     -snes_mf_operator : form matrix used to construct the preconditioner as set by the user,
566c4762a1bSJed Brown!                         but use matrix-free approx for Jacobian-vector
567c4762a1bSJed Brown!                         products within Newton-Krylov method
568c4762a1bSJed Brown!
569c4762a1bSJed Brown!  Note:  For the parallel case, vectors and matrices MUST be partitioned
570c4762a1bSJed Brown!     accordingly.  When using distributed arrays (DMDAs) to create vectors,
571c4762a1bSJed Brown!     the DMDAs determine the problem partitioning.  We must explicitly
572c4762a1bSJed Brown!     specify the local matrix dimensions upon its creation for compatibility
573c4762a1bSJed Brown!     with the vector distribution.  Thus, the generic MatCreate() routine
574c4762a1bSJed Brown!     is NOT sufficient when working with distributed arrays.
575c4762a1bSJed Brown!
576c4762a1bSJed Brown!     Note: Here we only approximately preallocate storage space for the
577c4762a1bSJed Brown!     Jacobian.  See the users manual for a discussion of better techniques
578c4762a1bSJed Brown!     for preallocating matrix memory.
579c4762a1bSJed Brown
580d8606c27SBarry Smith  PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-snes_mf', matrix_free, ierr))
581c4762a1bSJed Brown  if (.not. matrix_free) then
582d8606c27SBarry Smith    PetscCallA(DMSetMatType(da, MATAIJ, ierr))
583d8606c27SBarry Smith    PetscCallA(DMCreateMatrix(da, J, ierr))
584*2a8381b2SBarry Smith    PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, ctx, ierr))
585c4762a1bSJed Brown  end if
586c4762a1bSJed Brown
587c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
588c4762a1bSJed Brown!  Customize nonlinear solver; set runtime options
589c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
590c4762a1bSJed Brown!  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
591d8606c27SBarry Smith  PetscCallA(SNESSetDM(snes, da, ierr))
592d8606c27SBarry Smith  PetscCallA(SNESSetFromOptions(snes, ierr))
593c4762a1bSJed Brown
594c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
595c4762a1bSJed Brown!  Evaluate initial guess; then solve nonlinear system.
596c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
597c4762a1bSJed Brown!  Note: The user should initialize the vector, x, with the initial guess
598c4762a1bSJed Brown!  for the nonlinear solver prior to calling SNESSolve().  In particular,
599c4762a1bSJed Brown!  to employ an initial guess of zero, the user should explicitly set
600c4762a1bSJed Brown!  this vector to zero by calling VecSet().
601c4762a1bSJed Brown
602d8606c27SBarry Smith  PetscCallA(FormInitialGuess(snes, x, ierr))
603d8606c27SBarry Smith  PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
604d8606c27SBarry Smith  PetscCallA(SNESGetIterationNumber(snes, its, ierr))
605*2a8381b2SBarry Smith  if (ctx%rank == 0) then
606c4762a1bSJed Brown    write (6, 100) its
607c4762a1bSJed Brown  end if
608c4762a1bSJed Brown100 format('Number of SNES iterations = ', i5)
609c4762a1bSJed Brown
610c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
611c4762a1bSJed Brown!  Free work space.  All PETSc objects should be destroyed when they
612c4762a1bSJed Brown!  are no longer needed.
613c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
614d8606c27SBarry Smith  if (.not. matrix_free) PetscCallA(MatDestroy(J, ierr))
615d8606c27SBarry Smith  PetscCallA(VecDestroy(x, ierr))
616d8606c27SBarry Smith  PetscCallA(VecDestroy(r, ierr))
617d8606c27SBarry Smith  PetscCallA(SNESDestroy(snes, ierr))
618d8606c27SBarry Smith  PetscCallA(DMDestroy(da, ierr))
619c4762a1bSJed Brown
620d8606c27SBarry Smith  PetscCallA(PetscFinalize(ierr))
621c4762a1bSJed Brownend
622c4762a1bSJed Brown!
623c4762a1bSJed Brown!/*TEST
624c4762a1bSJed Brown!
625c4762a1bSJed Brown!   test:
626c4762a1bSJed Brown!      nsize: 4
627c4762a1bSJed Brown!      args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
628c4762a1bSJed Brown!      requires: !single
629c4762a1bSJed Brown!
630c4762a1bSJed Brown!   test:
631c4762a1bSJed Brown!      suffix: 2
632c4762a1bSJed Brown!      nsize: 4
633c4762a1bSJed Brown!      args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
634c4762a1bSJed Brown!      requires: !single
635c4762a1bSJed Brown!
636c4762a1bSJed Brown!   test:
637c4762a1bSJed Brown!      suffix: 3
638c4762a1bSJed Brown!      nsize: 3
639c4762a1bSJed Brown!      args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
640c4762a1bSJed Brown!      requires: !single
641c4762a1bSJed Brown!
642c4762a1bSJed Brown!   test:
643c4762a1bSJed Brown!      suffix: 4
644c4762a1bSJed Brown!      nsize: 3
645c4762a1bSJed Brown!      args: -snes_mf_operator -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
646c4762a1bSJed Brown!      requires: !single
647c4762a1bSJed Brown!
648c4762a1bSJed Brown!   test:
649c4762a1bSJed Brown!      suffix: 5
650c4762a1bSJed Brown!      requires: !single
651c4762a1bSJed Brown!
652c4762a1bSJed Brown!TEST*/
653