求3D POSSION 方程解算子程序, 用于网格光滑及其它相关数据光滑

叶轮机

求3D POSSION 方程解算子程序,改编后用于数据光滑

wendu

傲雪上有3D LAPLACE的源程序

chbo76

不如转贴一下啦。

harry

是peric的程序的求解器那部分吧，

firedragon

!
!%%%%%%%%%%%%%%%%%%%%%%%%%%% PROGRAM POISSON3 %%%%%%%%%%%%%%%%%%%%%%%%%%
!
   PROGRAM poisson3
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
!  Solves Poisson's equation on a rectangular domain with Dirichlet
!  boundary conditions by three methods:
!
!        1) Successive Over-Relaxation (SOR) method
!        2) Multi-grid method
!        3) Fast elliptic method
!
!  Author:  Dave Pruett
!           JMU Dept. of Mathematics and Statistics
!           08/15/00
!
!  Notes:  See Golub and Ortega, 1992
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
      USE data_types
      USE right_hand_side
      USE successive_over_relaxation
      USE multi_grid_algorithm
      USE fast_poisson_solver
      USE visualization
      IMPLICIT NONE
      INTEGER :: log2m, m, log2n, n, i, j, iskip, jskip, it_max, method
      INTEGER :: count0, count1, count_rate, count_max
      REAL(PREC) :: time
      REAL(PREC) :: x_max, y_max, x, y, dx, dy, pi, sigma
      REAL(PREC) :: r_max, f_max, rel_err
      REAL(PREC), DIMENSION(:,, ALLOCATABLE :: p, del_p, f, r, exact
!
! only for MULTI_GRID
!
      INTEGER :: iterations, max_vcycles, ierr, idiag
      REAL(PREC) :: tolerance
!
!  only for DFASTP
!
      INTEGER :: iflag = 0
      REAL(PREC) :: a, b, c, d, coef
      REAL(PREC), DIMENSION(:,, ALLOCATABLE :: ra1, ra2
      COMPLEX(PREC), DIMENSION(:,, ALLOCATABLE :: ca1, ca2
      REAL(PREC), DIMENSION(:), ALLOCATABLE :: rv
      COMPLEX(PREC), DIMENSION(:), ALLOCATABLE :: cv
!
!  initialize and allocate storage
!
      WRITE (6,*) ' enter 1=SOR, 2=MULTI-GRID, 3=FAST DIRECT'
      READ (5,*) method
      WRITE (6,*) ' enter log2m, log2n'
      READ (5,*) log2m, log2n
      m = 2**log2m
      n = 2**log2n
      ALLOCATE (p(0:m,0:n), del_p(0:m,0:n), f(0:m,0:n), r(0:m,0:n))
      ALLOCATE (exact(0:m,0:n))
!
!  form the right-hand-side and exact solution for comparison
!
      pi = ACOS(-1.0)
      x_max = pi / 4.d0
      y_max = 1.0d0
      dx = x_max / m
      dy = y_max / n
      DO i = 0,m
         x = i*dx
         DO j = 0,n
            y = j*dy
            exact(i,j) = RHS(x,y)
            f(i,j) = -3.d0*exact(i,j)
         END DO
      END DO
!
!  load the boundary conditions
!
      p = 0.0d0
      p(0:m,0) = exact(0:m,0)         
      p(0:m,n) = exact(0:m,n)         
      p(0,0:n) = exact(0,0:n)         
      p(m,0:n) = exact(m,0:n)         
!
!  solve by the SOR method
!
      IF (method == 1) THEN
         WRITE (6,*) '  '
         WRITE (6,*) ' solution by SOR '
         WRITE (6,*) ' enter maximum number of iterations, it_max '
         READ (5,*) it_max
         WRITE (6,*) ' enter relative error tolerance, rel_err '
         READ (5,*) rel_err
         WRITE (6,*) ' enter relaxation parameter, 0 < sigma < 2 '
         READ (5,*) sigma
         WRITE (6,*) ' enter 2 for movie+history; 1 history only; 0 otherwise '
         READ (5,*) idiag
         CALL SYSTEM_CLOCK (count0, count_rate, count_max)
         f = dx * dx * f     ! scale the RHS
         CALL sor (m, m, n, it_max, dx, dy, sigma, p, f, del_p, r, &
                   rel_err, idiag)
         CALL SYSTEM_CLOCK (count1, count_rate, count_max)
         time = REAL(count1 - count0) / count_rate
         WRITE (6,*) ' counts = ', count1 - count0
         WRITE (6,*) ' processor time for solution by SOR = ', time
!
!  solve by the MULTI-GRID method
!
      ELSE IF (method == 2) THEN
         WRITE (6,*) '  '
         WRITE (6,*) ' solution by MULTI-GRID'
         WRITE (6,*) ' enter maximum number of v-cycles, max_vcycles'
         READ (5,*) max_vcycles
         WRITE (6,*) ' enter iterations per v-cycle, iterations '
         READ (5,*) iterations
         WRITE (6,*) ' enter convergence tolerance, tolerance'
         READ (5,*) tolerance
         WRITE (6,*) ' enter 2 for movie+history; 1 history only; 0 otherwise '
         READ (5,*) idiag
         CALL SYSTEM_CLOCK (count0, count_rate, count_max)
         CALL multi_grid (log2m, log2n, dx, dy, p, f, iterations,  &
                          max_vcycles, tolerance, ierr, idiag)
         CALL SYSTEM_CLOCK (count1, count_rate, count_max)
         time = REAL(count1 - count0) / count_rate
         WRITE (6,*) ' counts = ', count1 - count0
         WRITE (6,*) ' processor time for solution by MULTI-GRID = ', time
!
!  solve by the FAST DIRECT method
!
      ELSE IF (method == 3) THEN
         WRITE (6,*) '  '
         WRITE (6,*) ' solution by FAST DIRECT '
         ALLOCATE (ra1(0:m,0:n), ra2(0:m,0:n))
         ALLOCATE (ca1(0:m,n/2), ca2(0:m,n/2))
         ALLOCATE (rv(0:m/2), cv(0:m/2))
         idiag = 0
         a = 0.0d0
         b = y_max
         c = 0.0d0
         d = x_max
         coef = 1.0
         CALL SYSTEM_CLOCK (count0, count_rate, count_max)
         CALL DFASTP (m+1,n,m,a,b,c,d,coef,p,f,ca1,ca2,cv, &
                      ra1,ra2,rv,iflag)                              
         WRITE (6,*) ' iflag = ', iflag
         IF (IFLAG /= n) THEN
            STOP ' DFASTP'
         END IF
         CALL SYSTEM_CLOCK (count1, count_rate, count_max)
         time = REAL(count1 - count0) / count_rate
         WRITE (6,*) ' counts = ', count1 - count0
         WRITE (6,*) ' processor time for solution by FAST DIRECT = ', time
         DEALLOCATE (ra1, ra2, ca1, ca2, rv, cv)
      ELSE
         WRITE (6,*) ' invalid choice for method: STOP'
         STOP ' method'
      END IF
!
!  output the solution
!
      iskip = m / 4
      jskip = n / 4
      WRITE (6,*) ' absolute error in solution'
      DO j = n,0,-jskip
         WRITE (6,'(6(1x,g12.6))') ((exact(i,j) - p(i,j)), &
                                     i=0,m,iskip)
      END DO
      IF (idiag > 1) THEN
         OPEN(UNIT=1,FILE='exact',STATUS='UNKNOWN')
         DO j = 0,n
            DO i = 0,m
               WRITE (1,'(1x,g13.7)') exact(i,j)
            END DO
         END DO
         CLOSE(1)
      END IF
      DEALLOCATE (p, del_p, f, r, exact)
      STOP
   END PROGRAM poisson3