source: src/main/main_drivers.F90 @ c28deba

Revision c28deba, 9.5 KB checked in by Oleg Batrashev <ogbash@…>, 3 years ago (diff)

Compact mesh plotting code.
  • 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    if (ismaster()) then ! MASTER
91       write(stream,*)
92       write(stream,*) 'master thread'
93       if (D_MSGLVL > 1) &
94            call MasterCtrlData_print()
95
96    else ! SLAVES
97       write(stream,'(a,i4,a)') 'slave [',myrank,'] thread'
98       if (D_MSGLVL > 1) &
99            call SharedCtrlData_print()
100    end if
101
102    ! =======================
103    ! Mesh and its Dual Graph
104    ! =======================
105    !
106    ! Create Mesh object
107
108    Msh = Mesh_New()
109
110    if (ismaster()) then
111       ! Initialise Mesh object
112       call Mesh_initFromFile(Msh, trim(mctls%info_file))
113    endif
114
115    ! Get from master Mesh's parameters: nell, ngf, mfrelt, nsd, nnode
116    call Mesh_paramsMPIBCAST(Msh)
117    if (D_MSGLVL > 1) &
118         call Mesh_printInfo(Msh)
119
120    ! Master reads in from files: feedom lists, coordinates, freedom map
121    if (ismaster()) then
122       call Mesh_readFromFile(Msh, &
123            fnFreelists = trim(mctls%freedom_lists_file), &
124            fnCoords    = trim(mctls%coords_file),        &
125            fnFreemap   = trim(mctls%freemap_file))
126    end if
127
128    ! Master non-blockingly sends mesh data: nfrelt, mhead
129    call Mesh_dataMPIISENDRECV(Msh, &
130         nfrelt   = .true., &
131         mhead    = .true.)
132    if (D_MSGLVL > 4) &
133         call Mesh_printElemFree(Msh)
134
135    ! For multi-variable problems which have more than one block
136    if (sctls%number_of_blocks > 1) then
137       if (ismaster()) then
138          call Mesh_readFromFile(Msh, &
139               fnFreemask = trim(mctls%freedom_mask_file))
140       end if
141       call Mesh_dataMPIISENDRECV(Msh, freemask = .true.)
142    end if
143
144
145    ! Build dual graph (Graph object is a data field in Mesh class)
146    ! let all procs do this - compare whith broadcasting the one bult on master
147    call Mesh_buildGraphDual(Msh)
148
149    ! Partition mesh's dual graph
150    if (ismaster()) then
151       if (sctls%plotting == D_PLOT_YES) then
152          call Mesh_pl2D_mesh(Msh)
153       end if
154
155       ! Partition graph: D_PART_PMETIS, D_PART_KMETIS, D_PART_VKMETIS
156       call Mesh_partitionDual(Msh, nparts, D_PART_VKMETIS, part_opts)
157       if (sctls%plotting == D_PLOT_YES) then
158          call Mesh_pl2D_partitions(Msh)
159       end if
160    endif
161
162    ! Distribute elements to partitons map among slaves
163    call Mesh_dataMPIISENDRECV(Msh, eptnmap=.true.)
164
165    ! Calculate number of elements in partitions
166    call Mesh_calcElemsInParts(Msh)
167
168    ! Build global to local, local to global maps and
169    ! inner/interface mask for freedoms
170    ! (also finds local number of freedoms 'Mesh%nlf';
171    !  this hidden appearance of 'Mesh_findNLF()' helps
172    !  to speed up calculations a bit)
173    call Mesh_buldMapsNghbrsMasksNLF(Msh)
174
175    ! ===============================
176    ! Assemble system matrix and RHS
177    ! ===============================
178    A = SpMtx_New()
179    allocate(b(Msh%nlf))
180    b = 0.0_rk
181
182    if (ismaster()) then
183       if (numprocs>1.and.present(A_interf)) then
184          A_interf = SpMtx_New()
185          call ElemMtxs_readAndDistribute(Msh, trim(mctls%elemmat_rhs_file), A, b, A_interf)
186       else
187          call ElemMtxs_readAndDistribute(Msh, trim(mctls%elemmat_rhs_file), A, b)
188       end if
189    else
190       if (numprocs>1.and.present(A_interf)) then
191          A_interf = SpMtx_New()
192          call ElemMtxs_recvAndAssemble(Msh, A, b, A_interf)
193       else
194          call ElemMtxs_recvAndAssemble(Msh, A, b)
195       end if
196    end if
197
198    if (sctls%verbose>9) then
199      write(stream,'(/a)') 'System matrix:'
200      call SpMtx_printInfo(A)
201      if (A%nrows <= 25) then
202         call SpMtx_printMat(A)
203      else if ((A%nrows > 25).and.(A%nrows <= 100)) then
204         call SpMtx_printRaw(A)
205      end if
206    endif
207
208    ! ==================
209    ! Finish assemble local RHS
210    ! ==================
211    if (A%nrows <= 25) & ! if (D_MSGLVL > 2) &
212         call Vect_Print(b, 'RHS assembled (local) ')
213
214    ! initialise auxiliary data for manipulating vectors
215    call Vect_setIntfEnd(sum(Msh%inner_interf_fmask))
216    call Vect_buildDotMask(Msh)
217    if (D_MSGLVL > 4) &
218         call Vect_Print(dot_intf_fmask,'dot_intf_fmask ')
219    call Vect_buildDotMap()
220    if (D_MSGLVL > 4) &
221         call Vect_Print(dot_intf_fmap,'dot_intf_fmap ')
222
223    ! Free mesh graph, not needed anymore
224    call Mesh_destroyGraph(Msh)
225
226    ! Update inner node count
227    Msh%ninner=Msh%nlf
228
229  end subroutine parallelAssembleFromElemInput
230
231
232  !----------------------------------------------------------------
233  !> Distribute assembled matrix and RHS from master to slaves
234  !----------------------------------------------------------------
235  subroutine parallelDistributeAssembledInput(Msh, A, b, A_interf)
236    implicit none
237
238    type(Mesh),     intent(in out) :: Msh !< Mesh
239    type(SpMtx),    intent(in out) :: A !< System matrix
240    float(kind=rk), dimension(:), pointer :: b !< local RHS
241    type(SpMtx),intent(in out),optional :: A_interf !< matrix at interface
242
243    integer :: n
244
245    ! ======================
246    ! Read matrix from file
247    ! ======================
248    if (ismaster()) then
249      write(stream,'(a,a)') ' ##### Assembled input file: ##### ', &
250            mctls%assembled_mtx_file
251      call ReadInSparseAssembled(A,trim(mctls%assembled_mtx_file))
252      allocate(b(A%nrows))
253      if (len_trim(mctls%assembled_rhs_file)>0) then
254        write(stream,'(a,a)') ' ##### Assembled RHS file: ##### ', &
255              mctls%assembled_rhs_file
256        call Vect_ReadFromFile(b, trim(mctls%assembled_rhs_file), mctls%assembled_rhs_format)
257      else
258        b=1.0_rk
259      end if
260    endif
261
262    ! =====================
263    ! Build mesh structure/distribute
264    ! =====================
265    if (numprocs==1) then
266      n=sqrt(1.0_rk*A%nrows)
267      if (n*n /= A%nrows) then
268        write (stream,*) 'Not a Cartesian Mesh!!!'
269        Msh=Mesh_New()
270        Msh%ngf=A%nrows
271        Msh%nlf=A%nrows
272        Msh%ninner=Msh%ngf
273      else
274        write (stream,*) 'Cartesian Mesh!!!'
275        call Mesh_BuildSquare(Msh,n)
276        Msh%ninner=Msh%ngf
277      endif
278    else ! numprocs>1
279      Msh=Mesh_New()
280      call SpMtx_DistributeAssembled(A,b,A_interf,Msh)
281    endif
282
283  end subroutine parallelDistributeAssembledInput
284
285end module main_drivers
Note: See TracBrowser for help on using the repository browser.