ElmerFoamFSI  2.0
ElmerFoamFSI is fluid-solid interaction simulation application built up from OpenFOAM CFD and Elmer CSM coupled through the IMPACT multiphysics software integration infrastructure.
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Groups Pages
ElmerCSCParallel.F90
Go to the documentation of this file.
1 !/*****************************************************************************/
2 ! *
3 ! * Elmer, A Finite Element Software for Multiphysical Problems
4 ! *
5 ! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6 ! *
7 ! * This library is free software; you can redistribute it and/or
8 ! * modify it under the terms of the GNU Lesser General Public
9 ! * License as published by the Free Software Foundation; either
10 ! * version 2.1 of the License, or (at your option) any later version.
11 ! *
12 ! * This library is distributed in the hope that it will be useful,
13 ! * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 ! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 ! * Lesser General Public License for more details.
16 ! *
17 ! * You should have received a copy of the GNU Lesser General Public
18 ! * License along with this library (in file ../LGPL-2.1); if not, write
19 ! * to the Free Software Foundation, Inc., 51 Franklin Street,
20 ! * Fifth Floor, Boston, MA 02110-1301 USA
21 ! *
22 ! *****************************************************************************/
23 
25 
26 !------------------------------------------------------------------------------
28 !------------------------------------------------------------------------------
29 
31 
32 !------------------------------------------------------------------------------
34 !------------------------------------------------------------------------------
35 
36 
37 #include "../config.h"
38 
39 ! ******************************************************************************
40 
41 SUBROUTINE initializer(global, runs, verbIn)
42  USE testobject
43  USE types
44  USE generalutils
45  USE parallelutils
46  USE timemodule
47 
48  IMPLICIT NONE
49 
50  include 'comf90.h'
51 
52  TYPE(t_global), POINTER :: global
53  INTEGER :: runs, counter, nindex, verbin
54 
55  REAL(KIND=dp) :: ct, rt
56 ! INTEGER, PARAMETER :: Init=0
57  CHARACTER(LEN=MAX_NAME_LEN) :: datestr, toutput
58  INTEGER, POINTER :: tempnodelist(:)
59  LOGICAL :: nodepresent
60  INTEGER :: lastfsielement
61  INTEGER :: paneid
62 ! INTEGER, POINTER :: Conn2(:)
63 #if defined(MINGW32)
64  CALL set_stdio_bufs()
65 #endif
66 
67  WRITE(6,*) 'ElmerCSCParallel:Initializer: Starting ...'
68 
69  !Allocate space for storing the verbosity
70  ALLOCATE(global%verbosity(1))
71  !Initialize verbosity to 1
72  global%verbosity = verbin
73 
74  WRITE(*,*) 'verbosity = ', global%verbosity
75 
76  CALL envir( 'ELMERSOLVER_OUTPUT_TOTAL'//char(0), toutput, tlen )
77  silent = toutput(1:1)=='0' .OR. toutput(1:5)=='false'
78 
79  ct = cputime()
80  rt = realtime()
81 
82  IF ( .NOT. silent ) THEN
83  datestr = formatdate()
84  WRITE( *,'(A,A)' ) "ELMER SOLVER (v " // version // ") STARTED AT: ", trim(datestr)
85  CALL flush(6)
86  END IF
87 
88  IF( myverbosity > 3) WRITE( *, * ) 'Calling ElmerInitialize ... '
89 
90  CALL elmerinitialize(runs, global%ElmerComm)
91 
92  global%MyModel => currentmodel
93 
94  !Assigning stuff & printing out to check
95  global%SolverId = -1
96  DO i=1, currentmodel % NumberOfSolvers
97  solver => global % MyModel % Solvers(i)
98  solverparams => getsolverparams(solver)
99  equation = getstring(solverparams, 'Equation', gotit)
100  WRITE(*,*) 'Equation = ',equation
101  IF( equation .eq. 'nonlinear elasticity' ) THEN
102  global%SolverId = i
103  EXIT
104  ELSE IF( equation .eq. 'linear elasticity' ) THEN
105  global%SolverId = i
106  EXIT
107  END IF
108  END DO
109 
110  !Find out which BC tag is the FSI BC
111  global%FSIbcId = -1
112  DO t=1,currentmodel % NumberOfBCs
113  isfsi = listgetlogical(currentmodel % BCs(t) % Values,&
114  'FSI BC', gotit )
115  IF ( isfsi ) THEN
116  IF (myverbosity > 3) WRITE(*,*) 'BC', t, 'is FSI'
117  global%FSIbcId = t
118  END IF
119  END DO
120 
121  IF (global%FSIbcId == -1) THEN
122  WRITE(*,*) 'Not registering data no FSI BC!'
123  RETURN
124  END IF
125 
126  !Determine the number of elements & nodes on the
127  !FSI BC so we can allocate the arrays
128  mymesh => currentmodel % Meshes
129  global%nNodes = 0
130  global%nElem = 0
131 
132  !Allocate a temporay array to avoid repeating nodes
133  WRITE (*,*) "Number of Bulk elemenets = ", mymesh % NumberOfBulkElements
134  WRITE (*,*) "Number of Boundary elemenets = ", mymesh % NumberOfBoundaryElements
135 
136  t = mymesh % NumberOfBulkElements
137  !WRITE (*,*) "Regular Element DIM = ", MyMesh % Elements(t) % TYPE % DIMENSION
138  mycurrentelement => mymesh % Elements(t+1)
139  t = mycurrentelement % TYPE % DIMENSION
140  counter = (mymesh % NumberOfBoundaryElements)*(t + 1)
141  ALLOCATE(tempnodelist(counter))
142  tempnodelist = -1
143 
144  !Loop over the boundary elements
145  counter = 0
146  DO t = mymesh % NumberOfBulkElements+1, &
147  mymesh % NumberOfBulkElements + &
148  mymesh % NumberOfBoundaryElements
149 
150  mycurrentelement => mymesh % Elements(t)
151  bc_id = getbcid(mycurrentelement)
152  !WRITE (*,*) "Boundary Element DIM = ", MyCurrentElement % TYPE % DIMENSION
153 
154  IF ( bc_id == global%FSIbcId ) THEN
155  lastfsielement = t
156  IF( myverbosity > 3) THEN
157  WRITE(*,*) 'Initializer element = ',t,' bc_id = ',bc_id
158  END IF
159  global%nElem = global%nElem + 1
160  DO nt = 1,mycurrentelement % TYPE % NumberOfNodes
161  mynodeindexes => mycurrentelement % NodeIndexes
162  IF( myverbosity > 3) THEN
163  WRITE(*,*) 'Node index = ', mynodeindexes(nt)
164  END IF
165  nodepresent = .false.
166  DO nk=1,counter
167  IF (tempnodelist(nk) == mynodeindexes(nt)) THEN
168  nodepresent = .true.
169  EXIT
170  END IF
171  END DO
172  IF ( nodepresent .eqv. .false. ) THEN
173  counter = counter + 1
174  IF( myverbosity > 4) THEN
175  WRITE(*,*) 'counter = ', counter
176  END IF
177  tempnodelist(counter) = mynodeindexes(nt)
178  END IF
179  END DO
180  END IF
181  END DO
182 
183  IF( myverbosity > 5) THEN
184  WRITE(*,*) 'SIZE(TempNodeList) =', SIZE(tempnodelist)
185  DO t = 1,SIZE(tempnodelist)
186  WRITE(*,*) tempnodelist(t)
187  END DO
188  END IF
189 
190  global%nNodes = counter
191 
192  IF( myverbosity > 4) THEN
193  WRITE(*,*) 'global%nElem = ',global%nElem
194  WRITE(*,*) 'global%nNodes = ',global%nNodes
195  END IF
196 
197  !Allocate space for storing the NodeDisplacements
198  !This location will be registered with IMPACT
199  IF( myverbosity > 3) WRITE(*,*) 'global%nNodes =', global%nNodes
200  ALLOCATE(global%NodeDisplacements(3*global%nNodes))
201  global%NodeDisplacements = 0.0
202  !Masoud
203  ALLOCATE(global%PreviousNodeDisplacements(3*global%nNodes))
204  global%PreviousNodeDisplacements = 0.d0
205  ALLOCATE(global%IsLoadFreeNode(global%nNodes))
206  global%IsLoadFreeNode = .false.
207  !Masoud End
208  ALLOCATE(global%MyToElmerNodes(global%nNodes))
209  IF( myverbosity > 3) THEN
210  WRITE(*,*) 'SIZE(NodeDisplacements)=',SIZE(global%NodeDisplacements)
211  END IF
212  ALLOCATE(currentmodel%ElmerToMyNodes(currentmodel % NumberOfNodes))
213  global%MyToElmerNodes = -1
214  currentmodel%ElmerToMyNodes = -1
215 
216  !Allocate space for storing the FaceLoads
217  !This location will be registered with IMPACT
218  ALLOCATE(global%FaceLoads(3*global%nElem))
219  global%FaceLoads = 0.0
220  !Allocate space for storing the FacePressures
221  !This location will be registered with IMPACT
222  ALLOCATE(global%FacePressures(global%nElem))
223  global%FacePressures = 0.0
224  !Allocate space for storing the NodePressures
225  !This location will be registered with IMPACT
226  ALLOCATE(global%NodePressures(global%nNodes))
227  global%NodePressures = 0
228  !Allocate space for storing the NodeLoads
229  !This location will be registered with IMPACT
230  ALLOCATE(global%NodeLoads(3*global%nNodes))
231  global%NodeLoads = 0.0
232  !Allocate space for storing the previous loads
233  !so linear interpolation can be done
234  ALLOCATE(global%PreviousLoads(3,global%nNodes))
235  global%PreviousLoads = 0.0
236  !Allocate the NodeLoads that gets passed to the user defined function
237  ALLOCATE(currentmodel%NodeLoadsPass(3,global%nNodes))
238  IF( myverbosity > 3) WRITE(*,*) 'SIZE(NodeLoads)=',SIZE(global%NodeLoads)
239  !Initialize the loads to 0.0
240  DO t = 1, global%nNodes
241  DO j =1,3
242 ! CurrentModel%NodeLoadsPass(j,t) = (t-1)*3.0 + j*1.0
243  currentmodel%NodeLoadsPass(j,t) = 0.0
244  global%PreviousLoads(j,t) = 0.0
245  END DO
246  END DO
247 
248  !Initialize the faceloads to 0.0
249 
250  !Allocate the mesh arrays now that their sizes are known
251  mycurrentelement => mymesh % Elements(lastfsielement)
252  IF (mycurrentelement % TYPE % DIMENSION == 1) THEN
253  IF (mycurrentelement % TYPE % NumberOfNodes == 2) THEN
254  global%MeshType = ':b2:'
255  global%nConn = 2
256  ELSE IF (mycurrentelement % TYPE % NumberOfNodes == 3) THEN
257  global%MeshType = ':b3:'
258  global%nConn = 3
259  ELSE
260  CALL fatal( ' ', 'FSI elements of 1 dimension must have 2 or 3 nodes!' )
261  END IF
262  ELSE IF (mycurrentelement % TYPE % DIMENSION == 2) THEN
263  IF( myverbosity > 3) THEN
264  WRITE(*,*) 'FSI BC element NumberOfNodes =',&
265  mycurrentelement % TYPE % NumberOfNodes
266  END IF
267  IF (mycurrentelement % TYPE % NumberOfNodes == 3) THEN
268  global%MeshType = ':t3:'
269  global%nConn = 3
270  ELSE IF (mycurrentelement % TYPE % NumberOfNodes == 6) THEN
271  global%MeshType = ':t6:'
272  global%nConn = 6
273  ELSE IF (mycurrentelement % TYPE % NumberOfNodes == 4) THEN
274  global%MeshType = ':q4:'
275  global%nConn = 4
276  ELSE IF (mycurrentelement % TYPE % NumberOfNodes == 8) THEN
277  global%MeshType = ':q8:'
278  global%nConn = 8
279  ELSE
280  CALL fatal( ' ', 'FSI elements of 2 dimensions must have 3 or 6 node&
281  triangles or 4 or 8 node quadrilaterals!' )
282  END IF
283  ELSE
284  CALL fatal( ' ', 'FSI elements cannot have dimension greater than 2!' )
285  END IF
286 
287  IF( myverbosity > 3) THEN
288  WRITE(*,*) 'global%MeshType = ',global%MeshType
289  WRITE(*,*) 'global%nConn =',global%nConn
290  END IF
291  ALLOCATE(global%Conn(global%nConn*global%nElem))
292  ALLOCATE(global%Coords(3*global%nNodes))
293 
294  !Populate the mesh arrays
295  !Loop over the boundary elements
296  counter = 0
297  DO t = mymesh % NumberOfBulkElements+1, &
298  mymesh % NumberOfBulkElements + &
299  mymesh % NumberOfBoundaryElements
300 
301  !WRITE(*,*) "ElmerCSCParallel element = ", t
302 
303  mycurrentelement => mymesh % Elements(t)
304  bc_id = getbcid(mycurrentelement)
305 
306  !For debugging - you can remove later
307  CALL getelementnodes(elementnodes,mycurrentelement)
308  !WRITE(*,*) 'Elmer elements:'
309  !DO nt = 1,MyCurrentElement % TYPE % NumberOfNodes
310  ! WRITE(*,*) nt,':',ElementNodes % x(nt),' ',ElementNodes % y(nt),' ',&
311  ! ElementNodes % z(nt)
312  !END DO
313  !end of debugging
314 
315  IF ( bc_id == global%FSIbcId ) THEN
316  counter = counter + 1
317  IF( myverbosity > 3) THEN
318  WRITE(*,*) 'element = ',t
319  END IF
320  CALL getelementnodes(elementnodes,mycurrentelement)
321 
322  global%FacePressures(counter) = 0.0
323  global%FaceLoads(3*(counter-1) + 1) = 0.0
324  global%FaceLoads(3*(counter-1) + 1) = 0.0
325  global%FaceLoads(3*(counter-1) + 1) = 0.0
326 
327  DO nt = 1,mycurrentelement % TYPE % NumberOfNodes
328  mynodeindexes => mycurrentelement % NodeIndexes
329 
330  IF( myverbosity > 3) THEN
331  WRITE(*,*) 'Node index = ', mynodeindexes(nt)
332  WRITE(*,*) elementnodes % x(nt), elementnodes % y(nt), &
333  elementnodes % z(nt)
334  END IF
335 
336  DO nk=1,global%nNodes
337  IF (tempnodelist(nk) == mynodeindexes(nt)) THEN
338  nindex = nk
339  EXIT
340  END IF
341  END DO
342 
343  !Conn2((counter-1)*2+nt) = 10
344  global%MyToElmerNodes(nindex) = mynodeindexes(nt)
345  currentmodel%ElmerToMyNodes(mynodeindexes(nt)) = nindex
346  global%Conn(global%nConn*(counter-1) + nt) = nindex
347  global%Coords(3*(nindex-1) + 1) = elementnodes % x(nt)
348  global%Coords(3*(nindex-1) + 2) = elementnodes % y(nt)
349  global%Coords(3*(nindex-1) + 3) = elementnodes % z(nt)
350  global%NodeLoads(3*(nindex-1) + 1) = 0.0
351  global%NodeLoads(3*(nindex-1) + 2) = 0.0
352  global%NodeLoads(3*(nindex-1) + 3) = 0.0
353  global%NodePressures(nindex) = 0.0
354  !Masoud: Marking LoadLessNodes
355  IF (abs(elementnodes % x(nt) - 0.6) < 10.0*epsilon(1.0)) THEN
356  !global%IsLoadFreeNode(nindex-1) = .TRUE.
357  END IF
358  !Masoud End
359  END DO
360  END IF
361  END DO
362 
363  !Masoud : printing IsLoadFreeNode
364  !WRITE(*,*) 'global%IsLoadFreeNode = ', global%IsLoadFreeNode
365  !Masoud End
366 
367  IF( myverbosity >= 5) THEN
368  WRITE(6,*) 'global%MyToElmerNodes = '
369  DO i=1,global%nNodes
370  WRITE(6,*) i,' = ',global%MyToElmerNodes(i)
371  END DO
372  WRITE(6,*) 'CurrentModel%ElmerToMyNodes = '
373  DO i=1,currentmodel % NumberOfNodes
374  WRITE(6,*) i,' = ',currentmodel%ElmerToMyNodes(i)
375  END DO
376 
377  WRITE(6,*)'-------------------------------------------'
378  WRITE(*,*) 'global%Conn = '
379  WRITE(6,*)'-------------------------------------------'
380  DO t = 1, global%nElem
381  WRITE(*,*) t,':',global%Conn(global%nConn*(t-1) + 1),&
382  global%Conn(global%nConn*(t-1) + 2)
383  END DO
384 
385  WRITE(6,*)
386 
387  WRITE(6,*)'-------------------------------------------'
388  WRITE(*,*) 'global%Coords = '
389  WRITE(6,*)'-------------------------------------------'
390  DO t = 1, global%nNodes
391  WRITE(*,*) t,':',global%Coords(3*(t-1)+1),global%Coords(3*(t-1)+2),&
392  global%Coords(3*(t-1)+3)
393  END DO
394  END IF
395 
396  WRITE(6,*)
397 
398  !Register the mesh
399  IF( myverbosity > 3) THEN
400  WRITE(*,*) 'window.mesh = ', &
401  trim(global%window_name)//'.'//trim(global%MeshType)
402  END IF
403  !Register data to pane starting from 100
404  paneid = 100 + global % procId
405  CALL com_set_size(trim(global%window_name)//'.nc',paneid,global%nNodes)
406  CALL com_set_array(trim(global%window_name)//'.nc',paneid,global%Coords,3)
407  CALL com_set_size(trim(global%window_name)//'.'//trim(global%MeshType),paneid,&
408  global%nElem)
409  CALL com_set_array(trim(global%window_name)//'.'//trim(global%MeshType),paneid,&
410  global%Conn, global%nConn)
411 
412  !CALL COM_SET_SIZE('Window1.:b2:real', 11,global%nElem)
413  !CALL COM_RESIZE_ARRAY('Window1.:b2:real',11)
414  !CALL COM_RESIZE_ARRAY(TRIM(global%window_name)//'.'//TRIM(global%MeshType),11)
415  !CALL COM_SET_ARRAY('Window1.:b2:real',11,Conn,2)
416 
417  !Set the displacments array with COM now that the mesh is registered
418  CALL com_new_dataitem(trim(global%window_name)//'.Displacements', 'n', com_double_precision, 3, 'm')
419  CALL com_set_array(trim(global%window_name)//'.Displacements',paneid,&
420  global%NodeDisplacements,3)
421  !Masoud
422  ! n: means node quantity, 3: number of columns, 11: window number (currently
423  ! default used everywhere), 'm': is units of the quanitity
424  CALL com_new_dataitem(trim(global%window_name)//'.PreviousDisplacements', 'n', com_double_precision, 3, 'm')
425  CALL com_set_array(trim(global%window_name)//'.PreviousDisplacements',paneid,&
426  global%PreviousNodeDisplacements,3)
427  !Masoud End
428 
429  !Set the loads array with COM now that the mesh is registered
430  CALL com_new_dataitem(trim(global%window_name)//'.Loads', 'n', com_double_precision, 3, '')
431  CALL com_set_array(trim(global%window_name)//'.Loads',paneid,&
432  global%NodeLoads,3)
433 
434  !Set the loads array with COM now that the mesh is registered
435  CALL com_new_dataitem(trim(global%window_name)//'.NodePressures',&
436  'n', com_double_precision, 1, '')
437  CALL com_set_array(trim(global%window_name)//'.NodePressures',paneid,&
438  global%NodePressures,1)
439 
440  !Set the pressures array with COM now that the mesh is registered
441  CALL com_new_dataitem(trim(global%window_name)//'.Pressures', 'e',&
442  com_double_precision, 1, '')
443  CALL com_set_array(trim(global%window_name)//'.Pressures',paneid,&
444  global%FacePressures,1)
445 
446  !Set the face loads array with COM now that the mesh is registered
447  CALL com_new_dataitem(trim(global%window_name)//'.FaceLoads', 'e',&
448  com_double_precision, 3, '')
449  CALL com_set_array(trim(global%window_name)//'.FaceLoads',paneid,&
450  global%FaceLoads,3)
451 
452  !Set the verbosity with COM
453  !CALL COM_NEW_DATAITEM(TRIM(global%window_name)//'.verbosity',&
454  ! 'w', COM_INTEGER, 1, '')
455  !CALL COM_SET_ARRAY(TRIM(global%window_name)//'.verbosity',11,&
456  ! global%verbosity,1)
457 
458  previoustime = 0.0
459 
460  IF( myverbosity > 4) THEN
461  WRITE(6,*)'-------------------------------------------'
462  WRITE(*,*) 'global%NodeDisplacements = '
463  WRITE(6,*)'-------------------------------------------'
464  DO i = 1,global%nNodes
465  DO j=1,3
466  global%NodeDisplacements(3*(i-1) + j) = 0.0
467  ENDDO
468  WRITE(6,*) (global%NodeDisplacements(3*(i-1) + j),j=1,3)
469  ENDDO
470 
471  WRITE(6,*)
472  END IF
473 
474 
475  WRITE(6,*) 'ElmerCSCParallel:Initializer: Done.....'
476 
477  nullify(tempnodelist)
478 
479 END SUBROUTINE initializer
480 
481 ! *****************************************************************************
482 
483 SUBROUTINE run(global, runs, tFinal)
484 
485  USE testobject
486  USE generalmodule
487  USE timemodule
488 
489  IMPLICIT NONE
490 
491  include 'comf90.h'
492 
493  TYPE(t_global), POINTER :: global
494  INTEGER :: runs
495  DOUBLE PRECISION, INTENT(IN) :: tfinal
496  DOUBLE PRECISION :: standardtimestep
497  DOUBLE PRECISION :: var1, var2, var3
498  DOUBLE PRECISION :: deltatime
499  !INTEGER, PARAMETER :: Initialize=0
500  !INTEGER :: TimeIntervals, CoupledMinIter, CoupledMaxIter, Timestep
501  !INTEGER, POINTER, SAVE :: OutputIntervals(:)
502  !LOGICAL :: LastSaved, Transient, Scanning
503 
504  INTERFACE
505  SUBROUTINE updatedisplacements(global,runs,tFinal)
506  USE testobject
507  USE generalmodule
508  USE timemodule
509 
510  TYPE(t_global), POINTER :: global
511  INTEGER :: runs
512  DOUBLE PRECISION :: tfinal
513 
514  END SUBROUTINE updatedisplacements
515  END INTERFACE
516 
517  currentmodel => global%MyModel
518 
519  finaltime = tfinal
520 
521  IF (previoustime >= finaltime) THEN
522  CALL fatal( ' ', 'Next time step passed to RUN is greater than &
523  previous time step!')
524  END IF
525 
526  IF (myverbosity > 3) WRITE(*,*) 'In RUN function'
527 
528  standardtimestep = timestepsizes(1,1)
529  deltatime = tfinal - stime(1)
530  timesteps(1) = floor(deltatime/standardtimestep)
531 
532  IF (myverbosity > 3) THEN
533  WRITE(*,*) 'tFinal = ', tfinal
534  WRITE(*,*) 'sTime(1) = ', stime(1)
535  WRITE(*,*) 'standardTimestep = ', standardtimestep
536  END IF
537 
538  var3 = deltatime - int(deltatime/standardtimestep)*standardtimestep
539 
540  IF ( abs(var3) < 1.0d-12 ) THEN
541  timeintervals = 1
542  IF (myverbosity > 3) WRITE(*,*) 'standardTimestep evenly divides time'
543  ELSE
544  timestepsizes(2,1) = deltatime - timesteps(1)*standardtimestep
545  timeintervals = 2
546  timesteps(2) = 1
547  IF (outputintervals(1) == 1) THEN
548  outputintervals(2) = 1
549  END IF
550  IF (myverbosity > 3) THEN
551  WRITE(*,*) 'Remainder for time'
552  WRITE(*,*) 'TimestepSizes(2,1) = ', timestepsizes(2,1)
553  END IF
554  END IF
555 
556  IF (myverbosity > 3) THEN
557  WRITE(*,*) 'In Run function'
558  WRITE(*,*) 'tFinal = ', tfinal
559  WRITE(*,*) 'sTime(1) = ', stime(1)
560  WRITE(*,*) 'Timesteps(1) = ', timesteps(1)
561  WRITE(*,*) 'standardTimestep = ', standardtimestep
562  !Printing out the Pressures as a check
563  WRITE(*,*) 'Pressures '
564  DO t = 1, global%nElem
565  WRITE(*,*) global%FacePressures(t)
566  END DO
567  !Printing out the FaceLoads as a check
568  WRITE(*,*) 'FaceLoads '
569  DO t = 1, global%nElem
570  WRITE(*,*) global%FaceLoads(3*(t-1) + 1), global%FaceLoads(3*(t-1) + 2),&
571  global%FaceLoads(3*(t-1) + 3)
572  END DO
573  !Print out the NodePressures as a check
574  WRITE(*,*) 'NodePressures '
575  DO t = 1, global%nNodes
576  WRITE(*,*) global%NodePressures(t)
577  END DO
578 
579  END IF
580 
581 !------------------------------------------------------------------------------
582 ! Here we actually start the simulation ....
583 ! First go trough timeintervals
584 !------------------------------------------------------------------------------
585  !WRITE(*,*) 'global%nNodes = ',global%nNodes
586  IF (myverbosity > 3) WRITE(*,*) 'Calling TimeStepper'
587  CALL timestepper(global,runs)
588  IF (myverbosity > 3) WRITE(*,*) 'Done with TimeStepper'
589 
590  IF (global%FSIbcId /= -1) THEN
591  IF (myverbosity > 3) WRITE(*,*) 'Calling UpdateDisplacements'
592  CALL updatedisplacements(global,runs, tfinal)
593  END IF
594 
595  !Setting PreviousLoads values to NodeLoads now that
596  !Run is finished
597  DO t = 1, global%nNodes
598  DO j =1,3
599  global%PreviousLoads(j,t) = global%NodeLoads(3*(t-1) + j)
600  END DO
601  END DO
602 
603  !Setting PreviousTime to FinalTime now that
604  !Run is finished
605  previoustime = finaltime
606 
607  WRITE(*,*) 'ElmerCSCParallel:Run: runs = ',runs
608  WRITE(*,*) 'ElmerCSCParallel:Run: CurrentModel%GetTestLoads = ',&
609  currentmodel%GetTestLoads
610  WRITE(*,*) 'ElmerCSCParallel:Run: CurrentModel%UDFUsed = ',&
611  currentmodel%UDFUsed
612  IF (runs == 1 .AND. currentmodel%GetTestLoads .eqv. .true. &
613  .AND. currentmodel%UDFUsed .eqv. .true.) THEN
614  runs = 2
615  END IF
616 
617 END SUBROUTINE run
618 
619 !**************************************************************************
620 
621 SUBROUTINE updatedisplacements(global,runs, tFinal)
622  USE testobject
623  USE generalmodule
624  USE timemodule
625 
626  IMPLICIT NONE
627 
628  include 'comf90.h'
629 
630  TYPE(t_global), POINTER :: global
631  INTEGER :: runs
632  DOUBLE PRECISION :: tfinal
633  INTEGER :: ncount, counter, ncountmax
634 
635  WRITE(*,*) 'ElmerCSCParallel:UpdateDisplacements:',&
636  ' Reporting displacements to the fluid solver'
637 
638  !access the displacements here
639  solver => currentmodel % Solvers(global%SolverId)
640 
641  !IF (TRIM(Equation) == "nonlinear elasticity") THEN
642  stresssol => solver % Variable
643  displacement => stresssol % Values
644  myperm => stresssol % Perm
645 
646  IF( myverbosity > 3) THEN
647  WRITE(*,*) 'The SIZE(Displacement) = ', SIZE(displacement)
648  WRITE(*,*) 'After SolverActivate Call, Displacement(1) = ',&
649  displacement(1)
650  DO t=1,SIZE(displacement)
651  WRITE(*,*) 'Displacement',t,'=',displacement(t)
652  END DO
653  END IF
654 
655  !Jess adding stuff to try accessing wetted surface
656  elementcount = getnofactive()
657  boundaryelementcount = getnofboundaryelements()
658 
659  mymesh => currentmodel % Meshes
660 
661  !Write out nodes, coordinates, and connectivities
662  !Store at registered address
663  ncount = 0
664  ncountmax = 0
665  counter = 0
666  IF( myverbosity > 3) THEN
667  WRITE(*,*) 'Number of bulk elements = ', mymesh % NumberOfBulkElements
668  WRITE(*,*) 'Number of boundary elements = ', mymesh % NumberOfBoundaryElements
669  END IF
670  DO t = mymesh % NumberOfBulkElements+1, &
671  mymesh % NumberOfBulkElements + &
672  mymesh % NumberOfBoundaryElements
673  mycurrentelement => mymesh % Elements(t)
674  bc_id = getbcid(mycurrentelement)
675  IF ( bc_id == global%FSIbcId ) THEN
676  counter = counter+1
677  IF( myverbosity > 3) THEN
678  WRITE(*,*) '*************************************'
679  WRITE(*,*) 'element ',counter,'on FSI boundary'
680  WRITE(*,*) 'Update t = ',t,' bc_id = ',bc_id
681  WRITE(*,*)
682  WRITE(*,*) 'Element no. of nodes:', &
683  mycurrentelement % TYPE % NumberOfNodes
684  END IF
685 
686  CALL getelementnodes(elementnodes,mycurrentelement)
687 
688  DO nt = 1,mycurrentelement % TYPE % NumberOfNodes
689  mynodeindexes => mycurrentelement % NodeIndexes
690 
691  IF( myverbosity > 3) THEN
692  WRITE(*,*) '*********************'
693  WRITE(*,*) 'Node index = ', mynodeindexes(nt)
694  WRITE(*,*) elementnodes % x(nt), elementnodes % y(nt), &
695  elementnodes % z(nt)
696  WRITE(*,*) 'MyPerm = ', myperm(mynodeindexes(nt))
697  END IF
698 
699  IF ( myperm(mynodeindexes(nt)) > 0 ) THEN
700  nk = (stresssol % DOFs)*(myperm(mynodeindexes(nt)) - 1)
701  IF( myverbosity > 3) WRITE(*,*) 'DOFs = ', stresssol % DOFs
702 
703  ncount = currentmodel%ElmerToMyNodes(mynodeindexes(nt))
704 
705  ! Masoud : Avoding displacement update duplicates and
706  ! implementing a correct displacement update
707 
708  IF ( ncount > ncountmax ) THEN
709  IF( myverbosity > 3) THEN
710  WRITE(6,*) 'Updating NodeDisplacement(',ncount,')'
711  END IF
712 
713  IF ( stresssol % DOFs == 1 ) THEN
714  IF( myverbosity > 3) WRITE(*,*) displacement(nk+1), 0.0, 0.0
715  global%NodeDisplacements(3*(ncount-1) + 1) = displacement(nk+1)
716  global%NodeDisplacements(3*(ncount-1) + 2) = 0.0d0
717  global%NodeDisplacements(3*(ncount-1) + 3) = 0.0d0
718  ELSE IF ( stresssol % DOFs == 2 ) THEN
719  IF( myverbosity > 3) WRITE(*,*) displacement(nk+1), displacement(nk+2)
720  global%NodeDisplacements(3*(ncount-1) + 1) = displacement(nk+1)
721  global%NodeDisplacements(3*(ncount-1) + 2) = displacement(nk+2)
722  global%NodeDisplacements(3*(ncount-1) + 3) = 0.0d0
723  ELSE IF ( stresssol % DOFs == 3 ) THEN
724  !Calculating displacement differences to pass to fluid solver
725  global%NodeDisplacements(3*(ncount-1) + 1) = displacement(nk+1)&
726  - global%PreviousNodeDisplacements(3*(ncount-1) + 1)
727  global%NodeDisplacements(3*(ncount-1) + 2) = displacement(nk+2)&
728  - global%PreviousNodeDisplacements(3*(ncount-1) + 2)
729  global%NodeDisplacements(3*(ncount-1) + 3) = displacement(nk+3)&
730  - global%PreviousNodeDisplacements(3*(ncount-1) + 3)
731  !Updating previous displacements to the current values for the
732  !next timestep
733  global%PreviousNodeDisplacements(3*(ncount-1) + 1) = displacement(nk+1)
734  global%PreviousNodeDisplacements(3*(ncount-1) + 2) = displacement(nk+2)
735  global%PreviousNodeDisplacements(3*(ncount-1) + 3) = displacement(nk+3)
736  !Priniting increament displacement values for the user
737  !WRITE(*,*) ' increament disps= ', global%NodeDisplacements(3*(nCount-1) + 1),&
738  !global%NodeDisplacements(3*(nCount-1) + 2), &
739  !global%NodeDisplacements(3*(nCount-1) + 3)
740  !WRITE(*,*) '----------------------'
741  ! Masoud : End
742  ELSE
743  WRITE(*,*) 'StressSol % DOFs = ', stresssol % DOFs
744  WRITE(*,*) 'DOFs are assumed to be <= 3'
745  CALL fatal( ' ', 'StressSol DOFs are greater than 3 ' )
746  END IF
747  ncountmax = max(ncount, ncountmax)
748  END IF
749  ! Masoud: End
750 
751  END IF
752 
753  END DO
754  END IF
755  END DO
756 
757  IF( myverbosity > 50) THEN
758  WRITE(6,*) 'tFinal = ', tfinal
759  WRITE(6,*) 'NodeDisplacements'
760  DO i = 1,global%nNodes
761  WRITE(6,*) (global%NodeDisplacements(3*(i-1) + j),j=1,3)
762  ENDDO
763  WRITE(6,*) 'Exiting UpdateDisplacements'
764  END IF
765 
766 END SUBROUTINE updatedisplacements
767 
768 ! ********************************************************************
769 
770 SUBROUTINE finalize(global, runs)
771  USE testobject
772  USE types
773  USE generalutils
774  USE parallelutils
775  USE generalmodule, only : firsttime
776 
777  IMPLICIT NONE
778 
779  include 'comf90.h'
780 
781  TYPE(t_global), POINTER :: global
782  INTEGER :: runs
783 
784  !INTEGER, PARAMETER :: Initialize=0
785  INTEGER :: tlen
786  LOGICAL :: silent
787  CHARACTER(LEN=MAX_NAME_LEN) :: datestr, toutput
788 
789  currentmodel => global%MyModel
790 
791  WRITE(*,*) 'ElmerCSCParallel:Finalize: Finishing simulation'
792 
793  CALL elmerfinalize(runs)
794 
795  CALL envir( 'ELMERSOLVER_OUTPUT_TOTAL'//char(0), toutput, tlen )
796  silent = toutput(1:1)=='0' .OR. toutput(1:5)=='false'
797 
798  IF ( .NOT. silent ) THEN
799  IF ( parenv % myPE == 0 ) THEN
800  !WRITE( *,'(a,F12.2,F12.2)' ) 'SOLVER TOTAL TIME(CPU,REAL): ', &
801  ! CPUTime()-CT, RealTime()-RT
802  datestr = formatdate()
803  WRITE( *,'(A,A)' ) 'ELMER SOLVER FINISHED AT: ', trim(datestr)
804  END IF
805  END IF
806 
807  WRITE(6,*) 'ElmerCSCParallel:Finalize: End Finalize. FirstTime = ',firsttime
808 
809 END SUBROUTINE finalize
810 
811 ! ********************************************************************
812 
813 SUBROUTINE elmercscparallel_load_module(name)
814 
815  USE testobject
816 
817  IMPLICIT NONE
818 
819  include "comf90.h"
820 
821  INTERFACE
822  SUBROUTINE com_set_pointer(attr,ptr,asso)
823  USE testobject
824  CHARACTER(*), INTENT(IN) :: attr
825  TYPE(t_global), POINTER :: ptr
826  EXTERNAL asso
827  END SUBROUTINE com_set_pointer
828 
829  SUBROUTINE initializer(global, runs, verbIn)
830  USE testobject
831  TYPE(t_global), POINTER :: global
832  INTEGER :: runs, verbin
833  END SUBROUTINE initializer
834 
835  SUBROUTINE run(global, runs, tFinal)
836  USE testobject
837  TYPE(t_global), POINTER :: global
838  INTEGER :: runs
839  DOUBLE PRECISION, INTENT(IN) :: tfinal
840  END SUBROUTINE run
841 
842  SUBROUTINE updatedisplacements(global, runs, tFinal)
843  USE testobject
844  TYPE(t_global), POINTER :: global
845  INTEGER :: runs
846  DOUBLE PRECISION :: tfinal
847  END SUBROUTINE updatedisplacements
848 
849  SUBROUTINE finalize(global, runs)
850  USE testobject
851  TYPE(t_global), POINTER :: global
852  INTEGER :: runs
853  END SUBROUTINE finalize
854 
855  END INTERFACE
856 
857  CHARACTER(*),intent(in) :: name
858  INTEGER :: com_types(7)
859  TYPE(t_global), POINTER :: glb
860  INTEGER :: ierror
861  INTEGER, target :: nproc, procid
862 
863  WRITE(*,'(A)') "Loading ElmerCSCParallel: "//trim(name)
864 
865 
866  ALLOCATE(glb)
867  glb%window_name = trim(name)
868  glb%other_window_handle = -1
869  glb%c_window_handle = -1
870  CALL com_new_window(trim(name), mpi_comm_null)
871 
872  glb%ElmerComm = com_get_default_communicator()
873  CALL mpi_comm_size(glb%ElmerComm, nproc, ierror)
874  glb%nProc = nproc
875  CALL mpi_comm_rank(glb%ElmerComm, procid, ierror)
876  glb%procId = procid
877 
878  CALL com_new_dataitem(trim(name)//'.global','w',com_f90pointer,1,'')
879  CALL com_allocate_array(trim(name)//'.global')
880  CALL com_new_dataitem(trim(name)//'.nProc','w',com_integer,1,'')
881  CALL com_set_size(trim(name)//'.nProc',0,1)
882  CALL com_set_array(name//'.nProc', 0, glb%nProc)
883 
884  com_types(1) = com_f90pointer
885  com_types(2) = com_integer
886  com_types(3) = com_integer
887 
888  CALL com_set_member_function(trim(name)//'.Initialize',initializer, &
889  trim(name)//'.global','bbi',com_types)
890 
891  com_types(3) = com_double_precision
892 
893  CALL com_set_member_function(trim(name)//'.Run',run, &
894  trim(name)//'.global','bbi',com_types)
895 
896  CALL com_set_member_function(trim(name)//'.Finalize',finalize, &
897  trim(name)//'.global','bb',com_types)
898 
899  CALL com_window_init_done(name)
900 
901  CALL com_set_pointer(name//'.global',glb,associate_pointer )
902 
903 END SUBROUTINE elmercscparallel_load_module
904 
905 
906 SUBROUTINE elmercscparallel_unload_module(name)
907  USE testobject
908  IMPLICIT NONE
909  include "comf90.h"
910  INTERFACE
911  SUBROUTINE com_get_pointer(attr,ptr,asso)
912  USE testobject
913  CHARACTER(*), INTENT(IN) :: attr
914  TYPE(t_global), POINTER :: ptr
915  EXTERNAL asso
916  END SUBROUTINE com_get_pointer
917  END INTERFACE
918  character(*),intent(in) :: name
919  TYPE(t_global), POINTER :: glb
920  INTEGER :: window_handle,other_window_handle,c_window_handle, owlen
921 
922  WRITE(*,'(A)') "Unloading ElmerCSCParallel: "//trim(name)
923  nullify(glb)
924  window_handle = com_get_window_handle(trim(name))
925  if(window_handle .gt. 0) then
926  CALL com_get_pointer(trim(name)//'.global',glb,associate_pointer)
927  IF(ASSOCIATED(glb).eqv..true.) THEN
928  WRITE(*,'(A)') 'Fortran module '//trim(glb%window_name)//' unloading name '//trim(name)
929  if(glb%other_window_handle .gt. 0) then
930  WRITE(*,*) 'Fortran module '//trim(glb%window_name)//&
931  ' unloading external Fortran module '//trim(glb%other_window_name)//'.'
932  other_window_handle = com_get_window_handle(trim(glb%other_window_name))
933  IF(other_window_handle .gt. 0) THEN
934  CALL com_unload_module("ElmerCSCParallel",trim(glb%other_window_name))
935  ENDIF
936  endif
937  if(glb%c_window_handle .gt. 0) then
938  WRITE(*,*) 'Fortran module '//trim(glb%window_name)//&
939  ' unloading external C module '//trim(glb%c_window_name)//'.'
940  c_window_handle = com_get_window_handle(trim(glb%c_window_name))
941  IF(c_window_handle .gt. 0) THEN
942  CALL com_unload_module("ElmerCSCParallel",trim(glb%c_window_name))
943  ENDIF
944  endif
945 ! DEALLOCATE(glb)
946  ENDIF
947  CALL com_delete_window(trim(name))
948  endif
949 END SUBROUTINE elmercscparallel_unload_module
950 
951 ! ******************************************************************************
952 
timemodule
ElmerLib Elmer library routines
Definition: TimeModule.F90:55
elmercscparallel_unload_module
subroutine elmercscparallel_unload_module(name)
Definition: ElmerCSCParallel.F90:906
timemodule::timestepper
subroutine timestepper(global, runs)
Definition: TimeModule.F90:184
elmerfinalize
subroutine elmerfinalize(runs)
ElmerLib Elmer library routines
Definition: Finalize.F90:55
finalize
subroutine finalize(global, runs)
Definition: ElmerCSC.F90:811
testobject::t_global
Definition: ModTestObject.F90:5
elmercscparallel_load_module
subroutine elmercscparallel_load_module(name)
Definition: ElmerCSCParallel.F90:813
elmerinitialize
subroutine elmerinitialize(runs)
ElmerLib Elmer library routines
Definition: Initialize.F90:55
updatedisplacements
subroutine updatedisplacements(global, runs, tFinal)
Definition: ElmerCSC.F90:605
testobject::associate_pointer
subroutine associate_pointer(attr, ptr)
Definition: ModTestObject.F90:30
generalmodule
ElmerLib Elmer library routines
Definition: GeneralModule.F90:55
initializer
subroutine initializer(global, runs, verbIn)
A caller for the Elmer main program.
Definition: ElmerCSC.F90:41
run
subroutine run(global, runs, tFinal)
Definition: ElmerCSC.F90:467