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