source: src/main/main_drivers.F90 @ 848affb

Revision 848affb, 10.8 KB checked in by oleg <oleg@…>, 5 years ago (diff)

Calculation of robust coarse spaces seem to be working for single
process.

git-svn-id:  svn://dougdevel.org/doug/trunk@235 146eac83-2b0f-0410-9363-9ebcc3d868b9
  • Property mode set to 100644
Line 
1! DOUG - Domain decomposition On Unstructured Grids
2! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
3! Department of Mathematics, University of Bath
4!
5! This library is free software; you can redistribute it and/or
6! modify it under the terms of the GNU Lesser General Public
7! License as published by the Free Software Foundation; either
8! version 2.1 of the License, or (at your option) any later version.
9!
10! This library is distributed in the hope that it will be useful,
11! but WITHOUT ANY WARRANTY; without even the implied warranty of
12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13! Lesser General Public License for more details.
14!
15! You should have received a copy of the GNU Lesser General Public
16! License along with this library; if not, write to the Free Software
17! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
18! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
19! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
20! mailto:info(at)dougdevel.org)
21
22module main_drivers
23
24  use doug_utils
25  use Graph_class
26  use Mesh_class
27  use ElemMtxs_mods
28  use SpMtx_mods
29  use Vect_mod
30  use DenseMtx_mod
31
32  implicit none
33
34#include<doug_config.h>
35
36#ifdef D_COMPLEX
37#define float complex
38#else
39#define float real
40#endif
41
42  public :: &
43       parallelDistributeInput, &
44       parallelAssembleFromElemInput ! TODO: make private, changing test_SpMtx_symmetry_at_pmvm first
45       
46contains
47
48  !----------------------------------------------------------------
49  !> Distributes data, chooses algorithm based on input type
50  !----------------------------------------------------------------
51  subroutine parallelDistributeInput(input_type, M, A, b, nparts, part_opts, A_interf)
52    implicit none
53
54    integer,        intent(in)     :: input_type !< Input Type
55    type(Mesh),     intent(in out) :: M !< Mesh
56    type(SpMtx),    intent(out) :: A !< System matrix
57    float(kind=rk), dimension(:), pointer :: b !< local RHS
58    ! Partitioning
59    integer, intent(in) :: nparts !< number of parts to partition a mesh
60    integer, dimension(6), intent(in) :: part_opts !< partition options (see METIS manual)
61    type(SpMtx),intent(in out),optional :: A_interf !< matrix at interface
62
63    select case (input_type)
64    case (DCTL_INPUT_TYPE_ELEMENTAL)
65       ! ELEMENTAL
66       call parallelAssembleFromElemInput(M,A,b,nparts,part_opts,A_interf)
67    case (DCTL_INPUT_TYPE_ASSEMBLED)
68       ! ASSEMBLED
69       call parallelDistributeAssembledInput(M,A,b,A_interf)
70    case default
71       call DOUG_abort('[DOUG main] : Unrecognised input type.', -1)
72    end select
73  end subroutine parallelDistributeInput
74
75  !----------------------------------------------------------------
76  !> Parallel assemble of system matrix and RHS from elemental input
77  !----------------------------------------------------------------
78  subroutine parallelAssembleFromElemInput(Msh, A, &
79               b, nparts, part_opts, A_interf)
80    implicit none
81
82    type(Mesh),     intent(in out) :: Msh !< Mesh
83    type(SpMtx),    intent(out) :: A !< System matrix
84    float(kind=rk), dimension(:), pointer :: b !< local RHS
85    ! Partitioning
86    integer, intent(in) :: nparts !< number of parts to partition a mesh
87    integer, dimension(6), intent(in) :: part_opts !< partition options (see METIS manual)
88    type(SpMtx),intent(in out),optional :: A_interf !< matrix at interface
89
90    ! =======================
91    ! Mesh and its Dual Graph
92    ! =======================
93    !
94    ! Create Mesh object
95    Msh = Mesh_New()
96
97    if (ismaster()) then ! MASTER
98       write(stream,*)
99       write(stream,*) 'master thread'
100
101       if (D_MSGLVL > 1) &
102            call MasterCtrlData_print()
103
104       ! Initialise Mesh object
105       call Mesh_initFromFile(Msh, trim(mctls%info_file))
106
107    else ! SLAVES
108       write(stream,'(a,i4,a)') 'slave [',myrank,'] thread'
109
110       if (D_MSGLVL > 1) &
111            call SharedCtrlData_print()
112    end if
113
114    ! Get from master Mesh's parameters: nell, ngf, mfrelt, nsd, nnode
115    call Mesh_paramsMPIBCAST(Msh)
116    if (D_MSGLVL > 1) &
117         call Mesh_printInfo(Msh)
118
119    ! Allocate data arrays (nfrelt, mhead, freemap, eptnmap) for mesh
120    call Mesh_allocate(Msh, &
121         nfrelt  =.true.,   &
122         mhead   =.true.,   &
123         freemap =.true.,   &
124         eptnmap  =.true.)
125
126    ! Master reads in from files: feedom lists, coordinates, freedom map
127    if (ismaster()) then
128       call Mesh_readFromFile(Msh, &
129            fnFreelists = trim(mctls%freedom_lists_file), &
130            fnCoords    = trim(mctls%coords_file),        &
131            fnFreemap   = trim(mctls%freemap_file))
132    end if
133
134    ! Master non-blockingly sends mesh data: nfrelt, mhead
135    call Mesh_dataMPIISENDRECV(Msh, &
136         nfrelt   = .true., &
137         mhead    = .true.)
138    if (D_MSGLVL > 4) &
139         call Mesh_printElemFree(Msh)
140
141    ! For multi-variable problems which have more than one block
142    if (sctls%number_of_blocks > 1) then
143       call Mesh_allocate(Msh, freemask=.true.)
144       if (ismaster()) then
145          call Mesh_readFromFile(Msh, &
146               fnFreemask = trim(mctls%freedom_mask_file))
147       end if
148       call Mesh_dataMPIISENDRECV(Msh, freemask = .true.)
149    end if
150
151
152    ! Build dual graph (Graph object is a data field in Mesh class)
153    ! let all procs do this - compare whith broadcasting the one bult on master
154    call Mesh_buildGraphDual(Msh)
155
156    ! Partition mesh's dual graph
157    if (ismaster()) then
158
159       !! Build dual graph (Graph object is a data field in Mesh class)
160       !call Mesh_buildGraphDual(Msh)
161
162       if (sctls%plotting == D_PLOT_YES) then
163
164          call Mesh_pl2D_pointCloud(Msh,D_PLPLOT_INIT)
165          ! Plots mesh's dual graph
166          call Mesh_pl2D_plotGraphDual(Msh,D_PLPLOT_END)
167          ! Mesh & its Dual Graph
168          call Mesh_pl2D_plotMesh(Msh, D_PLPLOT_INIT)
169          call Mesh_pl2D_plotGraphDual(Msh, D_PLPLOT_END)
170       end if
171
172       ! Partition graph: D_PART_PMETIS, D_PART_KMETIS, D_PART_VKMETIS
173       call Mesh_partitionDual(Msh, nparts, D_PART_VKMETIS, part_opts)
174
175       if (sctls%plotting == D_PLOT_YES) then
176          ! Draw colored partitoined graph
177          call Mesh_pl2D_plotGraphParted(Msh)
178
179          ! Plot partitions of the mesh
180          ! NB: Check for multivariable case! TODO
181          call Mesh_pl2D_Partition(Msh)
182          ! Partition with Dual Graph upon it
183          call Mesh_pl2D_Partition(Msh, D_PLPLOT_INIT)
184          call Mesh_pl2D_plotGraphDual(Msh, D_PLPLOT_CONT)
185          call Mesh_pl2D_pointCloud(Msh,D_PLPLOT_END)
186       end if
187
188       ! Destroy previously created graph (purely to save memory)
189       ! (If it is not killed here or somewere else Mesh_Destroy()
190       !  will kill it any way)
191       !call Mesh_destroyGraph(Msh)
192    else ! SLAVES
193       ! Number of partions the mesh was partitioned into
194       Msh%nparts = nparts
195       Msh%parted = .true.
196    end if
197
198    ! Distribute elements to partitons map among slaves
199    call Mesh_dataMPIISENDRECV(Msh, eptnmap=.true.)
200
201    ! Calculate number of elements in partitions
202    call Mesh_calcElemsInParts(Msh)
203
204    ! Build global to local, local to global maps and
205    ! inner/interface mask for freedoms
206    ! (also finds local number of freedoms 'Mesh%nlf';
207    !  this hidden appearance of 'Mesh_findNLF()' helps
208    !  to speed up calculations a bit)
209    call Mesh_buldMapsNghbrsMasksNLF(Msh)
210
211    ! ===============================
212    ! Assemble system matrix and RHS
213    ! ===============================
214    A = SpMtx_New()
215    allocate(b(Msh%nlf))
216    b = 0.0_rk
217
218    if (ismaster()) then
219       if (numprocs>1.and.present(A_interf)) then
220          A_interf = SpMtx_New()
221          call ElemMtxs_readAndDistribute(Msh, trim(mctls%elemmat_rhs_file), A, b, A_interf)
222       else
223          call ElemMtxs_readAndDistribute(Msh, trim(mctls%elemmat_rhs_file), A, b)
224       end if
225    else
226       if (numprocs>1.and.present(A_interf)) then
227          A_interf = SpMtx_New()
228          call ElemMtxs_recvAndAssemble(Msh, A, b, A_interf)
229       else
230          call ElemMtxs_recvAndAssemble(Msh, A, b)
231       end if
232    end if
233
234    if (sctls%verbose>9) then
235      write(stream,'(/a)') 'System matrix:'
236      call SpMtx_printInfo(A)
237      if (A%nrows <= 25) then
238         call SpMtx_printMat(A)
239      else if ((A%nrows > 25).and.(A%nrows <= 100)) then
240         call SpMtx_printRaw(A)
241      end if
242    endif
243
244    ! ==================
245    ! Finish assemble local RHS
246    ! ==================
247    if (A%nrows <= 25) & ! if (D_MSGLVL > 2) &
248         call Vect_Print(b, 'RHS assembled (local) ')
249
250    ! initialise auxiliary data for manipulating vectors
251    call Vect_setIntfEnd(sum(Msh%inner_interf_fmask))
252    call Vect_buildDotMask(Msh)
253    if (D_MSGLVL > 4) &
254         call Vect_Print(dot_intf_fmask,'dot_intf_fmask ')
255    call Vect_buildDotMap()
256    if (D_MSGLVL > 4) &
257         call Vect_Print(dot_intf_fmap,'dot_intf_fmap ')
258
259    ! Free mesh graph, not needed anymore
260    call Mesh_destroyGraph(Msh)
261
262    ! Update inner node count
263    Msh%ninner=Msh%nlf
264
265  end subroutine parallelAssembleFromElemInput
266
267
268  !----------------------------------------------------------------
269  !> Distribute assembled matrix and RHS from master to slaves
270  !----------------------------------------------------------------
271  subroutine parallelDistributeAssembledInput(Msh, A, b, A_interf)
272    implicit none
273
274    type(Mesh),     intent(in out) :: Msh !< Mesh
275    type(SpMtx),    intent(in out) :: A !< System matrix
276    float(kind=rk), dimension(:), pointer :: b !< local RHS
277    type(SpMtx),intent(in out),optional :: A_interf !< matrix at interface
278
279    integer :: n
280
281    ! ======================
282    ! Read matrix from file
283    ! ======================
284    if (ismaster()) then
285      write(stream,'(a,a)') ' ##### Assembled input file: ##### ', &
286            mctls%assembled_mtx_file
287      call ReadInSparseAssembled(A,trim(mctls%assembled_mtx_file))
288      allocate(b(A%nrows))
289      if (len_trim(mctls%assembled_rhs_file)>0) then
290        write(stream,'(a,a)') ' ##### Assembled RHS file: ##### ', &
291              mctls%assembled_rhs_file
292        call Vect_ReadFromFile(b, trim(mctls%assembled_rhs_file), mctls%assembled_rhs_format)
293      else
294        b=1.0_rk
295      end if
296    endif
297
298    ! =====================
299    ! Build mesh structure/distribute
300    ! =====================
301    if (numprocs==1) then
302      n=sqrt(1.0_rk*A%nrows)
303      if (n*n /= A%nrows) then
304        write (stream,*) 'Not a Cartesian Mesh!!!'
305        Msh=Mesh_New()
306        Msh%ngf=A%nrows
307        Msh%nlf=A%nrows
308        Msh%ninner=Msh%ngf
309      else
310        write (stream,*) 'Cartesian Mesh!!!'
311        call Mesh_BuildSquare(Msh,n)
312        Msh%ninner=Msh%ngf
313      endif
314    else ! numprocs>1
315      Msh=Mesh_New()
316      call SpMtx_DistributeAssembled(A,b,A_interf,Msh)
317    endif
318
319  end subroutine parallelDistributeAssembledInput
320
321end module main_drivers
Note: See TracBrowser for help on using the repository browser.