1c* ///////////////////////////////////////////////////////////////////////////
2c* @file    pded.f
3c* @author  Michael Holst
4c* @brief   PDE definition file.
5c* @version $Id: pdef.f,v 1.2 2010/08/12 05:53:09 fetk Exp $
6c* @attention
7c* @verbatim
8c*
9c* PMG -- Parallel algebraic MultiGrid
10c* Copyright (C) 1994-- Michael Holst.
11c*
12c* Michael Holst <mholst@math.ucsd.edu>
13c* University of California, San Diego
14c* Department of Mathematics, 5739 AP&M
15c* 9500 Gilman Drive, Dept. 0112
16c* La Jolla, CA 92093-0112 USA
17c* http://math.ucsd.edu/~mholst
18c*
19c* This file is part of PMG.
20c*
21c* PMG is free software; you can redistribute it and/or modify
22c* it under the terms of the GNU General Public License as published by
23c* the Free Software Foundation; either version 2 of the License, or
24c* (at your option) any later version.
25c*
26c* PMG is distributed in the hope that it will be useful,
27c* but WITHOUT ANY WARRANTY; without even the implied warranty of
28c* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
29c* GNU General Public License for more details.
30c*
31c* You should have received a copy of the GNU General Public License
32c* along with PMG; if not, write to the Free Software
33c* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
34c*
35c* Linking PMG statically or dynamically with other modules is making a
36c* combined work based on PMG. Thus, the terms and conditions of the GNU
37c* General Public License cover the whole combination.
38c*
39c* SPECIAL GPL EXCEPTION
40c* In addition, as a special exception, the copyright holders of PMG
41c* give you permission to combine the PMG program with free software
42c* programs and libraries that are released under the GNU LGPL or with
43c* code included in releases of ISIM, PMV, PyMOL, SMOL, VMD, and Vision.
44c* Such combined software may be linked with PMG and redistributed together
45c* in original or modified form as mere aggregation without requirement that
46c* the entire work be under the scope of the GNU General Public License.
47c* This special exception permission is also extended to any software listed
48c* in the SPECIAL GPL EXCEPTION clauses by the FEtk and APBS libraries.
49c*
50c* Note that people who make modified versions of PMG are not obligated
51c* to grant this special exception for their modified versions; it is
52c* their choice whether to do so. The GNU General Public License gives
53c* permission to release a modified version without this exception; this
54c* exception also makes it possible to release a modified version which
55c* carries forward this exception.
56c*
57c* @endverbatim
58c* ///////////////////////////////////////////////////////////////////////////
59
60c* ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
61c*
62c* DIFFERENTIAL EQUATION: Poisson's
63c*
64c* BOUNDARY CONDITIONS:
65c*
66c*     East   Face (xmin):  Dirichlet, homogeneous
67c*     West   Face (xmax):  Dirichlet, homogeneous
68c*     North  Face (ymin):  Dirichlet, homogeneous
69c*     South  Face (ymax):  Dirichlet, homogeneous
70c*     Top    Face (zmin):  Dirichlet, homogeneous
71c*     Bottom Face (zmax):  Dirichlet, homogeneous
72c*
73c* MESH:                  hx  = (xmax-xmin) / (nx-1)
74c*                        hy  = (ymax-ymin) / (ny-1)
75c*                        hz  = (zmax-zmin) / (nz-1)
76c*                        xi = xmin + (i-1) * hx,  i=1,...,nx
77c*                        yi = ymin + (j-1) * hy,  j=1,...,ny
78c*                        zi = zmin + (k-1) * hk,  k=1,...,nz
79c*
80c* CHOSEN TRUE SOLUTION:  u(x,y,z) = Sin(n pi x)*Sin(m pi y)*Sin(l pi z)
81c*
82c* author:  michael holst
83c* ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
84c*
85      function c_scal(coef,u,ipkey)
86c* *********************************************************************
87c* purpose:
88c*
89c*    define the nonlinearity (scalar version)
90c*
91c* author:  michael holst
92c* *********************************************************************
93      implicit         none
94      double precision coef,u,c_scal
95      integer          ipkey
96c*
97c*    *** create the nonlinear term ***
98      c_scal = coef * u
99c*
100c*    *** end it ***
101      return
102      end
103      function dc_scal(coef,u,ipkey)
104c* *********************************************************************
105c* purpose:
106c*
107c*    define the derivative of the nonlinearity (scalar version)
108c*
109c* author:  michael holst
110c* *********************************************************************
111      implicit         none
112      double precision coef,u,dc_scal
113      integer          ipkey
114c*
115c*    *** create the nonlinear term ***
116      dc_scal = coef
117c*
118c*    *** end it ***
119      return
120      end
121      subroutine c_vec(coef,uin,uout,nx,ny,nz,ipkey)
122c* *********************************************************************
123c* purpose:
124c*
125c*    define the nonlinearity (vector version)
126c*
127c* author:  michael holst
128c* *********************************************************************
129      implicit         none
130      integer          nx,ny,nz,ipkey
131      double precision uin(*),uout(*),coef(*)
132      integer          n,i,ii,nproc,ipara,ivect
133      parameter        (nproc=1)
134c*
135cmdir 0 0
136c*
137c*    *** find parallel loops (ipara), remainder (ivect) ***
138      n     = nx * ny * nz
139      ipara = n / nproc
140      ivect = mod(n,nproc)
141c*
142c*    *** do parallel loops ***
143cmdir 2 1
144      do 10 ii = 1, nproc
145cmdir 2 2
146         do 11 i = 1+(ipara*(ii-1)), ipara*ii
147            uout(i) = coef(i) * uin(i)
148 11      continue
149 10   continue
150c*
151c*    *** do vector loops ***
152cmdir 1 1
153      do 20 i = ipara*nproc+1, n
154         uout(i) = coef(i) * uin(i)
155 20   continue
156c*
157c*    *** end it ***
158      return
159      end
160      subroutine dc_vec(coef,uin,uout,nx,ny,nz,ipkey)
161c* *********************************************************************
162c* purpose:
163c*
164c*    define the derivative of the nonlinearity (vector version)
165c*
166c* author:  michael holst
167c* *********************************************************************
168      implicit         none
169      integer          nx,ny,nz,ipkey
170      double precision uin(*),uout(*),coef(*)
171      integer          n,i,ii,nproc,ipara,ivect
172      parameter        (nproc=1)
173c*
174cmdir 0 0
175c*
176c*    *** find parallel loops (ipara), remainder (ivect) ***
177      n     = nx * ny * nz
178      ipara = n / nproc
179      ivect = mod(n,nproc)
180c*
181c*    *** do parallel loops ***
182cmdir 2 1
183      do 10 ii = 1, nproc
184cmdir 2 2
185         do 11 i = 1+(ipara*(ii-1)), ipara*ii
186            uout(i) = coef(i)
187 11      continue
188 10   continue
189c*
190c*    *** do vector loops ***
191cmdir 1 1
192      do 20 i = ipara*nproc+1, n
193         uout(i) = coef(i)
194 20   continue
195c*
196c*    *** end it ***
197      return
198      end
199      subroutine fillco(iparm,rparm,nx,ny,nz,
200     2   xf,yf,zf,gxcf,gycf,gzcf,a1cf,a2cf,a3cf,ccf,fcf,tcf)
201c* *********************************************************************
202c* purpose:
203c*
204c*   this subroutine defines poisson's equation on the unit cube.
205c*   boundary conditions are zero dirichlet.
206c*
207c* details:
208c*
209c*   this routine sets up the coefficients of the following
210c*   three-dimensional, 2nd order linear elliptic partial
211c*   differential equation:
212c*
213c*      lu = f, u in omega
214c*       u = g, u on boundary of omega
215c*
216c*   where
217c*      omega = [xmin,xmax]x[ymin,ymax]x[zmin,zmax]
218c*
219c*   and l is taken to be in the following form:
220c*
221c*      lu = - \nabla \cdot (a \nabla u) + c u
222c*
223c*   the coefficients must be supplied as follows:
224c*
225c*      a1cf(,,,):  at x-half grid points, hence array is:  (nx-1)*ny*nz
226c*      a2cf(,,,):  at y-half grid points, hence array is:  nx*(ny-1)*nz
227c*      a3cf(,,,):  at z-half grid points, hence array is:  nx*ny*(nz-1)
228c*                  (we index all three as:  nx*ny*nz however)
229c*
230c*      ccf:        at grid points, hence array is:  nx*ny*nz
231c*      fcf:        at grid points, hence array is:  nx*ny*nz
232c*      u:          must be set to have appropriate boundary values
233c*
234c* author:  michael holst
235c* *********************************************************************
236      implicit         none
237      integer          iparm(*),nx,ny,nz
238      double precision rparm(*)
239      double precision a1cf(nx,ny,nz),a2cf(nx,ny,nz),a3cf(nx,ny,nz)
240      double precision ccf(nx,ny,nz),fcf(nx,ny,nz),tcf(nx,ny,nz)
241      double precision xf(nx),yf(ny),zf(nz)
242      double precision gxcf(ny,nz,*),gycf(nx,nz,*),gzcf(nx,ny,*)
243c*
244c*    *** variable declarations ***
245      integer          i,j,k,iinfo
246      double precision hx,hy,hz,xmin,xmax,ymin,ymax,zmin,zmax
247      double precision hxo2,hyo2,hzo2
248c*
249c*    *** some parameters ***
250      double precision zn,zm,zl,pi
251      parameter        (zn=1.0d0,zm=1.0d0,zl=1.0d0)
252      pi               = 4.0d0 * datan(1.0d0)
253c*
254cmdir 0 0
255c*
256c*    *** info messages ***
257      iinfo = iparm(12)
258c*
259c* *********************************************************************
260c* *** definition of the domain
261c* *********************************************************************
262c*
263c*    *** set the correct boundary ranges ***
264      xmin = 0.0d0
265      xmax = 1.0d0
266      ymin = 0.0d0
267      ymax = 1.0d0
268      zmin = 0.0d0
269      zmax = 1.0d0
270c*
271c*    *** store in rparm array ***
272      rparm(3)  = xmin
273      rparm(4)  = xmax
274      rparm(5)  = ymin
275      rparm(6)  = ymax
276      rparm(7)  = zmin
277      rparm(8)  = zmax
278c*
279c* *********************************************************************
280c* *** boundary value definitions
281c* *********************************************************************
282c*
283c*    *** determine mesh widths ***
284      hx  = (xmax-xmin) / dble(nx-1)
285      hy  = (ymax-ymin) / dble(ny-1)
286      hz  = (zmax-zmin) / dble(nz-1)
287      hxo2= hx / 2.0d0
288      hyo2= hy / 2.0d0
289      hzo2= hz / 2.0d0
290c*
291c* *********************************************************************
292c* *** pde coefficient value definitions
293c* *********************************************************************
294c*
295c*    *** fill coefficient arrays ***
296cmdir 3 1
297      do 10 k = 1, nz
298         zf(k) = zmin + dble(k-1) * hz
299cmdir 3 2
300         do 11 j = 1, ny
301            yf(j) = ymin + dble(j-1) * hy
302cmdir 3 3
303            do 12 i = 1, nx
304               xf(i) = xmin + dble(i-1) * hx
305c*
306c*             *** the diagonal tensor (2nd derivative) entries ***
307               a1cf(i,j,k) = (1.0d0)
308               a2cf(i,j,k) = (1.0d0)
309               a3cf(i,j,k) = (1.0d0)
310c*
311c*             *** the scalar (0th derivative) entry ***
312               ccf(i,j,k) = (0.0d0)
313c*
314c*             *** the rhs entry ***
315CZZZ           fcf(i,j,k) = ((zn*zn+zm*zm+zl*zl)*pi*pi
316CZZZ 2         *dsin(zn*pi*xf(i))*dsin(zm*pi*yf(j))*dsin(zl*pi*zf(k)))
317               fcf(i,j,k) = 1.0
318c*
319c*             *** the chosen true solution ***
320               tcf(i,j,k) =
321     2         (dsin(zn*pi*xf(i))*dsin(zm*pi*yf(j))*dsin(zl*pi*zf(k)))
322 12         continue
323 11      continue
324 10   continue
325c*
326c*    *** the (i=1) and (i=nx) boundaries (dirichlet) ***
327cmdir 2 1
328      do 50 k = 1, nz
329cmdir 2 2
330         do 51 j = 1, ny
331            gxcf(j,k,1) = tcf(1 ,j,k)
332            gxcf(j,k,2) = tcf(nx,j,k)
333            gxcf(j,k,3) = 0.0d0
334            gxcf(j,k,4) = 0.0d0
335 51      continue
336 50   continue
337c*
338c*    *** the (j=1) and (j=ny) boundaries (dirichlet) ***
339cmdir 2 1
340      do 60 k = 1, nz
341cmdir 2 2
342         do 61 i = 1, nx
343            gycf(i,k,1) = tcf(i,1 ,k)
344            gycf(i,k,2) = tcf(i,ny,k)
345            gycf(i,k,3) = 0.0d0
346            gycf(i,k,4) = 0.0d0
347 61      continue
348 60   continue
349c*
350c*    *** the (k=1) and (k=nz) boundaries (dirichlet) ***
351cmdir 2 1
352      do 70 j = 1, ny
353cmdir 2 2
354         do 71 i = 1, nx
355            gzcf(i,j,1) = tcf(i,j,1 )
356            gzcf(i,j,2) = tcf(i,j,nz)
357            gzcf(i,j,3) = 0.0d0
358            gzcf(i,j,4) = 0.0d0
359 71      continue
360 70   continue
361c*
362c*    *** end it ***
363      return
364      end
365
366