xref: /petsc/src/snes/tutorials/ex5f90.F90 (revision e7a95102f46630f317be643b805dc1c3f4655aeb)
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>
35c5e229c2SMartin Diehlmodule ex5f90module
36c4762a1bSJed Brown  use petscsnes
37dfbbaf82SBarry Smith  use petscdmda
38*e7a95102SMartin 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
47*e7a95102SMartin Diehl  interface
48*e7a95102SMartin Diehl    subroutine SNESSetApplicationContext(snes, ctx, ierr)
49*e7a95102SMartin Diehl      use petscsnes
50*e7a95102SMartin Diehl      import userctx
51*e7a95102SMartin Diehl      implicit none
52*e7a95102SMartin Diehl      SNES snes
53*e7a95102SMartin Diehl      type(userctx) ctx
54*e7a95102SMartin Diehl      PetscErrorCode ierr
55*e7a95102SMartin Diehl    end subroutine
56*e7a95102SMartin Diehl    subroutine SNESGetApplicationContext(snes, ctx, ierr)
57*e7a95102SMartin Diehl      use petscsnes
58*e7a95102SMartin Diehl      import userctx
59*e7a95102SMartin Diehl      implicit none
60*e7a95102SMartin Diehl      SNES snes
61*e7a95102SMartin Diehl      type(userctx), pointer :: ctx
62*e7a95102SMartin Diehl      PetscErrorCode ierr
63*e7a95102SMartin Diehl    end subroutine
64*e7a95102SMartin Diehl  end interface
65*e7a95102SMartin 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
136*e7a95102SMartin Diehl
137*e7a95102SMartin Diehl! ---------------------------------------------------------------------
138*e7a95102SMartin Diehl!
139*e7a95102SMartin Diehl!  FormInitialGuess - Forms initial approximation.
140*e7a95102SMartin Diehl!
141*e7a95102SMartin Diehl!  Input Parameters:
142*e7a95102SMartin Diehl!  X - vector
143*e7a95102SMartin Diehl!
144*e7a95102SMartin Diehl!  Output Parameter:
145*e7a95102SMartin Diehl!  X - vector
146*e7a95102SMartin Diehl!
147*e7a95102SMartin Diehl!  Notes:
148*e7a95102SMartin Diehl!  This routine serves as a wrapper for the lower-level routine
149*e7a95102SMartin Diehl!  "InitialGuessLocal", where the actual computations are
150*e7a95102SMartin Diehl!  done using the standard Fortran style of treating the local
151*e7a95102SMartin Diehl!  vector data as a multidimensional array over the local mesh.
152*e7a95102SMartin Diehl!  This routine merely handles ghost point scatters and accesses
153*e7a95102SMartin Diehl!  the local vector data via VecGetArray() and VecRestoreArray().
154*e7a95102SMartin Diehl!
155*e7a95102SMartin Diehl  subroutine FormInitialGuess(snes, X, ierr)
156*e7a95102SMartin Diehl!  Input/output variables:
157*e7a95102SMartin Diehl    SNES snes
158*e7a95102SMartin Diehl    type(userctx), pointer:: puser
159*e7a95102SMartin Diehl    Vec X
160*e7a95102SMartin Diehl    PetscErrorCode ierr
161*e7a95102SMartin Diehl    DM da
162*e7a95102SMartin Diehl
163*e7a95102SMartin Diehl!  Declarations for use with local arrays:
164*e7a95102SMartin Diehl    PetscScalar, pointer :: lx_v(:)
165*e7a95102SMartin Diehl
166*e7a95102SMartin Diehl    ierr = 0
167*e7a95102SMartin Diehl    PetscCallA(SNESGetDM(snes, da, ierr))
168*e7a95102SMartin Diehl    PetscCallA(SNESGetApplicationContext(snes, puser, ierr))
169*e7a95102SMartin Diehl!  Get a pointer to vector data.
170*e7a95102SMartin Diehl!    - For default PETSc vectors, VecGetArray() returns a pointer to
171*e7a95102SMartin Diehl!      the data array. Otherwise, the routine is implementation dependent.
172*e7a95102SMartin Diehl!    - You MUST call VecRestoreArray() when you no longer need access to
173*e7a95102SMartin Diehl!      the array.
174*e7a95102SMartin Diehl!    - Note that the interface to VecGetArray() differs from VecGetArray().
175*e7a95102SMartin Diehl
176*e7a95102SMartin Diehl    PetscCallA(VecGetArray(X, lx_v, ierr))
177*e7a95102SMartin Diehl
178*e7a95102SMartin Diehl!  Compute initial guess over the locally owned part of the grid
179*e7a95102SMartin Diehl    PetscCallA(InitialGuessLocal(puser, lx_v, ierr))
180*e7a95102SMartin Diehl
181*e7a95102SMartin Diehl!  Restore vector
182*e7a95102SMartin Diehl    PetscCallA(VecRestoreArray(X, lx_v, ierr))
183*e7a95102SMartin Diehl
184*e7a95102SMartin Diehl!  Insert values into global vector
185*e7a95102SMartin Diehl
186*e7a95102SMartin Diehl  end
187*e7a95102SMartin Diehl
188*e7a95102SMartin Diehl! ---------------------------------------------------------------------
189*e7a95102SMartin Diehl!
190*e7a95102SMartin Diehl!  InitialGuessLocal - Computes initial approximation, called by
191*e7a95102SMartin Diehl!  the higher level routine FormInitialGuess().
192*e7a95102SMartin Diehl!
193*e7a95102SMartin Diehl!  Input Parameter:
194*e7a95102SMartin Diehl!  x - local vector data
195*e7a95102SMartin Diehl!
196*e7a95102SMartin Diehl!  Output Parameters:
197*e7a95102SMartin Diehl!  x - local vector data
198*e7a95102SMartin Diehl!  ierr - error code
199*e7a95102SMartin Diehl!
200*e7a95102SMartin Diehl!  Notes:
201*e7a95102SMartin Diehl!  This routine uses standard Fortran-style computations over a 2-dim array.
202*e7a95102SMartin Diehl!
203*e7a95102SMartin Diehl  subroutine InitialGuessLocal(user, x, ierr)
204*e7a95102SMartin Diehl!  Input/output variables:
205*e7a95102SMartin Diehl    type(userctx) user
206*e7a95102SMartin Diehl    PetscScalar x(user%xs:user%xe, user%ys:user%ye)
207*e7a95102SMartin Diehl    PetscErrorCode ierr
208*e7a95102SMartin Diehl
209*e7a95102SMartin Diehl!  Local variables:
210*e7a95102SMartin Diehl    PetscInt i, j
211*e7a95102SMartin Diehl    PetscReal temp1, temp, hx, hy
212*e7a95102SMartin Diehl    PetscReal one
213*e7a95102SMartin Diehl
214*e7a95102SMartin Diehl!  Set parameters
215*e7a95102SMartin Diehl
216*e7a95102SMartin Diehl    ierr = 0
217*e7a95102SMartin Diehl    one = 1.0
218*e7a95102SMartin Diehl    hx = one/(user%mx - 1)
219*e7a95102SMartin Diehl    hy = one/(user%my - 1)
220*e7a95102SMartin Diehl    temp1 = user%lambda/(user%lambda + one)
221*e7a95102SMartin Diehl
222*e7a95102SMartin Diehl    do j = user%ys, user%ye
223*e7a95102SMartin Diehl      temp = min(j - 1, user%my - j)*hy
224*e7a95102SMartin Diehl      do i = user%xs, user%xe
225*e7a95102SMartin Diehl        if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then
226*e7a95102SMartin Diehl          x(i, j) = 0.0
227*e7a95102SMartin Diehl        else
228*e7a95102SMartin Diehl          x(i, j) = temp1*sqrt(min(hx*min(i - 1, user%mx - i), temp))
229*e7a95102SMartin Diehl        end if
230*e7a95102SMartin Diehl      end do
231*e7a95102SMartin Diehl    end do
232*e7a95102SMartin Diehl
233*e7a95102SMartin Diehl  end
234*e7a95102SMartin Diehl
235*e7a95102SMartin Diehl! ---------------------------------------------------------------------
236*e7a95102SMartin Diehl!
237*e7a95102SMartin Diehl!  FormFunctionLocal - Computes nonlinear function, called by
238*e7a95102SMartin Diehl!  the higher level routine FormFunction().
239*e7a95102SMartin Diehl!
240*e7a95102SMartin Diehl!  Input Parameter:
241*e7a95102SMartin Diehl!  x - local vector data
242*e7a95102SMartin Diehl!
243*e7a95102SMartin Diehl!  Output Parameters:
244*e7a95102SMartin Diehl!  f - local vector data, f(x)
245*e7a95102SMartin Diehl!  ierr - error code
246*e7a95102SMartin Diehl!
247*e7a95102SMartin Diehl!  Notes:
248*e7a95102SMartin Diehl!  This routine uses standard Fortran-style computations over a 2-dim array.
249*e7a95102SMartin Diehl!
250*e7a95102SMartin Diehl  subroutine FormFunctionLocal(x, f, user, ierr)
251*e7a95102SMartin Diehl!  Input/output variables:
252*e7a95102SMartin Diehl    type(userctx) user
253*e7a95102SMartin Diehl    PetscScalar x(user%gxs:user%gxe, user%gys:user%gye)
254*e7a95102SMartin Diehl    PetscScalar f(user%xs:user%xe, user%ys:user%ye)
255*e7a95102SMartin Diehl    PetscErrorCode ierr
256*e7a95102SMartin Diehl
257*e7a95102SMartin Diehl!  Local variables:
258*e7a95102SMartin Diehl    PetscScalar two, one, hx, hy, hxdhy, hydhx, sc
259*e7a95102SMartin Diehl    PetscScalar u, uxx, uyy
260*e7a95102SMartin Diehl    PetscInt i, j
261*e7a95102SMartin Diehl
262*e7a95102SMartin Diehl    one = 1.0
263*e7a95102SMartin Diehl    two = 2.0
264*e7a95102SMartin Diehl    hx = one/(user%mx - 1)
265*e7a95102SMartin Diehl    hy = one/(user%my - 1)
266*e7a95102SMartin Diehl    sc = hx*hy*user%lambda
267*e7a95102SMartin Diehl    hxdhy = hx/hy
268*e7a95102SMartin Diehl    hydhx = hy/hx
269*e7a95102SMartin Diehl
270*e7a95102SMartin Diehl!  Compute function over the locally owned part of the grid
271*e7a95102SMartin Diehl
272*e7a95102SMartin Diehl    do j = user%ys, user%ye
273*e7a95102SMartin Diehl      do i = user%xs, user%xe
274*e7a95102SMartin Diehl        if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then
275*e7a95102SMartin Diehl          f(i, j) = x(i, j)
276*e7a95102SMartin Diehl        else
277*e7a95102SMartin Diehl          u = x(i, j)
278*e7a95102SMartin Diehl          uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j))
279*e7a95102SMartin Diehl          uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1))
280*e7a95102SMartin Diehl          f(i, j) = uxx + uyy - sc*exp(u)
281*e7a95102SMartin Diehl        end if
282*e7a95102SMartin Diehl      end do
283*e7a95102SMartin Diehl    end do
284*e7a95102SMartin Diehl
285*e7a95102SMartin Diehl  end
286*e7a95102SMartin Diehl
287*e7a95102SMartin Diehl! ---------------------------------------------------------------------
288*e7a95102SMartin Diehl!
289*e7a95102SMartin Diehl!  FormJacobian - Evaluates Jacobian matrix.
290*e7a95102SMartin Diehl!
291*e7a95102SMartin Diehl!  Input Parameters:
292*e7a95102SMartin Diehl!  snes     - the SNES context
293*e7a95102SMartin Diehl!  x        - input vector
294*e7a95102SMartin Diehl!  dummy    - optional user-defined context, as set by SNESSetJacobian()
295*e7a95102SMartin Diehl!             (not used here)
296*e7a95102SMartin Diehl!
297*e7a95102SMartin Diehl!  Output Parameters:
298*e7a95102SMartin Diehl!  jac      - Jacobian matrix
299*e7a95102SMartin Diehl!  jac_prec - optionally different matrix used to construct the preconditioner (not used here)
300*e7a95102SMartin Diehl!
301*e7a95102SMartin Diehl!  Notes:
302*e7a95102SMartin Diehl!  This routine serves as a wrapper for the lower-level routine
303*e7a95102SMartin Diehl!  "FormJacobianLocal", where the actual computations are
304*e7a95102SMartin Diehl!  done using the standard Fortran style of treating the local
305*e7a95102SMartin Diehl!  vector data as a multidimensional array over the local mesh.
306*e7a95102SMartin Diehl!  This routine merely accesses the local vector data via
307*e7a95102SMartin Diehl!  VecGetArray() and VecRestoreArray().
308*e7a95102SMartin Diehl!
309*e7a95102SMartin Diehl!  Notes:
310*e7a95102SMartin Diehl!  Due to grid point reordering with DMDAs, we must always work
311*e7a95102SMartin Diehl!  with the local grid points, and then transform them to the new
312*e7a95102SMartin Diehl!  global numbering with the "ltog" mapping
313*e7a95102SMartin Diehl!  We cannot work directly with the global numbers for the original
314*e7a95102SMartin Diehl!  uniprocessor grid!
315*e7a95102SMartin Diehl!
316*e7a95102SMartin Diehl!  Two methods are available for imposing this transformation
317*e7a95102SMartin Diehl!  when setting matrix entries:
318*e7a95102SMartin Diehl!    (A) MatSetValuesLocal(), using the local ordering (including
319*e7a95102SMartin Diehl!        ghost points!)
320*e7a95102SMartin Diehl!        - Set matrix entries using the local ordering
321*e7a95102SMartin Diehl!          by calling MatSetValuesLocal()
322*e7a95102SMartin Diehl!    (B) MatSetValues(), using the global ordering
323*e7a95102SMartin Diehl
324*e7a95102SMartin Diehl!        - Set matrix entries using the global ordering by calling
325*e7a95102SMartin Diehl!          MatSetValues()
326*e7a95102SMartin Diehl!  Option (A) seems cleaner/easier in many cases, and is the procedure
327*e7a95102SMartin Diehl!  used in this example.
328*e7a95102SMartin Diehl!
329*e7a95102SMartin Diehl  subroutine FormJacobian(snes, X, jac, jac_prec, user, ierr)
330*e7a95102SMartin Diehl!  Input/output variables:
331*e7a95102SMartin Diehl    SNES snes
332*e7a95102SMartin Diehl    Vec X
333*e7a95102SMartin Diehl    Mat jac, jac_prec
334*e7a95102SMartin Diehl    type(userctx) user
335*e7a95102SMartin Diehl    PetscErrorCode ierr
336*e7a95102SMartin Diehl    DM da
337*e7a95102SMartin Diehl
338*e7a95102SMartin Diehl!  Declarations for use with local arrays:
339*e7a95102SMartin Diehl    PetscScalar, pointer :: lx_v(:)
340*e7a95102SMartin Diehl    Vec localX
341*e7a95102SMartin Diehl
342*e7a95102SMartin Diehl!  Scatter ghost points to local vector, using the 2-step process
343*e7a95102SMartin Diehl!     DMGlobalToLocalBegin(), DMGlobalToLocalEnd()
344*e7a95102SMartin Diehl!  Computations can be done while messages are in transition,
345*e7a95102SMartin Diehl!  by placing code between these two statements.
346*e7a95102SMartin Diehl
347*e7a95102SMartin Diehl    PetscCallA(SNESGetDM(snes, da, ierr))
348*e7a95102SMartin Diehl    PetscCallA(DMGetLocalVector(da, localX, ierr))
349*e7a95102SMartin Diehl    PetscCallA(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr))
350*e7a95102SMartin Diehl    PetscCallA(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr))
351*e7a95102SMartin Diehl
352*e7a95102SMartin Diehl!  Get a pointer to vector data
353*e7a95102SMartin Diehl    PetscCallA(VecGetArray(localX, lx_v, ierr))
354*e7a95102SMartin Diehl
355*e7a95102SMartin Diehl!  Compute entries for the locally owned part of the Jacobian preconditioner.
356*e7a95102SMartin Diehl    PetscCallA(FormJacobianLocal(lx_v, jac_prec, user, ierr))
357*e7a95102SMartin Diehl
358*e7a95102SMartin Diehl!  Assemble matrix, using the 2-step process:
359*e7a95102SMartin Diehl!     MatAssemblyBegin(), MatAssemblyEnd()
360*e7a95102SMartin Diehl!  Computations can be done while messages are in transition,
361*e7a95102SMartin Diehl!  by placing code between these two statements.
362*e7a95102SMartin Diehl
363*e7a95102SMartin Diehl    PetscCallA(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
364*e7a95102SMartin Diehl    if (jac /= jac_prec) then
365*e7a95102SMartin Diehl      PetscCallA(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
366*e7a95102SMartin Diehl    end if
367*e7a95102SMartin Diehl    PetscCallA(VecRestoreArray(localX, lx_v, ierr))
368*e7a95102SMartin Diehl    PetscCallA(DMRestoreLocalVector(da, localX, ierr))
369*e7a95102SMartin Diehl    PetscCallA(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
370*e7a95102SMartin Diehl    if (jac /= jac_prec) then
371*e7a95102SMartin Diehl      PetscCallA(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr))
372*e7a95102SMartin Diehl    end if
373*e7a95102SMartin Diehl
374*e7a95102SMartin Diehl!  Tell the matrix we will never add a new nonzero location to the
375*e7a95102SMartin Diehl!  matrix. If we do it will generate an error.
376*e7a95102SMartin Diehl
377*e7a95102SMartin Diehl    PetscCallA(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE, ierr))
378*e7a95102SMartin Diehl
379*e7a95102SMartin Diehl  end
380*e7a95102SMartin Diehl
381*e7a95102SMartin Diehl! ---------------------------------------------------------------------
382*e7a95102SMartin Diehl!
383*e7a95102SMartin Diehl!  FormJacobianLocal - Computes Jacobian matrix used to compute the preconditioner,
384*e7a95102SMartin Diehl!  called by the higher level routine FormJacobian().
385*e7a95102SMartin Diehl!
386*e7a95102SMartin Diehl!  Input Parameters:
387*e7a95102SMartin Diehl!  x        - local vector data
388*e7a95102SMartin Diehl!
389*e7a95102SMartin Diehl!  Output Parameters:
390*e7a95102SMartin Diehl!  jac_prec - Jacobian matrix used to compute the preconditioner
391*e7a95102SMartin Diehl!  ierr     - error code
392*e7a95102SMartin Diehl!
393*e7a95102SMartin Diehl!  Notes:
394*e7a95102SMartin Diehl!  This routine uses standard Fortran-style computations over a 2-dim array.
395*e7a95102SMartin Diehl!
396*e7a95102SMartin Diehl!  Notes:
397*e7a95102SMartin Diehl!  Due to grid point reordering with DMDAs, we must always work
398*e7a95102SMartin Diehl!  with the local grid points, and then transform them to the new
399*e7a95102SMartin Diehl!  global numbering with the "ltog" mapping
400*e7a95102SMartin Diehl!  We cannot work directly with the global numbers for the original
401*e7a95102SMartin Diehl!  uniprocessor grid!
402*e7a95102SMartin Diehl!
403*e7a95102SMartin Diehl!  Two methods are available for imposing this transformation
404*e7a95102SMartin Diehl!  when setting matrix entries:
405*e7a95102SMartin Diehl!    (A) MatSetValuesLocal(), using the local ordering (including
406*e7a95102SMartin Diehl!        ghost points!)
407*e7a95102SMartin Diehl!        - Set matrix entries using the local ordering
408*e7a95102SMartin Diehl!          by calling MatSetValuesLocal()
409*e7a95102SMartin Diehl!    (B) MatSetValues(), using the global ordering
410*e7a95102SMartin Diehl!        - Then apply this map explicitly yourself
411*e7a95102SMartin Diehl!        - Set matrix entries using the global ordering by calling
412*e7a95102SMartin Diehl!          MatSetValues()
413*e7a95102SMartin Diehl!  Option (A) seems cleaner/easier in many cases, and is the procedure
414*e7a95102SMartin Diehl!  used in this example.
415*e7a95102SMartin Diehl!
416*e7a95102SMartin Diehl  subroutine FormJacobianLocal(x, jac_prec, user, ierr)
417*e7a95102SMartin Diehl!  Input/output variables:
418*e7a95102SMartin Diehl    type(userctx) user
419*e7a95102SMartin Diehl    PetscScalar x(user%gxs:user%gxe, user%gys:user%gye)
420*e7a95102SMartin Diehl    Mat jac_prec
421*e7a95102SMartin Diehl    PetscErrorCode ierr
422*e7a95102SMartin Diehl
423*e7a95102SMartin Diehl!  Local variables:
424*e7a95102SMartin Diehl    PetscInt row, col(5), i, j
425*e7a95102SMartin Diehl    PetscInt ione, ifive
426*e7a95102SMartin Diehl    PetscScalar two, one, hx, hy, hxdhy
427*e7a95102SMartin Diehl    PetscScalar hydhx, sc, v(5)
428*e7a95102SMartin Diehl
429*e7a95102SMartin Diehl!  Set parameters
430*e7a95102SMartin Diehl    ione = 1
431*e7a95102SMartin Diehl    ifive = 5
432*e7a95102SMartin Diehl    one = 1.0
433*e7a95102SMartin Diehl    two = 2.0
434*e7a95102SMartin Diehl    hx = one/(user%mx - 1)
435*e7a95102SMartin Diehl    hy = one/(user%my - 1)
436*e7a95102SMartin Diehl    sc = hx*hy
437*e7a95102SMartin Diehl    hxdhy = hx/hy
438*e7a95102SMartin Diehl    hydhx = hy/hx
439*e7a95102SMartin Diehl
440*e7a95102SMartin Diehl!  Compute entries for the locally owned part of the Jacobian.
441*e7a95102SMartin Diehl!   - Currently, all PETSc parallel matrix formats are partitioned by
442*e7a95102SMartin Diehl!     contiguous chunks of rows across the processors.
443*e7a95102SMartin Diehl!   - Each processor needs to insert only elements that it owns
444*e7a95102SMartin Diehl!     locally (but any non-local elements will be sent to the
445*e7a95102SMartin Diehl!     appropriate processor during matrix assembly).
446*e7a95102SMartin Diehl!   - Here, we set all entries for a particular row at once.
447*e7a95102SMartin Diehl!   - We can set matrix entries either using either
448*e7a95102SMartin Diehl!     MatSetValuesLocal() or MatSetValues(), as discussed above.
449*e7a95102SMartin Diehl!   - Note that MatSetValues() uses 0-based row and column numbers
450*e7a95102SMartin Diehl!     in Fortran as well as in C.
451*e7a95102SMartin Diehl
452*e7a95102SMartin Diehl    do j = user%ys, user%ye
453*e7a95102SMartin Diehl      row = (j - user%gys)*user%gxm + user%xs - user%gxs - 1
454*e7a95102SMartin Diehl      do i = user%xs, user%xe
455*e7a95102SMartin Diehl        row = row + 1
456*e7a95102SMartin Diehl!           boundary points
457*e7a95102SMartin Diehl        if (i == 1 .or. j == 1 .or. i == user%mx .or. j == user%my) then
458*e7a95102SMartin Diehl          col(1) = row
459*e7a95102SMartin Diehl          v(1) = one
460*e7a95102SMartin Diehl          PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ione, col, v, INSERT_VALUES, ierr))
461*e7a95102SMartin Diehl!           interior grid points
462*e7a95102SMartin Diehl        else
463*e7a95102SMartin Diehl          v(1) = -hxdhy
464*e7a95102SMartin Diehl          v(2) = -hydhx
465*e7a95102SMartin Diehl          v(3) = two*(hydhx + hxdhy) - sc*user%lambda*exp(x(i, j))
466*e7a95102SMartin Diehl          v(4) = -hydhx
467*e7a95102SMartin Diehl          v(5) = -hxdhy
468*e7a95102SMartin Diehl          col(1) = row - user%gxm
469*e7a95102SMartin Diehl          col(2) = row - 1
470*e7a95102SMartin Diehl          col(3) = row
471*e7a95102SMartin Diehl          col(4) = row + 1
472*e7a95102SMartin Diehl          col(5) = row + user%gxm
473*e7a95102SMartin Diehl          PetscCallA(MatSetValuesLocal(jac_prec, ione, [row], ifive, col, v, INSERT_VALUES, ierr))
474*e7a95102SMartin Diehl        end if
475*e7a95102SMartin Diehl      end do
476*e7a95102SMartin Diehl    end do
477*e7a95102SMartin Diehl
478*e7a95102SMartin Diehl  end
479*e7a95102SMartin Diehl
480dfbbaf82SBarry Smithend module ex5f90module
481c4762a1bSJed Brown
482c4762a1bSJed Brownprogram main
483dfbbaf82SBarry Smith  use ex5f90module
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