DOUG 0.2

geom.F90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00028 
00029 program main_geom
00030 
00031   use Distribution_mod
00032   use Partitioning_mod
00033   use Partitioning_full_mod
00034   use Mesh_class
00035   use SpMtx_mods
00036   use Vect_mod
00037   use DenseMtx_mod
00038   use solvers_mod
00039   use CoarseAllgathers
00040   use Preconditioner_mod
00041   use FinePreconditioner_complete_mod
00042   use CoarsePreconditioner_geometric_mod
00043 
00044   implicit none
00045 
00046 #include<doug_config.h>
00047 
00048 #ifdef D_COMPLEX
00049 #define float complex
00050 #else
00051 #define float real
00052 #endif
00053 
00054   float(kind=rk), dimension(:), pointer :: xl ! local solution vector
00055   float(kind=rk), dimension(:), pointer :: x  ! global solution on master
00056 
00057   ! Partitioning
00058   integer               :: nparts ! number of partitons to partition a mesh
00059   integer, dimension(6) :: part_opts = (/0,0,0,0,0,0/)
00060 
00061   type(ConvInf) :: resStat
00062 
00063   ! +
00064   real(kind=rk) :: t1, time
00065   integer :: i,it,ierr,ol
00066   real(kind=rk) :: res_norm_local,res_norm,cond_num
00067   real(kind=rk), dimension(:), pointer :: r, y
00068   float(kind=rk), dimension(:), pointer :: yc, gyc, ybuf
00069 
00070   type(Distribution) :: D !< mesh and matrix distribution
00071   type(Partitionings) :: P !< mesh partitionings
00072   type(FinePreconditioner) :: FP !< fine preconditioner
00073   type(CoarsePreconditioner) :: CP !< coarse level preconditioner
00074 
00075   !DEBUG
00076   integer :: k
00077   real(kind=xyzk) :: mi(3),ma(3)
00078 
00079   ! Init DOUG
00080   call DOUG_Init()
00081 
00082   time = MPI_WTime()
00083 
00084   t1 = MPI_WTime()
00085   ! Master participates in calculations as well
00086   nparts = numprocs
00087 
00088   D = Distribution_NewInit(sctls%input_type,nparts,part_opts)
00089 
00090   if(pstream/=0) write(pstream, "(I0,':distribute time:',F0.3)") myrank, MPI_WTIME()-t1
00091 
00092   ! create partitionings
00093   P = Partitionings_New()
00094   call Partitionings_full_InitCoarse(P,D)
00095 
00096   ! create subdomains
00097   if (sctls%overlap<0) then ! autom. overlap from smoothing
00098     ol = max(sctls%smoothers,0)
00099   else
00100     ol = sctls%overlap
00101   endif
00102   FP = FinePreconditioner_New(D)
00103   call FinePreconditioner_Init(FP, D, P, ol)
00104   call FinePreconditioner_complete_Init(FP)
00105 
00106   ! conversion from elemental form to assembled matrix wanted?
00107   if (mctls%dump_matrix_only.eqv..true.) then
00108      call SpMtx_writeMatrix(D%A)
00109      call DOUG_Finalize()
00110      stop
00111   end if
00112   ! Geometric coarse grid processing
00113   CP = CoarsePreconditioner_New()
00114   if (sctls%input_type==DISTRIBUTION_TYPE_ELEMENTAL .and. sctls%levels==2) then
00115     t1 = MPI_WTime()    
00116     call CoarsePreconditioner_geometric_Init(CP, D)
00117     if(pstream/=0) write(pstream, "(I0,':coarse time:',F0.3)") myrank, MPI_WTIME()-t1
00118   endif
00119 
00120   ! Solve the system
00121   allocate(xl(D%mesh%nlf)); xl = 0.0_rk
00122 
00123   select case(sctls%solver)
00124   case (DCTL_SOLVE_PCG)
00125      t1 = MPI_WTIME()
00126      call pcg_weigs(D,x=xl,finePrec=FP,coarsePrec=CP,it=it,cond_num=cond_num)
00127      write(stream,*) 'time spent in pcg():',MPI_WTIME()-t1
00128      if(pstream/=0) write(pstream, "(I0,':pcg time:',F0.3)") myrank, MPI_WTIME()-t1
00129 
00130   case default
00131      call DOUG_abort('[DOUG main] : Wrong solution method specified', -1)
00132   end select
00133 
00134   ! Calculate solution residual (in parallel)
00135   allocate(r(size(xl)), y(size(xl)))
00136   call Distribution_pmvm(D,y,xl)
00137   r = y - D%rhs
00138   res_norm_local = Vect_dot_product(r, r)
00139   deallocate(r, y)
00140 
00141   ! Calculate global residual
00142   call MPI_REDUCE(res_norm_local, res_norm, 1, MPI_fkind, MPI_SUM, D_MASTER, MPI_COMM_WORLD, ierr)
00143 
00144   ! Assemble result on master and write it to screen and/or file
00145   if (ismaster()) then
00146      allocate(x(D%mesh%ngf)); x = 0.0_rk
00147   end if
00148   call Vect_Gather(xl, x, D%mesh)
00149   if (ismaster().and.(size(x) <= 100).and.(D_MSGLVL > 0)) &
00150        call Vect_Print(x, 'solution ')
00151   if (ismaster()) then
00152        write(stream,*) 'dsqrt(res_norm) =',dsqrt(res_norm)
00153        call WriteSolutionToFile(x)
00154   endif
00155   if (ismaster()) then
00156      deallocate(x)
00157   end if
00158 
00159   ! Destroy objects
00160   call Mesh_Destroy(D%mesh)
00161   call SpMtx_Destroy(D%A)
00162 
00163   if (sctls%input_type==DISTRIBUTION_TYPE_ELEMENTAL .and. sctls%levels==2) then
00164     call SpMtx_Destroy(CP%AC)
00165     call SpMtx_Destroy(CP%R)
00166     call SendData_Destroy(CP%cdat%send)
00167   endif
00168 
00169   call ConvInf_Destroy(resStat)
00170   call Vect_cleanUp()
00171   deallocate(D%rhs, xl)
00172 
00173   if(pstream/=0) write(pstream, "(I0,':total time:',F0.3)") myrank, MPI_WTIME()-time
00174 
00175   call DOUG_Finalize()
00176 
00177 end program main_geom