CoCalc -- MainUtils.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/src/MainUtils.F90
⁵²¹³ views
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
!
24
!/******************************************************************************
25
! *
26
! *  Authors: Juha Ruokolainen, Ville Savolainen
27
! *  Email:   [email protected]
28
! *  Web:     http://www.csc.fi/elmer
29
! *  Address: CSC - IT Center for Science Ltd.
30
! *           Keilaranta 14
31
! *           02101 Espoo, Finland 
32
! *
33
! *  Original Date: 08 Jun 1997
34
! *
35
! *****************************************************************************/
36

37
!> \ingroup ElmerLib
38
!> \}
39

40
#include "../config.h"
41
!------------------------------------------------------------------------------
42
!>  Utility routines for the Elmer main program.
43
!------------------------------------------------------------------------------
44
MODULE MainUtils
45
!------------------------------------------------------------------------------
46
  USE BlockSolve
47
  USE IterSolve, ONLY : NumericalError
48
  USE LoadMod, ONLY : ExecLocalAssembly, ExecSolver
49
  USE ModelDescription, ONLY : GetProcAddr
50
  USE MatrixAssembly, ONLY : CreateChildMatrix, GlueLocalSubMatrix, MoveRow, &
51
      SetMatrixElement
52
  USE ElementDescription, ONLY : SwapRefElemNodes
53
  USE ElementUtils, ONLY : CreateOdeMatrix, CreateMatrix
54
  
55
  USE MeshUtils, ONLY : BackCoordinateTransformation, Colouring_deallocate, &
56
      CoordinateTransformation, CreateDiscontMesh, ElmerColouringToGraph, &
57
      ElmerGraphColour, ElmerMeshToDualGraph, Graph_deallocate, LoadMesh2, &
58
      MakePermUsingMask, MeshStabParams, ReleaseMesh, SetActivEelementsTable, &
59
      SetCurrentMesh, SetMeshMaxDOFs, SplitMeshEqual, TransferCoordAndTime, &
60
      UpdateSolverMesh
61
  
62
  USE SolverUtils, ONLY : CalculateEntityWeights, &
63
      CalculateNodalWeights, CheckStepSize, ComputeChange, &
64
      ComputeNorm, CreateIpPerm, GenerateProjectors, GetPassiveBoundary, &
65
      InitializeTimestep, InitializeToZero, InvalidateVariable, &
66
      MatrixVectorMultiply, UpdateDependentObjects, UpdateExportedVariables, &
67
      FinalizeLumpedMatrix
68

69
  USE DefUtils, ONLY : GetString, GetCReal, GetElementNOFNodes, GetLogical, &
70
      DefaultDirichletBCs, GetMesh, GetInteger, GetMatrix, GetElementNOFNodes, &
71
      GetBodyForce, GetBC, Default2ndOrderTime, DefaultFinishBulkAssembly, &
72
      DefaultUpdateMass, DefaultFinishBoundaryAssembly, DefaultInitialize, &
73
      DefaultUpdateDamp, DefaultFinishAssembly, Default1stOrderTime
74
!------------------------------------------------------------------------------
75
  IMPLICIT NONE
76
!------------------------------------------------------------------------------
77
    
78
    LOGICAL, PRIVATE :: isParallel=.FALSE.
79

80
CONTAINS
81

82
!------------------------------------------------------------------------------
83
!> Routine checks the feasibility of solver options.
84
!------------------------------------------------------------------------------
85
  SUBROUTINE CheckLinearSolverOptions( Solver ) 
86
!------------------------------------------------------------------------------
87
    TYPE(Solver_t) :: Solver
88
!------------------------------------------------------------------------------
89
    TYPE(ValueList_t), POINTER :: Params
90
    LOGICAL :: Found, Parallel
91
    CHARACTER(:), ALLOCATABLE :: str
92
!------------------------------------------------------------------------------
93

94
    Params => ListGetSolverParams(Solver)
95
    Parallel = Solver % Parallel
96
    
97
    str = ListGetString( Params,'Linear System Solver', Found )
98

99
    IF ( str == 'direct' ) THEN
100
      str = ListGetString( Params,'Linear System Direct Method', Found )
101
      
102
      IF( Found ) THEN        
103
        IF ( Parallel ) THEN
104
          IF ( str /= 'smumps' .AND. str /= 'cmumps' .AND. str /= 'dmumps' .AND. str/= 'zmumps' &
105
                .AND. str /= 'mumps' .AND. str /= 'cpardiso' .AND. str /= 'permon' ) THEN
106
            CALL Warn( 'CheckLinearSolverOptions', 'Only MUMPS and CPardiso direct solver' // &
107
                ' interface implemented in parallel, trying MUMPS!')
108
            str = 'mumps' 
109
            CALL ListAddString( Params,'Linear System Direct Method', str)
110
          END IF
111
        END IF
112
        
113
#if !defined (HAVE_UMFPACK) && defined (HAVE_MUMPS)
114
        IF ( str == 'umfpack' .OR. str == 'big umfpack' ) THEN
115
          CALL Warn( 'CheckLinearSolverOptions', 'UMFPACK solver not installed, using MUMPS instead!' )
116
          str = 'mumps'
117
          CALL ListAddString( Params,'Linear System Direct Method', str)
118
        END IF
119
#endif
120

121
        SELECT CASE( str )
122
        CASE('banded' )
123
        CASE( 'umfpack', 'big umfpack' )
124
#ifndef HAVE_UMFPACK
125
          CALL Fatal( 'CheckLinearSolverOptions', 'UMFPACK (and MUMPS) solver has not been installed.' )
126
#endif
127
        CASE( 'mumps', 'mumpslocal', 'dmumps', 'zmumps', 'smumps', 'cmumps' )
128
#ifndef HAVE_MUMPS
129
          CALL Fatal( 'CheckLinearSolverOptions', 'MUMPS solver has not been installed.' )
130
#endif
131
        CASE( 'superlu' )
132
#ifndef HAVE_SUPERLU
133
          CALL Fatal( 'CheckLinearSolverOptions', 'SuperLU solver has not been installed.' )
134
#endif
135
        CASE( 'pardiso' )
136
#if !defined(HAVE_PARDISO) && !defined(HAVE_MKL)
137
          CALL Fatal( 'CheckLinearSolverOptions', 'Pardiso solver has not been installed.' )
138
#endif
139
        CASE( 'cpardiso')
140
#if !defined(HAVE_CPARDISO) || !defined(HAVE_MKL)
141
        CALL Fatal( 'CheckLinearSolverOptions', ' Cluster Pardiso solver has not been installed.' )
142
#endif
143
        CASE( 'cholmod','spqr' )
144
#ifndef HAVE_CHOLMOD
145
          CALL Fatal( 'CheckLinearSolverOptions', 'Cholmod solver has not been installed.' )
146
#endif
147
       CASE( 'permon')
148
#ifndef HAVE_FETI4I
149
          CALL Fatal( 'CheckLinearSolverOptions', 'FETI4I solver has not been installed.' )
150
#endif
151
        CASE DEFAULT
152
          CALL Fatal( 'CheckLinearSolverOptions', 'Unknown direct solver method: ' // TRIM(str) )
153
        END SELECT
154
        
155
      ELSE
156
        IF ( .NOT. Parallel ) THEN
157
#ifdef HAVE_UMFPACK
158
          str = 'umfpack'
159
#else
160
          str = 'banded'
161
#endif
162
        ELSE
163
#ifdef HAVE_MUMPS
164
          str = 'mumps'
165
#else
166
          CALL Fatal( 'CheckLinearSolverOptions', 'There is no direct parallel solver available (MUMPS)')
167
#endif
168
        END IF
169
        CALL Info('CheckLinearSolverOptions',&
170
            'Setting > Linear System Direct Method < to:'//TRIM(str), Level=6)
171
        CALL ListAddString( Params,'Linear System Direct Method', str )
172
      END IF
173

174
    ELSE IF ( str == 'feti' ) THEN
175
      IF( ParEnv % PEs <= 1 ) THEN
176
        CALL Fatal('CheckLinearSolverOptions','Feti not usable in serial!')
177
      END IF
178

179
    ELSE
180
      IF (ListGetLogical( Params,  &
181
          'Linear System Use Hypre', Found )) THEN
182
#ifndef HAVE_HYPRE
183
        CALL Fatal('CheckLinearSolverOptions','Hypre requested but not compiled with!')
184
#endif
185
      END IF
186

187
      IF (ListGetLogical( Params,  &
188
          'Linear System Use Trilinos', Found )) THEN        
189
        CALL Fatal('CheckLinearSolverOptions','Trilinos implementation was obsolite and has been romoved!')
190
      END IF
191

192
      ! Naming replacement: 'cholesky' for 'symmetric ILU'
193
      str = ListGetString(Params, 'Linear System Preconditioning', Found)
194
      IF ( LEN_TRIM(str) >= 8 ) THEN
195
        IF( str(1:8) == 'cholesky') THEN
196
          CALL ListAddString(Params, 'Linear System Preconditioning', 'ilu'//TRIM(str(9:)) )
197
          CALL ListAddLogical(Params, 'Linear System Symmetric', .TRUE.)
198
          CALL ListAddLogical(Params, 'Linear System Symmetric ILU', .TRUE.)
199
        END IF
200
      END IF
201

202
    END IF
203
  
204
!------------------------------------------------------------------------------
205
  END SUBROUTINE CheckLinearSolverOptions
206
!------------------------------------------------------------------------------
207

208

209
!------------------------------------------------------------------------------
210
! This subroutine enables the scaling of keywords by geometric entities. 
211
! If any Real type keyword has an associated keyword with suffix 'Normalize by Area'
212
! or 'Normalize by Volume', or the real valued keyword has the prefix -dist
213
! the keyword will be divided with the area/volume.
214
!------------------------------------------------------------------------------
215
   SUBROUTINE SetNormalizedKeywords(Model, Mesh)
216
     TYPE(Model_t) :: Model
217
     TYPE(Mesh_t), POINTER :: Mesh
218

219
     INTEGER :: cnt
220
     LOGICAL :: Found
221

222
     cnt = Model % NumberOfDistTags 
223

224
     ! Initialize the count if not done before
225
     IF( cnt == -1 ) THEN
226
       cnt = ListTagCount(Model,.TRUE.)
227
       Model % NumberOfDistTags = cnt
228
     END IF
229
     
230
     ! If no tags nothing to do
231
     IF( cnt == 0 ) RETURN
232

233
     ! If the weights have been already computed for linear cases no need to redo
234
     IF( Mesh % EntityWeightsComputed ) THEN
235
       IF( .NOT. ListGetLogical( Model % Simulation,&
236
           'Update Keyword Normalization',Found ) ) RETURN
237
     END IF
238
     
239
     ! Calculate entity weights
240
     CALL CalculateEntityWeights( Model, Mesh ) 
241

242
     ! Tag count for entity normalization tags, 2nd and 3rd parameter neglected
243
     ! for this USE CASE!
244
     CALL ListSetParameters( Model, 0, 0.0_dp, .FALSE., Found ) 
245

246
   END SUBROUTINE SetNormalizedKeywords
247
!------------------------------------------------------------------------------
248

249

250
!------------------------------------------------------------------------------
251
! This subroutine enables the rotation of vector valued keywords.
252
!------------------------------------------------------------------------------
253
   SUBROUTINE SetRotatedProperties(Model, Mesh)
254
     TYPE(Model_t) :: Model
255
     TYPE(Mesh_t), POINTER :: Mesh
256

257
     LOGICAL :: AnyBC, AnyBodyForce, AnyMat, AnyBody
258
     INTEGER :: list_ind
259
     REAL(KIND=dp) :: RotMatrix(3,3), Angles(3,1)
260
     REAL(KIND=dp), POINTER :: Pmatrix(:,:)
261
     INTEGER, TARGET :: NodeIndex(1)
262
     INTEGER, POINTER :: NodeIndexes(:)
263
     REAL(KIND=dp) :: EntityWeight
264
     LOGICAL :: Found
265
     TYPE(ValueList_t), POINTER :: List
266
     LOGICAL :: Visited = .FALSE.
267
     CHARACTER(:), ALLOCATABLE :: Suffix
268

269
     SAVE Visited, Suffix, AnyBC, AnyBodyForce, AnyMat, AnyBody
270

271
     IF(.NOT. Visited ) THEN
272
       Suffix = 'Property Rotate'
273
       AnyBC = ListCheckSuffixAnyBC( Model, Suffix )
274
       AnyBodyForce = ListCheckSuffixAnyBodyForce( Model, Suffix )
275
       AnyMat = ListCheckSuffixAnyMaterial( Model, Suffix )
276
       AnyBody = ListCheckSuffixAnyBody( Model, Suffix )
277
       Visited = .TRUE.
278
     END IF
279
       
280
     IF( .NOT. ( AnyBC .OR. AnyBodyForce .OR. AnyMat .OR. AnyBody ) ) RETURN
281

282
     NodeIndex(1) = 1
283
     NodeIndexes => NodeIndex
284

285
     IF( AnyBody ) THEN
286
       DO list_ind = 1,Model % NumberOfBodies
287
         List => Model % Bodies(list_ind) % Values                  
288
         IF(.NOT. ListCheckSuffix( List, Suffix ) ) CYCLE
289

290
         CALL ListGetRealVector( List,'Property Rotate',Angles,1,NodeIndexes,Found )
291
         IF(.NOT. Found ) THEN
292
           CALL Fatal('SetRotatedProperties','> Property Rotate < suffix requires angles!')
293
         END IF
294
         
295
         RotMatrix = AnglesToRotationMatrix( Angles(1:3,1) )
296

297
         PMatrix => ListGetConstRealArray( List,'Property Rotation Matrix',Found )
298
         IF( Found ) THEN
299
           PMatrix = RotMatrix(3,3) 
300
         ELSE
301
           CALL ListAddConstRealArray( List,'Property Rotation Matrix', 3, 3, RotMatrix )
302
         END IF
303
       END DO
304
     END IF
305

306

307
     IF( AnyBodyForce ) THEN
308
       DO list_ind = 1,Model % NumberOfBodyForces
309
         List => Model % BodyForces(list_ind) % Values                  
310
         IF(.NOT. ListCheckSuffix( List, Suffix ) ) CYCLE
311

312
         CALL ListGetRealVector( List,'Property Rotate',Angles,1,NodeIndexes,Found )
313
         IF(.NOT. Found ) THEN
314
           CALL Fatal('SetRotatedProperties','> Property Rotate < suffix requires angles!')
315
         END IF
316
         
317
         RotMatrix = AnglesToRotationMatrix( Angles(1:3,1) )
318

319
         PMatrix => ListGetConstRealArray( List,'Property Rotation Matrix',Found )
320
         IF( Found ) THEN
321
           PMatrix = RotMatrix(3,3) 
322
         ELSE
323
           CALL ListAddConstRealArray( List,'Property Rotation Matrix', 3, 3, RotMatrix )
324
         END IF
325
       END DO
326
     END IF
327

328

329
     IF( AnyMat ) THEN
330
       DO list_ind = 1,Model % NumberOfMaterials
331
         List => Model % Materials(list_ind) % Values                  
332
         IF(.NOT. ListCheckSuffix( List, Suffix ) ) CYCLE
333

334
         CALL ListGetRealVector( List,'Property Rotate',Angles,1,NodeIndexes,Found )
335
         IF(.NOT. Found ) THEN
336
           CALL Fatal('SetRotatedProperties','> Property Rotate < suffix requires angles!')
337
         END IF
338
         
339
         RotMatrix = AnglesToRotationMatrix( Angles(1:3,1) )
340

341
         PMatrix => ListGetConstRealArray( List,'Property Rotation Matrix',Found )
342
         IF( Found ) THEN
343
           PMatrix = RotMatrix(3,3) 
344
         ELSE
345
           CALL ListAddConstRealArray( List,'Property Rotation Matrix', 3, 3, RotMatrix )
346
         END IF
347
       END DO
348
     END IF
349

350

351
     IF( AnyBC ) THEN
352
       DO list_ind = 1,Model % NumberOfBCs
353
         List => Model % BCs(list_ind) % Values                  
354
         IF(.NOT. ListCheckSuffix( List, Suffix ) ) CYCLE
355

356
         CALL ListGetRealVector( List,'Property Rotate',Angles,1,NodeIndexes,Found )
357
         IF(.NOT. Found ) THEN
358
           CALL Fatal('SetRotatedProperties','> Property Rotate < suffix requires angles!')
359
         END IF
360
         
361
         RotMatrix = AnglesToRotationMatrix( Angles(1:3,1) )
362

363
         PMatrix => ListGetConstRealArray( List,'Property Rotation Matrix',Found )
364
         IF( Found ) THEN
365
           PMatrix = RotMatrix(3,3) 
366
         ELSE
367
           CALL ListAddConstRealArray( List,'Property Rotation Matrix', 3, 3, RotMatrix )
368
         END IF
369
       END DO
370
     END IF
371

372
   CONTAINS 
373

374
     FUNCTION AnglesToRotationMatrix( Angles, RotateOrder ) RESULT ( RotMatrix )
375

376
       REAL(KIND=dp) :: Angles(3), RotMatrix(3,3)
377
       INTEGER, OPTIONAL :: RotateOrder(3)
378

379
       INTEGER :: i,j
380
       REAL(KIND=dp) :: TrfMatrix(3,3), IdentityMatrix(3,3), Alpha
381

382
       IdentityMatrix = 0.0_dp
383
       DO i=1,3
384
         IdentityMatrix(i,i) = 1.0_dp
385
       END DO
386

387
       RotMatrix = IdentityMatrix
388

389
       DO i=1,3
390
         j = i
391
         IF( PRESENT( RotateOrder ) ) j = RotateOrder(i)
392
         Alpha = Angles(j) * PI / 180.0_dp
393

394
         IF( ABS(Alpha) < TINY(Alpha) ) CYCLE
395
         TrfMatrix = IdentityMatrix
396

397
         SELECT CASE(j)
398
         CASE(1)
399
           TrfMatrix(2,2) =  COS(Alpha)
400
           TrfMatrix(2,3) = -SIN(Alpha)
401
           TrfMatrix(3,2) =  SIN(Alpha)
402
           TrfMatrix(3,3) =  COS(Alpha)
403
         CASE(2)
404
           TrfMatrix(1,1) =  COS(Alpha)
405
           TrfMatrix(1,3) = -SIN(Alpha)
406
           TrfMatrix(3,1) =  SIN(Alpha)
407
           TrfMatrix(3,3) =  COS(Alpha)
408
         CASE(3)
409
           TrfMatrix(1,1) =  COS(Alpha)
410
           TrfMatrix(1,2) = -SIN(Alpha)
411
           TrfMatrix(2,1) =  SIN(Alpha)
412
           TrfMatrix(2,2) =  COS(Alpha)
413
         END SELECT
414

415
         RotMatrix = MATMUL( RotMatrix, TrfMatrix )
416
       END DO
417

418
     END FUNCTION AnglesToRotationMatrix
419

420

421
   END SUBROUTINE SetRotatedProperties
422
!------------------------------------------------------------------------------
423

424

425

426

427
!------------------------------------------------------------------------------
428
!> Get calling address of the procedure and add it to the Solver structure.
429
!------------------------------------------------------------------------------
430
  SUBROUTINE AddSolverProcedure( Solver,PROCEDURE  )
431
!------------------------------------------------------------------------------
432
    TYPE(Solver_t) :: Solver
433
    EXTERNAL :: PROCEDURE
434
    INTEGER  :: PROCEDURE
435
!------------------------------------------------------------------------------
436
    INTEGER(KIND=AddrInt) :: AddrFunc
437
    EXTERNAL :: AddrFunc
438
!------------------------------------------------------------------------------
439
    Solver % PROCEDURE = AddrFunc( PROCEDURE )
440
!------------------------------------------------------------------------------
441
  END SUBROUTINE AddSolverProcedure
442
!------------------------------------------------------------------------------
443

444

445
 
446
!------------------------------------------------------------------------------
447
   SUBROUTINE SwapMesh(Model,Mesh,Name,ExtMesh)
448
!------------------------------------------------------------------------------
449
     TYPE(Model_t) :: Model
450
     TYPE(Mesh_t), POINTER :: Mesh
451
     CHARACTER(LEN=*), OPTIONAL :: Name     
452
     TYPE(Mesh_t), POINTER, OPTIONAL :: ExtMesh
453
!------------------------------------------------------------------------------
454
     TYPE(Mesh_t), POINTER :: Newmesh, tmesh
455
     INTEGER :: Def_Dofs(10,6), i,j,k
456
     LOGICAL :: Found, Transient
457
     TYPE(Solver_t), POINTER :: Solver
458
!------------------------------------------------------------------------------
459

460
  INTERFACE
461
    SUBROUTINE InterpolateMeshToMesh( OldMesh, NewMesh, OldVariables, &
462
        NewVariables, UseQuadrantTree, Projector, MaskName, UnfoundNodes )
463
      USE Lists
464
      USE SParIterComm
465
      USE Interpolation
466
      USE CoordinateSystems
467
      USE MeshUtils, ONLY: ReleaseMesh
468
      TYPE(Mesh_t), TARGET  :: OldMesh, NewMesh
469
      TYPE(Variable_t), POINTER, OPTIONAL :: OldVariables, NewVariables
470
      LOGICAL, OPTIONAL :: UseQuadrantTree
471
      TYPE(Projector_t), POINTER, OPTIONAL :: Projector
472
      CHARACTER(LEN=*),OPTIONAL :: MaskName
473
      LOGICAL, POINTER, OPTIONAL :: UnfoundNodes(:)
474
    END SUBROUTINE InterpolateMeshToMesh
475
  END INTERFACE
476

477

478
     Def_Dofs = -1;
479
     DO i=1,Model % NumberOfSolvers
480
       DO j=1,10
481
         DO k=1,6
482
           Def_Dofs(j,k) = MAX(Def_Dofs(j,k), MAXVAL(Model % Solvers(i) % Def_Dofs(j,:,k)))
483
         END DO
484
       END DO
485
     END DO
486

487
     IF( PRESENT( ExtMesh ) ) THEN
488
       NewMesh => ExtMesh
489
     ELSE              
490
       Newmesh => LoadMesh2( Model, Name, Name, &
491
           .FALSE., Parenv % PEs, ParEnv % myPE, Def_Dofs )
492
       NewMesh % Name = Name
493
     END IF
494

495
     IF(.NOT.ASSOCIATED(NewMesh)) RETURN
496

497
     NewMesh % Next => Mesh % Next
498
     IF(ASSOCIATED(Mesh,  Model % Meshes)) THEN
499
       Model % Meshes => Newmesh
500
     ELSE
501
       Tmesh => Model % Meshes
502
       DO WHILE(ASSOCIATED(Tmesh % next))
503
         IF(ASSOCIATED(Mesh, Tmesh % next)) THEN
504
           Tmesh % Next => Newmesh
505
           EXIT
506
         END IF
507
       END DO
508
     END IF
509

510
     CALL TransferCoordAndTime(Mesh,NewMesh)
511

512
     DO i=1,Model % NumberOfSolvers
513
       Solver => Model % Solvers(i)
514
       IF(ASSOCIATED(Solver % Mesh, Mesh)) Solver % Mesh => Newmesh
515
     END DO
516

517
     IF(Mesh % DiscontMesh) CALL CreateDiscontMesh(Model,Newmesh,.TRUE.)
518

519
     IF(ASSOCIATED(Model % Mesh, Mesh)) Model % Mesh => NewMesh
520
     IF(ASSOCIATED(Model % Variables, Mesh % Variables)) Model % Variables => NewMesh % Variables
521

522
     Mesh % Next => Null()
523

524
     Transient = ListGetString( Model % Simulation, 'Simulation Type' ) == 'transient'
525

526
     DO i=1,Model % NumberOfSolvers
527
       Solver => Model % Solvers(i)
528
       IF(ASSOCIATED(Solver % Mesh, NewMesh)) THEN
529
         CALL FreeMatrix(Solver % Matrix)
530
         Model % Solver => Solver
531

532
         CALL AddEquationBasics( Solver, ListGetString(Solver % Values, &
533
                  'Variable', Found), Transient )
534
         CALL AddEquationSolution( Solver, Transient )
535
         IF ( Transient .AND. Solver % PROCEDURE /= 0 ) CALL InitializeTimestep(Solver)
536
       END IF
537
     END DO
538

539
     IF( Transient ) THEN
540
      IF(ParEnv % PEs > 1) CALL ParallelActive(.TRUE.)
541
      !Free quadrant tree to ensure its rebuilt in InterpolateMeshToMesh (bug fix)
542
      CALL FreeQuadrantTree(Mesh % RootQuadrant)
543
      CALL InterpolateMeshToMesh( Mesh, NewMesh, Mesh % Variables, NewMesh % Variables )
544
     END IF
545

546

547
     CALL ReleaseMesh( Mesh )
548

549
     CALL MeshStabParams( Newmesh )
550
     NewMesh % Changed = .TRUE.
551
    
552
!------------------------------------------------------------------------------
553
   END SUBROUTINE SwapMesh
554
!------------------------------------------------------------------------------
555

556
           
557
   SUBROUTINE CheckAndCreateDGIndexes( Mesh, ActiveElem ) 
558
     TYPE(Mesh_t), POINTER :: Mesh
559
     LOGICAL, OPTIONAL :: ActiveElem(:)
560

561
     TYPE(Element_t), POINTER :: Element     
562
     INTEGER :: i,n,t,DGIndex
563
     LOGICAL :: Failed
564

565
     Failed = .FALSE.
566
     
567
     DO t=1,Mesh % NumberOfBulkElements 
568
       IF( PRESENT( ActiveElem ) ) THEN
569
         IF( .NOT. ActiveElem(t) ) CYCLE
570
       END IF
571
       Element => Mesh % Elements(t)
572
       n = Element % TYPE % NumberOfNodes         
573
       IF( .NOT. ASSOCIATED( Element % DGIndexes ) ) THEN
574
         Failed = .TRUE.
575
       ELSE IF( SIZE( Element % DGIndexes ) /= n ) THEN
576
         Failed = .TRUE.
577
       END IF
578
     END DO
579
     IF(.NOT. Failed) RETURN
580
     
581
     ! Number the bulk indexes such that each node gets a new index
582
     CALL Info('CheckAndCreateDGIndexes','Creating DG indexes!',Level=12)
583

584
     DGIndex = 0       
585
     DO t=1,Mesh % NumberOfBulkElements 
586
       Element => Mesh % Elements(t)
587
       n = Element % TYPE % NumberOfNodes         
588
       IF( .NOT. ASSOCIATED( Element % DGIndexes ) ) THEN
589
         ALLOCATE( Element % DGindexes( n ) )
590
       ELSE IF( SIZE( Element % DGIndexes ) /= n ) THEN
591
         DEALLOCATE( Element % DGindexes )
592
         ALLOCATE( Element % DGindexes( n ) )
593
       END IF
594
         
595
       DO i=1,n
596
         DGIndex = DGIndex + 1
597
         Element % DGIndexes(i) = DGIndex
598
       END DO
599
     END DO
600

601
     ! Number the bulk indexes such that each node gets a new index
602
     CALL Info('CheckAndCreateDGIndexes','Creating DG '//I2S(DgIndex)//' indexes',Level=6)
603
    
604
   END SUBROUTINE CheckAndCreateDGIndexes
605
     
606
   
607

608
   !> Create permutation for discontinuous Galerkin type of fields optionally
609
   !> with a given mask.
610
   !-----------------------------------------------------------------------------------
611
   SUBROUTINE CreateDGPerm( Solver, DGPerm, DGCount, MaskName, SecName )
612

613
     TYPE(Solver_t), POINTER :: Solver
614
     INTEGER, POINTER :: DGPerm(:)
615
     INTEGER :: DGCount
616
     CHARACTER(LEN=*), OPTIONAL :: MaskName, SecName 
617
     
618
     TYPE(Mesh_t), POINTER :: Mesh
619
     TYPE(Element_t), POINTER :: Element, Parent
620
     INTEGER :: t, n, i, j, k, DGIndex
621
     LOGICAL :: Found, ActiveElem, HaveSome
622
     TYPE(ValueList_t), POINTER :: BF
623
     CHARACTER(:), ALLOCATABLE :: EquationName
624
     
625

626
     CALL Info('CreateDGPerm','Creating permutation for DG variable',Level=12)
627
     
628
     EquationName = ListGetString( Solver % Values, 'Equation', Found)
629
     IF( .NOT. Found ) THEN
630
       CALL Fatal('CreateDGPerm','Equation not present!')
631
     END IF     
632
     
633
     Mesh => Solver % Mesh
634
     DGIndex = 0
635

636
     ! Check whether the DG indexes might already have been allocated
637
     HaveSome = .FALSE.
638
     DO t=1,Mesh % NumberOfBulkElements 
639
       Element => Mesh % Elements(t)
640
       IF( ASSOCIATED( Element % DGIndexes ) ) THEN
641
         HaveSome = .TRUE.
642
         EXIT
643
       END IF
644
     END DO
645

646
     ! If they are allocated somewhere check that they are good for this variable as well 
647
     IF( HaveSome ) THEN
648
       CALL Info('CreateDGPerm','There are at least some associated DG elements',Level=15)
649
       DO t=1,Mesh % NumberOfBulkElements + Mesh % NumberOFBoundaryElements
650
         Element => Mesh % Elements(t)         
651
         IF( Element % PartIndex == ParEnv % myPE ) THEN
652
           IF ( CheckElementEquation( CurrentModel, Element, EquationName ) ) THEN             
653

654
             IF( PRESENT( MaskName ) ) THEN
655
               BF => ListGetSection( Element, SecName )
656
               ActiveElem = ListGetLogicalGen( BF, MaskName )
657
             ELSE
658
               ActiveElem = .TRUE.
659
             END IF
660

661
             IF( ActiveElem ) THEN
662
               IF( .NOT. ASSOCIATED( Element % DGIndexes ) ) THEN
663
                 CALL Fatal('CreateDGPerm','Either all or none of DGIndexes should preexist!')               
664
               END IF
665
             END IF
666
           END IF
667
         END IF
668
       END DO
669

670
       ! Find out the maximum index for the size of the permutation
671
       DO t=1,Mesh % NumberOfBulkElements 
672
         Element => Mesh % Elements(t)
673
         DGIndex = MAX( DGIndex, MAXVAL( Element % DGIndexes ) )
674
       END DO
675
     END IF
676

677

678
     ! If they have not been allocated before the allocate DGIndexes for all bulk elements
679
     IF(.NOT. HaveSome ) THEN       
680
       CALL Info('CreateDGPerm','Creating DG indexes for bulk elements',Level=15)
681

682
       ! Number the bulk indexes such that each node gets a new index
683
       DGIndex = 0       
684
       DO t=1,Mesh % NumberOfBulkElements !+ Mesh % NumberOfBoundaryElements
685
         Element => Mesh % Elements(t)
686
         n = Element % TYPE % NumberOfNodes         
687
         ALLOCATE( Element % DGindexes( n ) )
688
         DO i=1, n
689
           DGIndex = DGIndex + 1
690
           Element % DGIndexes(i) = DGIndex
691
         END DO
692
       END DO
693
       
694
       ! Make boundary elements to inherit the bulk indexes
695
       ! We neglect this as for now since it seems this causes problems in deallocation later on...
696
#if 1
697
       DO t=Mesh % NumberOfBulkElements + 1, &
698
           Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
699
         Element => Mesh % Elements(t)
700
         n = Element % TYPE % NumberOfNodes
701
               
702
         IF( .NOT. ASSOCIATED( Element % BoundaryInfo ) ) CYCLE
703
         DO k=1,2
704
           IF( k == 1 ) THEN
705
             Parent => Element % BoundaryInfo % Left
706
           ELSE 
707
             Parent => Element % BoundaryInfo % Right
708
           END IF
709
           IF(.NOT. ASSOCIATED( Parent ) ) CYCLE           
710
           IF( .NOT. ASSOCIATED( Parent % DGIndexes ) ) CYCLE
711

712
           IF( .NOT. ASSOCIATED( Element % DGIndexes ) ) THEN
713
             ALLOCATE( Element % DGIndexes(n) ) 
714
             Element % DgIndexes = 0
715
           END IF
716
           DO i = 1, n
717
             IF( Element % DGIndexes(i) > 0 ) CYCLE
718
             DO j = 1, Parent % TYPE % NumberOfNodes
719
               IF( Element % NodeIndexes(i) == Parent % NodeIndexes(j) ) THEN
720
                 Element % DGIndexes(i) = Parent % DGIndexes(j)
721
                 EXIT
722
               END IF
723
             END DO
724
           END DO
725
           IF( ANY( Element % DGIndexes == 0 ) ) PRINT *,'t dg:',t,n,Element % DGIndexes
726
           EXIT
727
         END DO
728
       END DO
729
#endif
730
     END IF
731
     
732
     CALL Info('CreateDGPerm','Size of DgPerm table: '//I2S(DGIndex),Level=12)
733
     
734
     ALLOCATE( DGPerm( DGIndex ) ) 
735
     DGPerm = 0
736

737
     ! If we consider boundary element above do that also here
738
     DO t=1,Mesh % NumberOfBulkElements + Mesh % NumberOFBoundaryElements
739
       Element => Mesh % Elements(t)
740
            
741
       IF( Element % PartIndex == ParEnv % myPE ) THEN
742
         IF ( CheckElementEquation( CurrentModel, Element, EquationName ) ) THEN
743
           
744
           IF( PRESENT( MaskName ) ) THEN
745
             BF => ListGetSection( Element, SecName )
746
             ActiveElem = ListGetLogicalGen( BF, MaskName )
747
           ELSE
748
             ActiveElem = .TRUE.
749
           END IF
750
           
751
           IF( ActiveElem ) THEN
752
             n = Element % TYPE % NumberOfNodes
753
             DGPerm( Element % DGIndexes ) = 1 
754
           END IF
755
           
756
         END IF
757
       END IF
758
     END DO
759

760
     DGCount = 0
761
     DO i=1, DGIndex
762
       IF( DGPerm(i) > 0 ) THEN
763
         DGCount = DGCount + 1
764
         DGPerm(i) = DGCount
765
       END IF
766
     END DO
767
     
768
     CALL Info('CreateDGPerm','Created permutation for DG nodes: '//I2S(DgCount),Level=8)  
769
     
770
   END SUBROUTINE CreateDGPerm
771

772

773
   !> Simple nodal permutation without any mask or other complications.
774
   !> For the more complex version there is MakePermUsingMask.
775
   !> This simple one could be used when primary variable is DG field and we want
776
   !> to create a nodal exported variable.
777
   !-----------------------------------------------------------------------------------
778
   SUBROUTINE CreateNodalPerm( Solver, NodalPerm, nSize )
779

780
     TYPE(Solver_t), POINTER :: Solver
781
     INTEGER, POINTER :: NodalPerm(:)
782
     INTEGER :: nSize
783
     
784
     TYPE(Mesh_t), POINTER :: Mesh
785
     TYPE(Element_t), POINTER :: Element
786
     INTEGER :: t, n, i, j
787
     LOGICAL :: Found
788
     CHARACTER(:), ALLOCATABLE :: EquationName
789
     
790
     CALL Info('CreateNodalPerm','Creating simple permutation for nodal variable',Level=12)
791
     
792
     EquationName = ListGetString( Solver % Values, 'Equation', Found)
793
     IF( .NOT. Found ) THEN
794
       CALL Fatal('CreateNodalPerm','Equation not present!')
795
     END IF     
796
     
797
     Mesh => Solver % Mesh
798
     n = Mesh % NumberOfNodes 
799

800
     ALLOCATE( NodalPerm(n) ) 
801
     NodalPerm = 0 
802
     
803
     DO t=1,Mesh % NumberOfBulkElements + Mesh % NumberOFBoundaryElements
804
       Element => Mesh % Elements(t)         
805
       IF ( CheckElementEquation( CurrentModel, Element, EquationName ) ) THEN             
806
         NodalPerm( Element % NodeIndexes ) = 1
807
       END IF
808
     END DO
809

810
     j = 0
811
     DO i=1,n
812
       IF( NodalPerm(i) > 0 ) THEN
813
         j = j + 1
814
         NodalPerm(i) = j
815
       END IF
816
     END DO
817
     nSize = j
818
     
819
     CALL Info('CreateNodalPerm','Number of active nodes in NodalPerm: '//I2S(nSize),Level=12)
820
     
821
   END SUBROUTINE CreateNodalPerm
822

823

824
   
825
   !> Create permutation for fields on elements, optional using mask
826
   !-----------------------------------------------------------------
827
   SUBROUTINE CreateElementsPerm( Solver, Perm, nsize, MaskName, SecName ) 
828
     TYPE(Solver_t),POINTER :: Solver
829
     INTEGER, POINTER :: Perm(:)
830
     INTEGER :: nsize
831
     CHARACTER(LEN=*), OPTIONAL :: MaskName, SecName
832

833
     INTEGER :: t, n, m
834
     TYPE(Element_t), POINTER :: Element
835
     LOGICAL :: Found, ActiveElem
836
     TYPE(Mesh_t), POINTER :: Mesh
837
     TYPE(ValueList_t), POINTER :: BF
838
     CHARACTER(:), ALLOCATABLE :: EquationName
839

840
     CALL Info('CreateElementsPerm','Creating permutation for elemental fields',Level=8)
841

842
     EquationName = ListGetString( Solver % Values, 'Equation', Found)
843
     IF( .NOT. Found ) THEN
844
       CALL Fatal('CreateElementsPerm','Equation not present!')
845
     END IF
846

847
     Mesh => Solver % Mesh
848

849

850
     NULLIFY( Perm ) 
851
     n = Mesh % NumberOfBulkElements + Mesh % NumberOFBoundaryElements
852
     ALLOCATE( Perm(n) )
853
     Perm = 0
854

855

856
     m = 0
857
     DO t=1,n
858
       Element => Solver % Mesh % Elements(t)
859
       IF( Element % PartIndex == ParEnv % myPE ) THEN
860
         IF ( CheckElementEquation( CurrentModel, Element, EquationName ) ) THEN
861
           IF( PRESENT( MaskName ) ) THEN
862
             BF => ListGetSection( Element, SecName )
863
             ActiveElem = ListGetLogicalGen( BF, MaskName )
864
           ELSE
865
             ActiveElem = .TRUE.
866
           END IF
867
         END IF
868

869
         IF( ActiveElem ) THEN
870
           m = m + 1
871
           Perm(t) = m
872
         END IF
873
       END IF
874
     END DO
875

876
     CALL Info('CreateElementsPerm','Number of active elements in permutation: '//I2S(m),Level=8)
877

878
     nsize = m 
879

880
   END SUBROUTINE CreateElementsPerm
881

882

883
   
884
   !> Create permutation table that follows the degrees of freedom of the primary
885
   !> permutation but is masked by a flag on the body force section.
886
   !---------------------------------------------------------------------------------
887
   SUBROUTINE CreateMaskedPerm( Solver, FullPerm, MaskName, MaskPerm, nsize, SecName )
888

889
     TYPE(Solver_t), POINTER :: Solver
890
     INTEGER, POINTER :: FullPerm(:)
891
     CHARACTER(LEN=*) :: MaskName
892
     INTEGER, POINTER :: MaskPerm(:) 
893
     INTEGER :: nsize
894
     CHARACTER(LEN=*), OPTIONAL :: SecName
895
         
896
     TYPE(Mesh_t), POINTER :: Mesh
897
     TYPE(Element_t), POINTER :: Element
898
     INTEGER :: t, n, m
899
     TYPE(ValueList_t), POINTER :: BF
900
     LOGICAL :: Found, ActiveElem
901
     INTEGER, ALLOCATABLE :: Indexes(:)
902
     CHARACTER(:), ALLOCATABLE :: EquationName
903
     
904
     
905
     CALL Info('CreateMaskedPerm','Creating variable with mask: '//TRIM(MaskName),Level=8)
906
            
907
     Mesh => Solver % Mesh
908
     NULLIFY( MaskPerm ) 
909

910
     n = SIZE( FullPerm ) 
911
     ALLOCATE( MaskPerm( n ), Indexes(100) ) 
912

913
     MaskPerm = 0     
914

915
     DO t=1,Mesh % NumberOfBulkElements + Mesh % NumberOFBoundaryElements
916
       Element => Mesh % Elements(t)
917

918
       IF( ParEnv % PEs > 1 ) THEN
919
         IF( t <= Mesh % NumberOfBulkElements ) THEN
920
           IF( Element % PartIndex /= ParEnv % myPE ) CYCLE
921
         END IF
922
       END IF
923
       
924
       IF( PRESENT( SecName ) ) THEN
925
         BF => ListGetSection( Element, SecName ) 
926
       ELSE
927
         IF( t <= Mesh % NumberOfBulkElements ) THEN
928
           BF => ListGetSection( Element,'body force')
929
         ELSE
930
           BF => ListGetSection( Element,'boundary condition')
931
         END IF
932
       END IF
933

934
       
935
       IF(.NOT. ASSOCIATED( BF ) ) CYCLE
936

937
       ActiveElem = ListGetLogicalGen( BF, MaskName )
938

939
       IF( ActiveElem ) THEN
940
         MaskPerm( Element % NodeIndexes ) = FullPerm( Element % NodeIndexes )
941
       END IF
942

943
     END DO
944
       
945
     m = 0
946
     DO t=1,n
947
       IF( MaskPerm( t ) > 0 ) THEN
948
         m = m + 1
949
         MaskPerm( t ) = m
950
       END IF
951
     END DO
952
            
953
     nsize = m
954
     
955
     CALL Info('CreateMaskedPerm','Created masked permutation for dofs: '//I2S(nsize),Level=8)  
956
     
957
   END SUBROUTINE CreateMaskedPerm
958

959
   !------------------------------------------------------------------
960
   ! Check for special solvers, to be executed only 
961
   ! at a certain instances during the simulation:
962
   !------------------------------------------------------------------  
963
   SUBROUTINE AddExecWhenFlag(Solver)
964
     TYPE(Solver_t), POINTER :: Solver
965

966
     TYPE(ValueList_t), POINTER :: SolverParams
967
     LOGICAL :: Found
968
     CHARACTER(:), ALLOCATABLE :: str
969
     
970
     SolverParams => ListGetSolverParams(Solver)
971

972
     ! Default value     
973
     Solver % SolverExecWhen = SOLVER_EXEC_ALWAYS
974

975
     str = ListGetString( SolverParams, 'Exec Solver', Found )
976

977
     IF( Found ) THEN    
978
       SELECT CASE( TRIM(str) )
979
       CASE( 'never' )
980
         Solver % SolverExecWhen = SOLVER_EXEC_NEVER
981
       CASE( 'always' )
982
         Solver % SolverExecWhen = SOLVER_EXEC_ALWAYS
983
       CASE( 'after simulation', 'after all' )
984
         Solver % SolverExecWhen = SOLVER_EXEC_AFTER_ALL
985
       CASE( 'before simulation', 'before all' )
986
         Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_ALL
987
       CASE( 'before timestep' )
988
         Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_TIME
989
       CASE( 'after timestep' )
990
         Solver % SolverExecWhen = SOLVER_EXEC_AFTER_TIME
991
       CASE( 'before saving' )
992
         Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_SAVE
993
       CASE( 'after saving' )
994
         Solver % SolverExecWhen = SOLVER_EXEC_AFTER_SAVE
995
       CASE( 'predictor-corrector' )
996
         Solver % SolverExecWhen = SOLVER_EXEC_PREDCORR
997
       CASE( 'when created' )
998
         Solver % SolverExecWhen = SOLVER_EXEC_WHENCREATED
999
       CASE( 'after control' )
1000
         Solver % SolverExecWhen = SOLVER_EXEC_AFTER_CONTROL
1001
       CASE DEFAULT
1002
         CALL Fatal('AddExecWhenFlag','Unknown "exec solver" flag: '//TRIM(str))
1003
       END SELECT
1004
     ELSE      
1005
       IF ( ListGetLogical( SolverParams, 'Before All', Found ) ) THEN
1006
         Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_ALL
1007
       ELSE IF ( ListGetLogical( SolverParams, 'Before Simulation', Found ) ) THEN
1008
         Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_ALL
1009
       ELSE IF ( ListGetLogical( SolverParams, 'After All', Found ) ) THEN
1010
         Solver % SolverExecWhen = SOLVER_EXEC_AFTER_ALL
1011
       ELSE IF ( ListGetLogical( SolverParams, 'After Simulation', Found ) ) THEN
1012
         Solver % SolverExecWhen = SOLVER_EXEC_AFTER_ALL
1013
       ELSE IF ( ListGetLogical( SolverParams, 'Before Timestep', Found ) ) THEN
1014
         Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_TIME
1015
       ELSE IF ( ListGetLogical( SolverParams, 'After Timestep', Found ) ) THEN
1016
         Solver % SolverExecWhen = SOLVER_EXEC_AFTER_TIME
1017
       ELSE IF ( ListGetLogical( SolverParams, 'Before Saving', Found ) ) THEN
1018
         Solver % SolverExecWhen = SOLVER_EXEC_AHEAD_SAVE
1019
       ELSE IF ( ListGetLogical( SolverParams, 'After Saving', Found ) ) THEN
1020
         Solver % SolverExecWhen = SOLVER_EXEC_AFTER_SAVE
1021
       ELSE IF ( ListGetLogical( SolverParams, 'Predictor-Corrector', Found ) ) THEN
1022
         Solver % SolverExecWhen = SOLVER_EXEC_PREDCORR
1023
       ELSE IF ( ListGetLogical( SolverParams, 'When Created', Found ) ) THEN
1024
         Solver % SolverExecWhen = SOLVER_EXEC_WHENCREATED
1025
       ELSE IF ( ListGetLogical( SolverParams, 'After Control', Found ) ) THEN
1026
         Solver % SolverExecWhen = SOLVER_EXEC_AFTER_CONTROL
1027
       END IF
1028
     END IF
1029

1030

1031
   END SUBROUTINE AddExecWhenFlag
1032
     
1033

1034
   SUBROUTINE PrintProblemSize( Solver )
1035
     TYPE(Solver_t) :: Solver
1036

1037
     TYPE(Matrix_t), POINTER :: Matrix    
1038
     INTEGER :: i, nz, nx, no, ns
1039
     INTEGER :: nzpar(0:2), nxpar(0:2), nopar(0:2), nspar(0:2)
1040
     LOGICAL :: Parallel, HaveParInfo
1041
     CHARACTER(*), PARAMETER :: Caller="PrintProblemSize"
1042

1043
     Matrix => Solver % Matrix
1044
     IF(.NOT. ASSOCIATED(Matrix)) RETURN
1045
     
1046
     nx = Matrix % NumberOfRows
1047
     nz = SIZE(Matrix % Values)
1048

1049

1050
     Parallel = ParEnv % PEs > 1 .AND. .NOT. Solver % Mesh % SingleMesh 
1051
       
1052
     ! Is there any size info to print?
1053
     i = ParallelReduction( nx )
1054
     IF(i==0) THEN
1055
       CONTINUE
1056

1057
     ELSE IF( Parallel ) THEN
1058
       no = 0; ns = 0
1059

1060
       IF( Parallel ) HaveParInfo = ASSOCIATED(Matrix % ParallelInfo)
1061
       IF(HaveParInfo) THEN
1062
         HaveParInfo = ASSOCIATED(Matrix % ParallelInfo % NeighbourList)
1063
       END IF
1064

1065
       IF(HaveParInfo) THEN
1066
         DO i=1,nx
1067
           IF( Matrix % ParallelInfo % NeighbourList(i) % Neighbours(1) == ParEnv % MyPe) no=no+1
1068
           IF( SIZE(Matrix % ParallelInfo % NeighbourList(i) % Neighbours) > 1) ns=ns+1
1069
         END DO
1070
       END IF
1071

1072
       DO i=0,2
1073
         nxpar(i) = ParallelReduction(nx,i)
1074
         nzpar(i) = ParallelReduction(nz,i)
1075
         IF(HaveParInfo) THEN
1076
           nopar(i) = ParallelReduction(no,i)
1077
           nspar(i) = ParallelReduction(ns,i)
1078
         END IF
1079
       END DO
1080

1081
       CALL Info(Caller,'Parallel size for Solver:           SUM       MIN       MAX')
1082
       WRITE(Message,'(A,T30,3I10)') '  Total dofs: ',nxpar
1083
       CALL Info(Caller,Message,Level=3)
1084
       IF( HaveParInfo ) THEN
1085
         WRITE(Message,'(A,T30,3I10)') '  Owned dofs: ',nopar
1086
         CALL Info(Caller,Message,Level=3)
1087
         WRITE(Message,'(A,T30,3I10)') '  Shared dofs: ',nspar
1088
         CALL Info(Caller,Message,Level=3)
1089
       END IF
1090
       WRITE(Message,'(A,T30,3I10)') '  Matrix nnz: ',nzpar
1091
       CALL Info(Caller,Message,Level=3)
1092
     ELSE
1093
       CALL Info(Caller,'Serial size for Solver:')
1094
       WRITE(Message,'(A,T30,1I10)') '  Total dofs: ',nx
1095
       CALL Info(Caller,Message,Level=3)
1096
       WRITE(Message,'(A,T30,1I10)') '  Matrix nnz: ',nz
1097
       CALL Info(Caller,Message,Level=3)
1098
     END IF
1099
     
1100
   END SUBROUTINE PrintProblemSize
1101
   
1102
   
1103
!------------------------------------------------------------------------------
1104
!> Add the generic stuff related to each Solver. 
1105
!> A few solvers are for historical reasons given a special treatment. 
1106
!------------------------------------------------------------------------------
1107
  SUBROUTINE AddEquationBasics( Solver, Name, Transient )
1108
!------------------------------------------------------------------------------
1109
    USE CoordinateSystems
1110
    TYPE(Solver_t), POINTER :: Solver
1111
    LOGICAL :: Transient
1112
    CHARACTER(LEN=*) :: Name
1113
!------------------------------------------------------------------------------
1114
    REAL(KIND=dp), POINTER :: Solution(:)
1115

1116
    INTEGER, POINTER :: Perm(:)
1117

1118
    INTEGER(KIND=AddrInt) :: InitProc, AssProc
1119

1120
    INTEGER :: MaxDGDOFs, MaxNDOFs, MaxEDOFs, MaxFDOFs, MaxBDOFs, MaxDOFsPerNode
1121
    INTEGER :: i,j,k,l,NDeg,Nrows,nSize,n,m,DOFs,dim,MatrixFormat,istat,Maxdim, AllocStat, &
1122
        i1,i2,i3,iostat
1123

1124
    LOGICAL :: Found, Stat, BandwidthOptimize, EigAnal, ComplexFlag, &
1125
        MultigridActive, VariableOutput, GlobalBubbles, HarmonicAnal, MGAlgebraic, &
1126
        VariableGlobal, VariableIP, VariableElem, VariableDG, VariableNodal, &
1127
        DG, NoMatrix, IsAssemblySolver, IsCoupledSolver, IsBlockSolver, IsProcedure, &
1128
        IsStepsSolver, LegacySolver, UseMask, TransientVar, InheritVarType, DoIt, &
1129
        GotSecName, NoPerm
1130
    
1131
!    CHARACTER(LEN=MAX_NAME_LEN) :: var_name
1132
    CHARACTER(:), ALLOCATABLE :: proc_name, tmpname, mask_name,sec_name,eq,str,var_name
1133

1134
    TYPE(ValueList_t), POINTER :: SolverParams
1135
    TYPE(Mesh_t),   POINTER :: NewMesh,OldMesh
1136
    TYPE(Element_t), POINTER :: CurrentElement
1137
    TYPE(Matrix_t), POINTER :: NewMatrix, Oldmatrix
1138

1139
    TYPE(Variable_t), POINTER :: Var, PVar
1140
    TYPE(Variable_t), POINTER :: NewVariable
1141

1142
    REAL(KIND=dp) :: tt, InitValue
1143
    REAL(KIND=dp), POINTER :: Component(:)
1144

1145
    TYPE(Graph_t) :: DualGraph
1146
    TYPE(GraphColour_t) :: GraphColouring, BoundaryGraphColouring
1147
    LOGICAL :: ConsistentColours
1148

1149
    LOGICAL :: ThreadedStartup, MultiColourSolver, HavePerm, Parallel
1150
    INTEGER :: VariableType
1151

1152
    CHARACTER(*), PARAMETER :: Caller="AddEquationBasics"
1153

1154
    
1155
    ! Set pointer to the list of solver parameters
1156
    !------------------------------------------------------------------------------
1157
    SolverParams => ListGetSolverParams(Solver)
1158

1159
    !------------------------------------------------------------------------------
1160
    ! Check the historical solvers that may be built-in on some .sif files
1161
    ! Therefore some special strategies are used for them.
1162
    !------------------------------------------------------------------------------
1163
    IsProcedure = ListCheckPresent( SolverParams, 'Procedure' )
1164
    Dim = CoordinateSystemDimension()        
1165
    Dofs = 1
1166
    InitValue = 0.0_dp
1167
    LegacySolver = .TRUE.
1168

1169
    SELECT CASE( Name )       
1170
      !------------------------------------------------------------------------------
1171
      
1172
      !------------------------------------------------------------------------------
1173
    CASE('navier-stokes')
1174
      !------------------------------------------------------------------------------
1175
      dofs = dim 
1176
      IF ( CurrentCoordinateSystem() == CylindricSymmetric ) DOFs = DOFs + 1
1177
      IF( dofs == 3 ) THEN
1178
        var_name = 'Flow Solution[Velocity:3 Pressure:1]'
1179
      ELSE
1180
        var_name = 'Flow Solution[Velocity:2 Pressure:1]'
1181
      END IF
1182
      proc_name = 'FlowSolve FlowSolver'
1183
      InitValue = 1.0d-6
1184
      ! We don't want to use automated setting of dofs later so set this back to one!
1185
      dofs = 1
1186

1187
      !------------------------------------------------------------------------------
1188
    CASE('magnetic induction')
1189
      !------------------------------------------------------------------------------
1190
      var_name = 'Magnetic Field'
1191
      proc_name = 'MagneticSolve MagneticSolver'
1192
      dofs = 3
1193
      CALL ListAddString( SolverParams,&
1194
          NextFreeKeyword('Exported Variable',SolverParams),&
1195
          'Electric Current[Electric Current:3]')                  
1196
      
1197
      !------------------------------------------------------------------------------
1198
    CASE('stress analysis')
1199
      !------------------------------------------------------------------------------
1200
      dofs = dim
1201
      var_name = 'Displacement'
1202
      proc_name = 'StressSolve StressSolver'      
1203
            
1204
      !------------------------------------------------------------------------------
1205
    CASE('mesh update')
1206
      !------------------------------------------------------------------------------
1207
      dofs = dim
1208
      var_name = 'Mesh Update'
1209
      proc_name = 'MeshSolve MeshSolver'
1210

1211
      IF( Transient ) THEN
1212
        IF( Dofs == 2 ) THEN
1213
          CALL ListAddString( SolverParams,&
1214
              NextFreeKeyword('Exported Variable',SolverParams),&
1215
              '-dofs 2 Mesh Velocity')        
1216
        ELSE
1217
          CALL ListAddString( SolverParams,&
1218
              NextFreeKeyword('Exported Variable',SolverParams),&
1219
              '-dofs 3 Mesh Velocity')                  
1220
        END IF
1221
      END IF
1222

1223
      !------------------------------------------------------------------------------
1224
    CASE('heat equation')
1225
      !------------------------------------------------------------------------------
1226
      var_name = 'Temperature'
1227
      proc_name = 'HeatSolve HeatSolver'
1228
      
1229
      CALL ListAddNewLogical( SolverParams,'Radiation Solver',.TRUE.)
1230
      !------------------------------------------------------------------------------
1231

1232
    CASE DEFAULT
1233
      LegacySolver = .FALSE.
1234
      
1235
    END SELECT
1236

1237
    IF( LegacySolver ) THEN
1238
      CALL Info(Caller,'Setting up keywords internally for legacy solver: '&
1239
          //TRIM(Name),Level=10)
1240
      IF( .NOT. ListCheckPresent( SolverParams,'Variable') ) THEN
1241
        CALL ListAddString( SolverParams,'Variable',var_name )
1242
        IF( Dofs > 1 ) CALL ListAddInteger( SolverParams,'Variable Dofs',dofs )
1243
      END IF
1244
      IF( .NOT. IsProcedure ) THEN      
1245
        CALL ListAddString(SolverParams, 'Procedure', proc_name,.FALSE.)
1246
      END IF
1247
    END IF
1248

1249
    ! We should have the procedure 
1250
    proc_name = ListGetString( SolverParams, 'Procedure',IsProcedure)
1251
    IF( IsProcedure ) THEN
1252
      CALL Info(Caller,'Using procedure: '//TRIM(proc_name),Level=10)
1253
    ELSE
1254
      CALL Warn(Caller,'Solver '//I2S(Solver % SolverId)//' may require "Procedure" to operate!')
1255
    END IF
1256

1257

1258
#if 0
1259
    ! Here we can for testing purposes change all solver named xxx to yyy
1260
    ! and then rerun all the tests, for example. Don't activate unless for this kind
1261
    ! of development work.
1262
    i = INDEX( proc_name,'HeatSolver')
1263
    IF( i /= 0 ) THEN
1264
      proc_name = 'HeatSolveVec HeatSolver'
1265
      CALL ListAddString(SolverParams, 'Procedure', proc_name,.FALSE.)
1266
      CALL Info(Caller,'Using procedure: '//TRIM(proc_name),Level=6)
1267
    END IF
1268
#endif
1269

1270
    ! Mark whether the solver will be treated as parallel in future
1271
    ! Mainly related to the linear system solution since assembly is trivially parallel.
1272
    ! There are some special cases where we have same mesh for multiple instances of almost
1273
    ! same independent problem. 
1274
    Parallel = ( ParEnv % PEs > 1 )
1275
    IF( Parallel ) THEN
1276
      IF(Solver % Mesh % SingleMesh ) THEN      
1277
        Parallel = ListGetLogical( SolverParams,'Enforce Parallel',Found )
1278
        IF(.NOT. Found) THEN
1279
          Parallel = ListGetLogical( CurrentModel % Simulation,'Enforce Parallel',Found )
1280
        END IF
1281
      END IF
1282
    END IF
1283
    Solver % Parallel = Parallel 
1284

1285
    
1286
    ! If there is a matrix level Flux Corrected Transport and/or nonlinear timestepping
1287
    ! then you must use global matrices for time integration.
1288
    !----------------------------------------------------------------------------------
1289
    DoIt = .FALSE.
1290
    IF( ListGetLogical( SolverParams,'Linear System FCT',Found ) ) DoIt = .TRUE.
1291
    IF( ListGetLogical( SolverParams,'Nonlinear Timestepping',Found ) ) DoIt = .TRUE.
1292
    IF( ListGetLogical( SolverParams,'Apply Conforming BCs',Found ) ) DoIt = .TRUE.
1293
    
1294
    IF( DoIt ) THEN
1295
      CALL Info(Caller,'Enforcing use of global mass matrix needed by other features!')
1296
      CALL ListAddLogical( SolverParams,'Use Global Mass Matrix',.TRUE.)
1297
    END IF
1298

1299
    ! Compute the mesh dimension for this solver
1300
    !----------------------------------------------------------------------------
1301
    eq = ListGetString( SolverParams, 'Equation', Found )
1302
    IF( Found ) THEN
1303
      CALL Info(Caller,'Setting up solver: '//TRIM(eq),Level=8)
1304
    END IF
1305

1306
    IF ( Found ) THEN
1307
      MAXdim = 0
1308
      DO i=1,Solver % Mesh % NumberOfBulkElements+Solver % Mesh % NumberOFBoundaryElements
1309
        CurrentElement => Solver % Mesh % Elements(i)
1310
        IF ( CheckElementEquation( CurrentModel, CurrentElement, eq ) ) THEN
1311
          Maxdim = MAX( CurrentElement % TYPE % DIMENSION, Maxdim )
1312
        END IF
1313
      END DO
1314
      CALL ListAddInteger( SolverParams, 'Active Mesh Dimension', Maxdim )
1315
    END IF
1316

1317
    ! Check the solver for initialization
1318
    ! This is utilized only by some solvers.
1319
    !-----------------------------------------------------------------
1320
    IF( IsProcedure ) THEN
1321
      InitProc = GetProcAddr( TRIM(proc_name)//'_Init', abort=.FALSE. )
1322
      CALL Info(Caller,'Checking for _init solver',Level=12)
1323
      IF ( InitProc /= 0 ) THEN
1324
        CALL ExecSolver( InitProc, CurrentModel, Solver, &
1325
            Solver % dt, Transient )
1326
      END IF
1327
    END IF
1328

1329
    Solver % SolverMode = SOLVER_MODE_DEFAULT
1330
    IF( ListGetLogical( SolverParams, 'Auxiliary Solver', Found ) ) &
1331
        Solver % SolverMode = SOLVER_MODE_AUXILIARY
1332
    IF( ListGetLogical( SolverParams, 'Coupled Solver', Found ) ) &
1333
        Solver % SolverMode = SOLVER_MODE_COUPLED
1334
    IF( ListGetLogical( SolverParams, 'Block Solver', Found ) ) &
1335
        Solver % SolverMode = SOLVER_MODE_BLOCK
1336
    IF( ListGetLogical( SolverParams, 'Assembly Solver', Found ) ) &
1337
        Solver % SolverMode = SOLVER_MODE_ASSEMBLY
1338

1339
    IF( Solver % SolverMode == SOLVER_MODE_DEFAULT ) THEN
1340
      eq = ListGetString( Solver  % Values, 'Equation', Found )
1341
      IF(.NOT. Found ) THEN
1342
        Solver % SolverMode = SOLVER_MODE_AUXILIARY 
1343
      ELSE
1344
        AssProc = GetProcAddr( TRIM(proc_name)//'_bulk', abort=.FALSE. )
1345
        CALL Info(Caller,'Checking for _bulk solver',Level=12)
1346
        IF ( AssProc /= 0 ) THEN
1347
          CALL Info(Caller,'Solver will be be performed in steps',Level=8)
1348
          Solver % SolverMode = SOLVER_MODE_STEPS
1349
        END IF        
1350
      END IF
1351
    END IF
1352

1353
    IsCoupledSolver = ( Solver % SolverMode == SOLVER_MODE_COUPLED ) 
1354
    IsBlockSolver = ( Solver % SolverMode == SOLVER_MODE_BLOCK ) 
1355
    IsAssemblySolver = ( Solver % SolverMode == SOLVER_MODE_ASSEMBLY ) 
1356
    IsAssemblySolver = IsAssemblySolver .OR. &
1357
	( IsCoupledSolver .AND. .NOT. IsProcedure ) .OR. &
1358
        ( IsBlockSolver .AND. .NOT. IsProcedure ) 
1359
    IsStepsSolver = ( Solver % SolverMode == SOLVER_MODE_STEPS ) 
1360

1361
    ! Get the procedure that really runs the solver
1362
    !------------------------------------------------------------------------------
1363
    IF( IsProcedure ) THEN
1364
      IF( IsStepsSolver ) THEN
1365
        Solver % PROCEDURE = AssProc
1366
      ELSE
1367
        Solver % PROCEDURE = GetProcAddr(proc_name)
1368
      END IF
1369
    END IF
1370

1371
    ! Default order of equation
1372
    !--------------------------
1373
    Solver % Order = 1
1374
    Solver % TimeOrder = 1
1375
    
1376
    ! Set up time-stepping strategies for transient problems
1377
    !------------------------------------------------------------------------------
1378
    IF ( Transient ) THEN
1379
      str = ListGetString( SolverParams, 'Predictor Method',Found )
1380
      IF ( .NOT. Found ) THEN
1381
        str = ListGetString( CurrentModel % Simulation, 'Predictor Method',Found )
1382
        IF ( Found ) THEN
1383
          CALL ListAddString( SolverParams, 'Predictor Method', str )
1384
        END IF
1385
      END IF
1386

1387
      str = ListGetString( SolverParams, 'Corrector Method',Found )
1388
      IF ( .NOT. Found ) THEN
1389
        str = ListGetString( CurrentModel % Simulation, 'Corrector Method',Found )
1390
        IF ( Found ) THEN
1391
          CALL ListAddString( SolverParams, 'Corrector Method', str )
1392
        END IF
1393
      END IF
1394

1395
      IF ( Found ) THEN
1396
        CALL ReadPredCorrParams( CurrentModel, SolverParams )
1397
        CALL ListAddString( SolverParams, 'Timestepping Method', str )  
1398
      END IF
1399

1400
      str = ListGetString( SolverParams, 'Timestepping Method',Found )
1401
      IF ( .NOT. Found ) THEN
1402
        str = ListGetString( CurrentModel % Simulation, 'Timestepping Method',Found )
1403
        IF ( Found ) THEN
1404
          CALL ListAddString( SolverParams, 'Timestepping Method', str )
1405
        END IF
1406
      END IF
1407
      
1408
      IF ( Found ) THEN
1409
        IF (str=='bdf') THEN
1410
          Solver % Order = ListGetInteger( SolverParams, &
1411
              'BDF Order', Found, minv=1, maxv=5 )
1412
          IF ( .NOT. Found ) THEN
1413
            Solver % Order = ListGetInteger( CurrentModel % &
1414
                Simulation, 'BDF Order', Found, minv=1, maxv=5 )
1415
          END IF
1416
          IF ( .NOT.Found ) THEN
1417
            Solver % Order = 2
1418
            CALL Warn(Caller, 'BDF order defaulted to 2.' )
1419
          END IF
1420
        ELSE IF ( str=='runge-kutta') THEN
1421
          Solver % Order = ListGetInteger( CurrentModel % &
1422
              Simulation, 'Runge-Kutta Order', Found, minv=2, maxv=4 )
1423
          IF ( .NOT.Found ) Solver % Order = 2
1424
        ELSE IF( SEQL(str,'adams')) THEN
1425
          Solver % Order = 2          
1426
        END IF
1427
        CALL Info(Caller,'Time stepping method is: '//TRIM(str),Level=10)
1428
      ELSE
1429
        CALL Warn(Caller, '> Timestepping method < defaulted to > Implicit Euler <' )
1430
        CALL ListAddString( SolverParams, 'Timestepping Method', 'Implicit Euler' )
1431
      END IF
1432

1433
    END IF
1434

1435
    ! Initialize and get the variable 
1436
    !------------------------------------------------------------------------    
1437
    Solver % TimeOrder = 0
1438
    NULLIFY( Solver % Matrix )    
1439
    var_name = ListGetString( SolverParams, 'Variable', Found )
1440

1441
    NoPerm = .FALSE.
1442
    
1443
    IF(.NOT. Found ) THEN
1444
      ! Variable does not exist
1445
      !------------------------------------------------------
1446
      CALL Info(Caller,'Creating null variable with no name',Level=15)
1447

1448
      ALLOCATE( Solver % Variable )
1449
      Solver % Variable % Name = ''
1450
      Solver % Variable % NameLen = 0
1451
      Solver % Variable % Norm = 0.0d0
1452
      NULLIFY( Solver % Variable % Perm )
1453
      NULLIFY( Solver % Variable % Values )
1454
      
1455
      
1456
    ELSE IF( IsCoupledSolver .AND. .NOT. IsProcedure ) THEN
1457
      ! Coupled solver may inherit the matrix only if procedure is given
1458
      !-----------------------------------------------------------------
1459

1460
    ELSE IF( IsBlockSolver .AND. .NOT. IsProcedure ) THEN
1461
      ! Block solver may inherit the matrix only if procedure is given
1462
      !-----------------------------------------------------------------
1463
      
1464
    ELSE
1465
      CALL Info(Caller,'Treating variable string: '//TRIM(var_name),Level=15)
1466

1467
      ! It may be a normal field variable or a global (0D) variable
1468
      !------------------------------------------------------------------------
1469
      VariableGlobal = ListGetLogical( SolverParams, 'Variable Global', Found )
1470
      
1471
      VariableOutput = ListGetLogical( SolverParams, 'Variable Output', Found )
1472
      IF ( .NOT. Found ) VariableOutput = .TRUE.
1473

1474
      VariableIp = ListGetLogical( SolverParams, 'Variable IP', Found )
1475
      VariableElem = ListGetLogical( SolverParams, 'Variable Elemental', Found )                 
1476

1477
      
1478
      DOFs = ListGetInteger( SolverParams, 'Variable DOFs', Found, minv=1 )
1479
      IF ( .NOT. Found ) THEN
1480
        j = 0
1481
        DOFs = 0
1482
        DO WHILE( .TRUE. )
1483
          i = INDEX( var_name(j+1:), ':' ) + j
1484
          IF ( i<=j ) EXIT
1485
          READ( var_name(i+1:),'(i1)',IOSTAT=iostat) k
1486
          IF(iostat /= 0) THEN
1487
            CALL Fatal(Caller,'Could not read component count of variable!')
1488
          END IF
1489
          DOFs = DOFs + k
1490
          j = i + 1
1491
        END DO
1492
      END IF
1493
      
1494
      DO WHILE( var_name(1:1) == '-' )
1495
        k = 0
1496
        j = LEN_TRIM(var_name)
1497

1498
        IF ( SEQL(var_name, '-noperm ') ) THEN
1499
          NoPerm = .TRUE.
1500
          k = 8
1501

1502
        ELSE IF ( SEQL(var_name, '-nooutput ') ) THEN
1503
          VariableOutput = .FALSE.
1504
          k = 10
1505
        
1506
        ELSE IF ( SEQL(var_name, '-global ') ) THEN
1507
          VariableGlobal = .TRUE.
1508
          k = 8
1509
        
1510
        ELSE IF ( SEQL(var_name, '-ip ') ) THEN
1511
          VariableIp = .TRUE.
1512
          k = 4
1513

1514
        ELSE IF ( SEQL(var_name, '-elem ') ) THEN
1515
          VariableElem = .TRUE.
1516
          k = 6 
1517
               
1518
        ELSE IF ( SEQL(var_name, '-dofs ') ) THEN
1519
          READ( var_name(7:), *, IOSTAT=iostat) DOFs
1520
          IF(iostat /= 0) THEN
1521
            CALL Fatal(Caller,'Could not read number after -dofs of variable!')
1522
          END IF
1523
          i = 7
1524
          DO WHILE( var_name(i:i) /= ' '  )
1525
            i = i + 1
1526
            IF ( i > j ) EXIT
1527
          END DO
1528
          k = i
1529
        ELSE
1530
          CALL Fatal(Caller,'Do not know how to parse: '//TRIM(var_name))
1531
        END IF
1532

1533
        IF(k>0) THEN
1534
          var_name(1:j-k) = var_name(k+1:j)
1535
          DO i=j-k+1,j
1536
            var_name(i:i) = ' '
1537
          END DO
1538
        END IF
1539

1540
      END DO
1541
      IF ( DOFs == 0 ) DOFs = 1
1542
      
1543
      n = LEN_TRIM(var_name)
1544
      
1545
      ! If the variable is 'global' it has nothing to do with the mesh and it may be simply 
1546
      ! allocated. 
1547
      !------------------------------------------------------------------------------------
1548
      IF( VariableGlobal ) THEN
1549
        CALL Info(Caller,'Creating global variable: '//var_name(1:n),Level=8)
1550

1551
        Solver % SolverMode = SOLVER_MODE_GLOBAL
1552
        ALLOCATE( Solution( DOFs ) )
1553
        Solution = 0.0_dp
1554
        
1555
        CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver, &
1556
            var_name(1:n), DOFs, Solution )
1557
        Solver % Variable => VariableGet( Solver % Mesh % Variables, var_name(1:n) )
1558
        IF( DOFs > 1 ) THEN
1559
          DO i=1,DOFs
1560
            tmpname = ComponentName( var_name(1:n), i )
1561
            Component => Solution( i:i )
1562
            CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
1563
                tmpname, 1, Component )
1564
          END DO
1565
        END IF
1566

1567
        IF( ListGetLogical( SolverParams,'Ode Matrix',Found ) ) THEN
1568
          CALL Info(Caller,'Creating dense matrix for ODE: '&
1569
              //I2S(Dofs),Level=8)
1570
          ALLOCATE( Solver % Variable % Perm(1) )
1571
          Solver % Variable % Perm(1) = 1
1572
          Solver % Matrix => CreateOdeMatrix( CurrentModel, Solver, Dofs )
1573
        END IF
1574

1575
      ELSE        
1576
        CALL Info(Caller,'Creating standard variable: '//var_name(1:n),Level=8)
1577

1578
        ! If the variable is a field variable create a permutation and matrix related to it
1579
        !----------------------------------------------------------------------------------
1580
        eq = ListGetString( SolverParams, 'Equation', Found )
1581
        IF(.NOT. Found) THEN
1582
          CALL Fatal(Caller,'Variable exists but > Equation < is not defined in Solver ')
1583
        END IF
1584
        Found = .FALSE.
1585
        DO i=1, CurrentModel % NumberOfEquations
1586
          IF( ListGetLogical( CurrentModel % Equations(i) % Values, TRIM(eq), Stat)) THEN
1587
            Found = .TRUE.
1588
            EXIT
1589
          END IF
1590
        END DO
1591
        IF(.NOT. Found ) THEN
1592
          CALL Fatal(Caller,'Variable > '//var_name(1:n)//&
1593
                     ' < exists but it is not associated to any equation')
1594
        END IF
1595
        
1596
        ! Computate the size of the permutation vector
1597
        !-----------------------------------------------------------------------------------------
1598
        CALL Info(Caller,'Computing size of permutation vector',Level=12)
1599
        Ndeg = 0
1600

1601
        IF(.TRUE.) THEN
1602
          eq = ListGetString( Solver  % Values, 'Equation', Found )
1603
          MaxNDOFs  = 0
1604
          MaxDOFsPerNode = 0
1605
          MaxDGDOFs = 0
1606
          MaxBDOFs = 0
1607
          DO i=1,Solver % Mesh % NumberOFBulkElements
1608
            CurrentElement => Solver % Mesh % Elements(i)
1609
            MaxDGDOFs = MAX( MaxDGDOFs, CurrentElement % DGDOFs )
1610
            MaxBDOFs  = MAX( MaxBDOFs,  CurrentElement % BDOFs )
1611
            MaxNDOFs  = MAX( MaxNDOFs,  CurrentElement % NDOFs )
1612
            MaxDOFsPerNode =  MAX( MaxDOFsPerNode, CurrentElement % NDOFs / &
1613
                CurrentElement % TYPE % NumberOfNodes )
1614
          END DO
1615
          
1616
          MaxEDOFs = 0
1617
          DO i=1,Solver % Mesh % NumberOFEdges 
1618
            CurrentElement => Solver % Mesh % Edges(i)
1619
            MaxEDOFs  = MAX( MaxEDOFs,  CurrentElement % BDOFs )
1620
          END DO
1621
          
1622
          MaxFDOFs = 0
1623
          DO i=1,Solver % Mesh % NumberOFFaces 
1624
            CurrentElement => Solver % Mesh % Faces(i)
1625
            MaxFDOFs  = MAX( MaxFDOFs,  CurrentElement % BDOFs )
1626
          END DO
1627
          
1628
          GlobalBubbles = ListGetLogical( SolverParams, 'Bubbles in Global System', Found )
1629
          IF (.NOT.Found) GlobalBubbles = .TRUE.
1630

1631
          Ndeg = Ndeg + MaxDOFsPerNode * Solver % Mesh % NumberOfNodes 
1632
          IF ( GlobalBubbles ) THEN
1633
            Ndeg = Ndeg + MaxBDOFs * Solver % Mesh % NumberOfBulkElements
1634

1635
            DO i=1,Solver % Mesh % NumberOfBoundaryElements
1636
              j = i + Solver % Mesh % NumberOfBulkElements
1637
              IF ( Solver % Mesh % Elements(j) % Type % ElementCode >= 300) THEN
1638
                MaxFDOFs = MAX( maxFDOFs, Solver % Mesh % Elements(j) % BDOFs )
1639
              ELSE
1640
                MaxEDOFs = MAX( maxEDOFs, Solver % Mesh % Elements(j) % BDOFs )
1641
              END IF
1642
            END DO
1643
          END IF
1644

1645
          IF ( MaxEDOFs > 0 ) Ndeg = Ndeg + MaxEDOFs * Solver % Mesh % NumberOFEdges
1646
          IF ( MaxFDOFs > 0 ) Ndeg = Ndeg + MaxFDOFs * Solver % Mesh % NumberOFFaces
1647

1648
          DG = ListGetLogical( SolverParams, 'Discontinuous Galerkin', Found )
1649
          IF( DG ) THEN
1650
            Ndeg = MAX( NDeg, MaxDGDOFs * (Solver % Mesh % NumberOfBulkElements+ &
1651
                Solver % Mesh % NumberOfBoundaryElements) )
1652
          END IF
1653
        END IF
1654

1655
        IF( ListGetLogical( SolverParams,'Radiation Solver',Found ) ) THEN
1656
          ! We need to precompute view factors if they are included in CRS matrix
1657
          ! Benefit of doing it at later stage is that we may modify the geometry, for example. 
1658
          IF( .NOT. (ListGetLogical( SolverParams,'Radiosity Model',Found ) .OR. &
1659
              ListGetLogical( SolverParams,'Spectral Model',Found ) .OR. &
1660
              ListGetLogical( SolverParams,'Update Gebhart Factors',Found ) ) ) THEN            
1661
            ! If we need to update the Gebhart factors we may not create the matrix topology at this
1662
            ! stage as we don't know the emissivities yet.
1663
            CALL RadiationFactors( Solver, .TRUE., .FALSE.)
1664
          END IF
1665
        END IF
1666
        
1667
        BandwidthOptimize = ListGetLogical( SolverParams, &
1668
            'Optimize Bandwidth', Found )
1669
        IF ( .NOT. Found ) BandwidthOptimize = .TRUE.
1670
        CALL CheckLinearSolverOptions( Solver )
1671

1672
        CALL Info(Caller,'Maximum size of permutation vector is: '//I2S(Ndeg),Level=12)
1673
        ALLOCATE( Perm(Ndeg), STAT=AllocStat )
1674
        IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for Perm')
1675
        Perm = 0
1676
        MatrixFormat = MATRIX_CRS
1677

1678
        ThreadedStartup = ListGetLogical( SolverParams,'Multithreaded Startup',Found )
1679
        IF( ThreadedStartup ) THEN
1680
          CALL Info(Caller,'Using multithreaded startup',Level=6)
1681
        END IF
1682
        
1683
        MultiColourSolver = ListGetLogical( SolverParams,'MultiColour Solver',Found )
1684

1685
        IF( MultiColourSolver .OR. ThreadedStartup ) THEN
1686
          CALL Info(Caller,'Creating structures for mesh colouring',Level=8)
1687
          ConsistentColours = .FALSE.
1688
          IF ( ListGetLogical(SolverParams,'MultiColour Consistent', Found) ) THEN
1689
            CALL Info(Caller,'Creating consistent colouring',Level=8)
1690
            ConsistentColours = .TRUE.
1691
          END IF
1692

1693
          IF (Solver % Mesh % NumberOfBulkElements > 0) THEN 
1694
             ! Construct the dual graph from Elmer mesh
1695
             CALL ElmerMeshToDualGraph(Solver % Mesh, DualGraph)
1696
             ! Colour the dual graph
1697
             CALL ElmerGraphColour(DualGraph, GraphColouring, ConsistentColours)
1698
             ! Deallocate dual graph as it is no longer needed
1699
             CALL Graph_Deallocate(DualGraph)
1700
          
1701
             ! Construct colour lists
1702
             ALLOCATE( Solver % ColourIndexList, STAT=AllocStat )
1703
             IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for ColourIndexList')
1704
             CALL ElmerColouringToGraph(GraphColouring, Solver % ColourIndexList)
1705
             CALL Colouring_Deallocate(GraphColouring)
1706
          END IF
1707

1708
          IF (Solver % Mesh % NumberOfBoundaryElements > 0) THEN 
1709
             ! Construct the dual graph from Elmer boundary mesh
1710
             CALL ElmerMeshToDualGraph(Solver % Mesh, DualGraph, UseBoundaryMesh=.TRUE.)
1711
             ! Colour the dual graph of boundary mesh
1712
             CALL ElmerGraphColour(DualGraph, BoundaryGraphColouring, ConsistentColours)
1713
             ! Deallocate dual graph as it is no longer needed
1714
             CALL Graph_Deallocate(DualGraph)
1715
                       
1716
             ALLOCATE( Solver % BoundaryColourIndexList, STAT=AllocStat )
1717
             IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for BoundaryColourIndexList')
1718
             CALL ElmerColouringToGraph(BoundaryGraphColouring, Solver % BoundaryColourIndexList)
1719
             CALL Colouring_Deallocate(BoundaryGraphColouring)
1720
          END IF
1721

1722
          ! DEBUG/TODO: Add a Solver keyword for enabling the functionality
1723
          ! CALL CheckColourings(Solver)
1724
        END IF
1725
        
1726
        CALL Info(Caller,'Creating solver matrix topology',Level=12)
1727
        Solver % Matrix => CreateMatrix( CurrentModel, Solver, Solver % Mesh, &
1728
            Perm, DOFs, MatrixFormat, BandwidthOptimize, eq(1:LEN_TRIM(eq)), DG, &
1729
            GlobalBubbles=GlobalBubbles, ThreadedStartup=ThreadedStartup )
1730
        Solver % GlobalBubbles = GlobalBubbles
1731

1732
        Nrows = DOFs * Ndeg
1733
        IF (ASSOCIATED(Solver % Matrix)) THEN
1734
          Nrows = Solver % Matrix % NumberOfRows
1735
          CALL Info(Caller,'Number of rows in CRS matrix: '//I2S(Nrows),Level=12)
1736
        END IF
1737
        
1738
        ! Check if mesh colouring is needed by the solver
1739
        IF( .NOT. MultiColourSolver .AND. ThreadedStartup ) THEN
1740
           ! Deallocate mesh colouring
1741
           IF (ASSOCIATED(Solver % ColourIndexList)) THEN
1742
              CALL Graph_Deallocate(Solver % ColourIndexList)
1743
              DEALLOCATE(Solver % ColourIndexList)
1744
              Solver % ColourIndexList => NULL()
1745
           END IF
1746

1747
           ! Deallocate boundary mesh colouring
1748
           IF (ASSOCIATED(Solver % BoundaryColourIndexList)) THEN
1749
              CALL Graph_Deallocate(Solver % BoundaryColourIndexList)
1750
              DEALLOCATE(Solver % BoundaryColourIndexList)
1751
              Solver % BoundaryColourIndexList => NULL()
1752
           END IF
1753
        END IF
1754

1755
        ! Basically the solver could be matrix free but still the matrix
1756
        ! is used here temporarily since it is needed when making the 
1757
        ! permutation vector
1758
        !-----------------------------------------------------------------
1759
        IF( ListGetLogical( SolverParams, 'No Matrix', Found ) ) THEN
1760
          CALL Info(Caller,'No matrix needed any more, freeing structures!',Level=12)
1761
          Solver % SolverMode = SOLVER_MODE_MATRIXFREE
1762
          CALL FreeMatrix( Solver % Matrix )
1763
        END IF
1764
       
1765
        CALL Info(Caller,'Creating solver variable',Level=12)
1766
        ALLOCATE(Solution(Nrows),STAT=AllocStat)
1767
        IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for Solution')
1768

1769
        !$OMP PARALLEL DO
1770
        DO i=1,Nrows
1771
          Solution(i) = InitValue
1772
        END DO
1773
        !$OMP END PARALLEL DO
1774
        CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver, &
1775
            var_name(1:n), DOFs, Solution, Perm, Output=VariableOutput )          
1776
        Solver % Variable => VariableGet( Solver % Mesh % Variables, var_name(1:n) )        
1777
        Solver % Variable % PeriodicFlipActive = Solver % PeriodicFlipActive
1778
        
1779
        IF ( DOFs > 1 ) THEN
1780
          DO i=1,DOFs
1781
            tmpname = ComponentName( var_name(1:n), i )
1782
            Component => Solution( i:Nrows-DOFs+i:DOFs )
1783
            CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
1784
                tmpname, 1, Component, Perm, Output=VariableOutput )
1785
            PVar => VariableGet( Solver % Mesh % Variables, tmpname )
1786
            PVar % PeriodicFlipActive = Solver % PeriodicFlipActive
1787
          END DO
1788
        END IF
1789

1790
        ! Optionally save the limiters as a field. This is allocated here so that
1791
        ! if it used as a dependent variable it is allocated before being used.
1792
        IF( ListGetLogical( SolverParams,'Apply Limiter', Found ) ) THEN
1793
          IF( ListGetLogical( SolverParams,'Save Limiter',Found ) ) THEN      
1794
            CALL Info(Caller,'Adding "contact active" field for '//var_name(1:n))
1795
            CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, Solver,&
1796
                TRIM(GetVarName(Solver % Variable))//' Contact Active', &
1797
                dofs = Solver % Variable % Dofs, Perm = Solver % Variable % Perm )
1798
          END IF
1799
        END IF
1800
                    
1801
        IF (ASSOCIATED(Solver % Matrix)) Solver % Matrix % Comm = ELMER_COMM_WORLD
1802

1803
        IF ( DG ) THEN
1804
          Solver % Variable % TYPE = Variable_on_nodes_on_elements
1805
        END IF
1806

1807
        IF( ListGetLogical(SolverParams,'Hcurl Basis',Found ) ) THEN
1808
          Solver % Variable % Type = Variable_on_edges
1809
        END IF
1810
        
1811
      END IF
1812
      !------------------------------------------------------------------------------
1813
    END IF
1814

1815
    IF( ListGetLogical( SolverParams,'Create Integration Points Table',Found ) ) THEN
1816
      CALL CreateIpPerm( Solver ) 
1817
    END IF
1818
    
1819
    !------------------------------------------------------------------------------
1820
    ! Add the exported variables which are typically auxiliary variables derived
1821
    ! from the solution without their own matrix equation.  
1822
    !------------------------------------------------------------------------------
1823
    l = 0
1824
    DO WHILE( .TRUE. )
1825
      l = l + 1
1826
      str = ComponentName( 'exported variable', l )
1827
      var_name = ListGetString( SolverParams, str, Found )
1828
      
1829
      IF(.NOT. Found) EXIT
1830
      
1831
      CALL Info(Caller,'Creating exported variable: '//TRIM(var_name),Level=12)
1832

1833
      str = TRIM( ComponentName( 'exported variable', l ) ) // ' Output'
1834
      VariableOutput = ListGetLogical( SolverParams, str, Found )
1835
      IF ( .NOT. Found ) VariableOutput = .TRUE.
1836

1837
      str = TRIM( ComponentName( 'exported variable', l ) ) // ' Mask'
1838
      tmpname = ListGetString( SolverParams, str, UseMask )
1839
     
1840
      GotSecName = .FALSE.
1841
      IF( UseMask ) THEN
1842
        i3 = LEN_TRIM(tmpname) 
1843
        i1 = INDEX(tmpname,':')
1844
        IF( i1 > 1 ) THEN
1845
          i2 = i1 + 1
1846
          i1 = i1 - 1
1847
          IF( i1 > 0 ) THEN
1848
            DO WHILE( tmpname(i2:i2) == ' ')
1849
              i2 = i2 + 1 
1850
              IF(i2>i3) EXIT
1851
            END DO
1852
          END IF
1853
          sec_name = tmpname(1:i1)
1854
          mask_name = tmpname(i2:i3)
1855
          CALL Info(Caller,'masking with section: '//TRIM(sec_name),Level=12)
1856
          CALL Info(Caller,'masking with keyword: '//TRIM(mask_name),Level=12)
1857
          GotSecName = .TRUE.
1858
        ELSE          
1859
          sec_name = 'body force'
1860
          mask_name = tmpname(1:i3)
1861
        END IF
1862
      END IF
1863
            
1864
      str = TRIM( ComponentName( 'exported variable', l ) ) // ' DOFs'
1865
      DOFs = ListGetInteger( SolverParams, str, Found )
1866
      IF ( .NOT. Found ) THEN
1867
        j = 0
1868
        DOFs = 0
1869
        DO WHILE( .TRUE. )
1870
          i = INDEX( var_name(j+1:), ':' ) + j
1871
          IF ( i<=j ) EXIT
1872
          READ( var_name(i+1:),'(i1)',IOSTAT=iostat) k
1873
          IF(iostat /= 0) THEN
1874
            CALL Fatal(Caller,'Could not read component dofs for exported variable '//I2S(l))
1875
          END IF
1876
          DOFs = DOFs + k
1877
          j = i + 1
1878
        END DO
1879
      END IF
1880

1881
            
1882
      VariableOutput = .TRUE.
1883
      VariableGlobal = .FALSE.
1884
      VariableIp = .FALSE.      
1885
      VariableElem = .FALSE.
1886
      VariableDG = .FALSE.
1887
      VariableNodal = .FALSE.
1888
      VariableType = Solver % Variable % TYPE
1889
      InheritVarType = .FALSE.
1890
            
1891
      DO WHILE( var_name(1:1) == '-' )
1892

1893
        k = 0
1894
        j = LEN_TRIM( var_name )
1895

1896
        IF ( SEQL(var_name, '-nooutput ') ) THEN
1897
          VariableOutput = .FALSE.
1898
          k = 10
1899
          
1900
        ! Different types of variables: global, ip, elem, dg
1901
        ELSE IF ( SEQL(var_name, '-global ') ) THEN
1902
          VariableGlobal = .TRUE.
1903
          k = 8
1904
          
1905
        ELSE IF ( SEQL(var_name, '-ip ') ) THEN
1906
          VariableIp = .TRUE.
1907
          k = 4
1908

1909
        ELSE IF ( SEQL(var_name, '-elem ') ) THEN
1910
          VariableElem = .TRUE.
1911
          k = 6
1912

1913
        ELSE IF ( SEQL(var_name, '-nodal ') ) THEN
1914
          VariableNodal = .TRUE.
1915
          k = 7
1916
          
1917
        ELSE IF ( SEQL(var_name, '-dg ') ) THEN
1918
          VariableDG = .TRUE.
1919
          k = 4
1920

1921
        ELSE IF ( SEQL(var_name, '-dofs ') ) THEN
1922
          READ( var_name(7:), *, IOSTAT=iostat ) DOFs 
1923
          IF(iostat /= 0) THEN
1924
            CALL Fatal(Caller,'Could not read -dofs parameter for exported variable '//I2S(l))
1925
          END IF
1926
          !j = LEN_TRIM( var_name )
1927
          k = 7
1928
          DO WHILE( var_name(k:k) /= ' '  )
1929
            k = k + 1
1930
            IF ( k > j ) EXIT
1931
          END DO
1932
        ELSE
1933
          CALL Fatal(Caller,'Do not know how to parse: '//TRIM(var_name))          
1934
        END IF
1935

1936
        ! Remove the processed section. 
1937
        IF(k>0) THEN
1938
          var_name(1:j-k) = var_name(k+1:j)
1939
          DO i=j-k+1,j
1940
            var_name(i:i) = ' '
1941
          END DO
1942
        END IF
1943
        
1944
      END DO
1945
      IF ( DOFs == 0 ) DOFs = 1
1946

1947
      NewVariable => VariableGet( Solver % Mesh % Variables, Var_name )
1948
    
1949
      IF ( .NOT. ASSOCIATED(NewVariable) ) THEN
1950
        CALL Info(Caller,'Creating exported variable: '//TRIM(var_name),Level=12)
1951

1952
        IF( NoPerm ) THEN
1953
          IF(.NOT. (VariableNodal .OR. VariableElem .OR. VariableDG ) ) THEN
1954
            CALL Fatal(Caller,'Invalid type for Noperm variable!')
1955
          END IF
1956
        END IF
1957

1958
        
1959
        IF( VariableIp ) THEN
1960
          VariableType = Variable_on_gauss_points
1961

1962
          IF( UseMask ) THEN
1963
            NULLIFY( Perm ) 
1964
            CALL CreateIpPerm( Solver, Perm, mask_name, sec_name ) 
1965
            nsize = MAXVAL( Perm ) 
1966
          ELSE
1967
            ! Create a table showing the offset for IPs within elements
1968
            CALL CreateIpPerm( Solver )             
1969
            nSize = Solver % IpTable % IpCount
1970
            Perm => Solver % IpTable % IpOffset
1971
          END IF
1972
          nsize = nsize * DOFs
1973
            
1974
        ELSE IF( VariableElem ) THEN
1975
          VariableType = Variable_on_elements
1976

1977
          ! We need to call these earlier than otherwise
1978
          IF( UseMask ) THEN            
1979
            CALL CreateElementsPerm( Solver, Perm, nsize, Mask_Name, sec_name ) 
1980
          ELSE IF( NoPerm ) THEN
1981
            nsize = Solver % Mesh % NumberOfBulkElements
1982
            Perm => NULL()
1983
          ELSE
1984
            CALL SetActiveElementsTable( CurrentModel, Solver, k, CreateInv = .TRUE.  )             
1985
            nSize = Solver % NumberOfActiveElements          
1986
            Perm => Solver % InvActiveElements
1987
          END IF
1988
          nSize = nSize * Dofs
1989
          CALL ListAddInteger( Solver % Values, 'Active Mesh Dimension', k )
1990
            
1991
        ELSE IF( VariableDG ) THEN
1992
          VariableType = Variable_on_nodes_on_elements
1993

1994
          NULLIFY( Perm ) 
1995
          IF( UseMask ) THEN
1996
            CALL CreateDGPerm( Solver, Perm, nsize, mask_name, sec_name ) 
1997
          ELSE IF( NoPerm ) THEN
1998
            nsize = 0
1999
            DO j=1, Solver % Mesh % NumberOfBulkElements
2000
              nsize = nsize + Solver % Mesh % Elements(j) % TYPE % NumberOfNodes
2001
            END DO
2002
            Perm => NULL()
2003
          ELSE
2004
            CALL CreateDGPerm( Solver, Perm, nsize )
2005
          END IF
2006
          nsize = nsize * DOFs
2007

2008
        ELSE IF( VariableNodal ) THEN
2009
          VariableType = Variable_on_nodes
2010
          NULLIFY( Perm ) 
2011

2012
          IF( UseMask ) THEN
2013
            CALL MakePermUsingMask( CurrentModel, Solver, Solver % Mesh, Mask_Name, &
2014
                .TRUE., Perm, nsize )
2015
            nsize = DOFs * nsize
2016
          ELSE IF( NoPerm ) THEN
2017
            nsize = Solver % Mesh % NumberOfNodes
2018
            Perm => NULL()
2019
          ELSE
2020
            ! Just simple permutation
2021
            CALL CreateNodalPerm( Solver, Perm, nSize )
2022
          END IF
2023
          
2024
        ELSE IF( VariableGlobal ) THEN
2025
          VariableType = Variable_global
2026
          nSize = DOFs
2027
          NULLIFY( Perm )
2028

2029
        ELSE ! Follow the primary type
2030
          InheritVarType = .TRUE.
2031
          IF( UseMask ) THEN
2032
            NULLIFY( Perm )
2033
            IF( GotSecName ) THEN
2034
              CALL CreateMaskedPerm( Solver, Solver % Variable % Perm, Mask_Name, &
2035
                  Perm, nsize, sec_name )
2036
            ELSE
2037
              CALL CreateMaskedPerm( Solver, Solver % Variable % Perm, Mask_Name, &
2038
                  Perm, nsize )
2039
            END IF              
2040
            nsize = DOFs * nsize
2041
          ELSE
2042
            nSize = DOFs * SIZE(Solver % Variable % Values) / Solver % Variable % DOFs          
2043
            Perm => Solver % Variable % Perm
2044
          END IF
2045
        END IF
2046
        
2047
        ALLOCATE( Solution(nSize), STAT = AllocStat )
2048
        IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for Solution')
2049
        
2050
        Solution = 0.0d0
2051
        IF( ASSOCIATED(Perm) ) THEN
2052
          CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
2053
              TRIM(var_name), DOFs, Solution, Perm, &
2054
              Output=VariableOutput, TYPE=VariableType )
2055
        ELSE          
2056
          CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
2057
              TRIM(var_name), DOFs, Solution, TYPE=VariableType )
2058
        END IF
2059
        NewVariable => VariableGet( Solver % Mesh % Variables, Var_name )
2060
        IF(ASSOCIATED( NewVariable ) ) THEN
2061
          CALL Info(Caller,'Succesfully created variable: '&
2062
              //ComponentName(var_name),Level=12)                    
2063
        ELSE
2064
          CALL Warn(Caller,'Could not create variable: '//TRIM(var_name))
2065
        END IF
2066
        
2067
        IF( InheritVarType ) NewVariable % PeriodicFlipActive = Solver % PeriodicFlipActive
2068
        
2069
        str = TRIM(ComponentName(Var_name ))//' Transient'
2070
        TransientVar = ListGetLogical( SolverParams, str, Found )        
2071
        
2072
        IF(.NOT. Found ) THEN
2073
          str = TRIM( ComponentName(Var_name) )//' Calculate Velocity'
2074
          TransientVar = ListGetLogical( SolverParams, str, Found )        
2075
        END IF
2076
        IF(.NOT. Found ) THEN
2077
          str = TRIM( ComponentName( Var_name ) )//' Calculate Acceleration'
2078
          TransientVar = ListGetLogical( SolverParams, str, Found )        
2079
        END IF
2080
        
2081
        IF( TransientVar ) THEN
2082
          n = 2
2083
          CALL Info(Caller,'Allocating prevvalues of size 2 for exported variable',Level=12)
2084
          ALLOCATE(NewVariable % PrevValues(nsize,n))
2085
          NewVariable % PrevValues = 0.0_dp
2086

2087
          ! Create upon request the variables for the external fields
2088
          CALL CreateTimeDerivativeVariables( Solver, NewVariable )
2089
        END IF
2090
          
2091
        IF ( DOFs > 1 ) THEN
2092
          n = LEN_TRIM( var_name )
2093
          DO j=1,DOFs
2094
            tmpname = ComponentName( var_name(1:n), j )
2095
            Component => Solution( j::DOFs )
2096

2097
            IF( ASSOCIATED(Perm) ) THEN
2098
              CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
2099
                  TRIM(tmpname), 1, Component, Perm,  &
2100
                  Output=VariableOutput, TYPE = VariableType )
2101
            ELSE
2102
              CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
2103
                  TRIM(tmpname), 1, Component, TYPE = VariableType )
2104
            END IF
2105
            NewVariable => VariableGet( Solver % Mesh % Variables, tmpname )
2106
            IF(ASSOCIATED( NewVariable ) ) THEN
2107
              CALL Info(Caller,'Succesfully created variable: '//TRIM(tmpname),Level=12)          
2108
            ELSE
2109
              CALL Warn(Caller,'Could not create variable: '//TRIM(tmpname))
2110
            END IF
2111

2112
            IF( InheritVarType ) NewVariable % PeriodicFlipActive = Solver % PeriodicFlipActive
2113
          END DO
2114
        END IF
2115

2116
      END IF
2117
    END DO
2118

2119
    IF(Doit) THEN
2120
      IF( ASSOCIATED( Solver % Matrix ) ) THEN
2121
        ALLOCATE(Solver % Matrix % MassValues(SIZE(Solver % Matrix % Values)));
2122
        Solver % Matrix % MassValues=0._dp
2123
      END IF
2124
    END IF
2125
     
2126
    Solver % LinBeforeProc = 0
2127
    str = ListGetString( Solver % Values, 'Before Linsolve', Found )
2128
    IF ( Found ) Solver % LinBeforeProc = GetProcAddr( str )
2129

2130
    Solver % LinAfterProc = 0
2131
    str = ListGetString( Solver % Values, 'After Linsolve', Found )
2132
    IF ( Found ) Solver % LinAfterProc = GetProcAddr( str )
2133

2134
    IF( ASSOCIATED( Solver % Matrix ) ) THEN
2135
      Solver % Matrix % MatVecSubr = 0
2136
      str = ListGetString( Solver % Values, 'Matrix Vector Proc', Found )
2137
      IF ( Found ) Solver % Matrix % MatVecSubr = GetProcAddr( str )
2138
    END IF
2139

2140
  END SUBROUTINE AddEquationBasics
2141

2142

2143

2144
  ! Create 1st and 2nd derirvatives of the solution either for postprocessing, 
2145
  ! dependencies, or for higher order restarts. The optional argument is intended
2146
  ! for exported variables.
2147
  !---------------------------------------------------------------------------
2148
  SUBROUTINE CreateTimeDerivativeVariables( Solver, Var )
2149
    TYPE(Solver_t), POINTER :: Solver
2150
    TYPE(Variable_t), POINTER, OPTIONAL :: Var
2151

2152
    TYPE(Variable_t), POINTER :: pVar
2153
    REAL(KIND=dp), POINTER :: Component(:)
2154
    INTEGER :: k
2155
    LOGICAL :: Found, DoIt
2156
    INTEGER :: TimeOrder
2157
    CHARACTER(:), ALLOCATABLE :: str,kword
2158
    
2159
    IF( PRESENT( Var ) ) THEN
2160
      pVar => Var
2161
      kword = TRIM(ComponentName(Var % Name))//' Calculate Velocity'
2162
    ELSE
2163
      pVar => Solver % Variable
2164
      kword = 'Calculate Velocity'
2165
    END IF
2166
    DoIt = ListGetLogical( Solver % Values, kword, Found)
2167
    IF(.NOT. DoIt) THEN
2168
      IF( PRESENT( Var ) ) THEN
2169
        kword = TRIM(ComponentName(Var % Name))//' Nonlinear Calculate Velocity'
2170
      ELSE
2171
        kword = 'Nonlinear Calculate Velocity'
2172
      END IF
2173
    END IF
2174
    DoIt = ListGetLogical( Solver % Values, kword, Found)
2175
         
2176
    IF( DoIt ) THEN
2177
      ! For exported variable we assume 1st order scheme
2178
      IF( PRESENT( Var ) ) THEN
2179
        TimeOrder = 1        
2180
      ELSE
2181
        TimeOrder = Solver % TimeOrder
2182
      END IF
2183

2184
      k = INDEX(pVar % Name,'[')-1
2185
      IF( k > 0 ) THEN
2186
        str = pVar % Name(1:k)//' Velocity'
2187
      ELSE
2188
        str = TRIM(pVar % Name)//' Velocity'
2189
      END IF
2190
      
2191
      IF( TimeOrder < 1 ) THEN
2192
        CALL Warn('CreateTimeDerivativeVariables',&
2193
            'Velocity computation implemented only for time-dependent variables')
2194
      ELSE IF ( TimeOrder == 1 ) THEN
2195
        CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, Solver, &
2196
            str, pVar % Dofs, Perm = pVar % Perm, VarType = pVar % TYPE )
2197
      ELSE IF ( Solver % TimeOrder >= 2 ) THEN
2198
        Component => pVar % PrevValues(:,1)
2199
        CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, Solver, &
2200
            str, pVar % Dofs, Component, pVar % Perm, Secondary = .TRUE., VarType = pVar % Type )
2201
      END IF
2202
    END IF
2203
   
2204
    IF( PRESENT( Var ) ) THEN
2205
      kword = TRIM(ComponentName(Var % Name))//' Calculate Acceleration'
2206
    ELSE
2207
      kword = 'Calculate Acceleration'
2208
    END IF
2209
    DoIt = ListGetLogical( Solver % Values, kword, Found)
2210
   
2211
    IF( DoIt ) THEN
2212
      IF( TimeOrder < 1 ) THEN
2213
        CALL Warn('CreateTimeDerivativeVariables',&
2214
            'Acceleration computation implemented only for time-dependent variables')      
2215
      ELSE IF ( TimeOrder == 1 ) THEN
2216
        str = TRIM(ComponentName(pVar % Name))//' Acceleration'
2217
        CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, Solver, &
2218
            str, Solver % Variable % Dofs, Perm = Solver % Variable % Perm, VarType = pVar % Type )
2219
      ELSE IF ( TimeOrder >= 2 ) THEN
2220
        Component => Solver % Variable % PrevValues(:,2)
2221
        str = TRIM( ComponentName( pVar % Name ) ) // ' Acceleration'
2222
        CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, Solver, &
2223
            str, pVar % Dofs, Component, pVar % Perm, Secondary = .TRUE., VarType = pVar % Type )
2224
      END IF
2225
    END IF
2226

2227
    IF( PRESENT( Var ) ) THEN
2228
      kword = TRIM(ComponentName(Var % Name))//' Transient Restart'
2229
    ELSE
2230
      kword = 'Transient Restart'
2231
    END IF    
2232
    DoIt = ListGetLogical( Solver % Values, kword, Found)
2233
        
2234
    IF( DoIt ) THEN
2235
      IF( .NOT. ASSOCIATED( pVar % PrevValues ) ) THEN
2236
        CALL Warn('CreateTimeDerivativeVariables',&
2237
            'Transient restart requires PrevValues!')
2238
      ELSE 
2239
        DO k = 1, SIZE( pVar % PrevValues, 2 )
2240
          Component => pVar % PrevValues(:,k)
2241
          str = TRIM( pVar % Name ) !//' PrevValues'//I2S(k)          
2242
          CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, Solver, &
2243
              str, pVar % Dofs, Component, pVar % Perm, Secondary = .TRUE., &
2244
              VarType = pvar % TYPE, VarSuffix = 'PrevValues'//I2S(k))
2245
        END DO
2246
      END IF
2247
    END IF
2248

2249
  END SUBROUTINE CreateTimeDerivativeVariables
2250
      
2251

2252
!------------------------------------------------------------------------------  
2253
!> Add information that is typically only needed if there's a matrix equation
2254
!> to work with. This should be called only after both the solution vector and
2255
!> matrix have been created.
2256
!------------------------------------------------------------------------------
2257
  SUBROUTINE AddEquationSolution(Solver, Transient )
2258
!------------------------------------------------------------------------------
2259
    TYPE(Solver_t), POINTER :: Solver
2260
    LOGICAL :: Transient
2261
!------------------------------------------------------------------------------
2262
    TYPE(Variable_t), POINTER :: Var
2263
    CHARACTER(LEN=MAX_NAME_LEN) :: str, var_name, tmpname
2264
    INTEGER :: i,j,k,n,m,nrows,DOFs
2265
    REAL(KIND=dp), POINTER :: Solution(:)
2266
    INTEGER, POINTER :: Perm(:)
2267
    LOGICAL :: Found, Stat, ComplexFlag, HarmonicAnal, EigAnal, &
2268
         VariableOutput, MGAlgebraic,MultigridActive, HarmonicMode, DoIt
2269
    INTEGER :: MgLevels, AllocStat
2270
    REAL(KIND=dp), POINTER :: Component(:)
2271
    REAL(KIND=dp), POINTER :: freqv(:,:)
2272
    REAL(KIND=dp) :: Freq
2273
    TYPE(Mesh_t),   POINTER :: NewMesh,OldMesh
2274
    TYPE(Element_t), POINTER :: CurrentElement
2275
    TYPE(Matrix_t), POINTER :: NewMatrix, Oldmatrix, SaveMatrix
2276
    CHARACTER(*), PARAMETER :: Caller="AddEquationSolution"   
2277

2278
    !------------------------------------------------------------------------------
2279
    Solver % DoneTime = 0
2280
    IF ( .NOT. ASSOCIATED( Solver % Variable ) ) RETURN
2281
    IF ( .NOT. ASSOCIATED( Solver % Variable % Values ) ) RETURN
2282
    IF (SIZE(Solver % Variable % Values)==0 ) THEN
2283
        DEALLOCATE(Solver % Variable % Values)
2284
        Solver % Variable % Values => Null(); RETURN
2285
    END IF
2286
    !------------------------------------------------------------------------------
2287
	
2288
    !------------------------------------------------------------------------------
2289
    ! If soft limiters are applied then also loads must be computed
2290
    !------------------------------------------------------------------------------
2291
    IF( ListGetLogical( Solver % Values,'Calculate Boundary Fluxes',Found) ) THEN
2292
      CALL ListAddLogical( Solver % Values,'Calculate Loads',.TRUE.)
2293
    END IF
2294
	    
2295
    !------------------------------------------------------------------------------
2296
    ! Create the variable needed for the computation of nodal loads and
2297
    ! residual: r=b-Ax. The difference here is at what stage A and b are stored.     
2298
    !------------------------------------------------------------------------------
2299
    DO k=1,2
2300
      IF(k==1) THEN
2301
        str = 'loads'
2302
      ELSE
2303
        str = 'residual'
2304
      END IF
2305

2306
      IF ( ListGetLogical( Solver % Values,'Calculate '//TRIM(str), Found ) ) THEN
2307
        Var_name = GetVarName(Solver % Variable) // ' '//TRIM(str)
2308
        Var => VariableGet( Solver % Mesh % Variables, var_name, .TRUE. )
2309
        IF ( .NOT. ASSOCIATED(Var) ) THEN
2310
          ALLOCATE( Solution(SIZE(Solver % Variable % Values)), STAT = AllocStat )
2311
          IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for '//TRIM(str))
2312

2313
          DOFs = Solver % Variable % DOFs
2314
          Solution = 0.0d0
2315
          nrows = SIZE( Solution ) 
2316
          Perm => Solver % Variable % Perm
2317

2318
          VariableOutput = ListGetLogical( Solver % Values,'Save '//TRIM(str),Found )
2319
          IF( .NOT. Found ) VariableOutput = Solver % Variable % Output
2320

2321
          CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
2322
              var_name, Solver % Variable % DOFs, Solution, &
2323
              Solver % Variable % Perm, Output=VariableOutput, TYPE = Solver % Variable % TYPE )
2324

2325
          IF ( DOFs > 1 ) THEN
2326
            n = LEN_TRIM( Var_name )
2327
            DO j=1,DOFs
2328
              tmpname = ComponentName( var_name(1:n), j )
2329
              Component => Solution( j:nRows-DOFs+j:DOFs )
2330
              CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
2331
                  tmpname, 1, Component, Perm, Output=VariableOutput, TYPE = Solver % Variable % TYPE )
2332
            END DO
2333
          END IF
2334
          NULLIFY( Solution )
2335
        END IF
2336
      END IF
2337
    END DO
2338

2339
    !------------------------------------------------------------------------------
2340
    ! Optionally create variable for saving permutation vector.
2341
    ! This is mainly for debugging purposes and is therefore commented out.
2342
    !------------------------------------------------------------------------------
2343
#if 0
2344
    DoIt = .FALSE.
2345
    IF( ASSOCIATED( Solver % Variable ) ) THEN
2346
      DoIt = ListGetLogical( Solver % Values,'Save Global Dofs', Found ) 
2347
      IF( Solver % Variable % TYPE /= Variable_on_nodes_on_elements  ) THEN              
2348
        CALL Info(Caller,'Saving of global dof indexes only for DG variable!')
2349
        DoIt = .FALSE.
2350
      END IF
2351
    END IF
2352

2353
    IF( DoIt ) THEN
2354
      Var_name = GetVarName(Solver % Variable) // ' GDofs'      
2355
      Var => VariableGet( Solver % Mesh % Variables, var_name )
2356
      
2357
      IF ( .NOT. ASSOCIATED(Var) ) THEN
2358
        n = SIZE(Solver % Variable % Values) / Solver % Variable % Dofs
2359
        ALLOCATE( Solution(n), STAT = AllocStat )
2360
        Solution = 0.0_dp
2361
        IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error Gdofs')
2362
        
2363
        Solution = 0.0d0
2364
        Perm => Solver % Variable % Perm
2365
        
2366
        CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
2367
            var_name, 1, Solution, Solver % Variable % Perm, TYPE = Solver % Variable % TYPE )
2368
        
2369
        IF( ParEnv % PEs == 1 ) THEN
2370
          DO i = 1, SIZE( Solution )
2371
            Solution(i) = 1.0_dp * i
2372
          END DO
2373
        END IF        
2374
        NULLIFY( Solution )
2375
      END IF
2376
    END IF
2377
#endif    
2378
    
2379
    ! Create a additional variable for the limiters. For elasticity, for example
2380
    ! the variable will be the contact load. 
2381
    IF ( ListGetLogical( Solver % Values,'Apply Limiter', Found ) .OR. &
2382
        ListGetLogical( Solver % Values,'Apply Contact BCs',Found ) ) THEN
2383
      Var_name = GetVarName(Solver % Variable) // ' Contact Load'
2384
      Var => VariableGet( Solver % Mesh % Variables, var_name )
2385
      IF ( .NOT. ASSOCIATED(Var) ) THEN
2386
        ALLOCATE( Solution(SIZE(Solver % Variable % Values)), STAT = AllocStat)
2387
        IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for Contact Loads')
2388

2389
        DOFs = Solver % Variable % DOFs
2390
        Solution = 0.0d0
2391
        nrows = SIZE( Solution ) 
2392
        Perm => Solver % Variable % Perm
2393
        VariableOutput = Solver % Variable % Output
2394

2395
        CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
2396
            var_name, Solver % Variable % DOFs, Solution, &
2397
            Solver % Variable % Perm, Output=VariableOutput, Type = Solver % Variable % Type )
2398

2399
        IF ( DOFs > 1 ) THEN
2400
          n = LEN_TRIM( Var_name )
2401
          DO j=1,DOFs
2402
            tmpname = ComponentName( var_name(1:n), j )
2403
            Component => Solution( j:nRows-DOFs+j:DOFs )
2404
            CALL VariableAdd( Solver % Mesh % Variables, Solver % Mesh, Solver,&
2405
                tmpname, 1, Component, Perm, Output=VariableOutput, Type = Solver % Variable % Type )
2406
          END DO
2407
        END IF
2408
        NULLIFY( Solution )
2409
      END IF
2410
    END IF
2411

2412
    Solver % NOFEigenValues = 0
2413
    Solver % MultiGridLevel = 1
2414
    Solver % MultiGridTotal = 0
2415
    Solver % MultiGridSolver = .FALSE.
2416
    Solver % MultiGridEqualSplit = .FALSE.
2417

2418

2419
    HarmonicAnal = ListGetLogical( Solver % Values, 'Harmonic Analysis', Found )
2420

2421
    HarmonicMode = ListGetLogical( Solver % Values, 'Harmonic Mode', Found )
2422
    
2423
    IF ( ASSOCIATED( Solver % Matrix ) ) THEN
2424
      IF(.NOT. ASSOCIATED(Solver % Matrix % RHS)) THEN
2425
        ALLOCATE( Solver % Matrix % RHS(Solver % Matrix % NumberOFRows), STAT=AllocStat )
2426
        IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for Rhs')
2427
        Solver % Matrix % RHS = 0.0d0
2428
        
2429
        Solver % Matrix % RHS_im => NULL()
2430
        IF ( HarmonicAnal .OR. HarmonicMode ) THEN
2431
          ALLOCATE( Solver % Matrix % RHS_im(Solver % Matrix % NumberOFRows), STAT=AllocStat)
2432
          IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for Rhs_im')                  
2433
          Solver % Matrix % RHS_im = 0.0d0
2434
        END IF
2435
      END IF
2436
    END IF
2437
!------------------------------------------------------------------------------
2438

2439
    EigAnal = ListGetLogical( Solver % Values, 'Eigen Analysis', Found )
2440

2441
    IF ( Transient .AND. .NOT. EigAnal .AND. .NOT. HarmonicAnal ) THEN
2442
      k = ListGetInteger( Solver % Values, 'Time Derivative Order', Found, &
2443
          minv=0, maxv=2 )
2444
      Solver % TimeOrder = 1
2445
      IF ( Found ) Solver % TimeOrder = MIN(MAX(1,k),2)
2446
      
2447
      IF ( ASSOCIATED( Solver % Matrix ) ) THEN
2448
        ALLOCATE( Solver % Matrix % Force(Solver % Matrix % NumberOFRows, &
2449
            Solver % TimeOrder+1), STAT=AllocStat)
2450
        IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for Force')
2451
        Solver % Matrix % Force = 0.0d0
2452
      END IF
2453
      
2454
      IF ( .NOT. ASSOCIATED( Solver % Variable % PrevValues ) ) THEN
2455

2456
        n = SIZE(Solver % Variable % Values)
2457
        IF ( Solver % TimeOrder == 2 ) THEN
2458
          m = 7
2459
        ELSE 
2460
          m = MAX( Solver % Order, Solver % TimeOrder )
2461
        END IF
2462

2463
        IF( m > 0 ) THEN
2464
          ALLOCATE(Solver % Variable % PrevValues( n, m ), STAT=AllocStat )
2465
          IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for PrevValues')
2466
          Solver % Variable % PrevValues = 0.0d0
2467
          
2468
          IF ( Solver % Variable % DOFs > 1 ) THEN
2469
            IF ( GetVarName( Solver % Variable ) == 'flow solution' ) THEN
2470
              DO k=1,Solver % Variable % DOFs-1
2471
                str = 'Velocity ' // CHAR(k+ICHAR('0'))
2472
                Var => VariableGet( Solver % Mesh % Variables, str, .TRUE. )
2473
                IF(.NOT. ASSOCIATED(Var)) THEN
2474
                  CALL Fatal(Caller,&
2475
                      "Failed to get variable: "//TRIM(str)//". Try specifying Variable =&
2476
                      & Flow Solution[velocity:DIM pressure:1] in .sif")
2477
                END IF
2478
                Var % PrevValues =>  &
2479
                    Solver % Variable % PrevValues(k::Solver % Variable % DOFs,:)
2480
              END DO
2481
              Var => VariableGet( Solver % Mesh % Variables, 'Pressure', .TRUE. )
2482
              IF(.NOT. ASSOCIATED(Var)) THEN
2483
                CALL Fatal(Caller,"Failed to get variable: &
2484
                    &pressure. Try specifying Variable = Flow Solution &
2485
                    &[velocity:DIM pressure:1] in .SIF")
2486
              END IF
2487
              Var % PrevValues =>  &
2488
                  Solver % Variable % PrevValues(k::Solver % Variable % DOFs,:)
2489
            ELSE
2490
              DO k=1,Solver % Variable % DOFs
2491
                str = ComponentName( Solver % Variable % Name, k ) 
2492
                Var => VariableGet( Solver % Mesh % Variables, str, .TRUE. )
2493
                IF( ASSOCIATED( Var ) ) THEN
2494
                  Var % PrevValues =>  &
2495
                      Solver % Variable % PrevValues(k::Solver % Variable % DOFs,:)
2496
                END IF
2497
              END DO
2498
            END IF
2499
          END IF
2500
        END IF
2501
      END IF
2502

2503
      ! Create 1st and 2nd derirvatives of the solution either for postprocessing, 
2504
      ! dependencies, or for higher order restarts
2505
      CALL CreateTimeDerivativeVariables( Solver )
2506
      
2507
    ELSE
2508
      Solver % TimeOrder = 0
2509

2510
      DoIt = ListGetLogical( Solver % Values,'Calculate Derivative',Found)
2511
      IF( .NOT. Found ) THEN
2512
        DoIt = ListGetLogical( Solver % Values,'Nonlinear Calculate Derivative',Found)
2513
      END IF
2514
      IF( DoIt ) THEN
2515
        str = GetVarname(Solver % Variable) // ' Derivative'
2516
        CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, Solver, &
2517
            str, Solver % Variable % Dofs, Perm = Solver % Variable % Perm, &
2518
            VarType = Solver % Variable % TYPE )
2519
        Var => VariableGet( Solver % Mesh % Variables, str, .TRUE. )
2520
        Var % Values = 0.0_dp
2521
      END IF
2522

2523
      IF ( EigAnal ) THEN
2524
        n = ListGetInteger( Solver % Values, 'Eigen System Values', Found )
2525
        IF ( Found .AND. n > 0 ) THEN
2526
          IF (ListGetLogical(Solver % Values, 'Linear System Skip Complex', Found)) THEN
2527
            ComplexFlag = .FALSE.
2528
          ELSE
2529
            ComplexFlag = ListGetLogical(Solver % Values, 'Linear System Complex', Found)
2530
          END IF
2531
          Solver % NOFEigenValues = n
2532
          IF ( .NOT. ASSOCIATED( Solver % Variable % EigenValues ) ) THEN
2533
            ALLOCATE( Solver % Variable % EigenValues(n) )
2534
            IF (ComplexFlag) THEN
2535
              ALLOCATE( Solver % Variable % EigenVectors(n, &
2536
                  SIZE( Solver % Variable % Values )/2 ), STAT=AllocStat)
2537
            ELSE
2538
              ALLOCATE( Solver % Variable % EigenVectors(n, &
2539
                  SIZE( Solver % Variable % Values ) ), STAT=AllocStat)
2540
            END IF
2541
            IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for EigenValues')
2542
           
2543
            Solver % Variable % EigenValues  = 0.0d0
2544
            Solver % Variable % EigenVectors = 0.0d0
2545

2546
            IF( .NOT. ComplexFlag .AND. Solver % Variable % DOFs > 1 ) THEN
2547
              CALL Info(Caller,'Repointing '//I2S(Solver % Variable % DOFs)//&
2548
                  ' eigenvalue components for: '//TRIM(Solver % Variable % Name))
2549
              
2550
              DO k=1,Solver % Variable % DOFs
2551
                str = ComponentName( Solver % Variable % Name, k )
2552
                Var => VariableGet( Solver % Mesh % Variables, str, .TRUE. )
2553
                
2554
                IF( ASSOCIATED( Var ) ) THEN
2555
                  CALL Info(Caller,'Eigenvalue component '&
2556
                      //I2S(k)//': '//TRIM(str))
2557
                  Var % EigenValues => Solver % Variable % EigenValues
2558
                  Var % EigenVectors =>  & 
2559
                      Solver % Variable % EigenVectors(:,k::Solver % Variable % DOFs )
2560
                END IF
2561
              END DO
2562
            END IF
2563
          END IF
2564
            
2565
          ALLOCATE( Solver % Matrix % MassValues(SIZE(Solver % Matrix % Values)), STAT=AllocStat)
2566
          IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for MassValues')
2567
          
2568
          Solver % Matrix % MassValues = 0.0d0
2569
        END IF
2570
      ELSE IF ( HarmonicAnal ) THEN
2571

2572
        n = ListGetInteger( Solver % Values,'Harmonic System Values',Found )
2573
        IF( n > 1 ) THEN
2574
          freqv => ListGetConstRealArray( Solver % Values, 'Frequency', Found )
2575
          IF(Found ) THEN
2576
            IF( SIZE( Freqv,1) < n ) THEN
2577
              CALL Fatal( Caller, 'Frequency must be at least same size as > Harmonic System Values <')
2578
            END IF
2579
          ELSE
2580
            CALL Fatal( Caller, '> Frequency < must be given for harmonic analysis.' )
2581
          END IF
2582
        ELSE
2583
          n = 1
2584
        END IF
2585

2586
        Solver % NOFEigenValues = n
2587
        IF ( .NOT. ASSOCIATED( Solver % Variable % EigenValues ) ) THEN
2588
          ALLOCATE( Solver % Variable % EigenValues(n) )
2589
          ALLOCATE( Solver % Variable % EigenVectors(n, &
2590
              SIZE( Solver % Variable % Values ) ), STAT=AllocStat)
2591
          IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for EigenValues')
2592
           
2593
          Solver % Variable % EigenValues  = 0.0d0
2594
          Solver % Variable % EigenVectors = 0.0d0
2595
          
2596
          DO k=1,Solver % Variable % DOFs
2597
            str = ComponentName( Solver % Variable % Name, k )
2598
            Var => VariableGet( Solver % Mesh % Variables, str, .TRUE. )
2599
            IF ( ASSOCIATED( Var ) ) THEN
2600
              Var % EigenValues => Solver % Variable % EigenValues
2601
              Var % EigenVectors => &
2602
                  Solver % Variable % EigenVectors(:,k::Solver % Variable % DOFs)
2603
            END IF
2604
          END DO
2605
        END IF
2606
        
2607
        ALLOCATE( Solver % Matrix % MassValues(SIZE(Solver % Matrix % Values)), STAT=AllocStat)
2608
        IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for MassValues')
2609
        Solver % Matrix % MassValues = 0.0d0
2610

2611
      ELSE IF( HarmonicMode ) THEN
2612

2613
        ALLOCATE( Solver % Matrix % MassValues(SIZE(Solver % Matrix % Values)), STAT=AllocStat)
2614
        IF( AllocStat /= 0 ) CALL Fatal(Caller,'Allocation error for MassValues')
2615
        Solver % Matrix % MassValues = 0.0d0        
2616

2617
      END IF
2618
    END IF
2619

2620

2621
!------------------------------------------------------------------------------
2622

2623
    IF ( ASSOCIATED( Solver % Matrix ) ) THEN
2624
      Solver % Matrix % Symmetric = ListGetLogical( Solver % Values, &
2625
          'Linear System Symmetric', Found )
2626
      
2627
      Solver % Matrix % Lumped = ListGetLogical( Solver % Values, &
2628
          'Lumped Mass Matrix', Found )
2629
      
2630
      MultigridActive = &
2631
          ListGetString( Solver % Values, 'Linear System Solver', Found ) == 'multigrid' .OR. &
2632
          ListGetString( Solver % Values, 'Linear System Preconditioning', Found ) == 'multigrid'
2633

2634

2635
!      Check for multigrid solver:
2636
!      ---------------------------
2637
       IF ( MultigridActive ) THEN
2638

2639
         ! Multigrid may be either solver or preconditioner, it it solver?
2640
         Solver % MultiGridSolver = ListGetString(Solver % Values, &
2641
             'Linear System Solver', Found ) == 'multigrid'
2642
 
2643
         ! There are four different methods: algrabraic, cluster, p and geometric
2644
         str = ListGetString( Solver % Values,'MG Method',Found) 
2645
         IF( Found ) THEN
2646
           MGAlgebraic = ( str == 'algebraic' ) .OR. ( str == 'cluster') .OR. ( str == 'p')
2647
         ELSE    
2648
           MGAlgebraic = ListGetLogical( Solver % Values, 'MG Algebraic', Found ) &
2649
               .OR. ListGetLogical( Solver % Values, 'MG Cluster', Found ) &
2650
               .OR. ListGetLogical( Solver % Values, 'MG Pelement', Found ) 
2651
         END IF
2652
         
2653
         
2654
         MgLevels = ListGetInteger( Solver % Values, &
2655
             'MG Levels', Found, minv=1 )         
2656
         IF ( .NOT. Found ) THEN
2657
           MgLevels = ListGetInteger( Solver % Values, &
2658
               'Multigrid Levels', Found, minv=1 )
2659
         END IF
2660
         IF ( .NOT. Found ) THEN
2661
           IF( MGAlgebraic ) THEN 
2662
             MgLevels = 10
2663
           ELSE
2664
             CALL Fatal(Caller,'> MG Levels < must be defined for geometric multigrid!')
2665
           END IF
2666
         END IF
2667
         Solver % MultiGridTotal = MgLevels
2668
 
2669

2670
!         In case of geometric multigrid make the hierarchical meshes
2671
!         -----------------------------------------------------------
2672
         IF(.NOT. MGAlgebraic ) THEN 
2673

2674
!         Check if h/2 splitting of mesh requested:
2675
!         ------------------------------------------
2676
           Solver % MultiGridEqualSplit = ListGetLogical( &
2677
               Solver % Values, 'MG Equal Split', Found )
2678
         
2679
           IF ( Solver % MultiGridEqualSplit ) THEN
2680
             CALL ParallelInitMatrix( Solver, Solver % Matrix )
2681
             Solver % MultiGridLevel = 1
2682
             
2683
             DO WHILE( Solver % MultiGridLevel < Solver % MultiGridTotal )
2684
               IF ( ASSOCIATED( Solver % Mesh % Child ) ) THEN
2685
                 NewMesh => Solver % Mesh % Child
2686
                 
2687
                 OldMesh   => Solver % Mesh
2688
                 OldMatrix => Solver % Matrix
2689
                 
2690
                 CALL UpdateSolverMesh( Solver, NewMesh )
2691
                 Solver % Mesh % Changed = .FALSE.
2692
               ELSE
2693
                 NewMesh => SplitMeshEqual( Solver % Mesh )
2694
                 CALL SetMeshMaxDofs(NewMesh)
2695
                 NewMesh % Next => CurrentModel % Meshes
2696
                 CurrentModel % Meshes => NewMesh
2697
                 
2698
                 OldMesh   => Solver % Mesh
2699
                 OldMatrix => Solver % Matrix
2700
                 
2701
                 CALL UpdateSolverMesh( Solver, NewMesh )
2702
                 Solver % Mesh % Changed = .FALSE.
2703
                 
2704
                 NewMesh % Parent => OldMesh
2705
                 OldMesh % Child  => NewMesh
2706
                 NewMesh % Name = OldMesh % Name
2707
               END IF
2708

2709
               Newmesh % OutputActive = .TRUE.
2710
               OldMesh % OutputActive = .FALSE.
2711
               
2712
               NewMatrix => Solver % Matrix
2713
               NewMatrix % Parent => OldMatrix
2714
               OldMatrix % Child  => NewMatrix
2715
               CALL ParallelInitMatrix( Solver, Solver % Matrix )
2716
               Solver % MultiGridLevel = Solver % MultiGridLevel + 1
2717
             END DO
2718
           ELSE
2719
             CALL ParallelInitMatrix( Solver, Solver % Matrix )
2720
             OldMesh   => Solver % Mesh
2721
             Var => Solver % Variable
2722
             SaveMatrix => Solver % Matrix
2723
             
2724
             Solver % MultiGridLevel = 1
2725
             DO WHILE( Solver % MultiGridLevel < Solver % MultigridTotal )
2726
               IF ( ASSOCIATED(Solver % Mesh % Parent) ) THEN
2727
                 NewMesh => Solver % Mesh % Parent
2728
                 OldMatrix => Solver % Matrix
2729
                 CALL UpdateSolverMesh( Solver, NewMesh )
2730
                 Solver % Mesh % Changed = .FALSE.
2731
                 NewMatrix => Solver % Matrix
2732
                 NewMatrix % Child => OldMatrix
2733
                 OldMatrix % Parent  => NewMatrix
2734
                 CALL ParallelInitMatrix(Solver, Solver % Matrix )
2735
               END IF
2736
               Solver % MultiGridLevel = Solver % MultiGridLevel+1
2737
             END DO
2738
             
2739
             Solver % Mesh => OldMesh
2740
             Solver % Variable => Var
2741
             Solver % Matrix => SaveMatrix
2742
             CALL SetCurrentMesh(CurrentModel,Solver % Mesh)
2743
           END IF
2744

2745
           CALL MeshStabParams( Solver % Mesh )
2746
         END IF
2747

2748
         Solver % MultiGridLevel = Solver % MultigridTotal
2749

2750
        END IF
2751

2752
       ! Set the default verbosity of the iterative solvers accordingly with the 
2753
       ! global verbosity.
2754
       !-----------------------------------------------------------------------------
2755
       IF( .NOT. ListCheckPresent( Solver % Values,'Linear System Residual Output') ) THEN
2756
         k = 1
2757
         IF( .NOT. OutputLevelMask(4) ) THEN
2758
           k = 0
2759
         ELSE IF( .NOT. OutputLevelMask(5) ) THEN
2760
           k = 10
2761
         END IF
2762
         IF( k /= 1 ) THEN
2763
           CALL ListAddInteger( Solver % Values,'Linear System Residual Output',k)
2764
         END IF
2765
       END IF
2766
     END IF
2767

2768

2769
     DoIt = ASSOCIATED( Solver % Matrix ) 
2770
     IF(DoIt) THEN
2771
       DoIt = InfoActive(20) .OR. ListGetLogical( CurrentModel % Simulation,'Size Info',Found)
2772
     END IF
2773
     IF(DoIt) CALL PrintProblemSize( Solver )
2774
         
2775
!------------------------------------------------------------------------------
2776
   END SUBROUTINE AddEquationSolution
2777
!------------------------------------------------------------------------------
2778

2779
    
2780

2781
!------------------------------------------------------------------------------
2782
!> Generate a similar solver instance as for the parent solver.
2783
!> The number of dofs may vary but the basis functions and permutation is reused.
2784
!> If also the number of dofs is the same also matrix topology is reused.
2785
!------------------------------------------------------------------------------
2786
   FUNCTION CreateChildSolver( ParentSolver, ChildVarName, ChildDofs, ChildPrefix, NoReuse) &
2787
       RESULT ( ChildSolver )
2788
     TYPE(Solver_t) :: ParentSolver
2789
     CHARACTER(LEN=*) :: ChildVarName
2790
     INTEGER, OPTIONAL :: ChildDofs
2791
     CHARACTER(LEN=*), OPTIONAL :: ChildPrefix
2792
     TYPE(Solver_t), POINTER :: ChildSolver
2793
     LOGICAL, OPTIONAL :: NoReuse
2794
     
2795
     INTEGER :: ParentDofs 
2796
     TYPE(Solver_t), POINTER :: Solver
2797
     REAL(KIND=dp), POINTER :: ChildVarValues(:)
2798
     INTEGER, POINTER :: ChildVarPerm(:)
2799
     TYPE(Variable_t), POINTER :: ChildVar
2800
     TYPE(Matrix_t), POINTER :: ChildMat, ParentMat
2801
     INTEGER :: n,m,dofs, i,j,k,l,ii, jj, nn
2802
     LOGICAL :: Found, Lvalue
2803

2804
     ParentDofs = ParentSolver % Variable % Dofs
2805
     IF( PRESENT( ChildDofs ) ) THEN
2806
       Dofs = ChildDofs
2807
     ELSE
2808
       Dofs = ParentDofs
2809
     END IF
2810

2811
     CALL Info('CreateChildSolver','Creating solver of size '//I2S(Dofs)//' for variable: &
2812
         '//TRIM(ChildVarName),Level=6)
2813

2814
     NULLIFY( Solver ) 
2815
     ALLOCATE( Solver )
2816
     ChildSolver => Solver
2817

2818
     Solver % Values => Null()
2819
     CALL ListAddString(Solver % Values,'Equation',TRIM(ChildVarName)//' solver' )
2820

2821
     IF( PRESENT( ChildPrefix ) ) THEN
2822
       CALL Info('CreateChildSolver','Copying keywords with prefix: '//TRIM(ChildPrefix),Level=8)
2823
       CALL ListCopyPrefixedKeywords( ParentSolver % Values, Solver % Values, &
2824
           ChildPrefix )
2825
     ELSE
2826
       CALL Info('CreateChildSolver','Copying all keywords',Level=8)
2827
       CALL ListCopyAllKeywords( ParentSolver % Values, Solver % Values )
2828
     END IF
2829

2830
     IF( .NOT. ASSOCIATED( ParentSolver % Mesh ) ) THEN
2831
       CALL Fatal('CreateChildSolver','Parent solver is missing mesh!')
2832
     END IF
2833
     Solver % Mesh => ParentSolver % Mesh
2834
     i = SIZE(ParentSolver % Def_Dofs,1)
2835
     j = SIZE(ParentSolver % Def_Dofs,2)
2836
     k = SIZE(ParentSolver % Def_Dofs,3)
2837
     ALLOCATE(Solver % Def_Dofs(i,j,k))
2838
     Solver % Def_Dofs = ParentSolver % Def_Dofs
2839

2840
     IF( .NOT. ASSOCIATED( ParentSolver % Variable ) ) THEN
2841
       CALL Fatal('CreateChildSolver','Parent solver is missing variable!')
2842
     END IF
2843

2844
     n = ( SIZE( ParentSolver % Variable % Values ) ) / ParentDofs    
2845
     
2846
     ALLOCATE( ChildVarValues( n * Dofs ) )
2847
     ChildVarValues = 0.0_dp
2848
     ChildVarPerm => ParentSolver % Variable % Perm
2849

2850
     Lvalue = ListGetLogical( Solver % Values,'Variable Output', Found )
2851
     IF(.NOT. Found ) Lvalue = ParentSolver % Variable % Output 
2852
          
2853
     CALL VariableAddVector( Solver % Mesh % Variables, Solver % Mesh, &
2854
         Solver, ChildVarName, Dofs, ChildVarValues, ChildVarPerm, Lvalue )
2855
     
2856

2857
     ChildVar => VariableGet( Solver % Mesh % Variables, ChildVarName )      
2858
     IF(.NOT. ASSOCIATED( ChildVar ) ) THEN
2859
       CALL Fatal('CreateChildSolver','Could not generate child variable!')
2860
     END IF
2861

2862
     ChildVar % Solver => Solver
2863
     
2864
     ChildVar % TYPE = ParentSolver % Variable % TYPE
2865
     Solver % Variable => ChildVar
2866

2867
     
2868
     CALL Info('CreateChildSolver','Creating matrix for solver variable',Level=8)    
2869
     Solver % Matrix => AllocateMatrix()
2870
     ChildMat => Solver % Matrix
2871

2872
     ParentMat => ParentSolver % Matrix
2873
     IF( .NOT. ASSOCIATED( ParentMat ) ) THEN
2874
       CALL Warn('CreateChildSolver','Parent matrix needed for child matrix!')
2875
       Solver % Matrix => NULL()
2876
     ELSE
2877
       ChildMat => CreateChildMatrix( ParentMat, ParentDofs, Dofs, Dofs, .TRUE., NoReuse = NoReuse )
2878
       ChildMat % Solver => Solver
2879
       Solver % Matrix => ChildMat
2880
     END IF
2881

2882
     ChildMat % COMPLEX = ListGetLogical( Solver % Values,'Linear System Complex',Found )
2883

2884
     IF(.NOT. Found ) THEN
2885
       IF( MODULO( ChildDofs, 2 ) == 0 ) THEN
2886
         ChildMat % COMPLEX = ParentMat % COMPLEX
2887
        ELSE
2888
         ChildMat % COMPLEX = .FALSE.
2889
        END IF
2890
     END IF
2891

2892
     Lvalue = ListGetLogical( Solver % Values,'Bubbles in Global System',Found )
2893
     IF( Found ) THEN
2894
       CALL ListAddNewLogical( ChildSolver % Values,'Bubbles in Global System',Lvalue )
2895
     END IF
2896
     
2897
     IF( ASSOCIATED( ParentSolver % ActiveElements ) ) THEN
2898
       Solver % ActiveElements => ParentSolver % ActiveElements
2899
       Solver % NumberOfActiveElements = ParentSolver % NumberOfActiveElements
2900
     END IF
2901

2902
     IF( ASSOCIATED( ParentSolver % ColourIndexList ) ) THEN
2903
       Solver % ColourIndexList => ParentSolver % ColourIndexList
2904
     END IF
2905
       
2906
     Solver % TimeOrder = ListGetInteger( Solver % values,'Time Derivative Order',Found )
2907
     IF(.NOT. Found ) Solver % TimeOrder = ParentSolver % TimeOrder
2908
     
2909
     Solver % MultigridTotal = 0
2910
     Solver % SolverExecWhen = SOLVER_EXEC_NEVER
2911
     Solver % LinBeforeProc = 0
2912
     Solver % LinAfterProc = 0
2913

2914
     IF ( Parenv  % PEs >1 ) THEN
2915
       CALL ParallelInitMatrix( Solver, Solver % Matrix )
2916
     END IF
2917
     
2918
     CALL Info('CreateChildSolver','All done for now!',Level=8)    
2919
   END FUNCTION CreateChildSolver
2920
!------------------------------------------------------------------------------
2921

2922
   
2923
   
2924
!------------------------------------------------------------------------------
2925
!> Solve the equations one-by-one. 
2926
!------------------------------------------------------------------------------
2927
  SUBROUTINE SolveEquations( Model, dt, TransientSimulation, &
2928
      CoupledMinIter, CoupledMaxIter, SteadyStateReached, &
2929
      RealTimestep, BeforeTime, AtTime, AfterTime )
2930
!------------------------------------------------------------------------------
2931
    TYPE(Model_t) :: Model
2932
    REAL(KIND=dp) :: dt
2933
    INTEGER, OPTIONAL :: RealTimestep
2934
    INTEGER :: CoupledMinIter, CoupledMaxIter
2935
    LOGICAL :: TransientSimulation, SteadyStateReached, TransientSolver
2936
    LOGICAL, OPTIONAL :: BeforeTime, AtTime, AfterTime
2937
!------------------------------------------------------------------------------
2938
    REAL(KIND=dp) :: RelativeChange, Tolerance, PrevDT = 0.0d0, Relaxation, SSCond
2939
    INTEGER :: i,j,k,l,n,ierr,istat,Visited=0, RKorder=0, nSolvers
2940
    LOGICAL :: Found, Stat, AbsNorm, Scanning, Convergence, RungeKutta, MeActive, &
2941
        NeedSol, CalculateDerivative, TestConvergence=.FALSE., TestDivergence, DivergenceExit, &
2942
        ExecSlot, CoupledAbort
2943
    LOGICAL, ALLOCATABLE :: DoneThis(:), AfterConverged(:)
2944
    TYPE(Solver_t), POINTER :: Solver
2945
    TYPE(Mesh_t),   POINTER :: Mesh
2946
    CHARACTER(LEN=max_name_len) :: When
2947
    TYPE(Variable_t), POINTER :: IterV, TimestepV
2948
    REAL(KIND=dp), POINTER :: steadyIt,nonlnIt
2949
    REAL(KIND=dp), POINTER :: k1(:), k2(:), k3(:), k4(:), sTime
2950

2951
    TYPE(Variable_t), POINTER :: TimeVar
2952
    REAL(KIND=dp) :: RK2_err
2953
    REAL(KIND=dp), ALLOCATABLE :: RK2_ErrorEstimate(:)
2954
    TYPE RungeKutta_t
2955
      REAL(KIND=dp), ALLOCATABLE :: k1(:),k2(:),k3(:),k4(:)
2956
    END TYPE RungeKutta_t
2957
    TYPE(RungeKutta_t), ALLOCATABLE, TARGET :: RKCoeff(:)
2958

2959
    TYPE(ParEnv_t) :: ParEnv_Save
2960
!------------------------------------------------------------------------------
2961

2962
    ParEnv_Save = ParEnv_Common
2963

2964
!------------------------------------------------------------------------------
2965
!   Initialize equation solvers for new timestep
2966
!------------------------------------------------------------------------------
2967
    nSolvers = Model % NumberOfSolvers
2968

2969
    IF(TransientSimulation) THEN
2970
       CoupledMinIter = ListGetInteger( Model % Simulation, &
2971
            'Steady State Min Iterations', Found )
2972

2973
       CoupledMaxIter = ListGetInteger( Model % Simulation, &
2974
            'Steady State Max Iterations', Found, minv=1 )
2975
       IF ( .NOT. Found ) CoupledMaxIter = 1
2976
    END IF
2977

2978
    Scanning = &
2979
      ListGetString( CurrentModel % Simulation, 'Simulation Type', Found ) == 'scanning'
2980

2981
    TestDivergence = .FALSE.
2982
    DivergenceExit = .FALSE.
2983
    
2984
    IF ( TransientSimulation ) THEN
2985
      DO k=1,nSolvers
2986
        Solver => Model % Solvers(k)
2987
        IF ( Solver % PROCEDURE /= 0 ) THEN
2988
          CALL InitializeTimestep(Solver)
2989
         END IF
2990
      END DO
2991
    END IF
2992

2993
    IF( TransientSimulation .OR. Scanning ) THEN
2994
       IterV => VariableGet(Model % Solvers(1) % Mesh % Variables, 'coupled iter')
2995
       steadyIt => IterV % Values(1)
2996
    END IF
2997
!------------------------------------------------------------------------------
2998

2999
    IF( PRESENT( BeforeTime ) ) THEN
3000
      ExecSlot = BeforeTime
3001
    ELSE
3002
      ExecSlot = .TRUE.
3003
    END IF
3004

3005
    IF( ExecSlot ) THEN
3006
      CALL Info('SolveEquations','Solvers before timestep',Level=12)
3007
      DO k=1,nSolvers
3008
        Solver => Model % Solvers(k)
3009
        IF ( Solver % PROCEDURE==0 ) CYCLE
3010
        IF ( Solver % SolverExecWhen == SOLVER_EXEC_AHEAD_TIME .OR. &
3011
            Solver % SolverExecWhen == SOLVER_EXEC_PREDCORR ) THEN
3012

3013
          IF( PRESENT( RealTimeStep ) ) THEN
3014
            IF( RealTimeStep == 1 ) THEN
3015
              IF( ListGetLogical( Solver % Values,'Skip First Timestep',Found ) ) CYCLE
3016
            END IF
3017
          END IF
3018

3019
          ! Use predictor method
3020
          IF( Solver % SolverExecWhen == SOLVER_EXEC_PREDCORR ) THEN
3021
            CALL Info('SolveEquations','Switching time-stepping method to predictor method',Level=7)
3022
            CALL ListAddString( Solver % Values, 'Timestepping Method', &
3023
                ListGetString( Solver % Values, 'Predictor Method',Found) )
3024
            IF(.NOT. Found ) THEN
3025
              CALL Fatal('SolveEquations','Predictor-corrector schemes require > Predictor Method <')
3026
            END IF
3027
            CALL ListAddLogical( Solver % Values,'Predictor Phase',.TRUE.)
3028
          END IF
3029
          
3030
          CALL SolverActivate( Model,Solver,dt,TransientSimulation )
3031
          CALL ParallelBarrier
3032
          
3033
          ! Use Corrector method
3034
          IF( Solver % SolverExecWhen == SOLVER_EXEC_PREDCORR ) THEN
3035
            CALL Info('SolveEquations','Switching time-stepping method to corrector method',Level=7)
3036
            CALL ListAddString( Solver % Values, 'Timestepping Method', &
3037
                ListGetString( Solver % Values, 'Corrector Method',Found) )
3038
            IF(.NOT. Found ) THEN
3039
              CALL Fatal('SolveEquations','Predictor-corrector schemes require > Corrector Method <')
3040
            END IF
3041
          END IF
3042
        END IF
3043
      END DO
3044
    END IF 
3045

3046
!------------------------------------------------------------------------------
3047
    IF( PRESENT( AtTime ) ) THEN
3048
      ExecSlot = AtTime
3049
    ELSE
3050
      ExecSlot = .TRUE.
3051
    END IF
3052

3053
    IF( ExecSlot ) THEN
3054
      TestDivergence = ListGetLogical( CurrentModel % Simulation, &
3055
          'Convergence Control Within Iterations', Found ) 
3056
      
3057
      ALLOCATE( DoneThis(nSolvers), AfterConverged(nSolvers) )
3058
      
3059
      DO i=1,nSolvers
3060
        Solver => Model % Solvers(i)
3061
        AfterConverged(i) = ListGetLogical( Solver % Values, &
3062
            'Coupled System After Others Converged', Found )
3063
      END DO
3064

3065
!------------------------------------------------------------------------------
3066
      CALL Info('SolveEquations','Solvers in main iteration loop',Level=12)
3067

3068
      TimeVar => VariableGet( Model % Variables, 'Time')
3069
      sTime => TimeVar % Values(1)
3070

3071
      RungeKutta = ListGetString( Model % Simulation, &
3072
          'Timestepping Method', Found ) == 'runge-kutta'
3073

3074
      IF( .NOT. RungeKutta ) THEN
3075
        ! Without Runge-Kutta the cycling over equations is pretty easy
3076
        CALL SolveCoupled()
3077
      ELSE
3078
        CALL Info('SolveEquations','Using Runge-Kutta time-stepping',Level=12)
3079

3080
        CALL Info('SolveEquations','Runge-Kutta predictor step',Level=12)
3081
        sTime = sTime - dt
3082
        CALL SolveCoupled()
3083

3084
        ! Perform Runge-Kutta steps for ru
3085
        !---------------------------------------------------------------
3086

3087
        DO i=1,nSolvers
3088
          Solver => Model % Solvers(i)
3089

3090
          IF ( Solver % PROCEDURE==0 ) CYCLE
3091
          IF ( .NOT. ASSOCIATED( Solver % Variable ) ) CYCLE
3092

3093
          RungeKutta = .FALSE.
3094
          IF ( TransientSimulation .AND. Solver % TimeOrder == 1 ) THEN
3095
            RungeKutta = ListGetString( Solver % Values, &
3096
                'Timestepping Method', Found ) == 'runge-kutta'
3097
          END IF
3098

3099
          IF ( .NOT. RungeKutta ) CYCLE
3100

3101
          CALL Info('SolveEquations','Solver '//I2S(i)//' is Runge-Kutta Solver',Level=12)
3102
          IF ( .NOT. ALLOCATED(RKCoeff) ) THEN
3103
            ALLOCATE(RKCoeff(nSolvers), RK2_ErrorEstimate(nSolvers))
3104
          END IF
3105

3106
          n = SIZE(Solver % Variable % Values)
3107
          ALLOCATE(RKCoeff(i) % k1(n), RKCoeff(i) % k2(n), RKCoeff(i) % k3(n), &
3108
              RKCoeff(i) % k4(n) )
3109

3110
          RKorder = Solver % Order
3111
          k1 => RKCoeff(i) % k1
3112
          k1 = Solver % Variable % Values-Solver % Variable % PrevValues(:,1)
3113
          Solver % Variable % Values = Solver % Variable % PrevValues(:,1) + 2*k1/RKorder
3114
        END DO
3115

3116
        IF ( .NOT. ALLOCATED(RKCoeff) ) THEN        
3117
          CALL Fatal('SolveEquations','No Runge-Kutta after all in any Solver?')
3118
        END IF
3119

3120

3121
        CALL Info('SolveEquations','Using Runge-Kutta Order: '//I2S(RKOrder),Level=12)
3122

3123
        IF(RKorder==4) THEN
3124
          dt = dt / 2
3125
          sTime = sTime + dt
3126
        ELSE
3127
          sTime = sTime + dt
3128
        END IF
3129
        CALL SolveCoupled()
3130

3131
        RK2_ErrorEstimate = 0._dp
3132
        RK2_err = 0.0_dp
3133
        DO i=1,nSolvers
3134
          Solver => Model % Solvers(i)
3135
          IF(.NOT. ASSOCIATED(Solver % Variable)) CYCLE
3136
          IF(.NOT. ASSOCIATED(Solver % Variable % PrevValues) ) CYCLE
3137
          IF(Solver % Variable % Norm < 1.0d-20) CYCLE 
3138
          
3139
          IF ( .NOT. ALLOCATED(RKCoeff(i) % k1)) CYCLE
3140

3141
          k1 => RKCoeff(i) % k1
3142
          k2 => RKCoeff(i) % k2
3143
          k2 = Solver % Variable % Values - Solver % Variable % PrevValues(:,1)
3144
          SELECT CASE(RKorder)
3145
          CASE(2)
3146
            Solver % Variable % Values = Solver % Variable % PrevValues(:,1) + (k1+k2)/2
3147
            RK2_Errorestimate(i) = SUM(((k2-k1)/2)**2)
3148
            RK2_ErrorEstimate(i) = SQRT( ParallelReduction(RK2_ErrorEstimate(i)) ) / &
3149
                Solver % Variable % Norm
3150
            RK2_err = MAX(RK2_err, RK2_ErrorEstimate(i) )
3151
          CASE(4)
3152
            Solver % Variable % Values = Solver % Variable % PrevValues(:,1) + k2
3153
            k2 = 2*k2
3154
          END SELECT
3155
        END DO
3156

3157
        ! Provide error measure for adaptive timestepping = ||RK2 (Heun) - Explicit Euler|| !
3158
        IF ( RKorder == 2 ) THEN
3159
          CALL ListAddConstReal( Model % Simulation, 'Adaptive Error Measure', RK2_err )
3160
        END IF
3161

3162

3163
        ! For 4th order R-K we don't have error estimate but we have more steps to do
3164
        IF( RKOrder == 4 ) THEN
3165
          CALL SolveCoupled()
3166

3167
          DO i=1,nSolvers
3168
            Solver => Model % Solvers(i)
3169
            IF ( .NOT. ALLOCATED(RKCoeff(i) % k1)) CYCLE
3170

3171
            k3 => RKCoeff(i) % k3
3172
            k3 = 2*(Solver % Variable % Values - Solver % Variable % PrevValues(:,1))
3173
            Solver % Variable % Values = Solver % Variable % PrevValues(:,1) + k3
3174
          END DO
3175

3176
          sTime = sTime + dt
3177
          dt = 2 * dt
3178
          CALL SolveCoupled()
3179

3180
          DO i=1,nSolvers
3181
            Solver => Model % Solvers(i)
3182
            IF ( .NOT. ALLOCATED(RKCoeff(i) % k1)) CYCLE
3183

3184
            k1 => RKCoeff(i) % k1
3185
            k2 => RKCoeff(i) % k2
3186
            k3 => RKCoeff(i) % k3
3187
            k4 => RKCoeff(i) % k4
3188
            k4 = Solver % Variable % Values - Solver % Variable % PrevValues(:,1)
3189

3190
            Solver % Variable % Values = Solver % Variable % PrevValues(:,1) + &
3191
                ( k1 + 2*k2 + 2*k3 + k4 ) / 6
3192
          END DO
3193
        END IF
3194

3195
        DO i=1,nSolvers
3196
          Solver => Model % Solvers(i)
3197
          IF ( ALLOCATED(RKCoeff(i) % k1) ) THEN
3198
            DEALLOCATE(RKCoeff(i) % k1, RKCoeff(i) % k2, RKCoeff(i) % k3, RKCoeff(i) % k4)
3199
          END IF
3200

3201
          IF( ASSOCIATED( Solver % Variable ) ) THEN
3202
            IF(ASSOCIATED(Solver % Variable % Values)) THEN
3203
              n = SIZE(Solver % Variable % Values)
3204
              IF(n>0) &
3205
                Solver % Variable % Norm = ComputeNorm( Solver, n, Solver % Variable % Values)
3206
            END IF
3207
          END IF
3208
        END DO
3209
        DEALLOCATE(RKCoeff)
3210
      END IF
3211

3212
      IF ( .NOT.TransientSimulation ) SteadyStateReached = ALL(DoneThis)
3213
      DEALLOCATE( DoneThis, AfterConverged )
3214
    END IF
3215
!------------------------------------------------------------------------------      
3216

3217
    
3218
!------------------------------------------------------------------------------
3219
    IF( PRESENT( AfterTime ) ) THEN
3220
      ExecSlot = AfterTime
3221
    ELSE
3222
      ExecSlot = .TRUE.
3223
    END IF
3224

3225
    IF( ExecSlot ) THEN
3226
      CALL Info('SolveEquations','Solvers after timestep',Level=12)
3227
      DO k=1,nSolvers
3228
        Solver => Model % Solvers(k)
3229
        IF ( Solver % PROCEDURE==0 ) CYCLE
3230
        IF ( Solver % SolverExecWhen == SOLVER_EXEC_AFTER_TIME .OR. &
3231
            Solver % SolverExecWhen == SOLVER_EXEC_PREDCORR ) THEN
3232

3233
          IF( PRESENT( RealTimeStep ) ) THEN
3234
            IF( RealTimeStep == 1 ) THEN
3235
              IF( ListGetLogical( Solver % Values,'Skip First Timestep',Found ) ) CYCLE
3236
            END IF
3237
          END IF
3238

3239
          IF( Solver % SolverExecWhen == SOLVER_EXEC_PREDCORR ) THEN
3240
            CALL InitializeTimestep(Solver)
3241
            CALL ListAddLogical( Solver % Values,'Predictor Phase',.FALSE.)
3242
          END IF
3243

3244
          CALL SolverActivate( Model,Solver,dt,TransientSimulation )
3245
          CALL ParallelBarrier
3246
        END IF
3247
      END DO
3248
    END IF
3249

3250
    ParEnv_Common = ParEnv_Save
3251
    ParEnv => ParEnv_Common
3252
    IF(ParEnv % PEs>1) THEN
3253
      IF(.NOT.ASSOCIATED(ParEnv % Active)) ALLOCATE(ParEnv % Active(ParEnv % PEs))
3254
      ParEnv % Active = .TRUE.
3255
      ParEnv % ActiveComm = ELMER_COMM_WORLD
3256
    END IF
3257

3258
CONTAINS
3259

3260
    SUBROUTINE SolveCoupled()
3261

3262
    TYPE(Mesh_t), POINTER, SAVE :: PrevMesh
3263
    INTEGER, SAVE :: PrevMeshNoNodes
3264
    LOGICAL, SAVE :: FirstTime=.TRUE.
3265

3266
    IF(FirstTime) THEN
3267
      PrevMesh => Model % Mesh
3268
      PrevMeshNoNodes = Model % Mesh % NumberOfNodes
3269
      FirstTime = .FALSE.
3270
    END IF
3271

3272
!------------------------------------------------------------------------------
3273

3274
    CALL Info('SolveEquations','Performing set of solvers in sequence',Level=12)
3275
    
3276
     DO i=1,CoupledMaxIter
3277
       IF ( TransientSimulation .OR. Scanning ) THEN
3278
         IF( CoupledMaxIter > 1 ) THEN
3279
           CALL Info( 'SolveEquations', '-------------------------------------', Level=3 )
3280
           WRITE( Message, * ) 'Coupled system iteration: ', i
3281
           CALL Info( 'SolveEquations', Message, Level=3 )
3282
           CALL Info( 'SolveEquations', '-------------------------------------', Level=3 )         
3283
         END IF
3284
         steadyIt = i
3285
       END IF
3286
        
3287
       IF( GetNamespaceCheck() ) CALL ListPushNamespace('coupled '//i2s(i)//': ')
3288

3289
       DoneThis = .FALSE.
3290

3291
!      Initialize the mesh output flag to FALSE here, reactivated
3292
!      later for meshes connected to active solvers.
3293
!      ----------------------------------------------------------
3294
       Mesh => Model % Meshes
3295
       DO WHILE( ASSOCIATED( Mesh ) )
3296
         Mesh % OutputActive = .FALSE.
3297
         Mesh => Mesh % Next
3298
       END DO
3299

3300
!------------------------------------------------------------------------------
3301
!      Go through number of solvers (heat,laminar or turbulent flow, etc...)
3302
!------------------------------------------------------------------------------
3303
       DO k=1,Model % NumberOfSolvers
3304
!------------------------------------------------------------------------------
3305
          Solver => Model % Solvers(k)
3306

3307
          IF ( Solver % PROCEDURE == 0 ) THEN
3308
            IF( .NOT. ( Solver % SolverMode == SOLVER_MODE_COUPLED .OR. &
3309
              Solver % SolverMode == SOLVER_MODE_ASSEMBLY .OR. &
3310
              Solver % SolverMode == SOLVER_MODE_BLOCK ) ) THEN
3311

3312
              CALL Warn('SolveEquations','No routine related to solver!')
3313
              DoneThis(k) = .TRUE.
3314
              CYCLE
3315
            END IF
3316
          END IF
3317

3318
          IF( PRESENT( RealTimeStep ) ) THEN
3319
            IF( RealTimeStep == 1 ) THEN
3320
              IF( ListGetLogical( Solver % Values,'Skip First Timestep',Found ) ) THEN
3321
                DoneThis(k) = .TRUE.
3322
                CYCLE
3323
              END IF
3324
            END IF
3325
          END IF
3326
            
3327
          IF ( Solver % SolverExecWhen /= SOLVER_EXEC_ALWAYS ) THEN
3328
            DoneThis(k) = .TRUE.
3329
            CYCLE
3330
          END IF
3331

3332
          IF ( AfterConverged(k) .AND. .NOT. ALL(AfterConverged .OR. DoneThis) ) CYCLE
3333
!------------------------------------------------------------------------------
3334

3335
          n = 0
3336
          IF ( ASSOCIATED(Solver % Variable) ) THEN
3337
            IF ( ASSOCIATED(Solver % Variable % Values) ) &
3338
                n = SIZE(Solver % Variable % Values)
3339
            Solver % Variable % PrevNorm = Solver % Variable % Norm
3340
          END IF
3341

3342
          ! There are some operations that require that the previous steady state values
3343
          ! are present. Check for these operations.
3344
          !------------------------------------------------------------------------------
3345
          CalculateDerivative = ListGetLogical( Solver % Values, &
3346
                'Calculate Derivative', Stat )
3347
          IF( CalculateDerivative .AND. .NOT. Scanning ) THEN
3348
            CALL Fatal('SolveEquations','> Calculate Derivative < should be active only for scanning!')
3349
          END IF
3350

3351
          NeedSol = CalculateDerivative
3352
          IF(.NOT. NeedSol ) THEN
3353
            NeedSol = ( ListGetString( Solver % Values, &
3354
                'Steady State Convergence Measure', Stat ) /= 'norm')  
3355
            NeedSol = NeedSol .AND. Stat
3356
          END IF
3357
          IF(.NOT. NeedSol ) THEN
3358
            NeedSol = ListCheckPresent( Solver % Values, &
3359
                'Steady State Relaxation Factor')
3360
          END IF
3361

3362
          IF ( NeedSol .AND. n > 0 ) THEN
3363
            Stat = ASSOCIATED(Solver % Variable % SteadyValues)
3364
            IF(Stat) THEN
3365
              IF(SIZE(Solver % Variable % SteadyValues) /= n) THEN
3366
                DEALLOCATE(Solver % Variable % SteadyValues)
3367
                Stat = .FALSE.
3368
              END IF
3369
            END IF
3370
            IF(.NOT. Stat) THEN
3371
              ALLOCATE( Solver % Variable % SteadyValues(n), STAT=istat ) 
3372
              IF ( istat /= 0 ) CALL Fatal( 'SolveEquations', 'Memory allocation error.' )
3373
            END IF
3374
            Solver % Variable % SteadyValues(1:n) = Solver % Variable % Values(1:n)
3375
          END IF
3376

3377
          ! The solver may conditionally be turned into steady state 
3378
          ! if the > Steady State Condition < is set positive.
3379
          !------------------------------------------------------------------------
3380
          TransientSolver = TransientSimulation
3381
          IF( TransientSolver ) THEN
3382
            SSCond = ListGetCReal( Solver % Values,'Steady State Condition',Found )
3383
            IF( Found .AND. SSCond > 0.0_dp ) THEN
3384
              TransientSolver = .FALSE.
3385
              CALL Info('SolveEquations','Running solver in steady state',Level=6)
3386
            END IF
3387
          END IF
3388

3389
          CALL SolverActivate(Model,Solver,dt,TransientSolver)
3390
          
3391
          IF( TestDivergence ) THEN
3392
            DivergenceExit = ( Solver % Variable % NonlinConverged > 1 )
3393
            EXIT
3394
          END IF
3395
            
3396
          IF ( ASSOCIATED(Solver % Variable) ) THEN
3397
            IF ( ASSOCIATED(Solver % Variable % Values) ) &
3398
                n = SIZE(Solver % Variable % Values)
3399
          END IF
3400
!------------------------------------------------------------------------------
3401
!         check for coupled system convergence
3402
!------------------------------------------------------------------------------
3403

3404
           IF ( Scanning .OR. TransientSimulation ) THEN             
3405
             IF( CoupledMaxIter == 1 ) THEN
3406
               TestConvergence = .FALSE.
3407
               ! This means that the nonlinear system norm has not been computed
3408
               IF( Solver % Variable % NonlinConverged < 0 )  THEN
3409
                 TestConvergence = ListCheckPresent( Solver % Values,'Reference Norm' )
3410
               END IF
3411
               IF(.NOT. TestConvergence) THEN
3412
                 TestConvergence = ListCheckPresent( Solver % Values,'Steady State Relaxation Factor')
3413
               END IF
3414
             ELSE
3415
               TestConvergence = ( i >= CoupledMinIter )
3416
             END IF
3417
           ELSE    ! Steady-state
3418
             TestConvergence = .TRUE.
3419
           END IF
3420
           
3421
           IF( TestConvergence .OR. CalculateDerivative ) THEN
3422
             IF ( ParEnv % PEs > 1 ) THEN
3423
               IF ( ParEnv % Active(ParEnv % MyPE+1) ) THEN
3424
                 CALL ComputeChange(Solver,.TRUE., n)
3425
               ELSE
3426
                 IF(.NOT.ASSOCIATED(Solver % Variable)) ALLOCATE(Solver % Variable)
3427
                 Solver % Variable % SteadyConverged = 1
3428
               END IF
3429
             ELSE
3430
               CALL ComputeChange(Solver,.TRUE.)
3431
             END IF
3432
           END IF
3433

3434
           ! The ComputeChange subroutine sets a flag to zero if not yet
3435
           ! converged (otherwise -1/1)
3436
           !------------------------------------------------------------
3437
           IF( TestConvergence ) THEN
3438
             DoneThis(k) = ( Solver % Variable % SteadyConverged /= 0 ) 
3439
           END IF
3440

3441
           IF( Solver % Mesh % AdaptiveFinished .AND. .NOT. DoneThis(k)) THEN
3442
             CALL Info('SolveEquations','Overriding convergence due to Adaptive Meshing Finished!')
3443
             DoneThis(k) = .TRUE.
3444
           END IF                    
3445
           
3446
           CALL ParallelAllReduceAnd( DoneThis(k) )
3447
           IF( ALL(DoneThis) ) EXIT
3448
!------------------------------------------------------------------------------
3449
         END DO
3450
!------------------------------------------------------------------------------
3451
         CALL ListPopNamespace()
3452

3453
         !Check if the mesh changed - should do this elsewhere too?
3454
         IF(ASSOCIATED(Model % Mesh, PrevMesh) .AND. &
3455
              Model % Mesh % NumberOfNodes == PrevMeshNoNodes) THEN
3456
           Model % Mesh % Changed = .FALSE.
3457
         ELSE
3458
           PrevMesh => Model % Mesh
3459
           PrevMeshNoNodes = Model % Mesh % NumberOfNodes
3460
           Model % Mesh % Changed = .TRUE.
3461
         END IF
3462

3463
         IF( DivergenceExit ) EXIT
3464

3465
         IF ( ALL(DoneThis) ) EXIT
3466
      END DO
3467
      
3468
      IF( TestConvergence .AND. CoupledMaxIter > 1 ) THEN
3469
        IF ( TransientSimulation .AND. .NOT. ALL(DoneThis) ) THEN
3470
          CALL Info( 'SolveEquations','Coupled system iteration: '//I2S(MIN(i,CoupledMaxIter)),Level=4)
3471
          CoupledAbort = ListGetLogical( Model % Simulation,  &
3472
              'Coupled System Abort Not Converged', Found )
3473
          CALL NumericalError('SolveEquations','Coupled system did not converge',CoupledAbort)
3474
        END IF
3475
      END IF
3476
        
3477
    END SUBROUTINE SolveCoupled
3478
!------------------------------------------------------------------------------
3479
  END SUBROUTINE SolveEquations
3480
!------------------------------------------------------------------------------
3481
  
3482

3483
!------------------------------------------------------------------------------
3484
!> This is a line of monolithic solvers where different physics and constraints 
3485
!> are assemblied to the same matrix.
3486
!> Provide assembly loop and solution of linear and nonlinear systems
3487
!> This routine uses minimalistic assembly routines to create the 
3488
!> matrices. Often the results to less labour in coding which may be 
3489
!> comprimized by less flexibility.
3490
! 
3491
! There are currently two modes
3492
! IsCoupledSolver:  the variable consists of several pieces and the matrix equation
3493
!                   is created within this solver.
3494
! IsAssemblySolver: only one assembly routine is used and the matrix equation may
3495
!                   be formed externally using standard procedures.
3496
! TODO:
3497
! non-nodal elements (internal allocation & right elementtype)
3498
! solution (and perm) vectors could be reused by setting pointers to correct positions
3499
! check the time-dependency
3500
! check for eigenmode analysis
3501
! different lumped functionalities (energy,...)
3502
! bandwidth optimization for coupled system
3503
!------------------------------------------------------------------------------
3504
  SUBROUTINE CoupledSolver( Model, Solver, dt, Transient )
3505
!------------------------------------------------------------------------------    
3506
    IMPLICIT NONE
3507
!------------------------------------------------------------------------------
3508
    TYPE(Solver_t), TARGET :: Solver
3509
    TYPE(Model_t) :: Model
3510
    REAL(KIND=dp) :: dt
3511
    LOGICAL :: Transient
3512
!------------------------------------------------------------------------------
3513
! Local variables
3514
!------------------------------------------------------------------------------
3515
    TYPE(Solver_t), POINTER :: PSolver, PSolver2
3516
    TYPE(Variable_t), POINTER :: Var
3517
    TYPE(Element_t),POINTER :: Element
3518
    INTEGER :: i,j,k,l,t,n,nd,NoIterations,iter,istat
3519
    REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), DAMP(:,:), MASS(:,:), FORCE(:)
3520
    LOGICAL :: GotIt, GotIt2, GotProc, BulkMode, ThisConstraint, &
3521
        IsCoupledSolver, IsAssemblySolver, IsProcedure, IsListMatrix
3522
    INTEGER :: Row, Col, ColDofs, RowDofs, ColInd0, RowInd0, MaxDofs, &
3523
         ColVar, RowVar, Nrow, Ncol, NoVar, NoCons, TotSize, ConDofs, OffSet(20), &
3524
         VarSizes(20),VarDofs(20)
3525
    INTEGER :: ElementsFirst, ElementsLast
3526
    INTEGER, POINTER :: VarPerm(:), ColPerm(:), RowPerm(:), ColInds(:), RowInds(:), DirPerm(:)
3527
    REAL(KIND=dp), POINTER :: ConsValues(:)
3528
    REAL(KIND=dp) :: NonlinearTol, Norm, PrevNorm, ConsValue, ConsCoeff, ConsVolume
3529
    CHARACTER(LEN=max_name_len) :: str, VarName, ColName, RowName, ConsType
3530
    LOGICAL :: Coupling, Equality, AssemblySymmetric, AssemblyAntisymmetric, &
3531
        AllDirFlag, Robust, ReduceStep
3532
    INTEGER, POINTER :: Rows(:),Cols(:),Indexes(:),AllPerm(:)
3533
    TYPE(ListMatrix_t), POINTER :: Alist(:) => NULL()
3534
    REAL(KIND=dp), POINTER :: AllValues(:)
3535
    REAL(KIND=dp), POINTER CONTIG :: ForceVector(:)
3536
    LOGICAL, POINTER :: AllDir(:)
3537
    TYPE (Matrix_t), POINTER :: Amat
3538
    REAL(KIND=dp), POINTER :: Component(:)
3539
    INTEGER, POINTER :: VarInds(:)
3540
    TYPE(Mesh_t), POINTER :: Mesh
3541
    TYPE(ValueList_t), POINTER :: SolverParams
3542

3543
    SolverParams => ListGetSolverParams(Solver)
3544

3545
    IsCoupledSolver = .FALSE.
3546
    IsAssemblySolver = .FALSE.
3547

3548
    IsProcedure = ListCheckPresent( SolverParams,'Procedure')
3549

3550
    CALL Info('CoupledSolver','---------------------------------------',Level=5)
3551
    IF( IsProcedure ) THEN
3552
      CALL Info('CoupledSolver','Inheriting matrix from procedure',Level=5)     
3553
    ELSE IF( Solver % SolverMode == SOLVER_MODE_COUPLED ) THEN
3554
      CALL Info('CoupledSolver','Solving a system of equations',Level=5)
3555
      IsCoupledSolver = .TRUE.
3556
    ELSE IF( Solver % SolverMode == SOLVER_MODE_ASSEMBLY ) THEN
3557
      CALL Info('CoupledSolver','Solving one equation',Level=5)      
3558
      IsAssemblySolver = .TRUE.
3559
    ELSE
3560
      CALL Fatal('CoupledSolver','You should maybe not be here?')
3561
    END IF
3562
    CALL Info('CoupledSolver','---------------------------------------',Level=5)
3563

3564

3565
    Mesh => GetMesh()
3566
    PSolver => Solver
3567

3568
    !------------------------------------------------------------------------------
3569
    ! Check out which variables the coupled model includes
3570
    ! and compure size information related to the new coupled dof.
3571
    !------------------------------------------------------------------------------
3572
    Offset = 0
3573
    VarSizes = 0    
3574
    NoVar = 0
3575
    NoCons = 0
3576
    VarDofs = 0
3577

3578
    AllDirFlag = .FALSE.
3579
    IF( IsProcedure .OR. IsAssemblySolver ) THEN
3580
      CALL Info('CoupledSolver','Using existing variable and matrix',Level=8)
3581
      IF(.NOT. ASSOCIATED(Solver % Matrix)) THEN
3582
        CALL Fatal('CoupledSolver','In legacy solver mode the Matrix should exist!')
3583
      END IF
3584
      IF(.NOT. ASSOCIATED(Solver % Variable)) THEN
3585
        CALL Fatal('CoupledSolver','In legacy solver mode the Variable should exist!')
3586
      END IF
3587
      NoVar = 1
3588
      Var => Solver % Variable
3589
      VarDofs(1) = Var % Dofs
3590
      VarSizes(1) = SIZE( Var % Values )
3591
      TotSize = VarSizes(1)
3592
      MaxDofs = VarDofs(1)
3593
    ELSE
3594

3595
      DO i = 1,100
3596
        WRITE (str,'(A,I0)') 'Variable ',i
3597
        VarName = ListGetString( SolverParams, TRIM(str), GotIt )
3598
        IF(.NOT. GotIt) EXIT
3599
        Var => VariableGet( Mesh % Variables, TRIM(VarName) )
3600
        
3601
        IF(.NOT. ASSOCIATED( Var )) THEN
3602
          CALL Info('CoupledSolver','Variable '//TRIM(VarName)//' does not exist, creating',Level=10)
3603
          PSolver => Solver
3604
          Var => CreateBlockVariable(PSolver, i, VarName )
3605
        END IF
3606

3607
        NoVar = NoVar + 1
3608
        VarDofs(NoVar) = Var % Dofs
3609
        VarSizes(NoVar) = SIZE( Var % Values )
3610
        Offset(NoVar+1) = Offset(NoVar) + VarSizes(NoVar)
3611
      END DO
3612

3613
      !------------------------------------------------------------------------------------
3614
      ! Here is a hack for taking constraints into account where dofs are created on-the-fly
3615
      ! might be better to create special solvers somewhere else.
3616
      ! The size of constraint is deduced directly from its type.
3617
      !-------------------------------------------------------------------------------------
3618
      DO i = 1,100
3619
        WRITE (str,'(A,I0)') 'Constraint ',i
3620
        VarName = ListGetString( SolverParams, TRIM(str), GotIt )
3621
        IF(.NOT. GotIt) EXIT
3622
        
3623
        NoCons = NoCons + 1      
3624
        
3625
        WRITE (str,'(A,I0,A)') 'Constraint ',i,' Type'
3626
        ConsType = ListGetString( SolverParams,TRIM(str))
3627
        
3628
        SELECT CASE( ConsType )
3629
          
3630
        CASE('integral')
3631
          ConDofs = 1
3632
          
3633
        CASE('floating')
3634
          ConDofs = 1
3635
          AllDirFlag = .TRUE.
3636
          
3637
        CASE('equality')
3638
          ConDofs = 0
3639
          AllDirFlag = .TRUE.
3640
          
3641
        CASE DEFAULT          
3642
          CALL Warn('CoupledSolver','Coupled constraint does not really work yet')
3643
          WRITE (str,'(A,I0,A)') 'Constraint ',i,' DOFs'
3644
          ConDofs = ListGetInteger( SolverParams, str)
3645

3646
        END SELECT
3647

3648
        ! Create the constrained variables for possible other use
3649
        !-----------------------------------------------------------
3650
        IF( ConDofs > 0 ) THEN
3651
          Var => VariableGet( Mesh % Variables, VarName )
3652
          IF( .NOT. ASSOCIATED(Var) ) THEN
3653
            CALL Info('CoupledSolver','Constraint '//TRIM(VarName)//' does not exist, creating',Level=10)
3654
            ALLOCATE(ConsValues(ConDofs))
3655
            CALL VariableAdd( Mesh % Variables, Mesh, Solver, &
3656
                VarName, ConDofs, ConsValues, Output = .FALSE. )          
3657
            Var => VariableGet( Mesh % Variables, VarName )         
3658
          END IF
3659
        END IF
3660
        
3661
        j = NoVar + NoCons
3662
        VarDofs(j) = ConDofs
3663
        VarSizes(j) = 1
3664
        Offset(j+1) = Offset(j) + VarSizes(j)
3665
      END DO
3666

3667
      DO j=1,NoVar+NoCons
3668
        WRITE(Message,'(A,I0,A,T35,I0)') 'Permutation offset ',j,': ',OffSet(j)
3669
        CALL Info('CoupledSolver',Message,Level=8)
3670
      END DO
3671
      
3672
      TotSize = SUM( VarSizes )
3673
      MaxDofs = MAXVAL( VarDofs )
3674
      
3675
      WRITE(Message,'(A,T35,I0)') 'Number of coupled variables: ',NoVar
3676
      CALL Info('CoupledSolver',Message,Level=8)
3677
      
3678
      WRITE(Message,'(A,T35,I0)') 'Number of constraints: ',NoCons
3679
      CALL Info('CoupledSolver',Message,Level=8)
3680
      
3681
      WRITE(Message,'(A,T35,I0)') 'Size of coupled system: ',TotSize
3682
      CALL Info('CoupledSolver',Message,Level=8)
3683
      
3684
      ! For the 1st time the matrix should be a list, later CRS
3685
      !--------------------------------------------------------
3686
      IF(.NOT. ASSOCIATED(Solver % Matrix)) THEN
3687
        Amat => AllocateMatrix()
3688
        Amat % ListMatrix => Alist
3689
        Amat % FORMAT = MATRIX_LIST      
3690
        Solver % Matrix => Amat
3691
      END IF
3692

3693
      ! This actual variable is not saved, so a dummy name suffices
3694
      !------------------------------------------------------------
3695
      VarName = ListGetString( SolverParams,'Variable', GotIt )
3696
      IF(.NOT. GotIt) THEN
3697
        CALL Info('CoupledSolver','New coupled variable added: Alldofs', Level=6)
3698
        VarName = 'alldofs'
3699
      END IF
3700

3701
      ! If the variable hasn't been created do it now
3702
      !----------------------------------------------
3703
      Var => VariableGet(Mesh % Variables, VarName )
3704
      IF(.NOT. ASSOCIATED(Var) ) THEN
3705

3706
        ALLOCATE(AllPerm(TotSize),AllValues(TotSize),ForceVector(TotSize))
3707
        DO i=1,TotSize
3708
          AllPerm(i) = i
3709
          AllValues(i) = 0.0_dp
3710
        END DO
3711
        CALL VariableAdd( Mesh % Variables, Mesh, Solver, &
3712
            VarName,1,AllValues,AllPerm,Output=.FALSE.)        
3713
        Amat % Rhs => ForceVector
3714
        Solver % Variable => VariableGet(Mesh % Variables, VarName )
3715
        Var => Solver % Variable
3716
       
3717
        ! Map the original vectors to the monolithic vector if requested
3718
        !----------------------------------------------------------------------------------
3719
        IF( ListGetLogical( SolverParams,'Coupled Initial Guess',GotIt)) THEN
3720
          CALL SingleToCoupledVector()
3721
        END IF
3722
      END IF
3723

3724
    END IF  ! IsProcedure
3725
      
3726

3727
!------------------------------------------------------------------------------
3728
! Do some initial stuff common for both modes
3729
!------------------------------------------------------------------------------
3730
  
3731
    N = Mesh % MaxElementDOFs
3732
    
3733
    ALLOCATE( FORCE( MaxDofs*N ),      &
3734
        STIFF( MaxDofs*N, MaxDofs*N ), &
3735
        DAMP( MaxDofs*N, MaxDofs*N ),  &
3736
        MASS(  MaxDofs*N, MaxDofs*N ), &
3737
        ColInds( N ), RowInds( N ),    &
3738
        Indexes( n ), STAT=istat )
3739
    IF ( istat /= 0 ) CALL FATAL('CoupledSolver','Memory allocation error')
3740
    
3741
    NoIterations = GetInteger( SolverParams,'Nonlinear System Max Iterations',GotIt)
3742
    IF(.NOT. GotIt) NoIterations = 1
3743
    NonlinearTol = GetCReal( SolverParams,'Nonlinear System Convergence Tolerance',gotIt)
3744
    Robust = ListGetLogical(SolverParams,'Nonlinear System Linesearch',GotIt)
3745
    Amat => GetMatrix()
3746
    ForceVector => Amat % Rhs
3747
    IF( AllDirFlag ) ALLOCATE( AllDir(TotSize) ) 
3748
   
3749
!------------------------------------------------------------------------------
3750
! Iterate over any nonlinearity of material or source
3751
!------------------------------------------------------------------------------
3752
    CALL Info('CoupledSolver','-------------------------------------------------',Level=5)
3753

3754
    IsListMatrix = ( Amat % FORMAT == MATRIX_LIST )
3755
    
3756
    DO iter = 1,NoIterations
3757

3758
      WRITE(Message,'(A,T35,I5)') 'Coupled iteration:',iter
3759
      CALL Info('CoupledSolver',Message,Level=5)
3760
   
3761
100   IF( IsProcedure ) THEN
3762

3763
        ! Perform the normal solver procedure with just one iteration 
3764
        ! linear system solver being disabled i.e. just do the assembly.
3765
        ! Also the possible internal linesearch is disabled.
3766
        !---------------------------------------------------------------------
3767
        CALL ListAddInteger( SolverParams,'Nonlinear System Max Iterations',1)
3768
        CALL ListAddLogical( SolverParams,'Linear System Solver Disabled',.TRUE.)
3769
        CALL ListAddLogical( SolverParams,'Skip Compute Nonlinear Change',.TRUE.)
3770
        IF( Robust ) THEN
3771
          CALL ListRemove( SolverParams,'Nonlinear System Linesearch')
3772
        END IF
3773

3774
        CALL SingleSolver( Model, PSolver, dt, Transient )
3775

3776
        CALL ListAddInteger( SolverParams,'Nonlinear System Max Iterations',NoIterations)
3777
        CALL ListRemove( SolverParams,'Linear System Solver Disabled')
3778
        IF(Robust ) THEN
3779
          CALL ListAddLogical(SolverParams,'Nonlinear System Linesearch',.TRUE.)
3780
        ELSE
3781
          CALL ListAddLogical(SolverParams,'Skip Compute Nonlinear Change',.FALSE.)
3782
        END IF          
3783
     
3784
      ELSE
3785
        IF(.NOT. IsListMatrix ) CALL DefaultInitialize()
3786
        IF( AllDirFlag ) AllDir = .FALSE.
3787

3788
!----------------------------------------------------------------
3789
! Here we use the same assembly for coupled system and block system
3790
! approach.
3791
!----------------------------------------------------------------
3792
        DO RowVar = 1,NoVar 
3793
          DO ColVar = 1,NoVar          
3794
            CALL BlockSystemAssembly(PSolver,dt,Transient,RowVar,ColVar,&
3795
                    Offset(RowVar),Offset(ColVar))
3796
          END DO
3797
        END DO
3798
    
3799
!-------------------------------------------------------------------------
3800
! Tailored assembly routines for setting the constraints to the same matrix
3801
!-------------------------------------------------------------------------
3802

3803
        IF( IsCoupledSolver .AND. NoCons > 0 ) THEN
3804
          CALL CoupledConstraintAssembly()
3805
        END IF
3806

3807
!----------------------------------------------------------------------
3808
! The CRS matrix may be created only when the matrix structure is known
3809
! and some initialization may be done only when the matrix exists
3810
!----------------------------------------------------------------------
3811
        IF( IsListMatrix ) THEN
3812
          CALL List_ToCRSMatrix(Amat)
3813
          IsListMatrix = .FALSE.
3814
          CALL AddEquationSolution(PSolver, Transient )
3815
          GOTO 100  
3816
        END IF
3817
        CALL DefaultFinishAssembly()
3818
      
3819
!------------------------------------------------------------------------------
3820
!    Do the Dirichlet conditions using offset
3821
!------------------------------------------------------------------------------          
3822
        IF( IsCoupledSolver ) THEN
3823
          CALL CoupledSystemDirichlet()
3824
        ELSE      
3825
          CALL DefaultDirichletBCs()
3826
        END IF
3827
      END IF
3828

3829
!------------------------------------------------------------------------------
3830
!    Check the stepsize of nonlinear iteration using the Armijo-GoldStein 
3831
!    criterion for the stepsize.
3832
!------------------------------------------------------------------------------          
3833
      IF( Robust  ) THEN   
3834
        IF( CheckStepSize( Solver, iter == 1) ) THEN
3835
          IF( IsCoupledSolver ) CALL CoupledToSingleVector()
3836
          GOTO 100
3837
        ELSE
3838
          IF(Solver % Variable % NonlinConverged > 0) EXIT
3839
        END IF
3840
      END IF
3841

3842
!------------------------------------------------------------------------------
3843
!   Finally solve the system
3844
!------------------------------------------------------------------------------          
3845
      Norm = DefaultSolve()     
3846
      CALL Info('CoupledSolver','-------------------------------------------------',Level=5)
3847

3848
      IF( .NOT. ( IsAssemblySolver .OR. IsProcedure ) ) THEN
3849
        CALL CoupledToSingleVector()
3850
      END IF
3851

3852
      ! The non-robust solver will use one assembly routine less 
3853
      IF(.NOT. Robust) THEN
3854
        IF(Solver % Variable % NonlinChange < NonlinearTol) EXIT
3855
      END IF
3856
    END DO
3857

3858
    
3859
    DEALLOCATE( FORCE, STIFF, DAMP, MASS, ColInds, RowInds, Indexes )
3860
    IF( AllDirFlag ) DEALLOCATE( AllDir ) 
3861
    
3862
    CALL Info('CoupledSolver','All done',Level=8)
3863
    CALL Info('CoupledSolver','-------------------------------------------------',Level=5)
3864

3865

3866
  CONTAINS 
3867

3868
!------------------------------------------------------------------------------
3869
!> Integration routine for integral type of constraints i.e. body or boundary
3870
!> integral of some dof is known a priori.
3871
!------------------------------------------------------------------------------
3872
   SUBROUTINE IntegralConstraint( Mass, Damp, Stiff, Force, Element, n )
3873
!------------------------------------------------------------------------------
3874
     
3875
      IMPLICIT NONE
3876
!------------------------------------------------------------------------------
3877
      REAL(KIND=dp) :: Stiff(:,:), Damp(:,:), Mass(:,:), Force(:)
3878
      TYPE(Element_t), POINTER :: Element
3879
      INTEGER :: n
3880
!------------------------------------------------------------------------------
3881
!    Local variables
3882
!------------------------------------------------------------------------------
3883
      TYPE(GaussIntegrationPoints_t) :: IntegStuff
3884
      REAL(KIND=dp) :: Basis(n)
3885
      REAL(KIND=dp) :: Weight, DetJ
3886
      INTEGER :: i,j,t,p,q
3887
      LOGICAL :: Found
3888
      
3889
      TYPE(Nodes_t),SAVE :: Nodes
3890
      
3891
      CALL GetElementNodes( Nodes ) 
3892
      IntegStuff = GaussPoints( Element )
3893
      
3894
      DO t=1,IntegStuff % n
3895
        Found = ElementInfo( Element, Nodes, IntegStuff % u(t), &
3896
            IntegStuff % v(t), IntegStuff % w(t), detJ, Basis )     
3897
        Weight = IntegStuff % s(t) * detJ
3898
        DO p=1,n
3899
          STIFF(1,p) = STIFF(1,p) + Weight * Basis(p)
3900
          FORCE(p) = FORCE(p) + ConsValue * Weight * Basis(p)
3901
        END DO
3902
        ConsVolume = ConsVolume + Weight 
3903
      END DO
3904

3905
!------------------------------------------------------------------------------
3906
    END SUBROUTINE IntegralConstraint
3907
!------------------------------------------------------------------------------
3908
 
3909
   
3910
    !------------------------------------------------------------------------------          
3911
    !> Subroutine sets the Dirichlet conditions for the coupled system 
3912
    !> using the single system vector names and sizes.
3913
    !----------------------------------------------------------------------------------
3914
    SUBROUTINE CoupledSystemDirichlet()
3915
      CALL Info( 'CoupledSolver', 'Setting coupled system Dirichlet conditions', Level=6 )
3916
      
3917
      DO i = 1,NoVar
3918
        WRITE (str,'(A,I0)') 'Variable ',i
3919
        VarName = ListGetString( SolverParams, TRIM(str), GotIt )
3920
        IF(.NOT. GotIt) EXIT
3921
        
3922
        Var => VariableGet( Mesh % Variables, TRIM(VarName) )
3923
        IF (.NOT. ASSOCIATED(Var)) EXIT 
3924
        
3925
        CALL DefaultDirichletBCs( Ux=Var, UOffset=Offset(i) )
3926
      END DO
3927
 
3928
    END SUBROUTINE CoupledSystemDirichlet
3929

3930

3931
    !------------------------------------------------------------------------------          
3932
    ! Subroutine copies results from the single system vectors to a coupled system
3933
    ! initial guess.
3934
    !----------------------------------------------------------------------------------
3935
    SUBROUTINE SingleToCoupledVector()
3936
      CALL Info('CoupledSolver','Copying an initial guess',Level=8)
3937
      
3938
      DO i = 1,NoVar + NoCons
3939
        IF( VarDofs(i) == 0 ) CYCLE
3940
        
3941
        IF( i <= NoVar ) THEN
3942
          WRITE (str,'(A,I0)') 'Variable ',i
3943
        ELSE
3944
          WRITE (str,'(A,I0)') 'Constraint ',i-NoVar        
3945
        END IF
3946
        VarName = ListGetString( SolverParams, TRIM(str), GotIt )
3947
        Var => VariableGet( Mesh % Variables, TRIM(VarName) )
3948
        
3949
        ! CALL Info('CoupledSolver','Copying from variable: '//TRIM(VarName),Level=5)
3950
        DO j=1,SIZE(Var % Values)
3951
          Solver % Variable % Values(Offset(i)+j) = Var % Values(j)
3952
        END DO
3953
      END DO
3954

3955
    END SUBROUTINE SingleToCoupledVector
3956

3957

3958
    !------------------------------------------------------------------------------          
3959
    !> Subroutine copies results from the coupled system vector back to the 
3960
    !> original vectors.
3961
    !----------------------------------------------------------------------------------
3962
    SUBROUTINE CoupledToSingleVector()
3963
      CALL Info('CoupledSolver','Copying results into original variables',Level=8)
3964
      DO i = 1,NoVar + NoCons
3965
        IF( VarDofs(i) == 0 ) CYCLE
3966
        
3967
        IF( i <= NoVar ) THEN
3968
          WRITE (str,'(A,I0)') 'Variable ',i
3969
        ELSE
3970
          WRITE (str,'(A,I0)') 'Constraint ',i-NoVar        
3971
        END IF
3972
        VarName = ListGetString( SolverParams, TRIM(str), GotIt )
3973
        Var => VariableGet( Mesh % Variables, TRIM(VarName) )
3974
        
3975
        ! CALL Info('CoupledSolver','Copying to variable: '//TRIM(VarName),Level=5)
3976
        DO j=1,SIZE(Var % Values)
3977
          Var % Values(j) = Solver % Variable % Values(Offset(i)+j)
3978
        END DO
3979

3980
        CALL InvalidateVariable( Model % Meshes, Solver % Mesh, VarName )
3981
      END DO
3982

3983
    END SUBROUTINE CoupledToSingleVector
3984

3985

3986
    !-----------------------------------------------------------------------------------
3987
    !> Perform constraints assembly
3988
    !> Some constraints are assembled by integration while others are done by elimination.
3989
    !> It is important that constraints are applied after the normal assembly process. 
3990
    !-----------------------------------------------------------------------------------
3991
    SUBROUTINE CoupledConstraintAssembly()
3992
      !-----------------------------------------------------------------------------------
3993
      INTEGER :: TargetDof, TmpInds(4) 
3994

3995
      CALL Info('CoupledSolver','Starting constraint assembly',Level=8)
3996
      
3997
      BulkMode = .TRUE.
3998
200   IF(BulkMode) THEN
3999
        ! CALL Info('CoupledSolver','Starting constraint bulk assembly',Level=5)
4000
        ElementsFirst = 1
4001
        ElementsLast = Mesh % NumberOfBulkElements 
4002
      ELSE
4003
        ! CALL Info('CoupledSolver','Starting boundary assembly',Level=5)
4004
        ElementsFirst = Mesh % NumberOfBulkElements + 1
4005
        ElementsLast =  Mesh % NumberOfBulkElements + &
4006
              Mesh % NumberOfBoundaryElements
4007
      END IF
4008
      
4009
      ! Variables over rows
4010
      !-------------------------------------------
4011
      DO RowVar = 1,NoCons        
4012

4013
        RowDofs = VarDofs(NoVar + RowVar) 
4014
        RowInd0 = Offset(NoVar + RowVar)
4015

4016
        AssemblySymmetric = .FALSE.
4017
        AssemblyAntiSymmetric = .FALSE.
4018

4019
        WRITE (str,'(A,I0)') 'Constraint ',RowVar
4020
        RowName = ListGetString( Solver % Values, TRIM(str), GotIt )
4021

4022
        WRITE (str,'(A,I0,A)') 'Constraint ',RowVar,' Variables'
4023
        VarInds => ListGetIntegerArray( Solver % Values, TRIM(str) )
4024

4025
        IF( ASSOCIATED(VarInds)) THEN
4026
          ColVar =  VarInds(1) 
4027
        ELSE        
4028
          CALL Fatal('CoupledSolver','Cannot continue without pointer to variables')
4029
        END IF
4030
        
4031
        WRITE (str,'(A,I0,A)') 'Constraint ',RowVar,' Components'
4032
        TargetDof = ListGetInteger( Solver % Values, TRIM(str), GotIt )
4033
      
4034
        WRITE (str,'(A,I0,A)') 'Constraint ',RowVar,' Value'
4035
        ConsValue = GetCReal(Solver % Values,TRIM(str), GotIt)
4036

4037
        WRITE (str,'(A,I0,A)') 'Constraint ',RowVar,' Coeff'
4038
        ConsCoeff = GetCReal(Solver % Values,TRIM(str), GotIt)
4039
        IF(.NOT. GotIt) ConsCoeff = 1.0_dp
4040

4041

4042
        IF ( ConsType == 'equality' ) THEN
4043
          WRITE (str,'(A,I0)') 'Variable ',VarInds(2)
4044
          VarName = ListGetString( SolverParams, TRIM(str) )
4045
          Var => VariableGet( Mesh % Variables, TRIM(VarName) ) 
4046
          RowPerm => Var % Perm
4047
          RowDofs = VarDofs(VarInds(2))
4048
          RowInd0 = Offset(VarInds(2)) 
4049
        ELSE
4050
          Nrow = 1
4051
          RowInds(1:Nrow) = 1         
4052
          Var => VariableGet( Mesh % Variables, TRIM(RowName) )      
4053
          
4054
          ! Add the diagonal entry expected by some subroutines
4055
          !-------------------------------------------------------          
4056
          DO i=1,Nrow
4057
            DO j=1,RowDofs
4058
              Row = RowInd0 + RowDofs * (RowInds(i)-1) + j
4059
              CALL AddToMatrixElement( Amat, Row, Row, 0.0_dp )
4060
            END DO
4061
          END DO
4062
        END IF
4063
         
4064
        IF( ConsType == 'integral') THEN
4065
          AssemblySymmetric = .TRUE.
4066
          ConsVolume = 0.0_dp
4067
        END IF
4068

4069
        WRITE (str,'(A,I0)') 'Variable ',ColVar
4070
        ColName = ListGetString( Solver % Values, TRIM(str), GotIt )
4071
        Var => VariableGet( Mesh % Variables, TRIM(ColName) )
4072
          
4073
        ! Constrain only the target variable
4074
        !-----------------------------------------------
4075
        
4076
        ColPerm => Var % Perm
4077
        ColDofs = Var % Dofs
4078
        ColInd0 = Offset(ColVar)
4079
          
4080
        ! The assembly loop for a submatrix starts here
4081
        !------------------------------------------------------------------------------             
4082
        DO t=ElementsFirst,ElementsLast
4083

4084
          Element => Mesh % Elements(t)
4085
          Model % CurrentElement => Element
4086
          
4087
          ! How to treat non-nodal elements must be rethought 
4088
          ! nd = GetElementNOFDOFs( Element, Solver )                  
4089
          !-----------------------------------------------------------------
4090
          n  = GetElementNOFNodes()
4091
          nd = GetElementDOFs(Indexes)
4092
                        
4093
          ! Set the permutations for equality constraint
4094
          !----------------------------------------------------
4095
          IF( ConsType == 'equality') THEN
4096
            Nrow = nd
4097
            RowInds(1:Nrow) = RowPerm(Indexes(1:nrow))
4098
            IF(.NOT. ALL(RowInds(1:Nrow) > 0)) CYCLE
4099
          END IF
4100
            
4101
          Ncol = nd
4102
          ColInds(1:n) = ColPerm(Indexes(1:n))
4103
          IF(.NOT. ALL(ColInds(1:Ncol) > 0)) CYCLE                 
4104

4105
            
4106
          ! Check where constraints are active, both bodies and BCs
4107
          !-------------------------------------------------------------------                         
4108
          Coupling = .FALSE.              
4109
          IF( BulkMode ) THEN
4110
            ! Check coupling to bodies using Body/Equation section
4111
            Coupling = GetLogical( GetBodyForce(), RowName, gotIt)
4112
          ELSE 
4113
            Coupling = GetLogical( GetBC(), RowName, gotIt)
4114
          END IF
4115
          IF( .NOT. Coupling ) CYCLE
4116
              
4117
          ! These two constraints are based on moving already assembled information
4118
          ! rather than assembling new information. If the row is already treated cycle.
4119
          ! The constraints have been tested only with one-component cases.
4120
          !-----------------------------------------------------------------------------
4121
          IF( ConsType == 'floating' .OR. ConsType == 'equality') THEN
4122

4123
            DO i=1,Ncol                
4124
              DO k=0,ColDofs-1
4125
                
4126
                ! Note that for this type the column is rather also a row
4127
                !--------------------------------------------------------
4128
                Col  = ColInd0 + ColDofs * ColInds(i) - k                  
4129
                IF( AllDir(Col) ) CYCLE
4130
                AllDir(Col) = .TRUE.
4131
                
4132
                IF( ConsType == 'floating') THEN
4133
                  Row = RowInd0 + 1
4134
                ELSE IF( ConsType == 'equality') THEN
4135
                  Row = RowInd0 + RowDofs * RowInds(i) - k
4136
                END IF
4137

4138
                IF( IsListMatrix ) THEN
4139
                  CALL MoveRow( Amat, Col, Row )
4140
                  CALL SetMatrixElement( Amat,Col,Col,0.0_dp )
4141
                  CALL SetMatrixElement( Amat,Col,Row,0.0_dp )
4142
                  CYCLE
4143
                END IF
4144

4145
                CALL MoveRow( Amat, Col, Row, ConsCoeff ) 
4146
                ForceVector(Row) = ForceVector(Row) + ForceVector(Col)
4147
                ForceVector(Col) = 0.0_dp
4148
                
4149
                ForceVector(Row) = ForceVector(Row) + ConsValue
4150
                CALL SetMatrixElement( Amat,Col,Col,1.0_dp )
4151
                CALL SetMatrixElement( Amat,Col,Row,-ConsCoeff)
4152
              END DO
4153
            END DO
4154
            CYCLE
4155
          END IF
4156

4157

4158
          ! Do the assembly, now really (active only for some constraints)
4159
          !----------------------------------------------------------------
4160
            
4161
          STIFF = 0.0_dp
4162
          DAMP = 0.0_dp
4163
          MASS = 0.0_dp
4164
          FORCE = 0.0_dp
4165
          
4166
          CALL IntegralConstraint( MASS, DAMP, STIFF, FORCE, Element, Ncol )
4167
          
4168
          IF ( Transient ) THEN
4169
            IF( Solver % TimeOrder == 1 ) THEN
4170
              CALL Default1stOrderTime( MASS,STIFF,FORCE)
4171
            ELSE IF( Solver % TimeOrder == 2) THEN
4172
              CALL Default2ndOrderTime( MASS,DAMP,STIFF,FORCE )
4173
            END IF
4174
          END IF
4175

4176
           
4177
          ! Assemble the matrix with offset
4178
          ! Because we want to have constraints component-wise
4179
          ! There is somewhat dirty hack for calling the gluematrix
4180
          !---------------------------------------------------------
4181
          DO k=1,ColDofs
4182
            IF( TargetDof /= 0 .AND. TargetDof /= k) CYCLE
4183
            TmpInds(1:Ncol) = ColDofs * (ColInds(1:Ncol)-1) + k
4184

4185
            IF( IsListMatrix ) THEN            
4186
              CALL GlueLocalSubMatrix( Amat, &
4187
                  RowInd0,ColInd0,Nrow,Ncol,RowInds,TmpInds,&
4188
                  RowDofs,1,STIFF )
4189
              IF( AssemblySymmetric .OR. AssemblyAntisymmetric ) THEN
4190
                CALL GlueLocalSubMatrix( Amat, &
4191
                    ColInd0,RowInd0,Ncol,Nrow,TmpInds,RowInds,&
4192
                    1,RowDofs,STIFF )               
4193
              END IF
4194
              CYCLE
4195
            END IF
4196

4197
            CALL GlueLocalSubMatrix( Amat, &
4198
                RowInd0,ColInd0,Nrow,Ncol,RowInds,TmpInds,&
4199
                RowDofs,1,STIFF )
4200

4201
            ! For some constraints assemble also the transpose
4202
            !--------------------------------------------------
4203
            IF( AssemblySymmetric ) THEN
4204
              CALL GlueLocalSubMatrix( Amat, &
4205
                  ColInd0,RowInd0,Ncol,Nrow,TmpInds,RowInds,&
4206
                  1,RowDofs,TRANSPOSE(STIFF) )               
4207
            ELSE IF( AssemblyAntisymmetric ) THEN
4208
              CALL GlueLocalSubMatrix( Amat, &
4209
                  ColInd0,RowInd0,Ncol,Nrow,TmpInds,RowInds,&
4210
                  1,RowDofs,-TRANSPOSE(STIFF) )               
4211
            END IF
4212
            
4213
            ! Assemble the r.h.s with offset
4214
            !-----------------------------------------------
4215
            DO i=1,Nrow
4216
              DO j=1,RowDofs
4217
                Row = RowInd0 + RowDofs * (RowInds(i)-1) + j
4218
                ForceVector(Row) = ForceVector(Row) + &
4219
                    FORCE(RowDofs*(i-1)+j)
4220
              END DO
4221
            END DO
4222
          
4223
          END DO          
4224
        END DO
4225
      
4226
        ! For constraints do some special setting
4227
        !---------------------------------------------------
4228
        IF( ConsType == 'integral' ) THEN
4229
          ! PRINT *,'Integral constraint sum of weights: ',ConsVolume
4230
          DO i=1,Nrow
4231
            DO j=1,RowDofs
4232
              Row = RowInd0 + RowDofs * (RowInds(i)-1) + j
4233
              ForceVector(Row) = ConsValue
4234
            END DO
4235
          END DO
4236
        END IF
4237
      END DO
4238
      
4239
      IF(BulkMode) THEN
4240
        ! CALL Info( 'CoupledSolver', 'Bulk assembly done for constraints', Level=4 )
4241
        BulkMode = .FALSE.
4242
        GOTO 200
4243
      ELSE 
4244
        ! CALL Info( 'CoupledSolver', 'Boundary assembly done for constraints', Level=4 )
4245
      END IF
4246

4247
    END SUBROUTINE CoupledConstraintAssembly
4248

4249
!------------------------------------------------------------------------------
4250
  END SUBROUTINE CoupledSolver
4251
!------------------------------------------------------------------------------
4252
 
4253

4254

4255

4256
!------------------------------------------------------------------------------
4257
!> This is a line of solvers where a matrix of matrices and a vector of vectors 
4258
!> are created to allow different kinds of block strategies for the solvers.
4259
!> This strategy has optimal memory consumption even if block strategies are 
4260
!> employed on the linear system level. 
4261
!------------------------------------------------------------------------------
4262
  SUBROUTINE BlockSolver( Model, Solver, dt, Transient )
4263
!------------------------------------------------------------------------------
4264
    IMPLICIT NONE
4265
 !------------------------------------------------------------------------------
4266
    TYPE(Solver_t), TARGET :: Solver
4267
    TYPE(Model_t) :: Model
4268
    REAL(KIND=dp) :: dt
4269
    LOGICAL :: Transient
4270
!------------------------------------------------------------------------------
4271
! Local variables
4272
!------------------------------------------------------------------------------
4273

4274
    TYPE(Solver_t), POINTER :: PSolver
4275
    TYPE(Variable_t), POINTER :: Var, SolverVar
4276
    INTEGER :: i,j,k,l,n,nd,NonLinIter,tests,NoTests,iter
4277
    LOGICAL :: GotIt, GotIt2, BlockPrec, BlockGS
4278
    INTEGER :: ColVar, RowVar, NoVar, BlockDofs
4279

4280
    REAL(KIND=dp) :: NonlinearTol, Norm, PrevNorm, Residual, PrevResidual, &
4281
        TotNorm, MaxChange, alpha, beta, omega, rho, oldrho, s, r, PrevTotNorm, &
4282
        Coeff
4283
    CHARACTER(LEN=max_name_len) :: str, VarName, ColName, RowName
4284
    LOGICAL :: Robust, LinearSearch, AcceptStep, IsProcedure, ScaleSystem,&
4285
        ReuseMatrix, LS, InitDone
4286
    INTEGER, POINTER :: VarPerm(:)
4287
    
4288
    TYPE (Matrix_t), POINTER :: Amat, SolverMatrix
4289
    TYPE(Mesh_t), POINTER :: Mesh
4290
    TYPE(ValueList_t), POINTER :: SolverParams
4291

4292
    CALL Info('BlockSolver','---------------------------------------',Level=5)
4293
    IF( Solver % SolverMode /= SOLVER_MODE_BLOCK ) THEN
4294
      CALL Fatal('BlockSolver','You should maybe not be here?')
4295
    ELSE
4296
      CALL Info('BlockSolver','Solving system of equations utilizing block strategies',Level=5)
4297
    END IF
4298
    CALL Info('BlockSolver','---------------------------------------',Level=5)
4299

4300
    SolverParams => ListGetSolverParams(Solver)
4301
    Mesh => Solver % Mesh
4302
    PSolver => Solver
4303

4304
    SolverRef => Solver
4305

4306
    isParallel = ParEnv % PEs > 1
4307
    
4308
         
4309
    ! Determine some parameters related to the block strategy
4310
    !------------------------------------------------------------------------------
4311
    BlockPrec = GetLogical( SolverParams,'Block Preconditioner',GotIt)
4312
    IF(.NOT. GotIt ) BlockPrec = .TRUE.
4313
    BlockGS = GetLogical( SolverParams,'Block Gauss-Seidel',GotIt)    
4314
    IsProcedure = ListCheckPresent( SolverParams,'Procedure')
4315

4316
    ! Initialize the matrix
4317
    ! 1) Matrix is inherited from a monolithic matrix equation, or
4318
    ! 2) Matrix is assemblied using the compact assembly process
4319
    !------------------------------------------------------------------------------
4320
    IF( IsProcedure ) THEN 
4321
      BlockDofs = Solver % Variable % Dofs      
4322

4323
      CALL BlockInitMatrix( Solver, TotMatrix, BlockDofs )
4324

4325
      NoVar = TotMatrix % NoVar
4326

4327
      TotMatrix % Solver => Solver
4328
      SolverMatrix => Solver % Matrix
4329
      SolverVar => Solver % Variable
4330

4331
    ELSE
4332
      TotMatrix => Solver % BlockMatrix
4333
      IF( .NOT.ASSOCIATED(TotMatrix)) THEN
4334
        DO i = 1,9
4335
          WRITE (str,'(A,I0)') 'Variable ',i
4336
          IF(.NOT. ListCheckPresent( SolverParams, TRIM(str)) ) EXIT
4337
        END DO
4338
        NoVar = i-1
4339
        
4340
        CALL BlockInitMatrix( Solver, TotMatrix, NoVar )
4341
        TotMatrix % Solver => Solver
4342
        
4343
        WRITE(Message,'(A,I0)') 'Number of coupled variables: ',NoVar
4344
        CALL Info('BlockSolver',Message)
4345
        IF( NoVar == 0 ) CALL Fatal('BlockSolver','No variables, nothing to do!')
4346

4347
        ! Create the matrix structures using the list structure as 
4348
        !------------------------------------------------------------------------------    
4349
        CALL Info('BlockSolver','Creating matrix structured using list matrices',Level=12)
4350
        DO RowVar=1,NoVar
4351
          Solver % Variable => TotMatrix % SubVector(RowVar) % Var
4352
          
4353
          DO ColVar=1,NoVar            
4354
            Solver % Variable => TotMatrix % SubVector(ColVar) % Var        
4355
            Solver % Matrix => TotMatrix % Submatrix(RowVar,ColVar) % Mat
4356
            
4357
            Amat => Solver % Matrix        
4358
            Amat % ListMatrix => NULL()
4359
            Amat % FORMAT = MATRIX_LIST      
4360
            Amat % NumberOfRows = 0
4361
            
4362
            CALL BlockSystemAssembly(PSolver,dt,Transient,RowVar,ColVar)
4363
            
4364
            IF( .NOT. ASSOCIATED( Amat % ListMatrix ) ) THEN
4365
              Amat % FORMAT = MATRIX_CRS
4366
            ELSE
4367
              CALL List_ToCRSMatrix(Amat)
4368
              CALL AddEquationSolution(PSolver, Transient )            
4369
            END IF
4370

4371
          END DO
4372
        END DO
4373
      END IF
4374

4375
      ! The variable and solver pointer should direct somewhere as they are used to 
4376
      ! monitor the steady-state solution, for example. Therefore the user may
4377
      ! choose the component if not preferring one. 
4378
      !---------------------------------------------------------------------------
4379
      RowVar = ListGetInteger( SolverParams,'Primary Variable',GotIt)
4380
      IF(.NOT. GotIt) RowVar = 1
4381
      SolverVar => TotMatrix % SubVector(RowVar) % Var
4382
      SolverMatrix => TotMatrix % Submatrix(RowVar,RowVar) % Mat
4383
    END IF
4384

4385
    !------------------------------------------------------------------------------
4386
    ! This is the nonlinear loop that may be used either for true 
4387
    ! nonlinearities.
4388
    !------------------------------------------------------------------------------    
4389
    NonLinIter = GetInteger( SolverParams,'Nonlinear System Max Iterations',GotIt)
4390
    IF(.NOT. GotIt) NonLinIter = 1
4391
    NonlinearTol = GetCReal( SolverParams,'Nonlinear System Convergence Tolerance',gotIt)
4392

4393
    Robust = ListGetLogical(SolverParams,'Nonlinear System Linesearch',GotIt)
4394
    IF( Robust ) THEN
4395
      LinearSearch = ListGetLogical( SolverParams,'Nonlinear System Linesearch Linear')
4396
      NoTests = GetInteger( SolverParams,'Nonlinear System Linesearch Iterations',GotIt)
4397
      IF(.NOT. GotIt) NoTests = NonLinIter
4398
      CALL ListAddString(SolverParams,'Nonlinear System Convergence Measure','residual')
4399
      CALL ListAddLogical(SolverParams,'Skip Compute Nonlinear Change',.TRUE.)
4400
    END IF       
4401

4402
    CALL Info('BlockSolver','-------------------------------------------------',Level=6)
4403
    Residual = -1.0_dp
4404
    PrevResidual = -1.0_dp
4405
    
4406
    DO iter = 1,NonLinIter
4407
      
4408
      WRITE(Message,'(A,T35,I0)') 'Coupled iteration: ',iter
4409
      CALL Info('BlockSolver',Message,Level=6)
4410

4411
      tests = 0
4412

4413
      ! Assembly of matrices either using legacy or compact strategy
4414
      !------------------------------------------------------------------
4415
100   IF( IsProcedure ) THEN
4416

4417
        ! Perform the normal solver procedure with just one iteration 
4418
        ! linear system solver being disabled i.e. just do the assembly.
4419
        !---------------------------------------------------------------------
4420
        CALL ListAddInteger( SolverParams,'Nonlinear System Max Iterations',1)
4421
        CALL ListAddLogical( SolverParams,'Linear System Solver Disabled',.TRUE.)
4422

4423
        CALL SingleSolver( Model, PSolver, dt, Transient )
4424

4425
        ! Scaling of linear system. This could also be fused with the submatrix
4426
        ! picking procedure.
4427
        !---------------------------------------------------------------------
4428
        ScaleSystem = ListGetLogical( SolverParams,'Linear System Scaling',GotIt) 
4429
        IF(.NOT. GotIt) ScaleSystem = .TRUE.
4430
        
4431
        IF( ScaleSystem ) THEN
4432
	  CALL Info('BlockSolver','Applying scaling',Level=8)
4433
          CALL ScaleLinearSystem(Solver, SolverMatrix, &
4434
              SolverMatrix % Rhs, Solver % Variable % Values )
4435
          CALL ListAddLogical( Solver % Values,'Linear System Scaling',.FALSE.)
4436
        END IF
4437

4438
        CALL BlockPickMatrix( Solver, NoVar )
4439

4440
        CALL ListAddInteger( SolverParams,'Nonlinear System Max Iterations',NonLinIter)
4441
        CALL ListRemove( SolverParams,'Linear System Solver Disabled')
4442
      ELSE
4443
                
4444
        DO RowVar=1,NoVar          
4445
          Solver % Variable => TotMatrix % SubVector(RowVar) % Var          
4446
          DO ColVar=1,NoVar            
4447
            
4448
            Solver % Matrix => TotMatrix % Submatrix(RowVar,ColVar) % Mat
4449
            IF( Solver % Matrix % NumberOfRows == 0 ) CYCLE
4450
            Solver % Variable => TotMatrix % SubVector(ColVar) % Var
4451
            CALL InitializeToZero(Solver % Matrix, Solver % Matrix % rhs)
4452
            
4453
            CALL ListPushNameSpace('block:')
4454
            CALL ListPushNameSpace('block '//i2s(RowVar)//i2s(ColVar)//':')
4455
            CALL BlockSystemAssembly(PSolver,dt,Transient,RowVar,ColVar)
4456
            
4457
            ! Mainly sets the r.h.s. in transient case correctly
4458
            CALL DefaultFinishAssembly()                    
4459
            
4460
            CALL BlockSystemDirichlet(TotMatrix,RowVar,ColVar)
4461
            CALL ListPopNameSpace(); CALL ListPopNameSpace()
4462
          END DO
4463
        END DO
4464
      END IF
4465

4466
      ! The user may give a user defined preconditioner matrix
4467
      !-----------------------------------------------------------
4468
      CALL BlockPrecMatrix( Solver, NoVar ) 
4469
      
4470
      IF (isParallel) THEN
4471
        DO RowVar=1,NoVar
4472
          DO ColVar=1,NoVar
4473
            Amat => TotMatrix % SubMatrix(RowVar,ColVar) % Mat
4474
            Amat % Comm = ELMER_COMM_WORLD
4475
            Parenv % ActiveComm = Amat % Comm
4476
            Solver % Variable => TotMatrix % SubVector(ColVar) % Var
4477
            CALL ParallelInitMatrix(Solver,Amat)
4478

4479
            Amat % ParMatrix % ParEnv % ActiveComm = Amat % Comm
4480
            ParEnv => Amat % ParMatrix % ParEnv
4481
            CALL ParallelActive( .TRUE.)
4482
          END DO
4483
        END DO
4484
     END IF
4485

4486

4487
!------------------------------------------------------------------------------
4488
!    Check the stepsize of nonlinear iteration using the Armijo-GoldStein 
4489
!    criterion for the stepsize, if requested
4490
!------------------------------------------------------------------------------          
4491
      IF( Robust  ) THEN   
4492

4493
200     AcceptStep = CheckStepSizeBlock(TotMatrix,tests==0,PrevResidual,Residual)
4494
        tests = tests + 1
4495
        
4496
        IF( tests >  NoTests ) THEN
4497
          CALL Fatal('BlockSolver','Maximum number of linesearch steps exceeded')
4498
        END IF
4499
        
4500
        ! Update the reference residual only when new step is accepted
4501
        IF( iter == 1 ) THEN
4502
          PrevResidual = Residual
4503
        ELSE
4504
          IF( AcceptStep ) THEN 
4505
            PrevResidual = Residual
4506
            IF(Solver % Variable % NonlinChange < NonlinearTol) EXIT
4507
          ELSE
4508
            IF( LinearSearch ) THEN
4509
              GOTO 100
4510
            ELSE
4511
              GOTO 200
4512
            END IF
4513
          END IF
4514
        END IF
4515
      END IF
4516
      
4517

4518
      !------------------------------------------------------------------------------
4519
      ! Finally solve the system using 'outer: ' as the optional namespace
4520
      ! for the linear system setting.
4521
      !------------------------------------------------------------------------------          
4522
      
4523
      TotNorm = 0.0_dp
4524
      MaxChange = 0.0_dp
4525

4526
      CALL ListPushNameSpace('outer:')
4527

4528
      ! The case with one block is mainly for testing and developing features
4529
      ! related to nonlinearity and assembly.
4530
      !----------------------------------------------------------------------
4531
      IF( NoVar == 1 ) THEN
4532
	CALL Info('BlockSolver','Solving in standard manner',Level=8)
4533

4534
        Solver % Variable => TotMatrix % SubVector(1) % Var
4535
        Solver % Matrix => TotMatrix % Submatrix(1,1) % Mat
4536

4537
        TotNorm = DefaultSolve()
4538
        MaxChange = Solver % Variable % NonlinChange 
4539

4540
      ELSE IF( BlockPrec ) THEN
4541
	CALL Info('BlockSolver','Using block preconditioning strategy',Level=8)        
4542
        CALL BlockKrylovIter( Solver, MaxChange )
4543

4544
      ELSE
4545
        Solver % Variable => TotMatrix % SubVector(1) % Var
4546
        Solver % Matrix => TotMatrix % Submatrix(1,1) % Mat
4547

4548
	CALL Info('BlockSolver','Using block solution strategy',Level=8)
4549
        CALL BlockStandardIter( Solver, MaxChange )
4550
      END IF      
4551
      CALL ListPopNameSpace()
4552

4553
      ! For legacy matrices do the backmapping 
4554
      !------------------------------------------
4555
      IF( IsProcedure ) THEN
4556
        Solver % Matrix => SolverMatrix
4557
        Solver % variable => SolverVar
4558
        IF( ScaleSystem ) THEN
4559
          CALL BackScaleLinearSystem(Solver, SolverMatrix, &
4560
              SolverMatrix % Rhs, Solver % Variable % Values )
4561
          CALL ListAddLogical( Solver % Values,'Linear System Scaling',.TRUE.)
4562
        END IF
4563
      END IF
4564
      
4565
      IF(.NOT. Robust ) THEN
4566
        IF( MaxChange < NonlinearTol) EXIT
4567
      END IF
4568
    END DO
4569
       
4570
    ! Before returning set the pointers appropriately as saved before
4571
    !---------------------------------------------------------------------------
4572
    Solver % Variable => SolverVar
4573
    Solver % Matrix => SolverMatrix
4574

4575
    CALL Info('BlockSolver','All done',Level=8)
4576
    CALL Info('BlockSolver','-------------------------------------------------',Level=6)
4577

4578

4579
  CONTAINS 
4580

4581

4582
    !------------------------------------------------------------------------------          
4583
    !> Subroutine sets the Dirichlet conditions for the block system 
4584
    !> using the single system vector names and sizes. For that aim there is an 
4585
    !> additional flag that is used to detect that the matrix cannot have the digonal 
4586
    !> entry that is by default set to unity and the r.h.s. to the target value.
4587
    !> Now for off-diagonal matrices they will be both omitted. 
4588
    !> The routine assumes that the r.h.s. of the off diagonal matrices is zero.
4589
    !----------------------------------------------------------------------------------
4590
    SUBROUTINE BlockSystemDirichlet(BlockMatrix,NoRow,NoCol)
4591
      
4592
      TYPE(BlockMatrix_t) :: BlockMatrix
4593
      INTEGER :: NoRow, NoCol
4594

4595
      LOGICAL :: OffDiagonal
4596

4597
      CALL Info( 'BlockSolver', 'Setting block system Dirichlet conditions', Level=8 )
4598
      
4599
      Solver % Matrix => BlockMatrix % SubMatrix( NoRow, NoCol ) % Mat
4600
      IF( Solver % Matrix % NumberOfRows == 0 ) RETURN
4601
       
4602
      WRITE (str,'(A,I0)') 'Variable ',NoRow
4603
      VarName = ListGetString( SolverParams, TRIM(str), GotIt )
4604
      IF(.NOT. GotIt) RETURN
4605

4606
      Var => VariableGet( Mesh % Variables, TRIM(VarName) )
4607
      IF (.NOT. ASSOCIATED(Var)) RETURN
4608
 
4609
      Solver % Variable => Var
4610
      OffDiagonal = ( NoRow /= NoCol )
4611
      CALL DefaultDirichletBCs( Ux=Var, OffDiagonalMatrix = OffDiagonal )
4612
 
4613
    END SUBROUTINE BlockSystemDirichlet
4614

4615

4616
    !------------------------------------------------------------------------------          
4617
    !> Compute the block norm, optionally. 
4618
    !----------------------------------------------------------------------------------
4619
    FUNCTION ComputeBlockNorm(BlockMatrix, DiagonalOnly, MatrixOnly ) RESULT ( Residual )
4620
      
4621
      TYPE(BlockMatrix_t), TARGET :: BlockMatrix
4622
      LOGICAL, OPTIONAL :: DiagonalOnly, MatrixOnly
4623
      REAL(KIND=dp) :: Residual
4624

4625
      REAL(KIND=dp) :: rnorm, xnorm, bnorm, TotXnorm, TotRnorm, TotBnorm, TotRelNorm
4626
      REAL(KIND=dp), POINTER :: x(:),res(:),b(:),rtmp(:)
4627
      INTEGER :: n, NoRow,NoCol, NoVar
4628
      TYPE(Matrix_t), POINTER :: A
4629

4630

4631
      CALL Info('BlockSolver','Computing block matrix norm',Level=8)
4632
      
4633
      NoVar = BlockMatrix % NoVar
4634
      ALLOCATE( rtmp( BlockMatrix % MaxSize ), res( BlockMatrix % MaxSize) )
4635

4636
      DO NoRow = 1,NoVar 
4637
        Var => BlockMatrix % SubVector(NoRow) % Var
4638
        n = SIZE( Var % Values )
4639

4640
        x => Var % Values
4641
        xNorm = ComputeNorm(Solver, n, x)
4642

4643
        A => BlockMatrix % SubMatrix(NoRow,NoRow) % Mat
4644
        Solver % Matrix => A
4645
        b => A % rhs
4646
        
4647
        xNorm = ComputeNorm(Solver, n, x)
4648
        bNorm = ComputeNorm(Solver, n, b)
4649

4650
	res = 0.0_dp        
4651
        rtmp = 0.0_dp  
4652
        
4653
        DO NoCol = 1,NoVar           
4654
          IF( PRESENT( DiagonalOnly ) ) THEN
4655
            IF( DiagonalOnly .AND. NoCol /= NoRow ) CYCLE
4656
          END IF
4657
          
4658
          Var => BlockMatrix % SubVector(NoCol) % Var
4659
          x => Var % Values
4660

4661
          A => BlockMatrix % SubMatrix( NoRow, NoCol )  % Mat
4662
          IF( A % NumberOfRows == 0 ) CYCLE
4663
          b => A % rhs
4664
          
4665
          CALL MatrixVectorMultiply( A, x, rtmp)      
4666
          res = res + rtmp
4667
        END DO
4668
        
4669
        IF( PRESENT( MatrixOnly ) ) THEN
4670
          IF( .NOT. MatrixOnly ) THEN              
4671
            res = res - b
4672
          END IF
4673
        ELSE
4674
          res = res - b
4675
        END IF
4676
	
4677
        rNorm = ComputeNorm(Solver, n, res)
4678

4679
        ! PRINT *,'comp norms:',NoRow,Rnorm,Xnorm
4680

4681
        BlockMatrix % SubVector(NoRow) % rnorm = rnorm
4682
        BlockMatrix % SubVector(NoRow) % xnorm = xnorm
4683
        BlockMatrix % SubVector(NoRow) % bnorm = bnorm
4684
      END DO
4685

4686
      TotRnorm = SUM( BlockMatrix % SubVector(1:NoVar) % rnorm ** 2) 
4687
      TotXnorm = SUM( BlockMatrix % SubVector(1:NoVar) % xnorm ** 2) 
4688
      TotBnorm = SUM( BlockMatrix % SubVector(1:NoVar) % bnorm ** 2) 
4689

4690
      TotRnorm = SQRT( TotRnorm / NoVar )
4691
      TotXnorm = SQRT( TotXnorm / NoVar )
4692
      TotBnorm = SQRT( TotBnorm / NoVar )
4693

4694
      BlockMatrix % Rnorm = TotRNorm
4695
      BlockMatrix % Xnorm = TotXnorm
4696
      BlockMatrix % Bnorm = TotBnorm
4697

4698
      ! PRINT *,'tot norms:',TotRnorm,TotXnorm
4699
      DEALLOCATE( rtmp, res )
4700

4701
      Residual = BlockMatrix % Rnorm
4702

4703
    END FUNCTION ComputeBlockNorm
4704
    
4705

4706

4707

4708
!------------------------------------------------------------------------------
4709
    FUNCTION CheckStepSizeBlock(BlockMatrix,FirstTrial,PrevResidual,Residual) &
4710
        RESULT (Success) 
4711
!------------------------------------------------------------------------------
4712
      TYPE(BlockMatrix_t) :: BlockMatrix
4713
      REAL(KIND=dp) :: PrevResidual, Residual
4714
      LOGICAL :: FirstTrial,Success
4715
!------------------------------------------------------------------------------
4716
      INTEGER :: i,n,niter
4717
      TYPE(Matrix_t), POINTER :: A
4718
      REAL(KIND=dp), POINTER :: b(:), x(:), x0(:), r(:)
4719
      REAL(KIND=dp) :: Norm, PrevNorm, rNorm, bNorm, Relaxation, Alpha, Myy
4720
      TYPE(Variable_t), POINTER :: iterV
4721
      LOGICAL :: Stat
4722
      CHARACTER(LEN=MAX_NAME_LEN) :: SolverName
4723
      
4724
      
4725
      SAVE Alpha, Relaxation, Myy
4726
      
4727
      !--------------------------------------------------------------------------
4728
      ! This is the real residual r=b-Ax
4729
      !--------------------------------------------------------------------------
4730
      Residual = ComputeBlockNorm( BlockMatrix ) 
4731
      
4732
      
4733
      ! Negative (impossible) value may be used as indicator that's its the 1st step
4734
      !-----------------------------------------------------------------------------
4735
      IF( PrevResidual < 0.0 ) THEN
4736
        Success = .FALSE.
4737
        RETURN
4738
      END IF
4739
      
4740
      ! At the first step set the relaxation
4741
      !-----------------------------------------------------------------------------
4742
      IF( FirstTrial ) THEN
4743
        Alpha = 1.0_dp
4744
        Relaxation = ListGetConstReal( Solver % Values, &
4745
            'Nonlinear System Linesearch Factor', Stat )
4746
        IF(.NOT. Stat) Relaxation = 0.5_dp
4747
        Myy = ListGetConstReal( Solver % Values, &
4748
            'Nonlinear System Linesearch Limit', Stat )
4749
        IF(.NOT. Stat) Myy = 0.5_dp
4750
      END IF
4751
      
4752
      ! Armijo GoldStein Criterion for accepting stepsize
4753
      !-----------------------------------------------------------------
4754
      Success = ( PrevResidual - Residual > Myy * Alpha * PrevResidual)
4755
      
4756
      IF( Success ) THEN      
4757
        iterV => VariableGet( Solver % Mesh % Variables, 'nonlin iter' )
4758
        niter = NINT(iterV % Values(1))
4759
        
4760
        DO i=1,BlockMatrix % NoVar
4761
          Var => BlockMatrix % SubVector(i) % Var
4762
          b => BlockMatrix % SubMatrix(i,i) % Mat % rhs
4763
          x => Var % Values
4764
          n = SIZE( b )
4765
          
4766
          Solver % Variable => Var
4767
          bNorm = ComputeNorm(Solver, n, b)
4768
          
4769
          Var % NonlinChange = Residual / bNorm
4770
          
4771
          Norm = ComputeNorm(Solver, n, x)
4772
          Solver % Variable % Norm = Norm
4773
          
4774
          SolverName = ListGetString( Solver % Values, 'Equation',Stat)
4775
          IF(.NOT. Stat) SolverName = Solver % Variable % Name
4776
          
4777
          WRITE( Message, '(a,g15.8,g15.8,a,I2)') &
4778
              'NS (ITER='//i2s(niter)//') (NRM,RELC): (',Norm, Residual / bNorm,&
4779
              ' ) :: '// TRIM(SolverName),i
4780
          CALL Info( 'CheckStepSize', Message, Level=3 )       
4781
        END DO
4782
        iterV % Values(1) = niter + 1 
4783
      ELSE
4784
        
4785
        DO i=1,BlockMatrix % NoVar
4786
          Var => BlockMatrix % SubVector(i) % Var
4787
          IF(.NOT. ASSOCIATED(Var % NonlinValues)) &
4788
              CALL Fatal('CheckStepSize','Previous nonlinear solution is needed')       
4789
          x0 => Var % NonlinValues
4790

4791
          x => Var % Values
4792

4793
          ! PRINT *,'Before Range:',i,MINVAL(x),MAXVAL(x),MINVAL(x0),MAXVAL(x0)
4794

4795
          x = (1-Relaxation) * x0 + Relaxation * x
4796

4797
          ! PRINT *,'After Range:',i,MINVAL(x),MAXVAL(x),MINVAL(x0),MAXVAL(x0)
4798

4799
        END DO
4800
        
4801
        Alpha = Alpha * Relaxation
4802
        CALL Info( 'CheckStepSize','Step rejected, increasing relaxation', Level=5 )      
4803
        ! PRINT *,'Residual',Residual,PrevResidual,Alpha
4804

4805
      END IF
4806

4807
      RowVar = ListGetInteger( SolverParams,'Primary Variable',GotIt)
4808
      IF(.NOT. GotIt) RowVar = 1
4809
      Solver % Matrix => TotMatrix % Submatrix(RowVar,RowVar) % Mat
4810
      Solver % Variable => TotMatrix % SubVector(RowVar) % Var
4811

4812
      
4813
 !------------------------------------------------------------------------------
4814
    END FUNCTION CheckStepSizeBlock
4815
 !------------------------------------------------------------------------------
4816
       
4817
!------------------------------------------------------------------------------
4818
  END SUBROUTINE BlockSolver
4819
!------------------------------------------------------------------------------
4820
 
4821
!---------------------------------------------------
4822
!> Perform assembly for the block system linear
4823
!> system of equations.
4824
!---------------------------------------------------
4825
  SUBROUTINE BlockSystemAssembly(Solver,dt,Transient,RowVar,ColVar,&
4826
      RowIndOffset,ColIndOffset)
4827
!---------------------------------------------------
4828
    TYPE(Solver_t), POINTER :: Solver
4829
    REAL(KIND=dp) :: dt
4830
    LOGICAL :: Transient
4831
    INTEGER :: RowVar, ColVar
4832
    INTEGER, OPTIONAL :: RowIndOffset,ColIndOffset
4833

4834
    INTEGER :: RowInd0, ColInd0
4835
    REAL(KIND=dp) :: SymmCoeff
4836
    TYPE(Element_t), POINTER :: Element
4837
    TYPE(Variable_t), POINTER :: Var
4838
    TYPE(Matrix_t), POINTER :: Amat
4839
    TYPE(Mesh_t), POINTER :: Mesh
4840
    TYPE(Valuelist_t), POINTER :: SolverParams
4841
    INTEGER :: i,j,t,n,nd,istat
4842
    INTEGER :: ElementsFirst, ElementsLast, Row, Col, RowDofs, ColDofs, Ncol, Nrow
4843
    INTEGER :: AllocCols, AllocRows
4844
    INTEGER, POINTER :: ColPerm(:), RowPerm(:), Indexes(:),ColInds(:),RowInds(:) 
4845
    LOGICAL :: GotIt
4846
    REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), DAMP(:,:), MASS(:,:), FORCE(:)
4847
    REAL(KIND=dp), POINTER :: ForceVector(:)
4848
    CHARACTER(LEN=MAX_NAME_LEN) :: ProcName, RowName, ColName, str
4849
    INTEGER(KIND=AddrInt) :: ProcPntr    
4850
    LOGICAL :: BulkMode, AssemblySymmetric, AssemblyAntiSymmetric, IsListMatrix
4851
    LOGICAL :: AllocationsDone = .FALSE., Diagonal
4852
    CHARACTER(*), PARAMETER :: Caller="BlockSystemAssembly"   
4853

4854
    SAVE :: AllocationsDone, AllocCols, AllocRows, &
4855
        FORCE, STIFF, DAMP, MASS, ColInds, RowInds, indexes
4856

4857
    Mesh => Solver % Mesh
4858
    SolverParams => Solver % Values
4859
    Amat => Solver % Matrix
4860
    ForceVector => Amat % rhs
4861

4862
    RowInd0 = 0
4863
    ColInd0 = 0
4864
    IF( PRESENT( RowIndOffset ) ) RowInd0 = RowIndOffset
4865
    IF( PRESENT( ColIndOffset ) ) ColInd0 = ColIndOffset
4866

4867
    ! Row Variable
4868
    !-------------------------------------------
4869
    WRITE (str,'(A,I0)') 'Variable ',RowVar
4870
    RowName = ListGetString( SolverParams, TRIM(str), GotIt )
4871
    Var => VariableGet( Mesh % Variables, TRIM(RowName) ) 
4872
    IF(.NOT. ASSOCIATED( Var ) .AND. RowVar == 1 .AND. ColVar == 1 ) THEN
4873
      Var => Solver % Variable
4874
    END IF
4875
    IF( .NOT. ASSOCIATED( Var ) ) THEN
4876
      CALL Fatal(Caller,'Could not find variable: '//I2S(RowVar))
4877
    END IF
4878
    RowDofs = Var % Dofs
4879
    RowPerm => Var % Perm
4880
    IF( .NOT. ASSOCIATED( RowPerm ) ) THEN
4881
      CALL Fatal(Caller,'Could not find permutation: '//I2S(RowVar))
4882
    END IF
4883
    
4884
    ! Column variable
4885
    !------------------------------------------
4886
    IF( ColVar /= RowVar ) THEN
4887
      WRITE (str,'(A,I0)') 'Variable ',ColVar
4888
      ColName = ListGetString( SolverParams, TRIM(str), GotIt )
4889
      Var => VariableGet( Mesh % Variables, TRIM(ColName) )
4890
    END IF          
4891
    IF( .NOT. ASSOCIATED( Var ) ) THEN
4892
      CALL Fatal(Caller,'Could not find variable: '//I2S(ColVar))
4893
    END IF
4894
    ColDofs = Var % Dofs
4895
    ColPerm => Var % Perm
4896
    IF( .NOT. ASSOCIATED( ColPerm ) ) THEN
4897
      CALL Fatal(Caller,'Could not find permutation: '//I2S(ColVar))
4898
    END IF
4899

4900
    ! These could be user provided for each block
4901
    !-----------------------------------------
4902
    AssemblySymmetric = ListGetLogical( SolverParams,'Symmetric Assembly',GotIt)
4903
    AssemblyAntiSymmetric = ListGetLogical( SolverParams,'AntiSymmetric Assembly',GotIt)
4904

4905
    IsListMatrix = ( Solver % Matrix % FORMAT == MATRIX_LIST ) 
4906
    N = Mesh % MaxElementDOFs    
4907

4908
    IF( AllocationsDone ) THEN
4909
      IF( RowDofs * n > AllocRows .OR. ColDofs * n > AllocCols )  THEN
4910
        DEALLOCATE( FORCE, STIFF, DAMP, MASS, ColInds, RowInds )        
4911
        AllocationsDone = .FALSE.
4912
      END IF
4913
    END IF
4914
    
4915
    IF( .NOT. AllocationsDone ) THEN
4916
      AllocRows = RowDofs * n
4917
      AllocCols = ColDofs * n
4918
      ALLOCATE(Indexes(AllocRows))
4919

4920
      ALLOCATE( FORCE( AllocRows ),      &
4921
          STIFF( AllocRows, AllocCols ), &
4922
          DAMP( AllocRows, AllocCols ),  &
4923
          MASS( AllocRows, AllocCols ), &
4924
          ColInds( N ), RowInds( N ),    &
4925
          STAT=istat )
4926
      IF ( istat /= 0 ) CALL FATAL(Caller,'Memory allocation error')
4927
      AllocationsDone = .TRUE.
4928
      STIFF = 0.0_dp
4929
      DAMP = 0.0_dp
4930
      MASS = 0.0_dp
4931
      FORCE = 0.0_dp
4932
      ColInds = 0
4933
      RowInds = 0
4934
    END IF
4935
      
4936
    CALL Info(Caller,'Starting block system assembly',Level=8)
4937
    
4938
    BulkMode = .TRUE.
4939
    
4940
100 IF(BulkMode) THEN
4941
      ! CALL Info(Caller,'Starting bulk assembly',Level=5)
4942
      ElementsFirst = 1
4943
      ElementsLast = Mesh % NumberOFBulkElements
4944
    ELSE
4945
      ! CALL Info(Caller,'Starting boundary assembly',Level=5)
4946
      ElementsFirst = Mesh % NumberOFBulkElements+1
4947
      ElementsLast = Mesh % NumberOfBulkElements + &
4948
          Mesh % NumberOFBoundaryElements
4949
    END IF
4950
        
4951
      
4952
    ! Load the assembly procedure
4953
    !-----------------------------------------
4954
    IF( BulkMode ) THEN
4955
      WRITE (str,'(A,I0,I0)') 'Bulk Assembly Procedure ',RowVar,ColVar
4956
    ELSE
4957
      WRITE (str,'(A,I0,I0)') 'Boundary Assembly Procedure ',RowVar,ColVar
4958
    END IF    
4959
    ProcName = ListGetString( SolverParams, TRIM(str), GotIt )
4960

4961
    ! Test if the 11 block is not given with 'ij' indexes
4962
    !--------------------------------------------------------
4963
    IF(.NOT. GotIt .AND. RowVar == 1 .AND. ColVar == 1) THEN
4964
      IF( BulkMode ) THEN
4965
        WRITE (str,'(A)') 'Bulk Assembly Procedure'
4966
      ELSE
4967
        WRITE (str,'(A)') 'Boundary Assembly Procedure'
4968
      END IF
4969
      ProcName = ListGetString( SolverParams, TRIM(str), GotIt )
4970
      IF(.NOT. GotIt) THEN
4971
         CALL Fatal(Caller,'Bulk Assembly Procedure not given!')
4972
      END IF
4973
    END IF
4974

4975
    Diagonal = (RowVar == ColVar) 
4976
    
4977
    IF( .NOT. GotIt ) THEN
4978
      IF( BulkMode .AND. Diagonal ) THEN
4979
        CALL Warn(Caller,'Diagonal bulk entries should be assembled!')
4980
      END IF
4981
      RETURN
4982
    END IF
4983

4984
    ProcPntr = GetProcAddr( TRIM(ProcName), abort=.FALSE.)
4985
    IF ( ProcPntr == 0 ) THEN
4986
      CALL Fatal(Caller,'Assembly routine not found: '//TRIM(ProcName))
4987
    ELSE
4988
      CALL Info(Caller,'Using assembly routine: '//TRIM(ProcName),Level=8)
4989
    END IF
4990
    
4991
    ! These may be fetched within the assembly routine, if needed
4992
    CALL ListAddInteger( SolverParams,'Block Matrix Row',RowVar )
4993
    CALL ListAddInteger( SolverParams,'Block Matrix Column',ColVar )
4994

4995

4996
    ! The assembly loop for a submatrix starts here
4997
    !------------------------------------------------------------------------------             
4998
    DO t=ElementsFirst,ElementsLast
4999
      
5000
      Element => Mesh % Elements(t)
5001
      CurrentModel % CurrentElement => Element
5002
      
5003
      !-----------------------------------------------------------------
5004
      n  = GetElementNOFNodes()
5005
      nd = GetElementDOFs(Indexes)
5006
!     Indexes => Element % NodeIndexes
5007
!     n = Element % TYPE % NumberOfnodes
5008
!     nd = n
5009
      Nrow = nd
5010

5011
      RowInds(1:Nrow) = RowPerm(Indexes(1:nd))
5012
      IF(.NOT. ALL(RowInds(1:Nrow) > 0)) CYCLE
5013

5014
      Ncol = nd
5015
      ColInds(1:Ncol) = ColPerm(Indexes(1:nd))
5016
      IF(.NOT. ALL(ColInds(1:Ncol) > 0)) CYCLE                 
5017

5018
      ! Here just the matrix structure is set, no values
5019
      !---------------------------------------------------------------
5020
      IF( IsListMatrix ) THEN
5021
        CALL GlueLocalSubMatrix( Amat, &
5022
            RowInd0,ColInd0,Nrow,Ncol,RowInds,ColInds,&
5023
            RowDofs,ColDofs,STIFF )
5024
        IF( AssemblySymmetric .OR. AssemblyAntiSymmetric ) THEN
5025
          CALL GlueLocalSubMatrix( Amat, &
5026
              ColInd0,RowInd0,Ncol,Nrow,ColInds,RowInds,&
5027
              ColDofs,RowDofs,STIFF )          
5028
        END IF
5029
        CYCLE
5030
      END IF
5031
      
5032
      ! Do the assembly, now really
5033
      !---------------------------------------------------------------
5034
      STIFF = 0.0_dp
5035
      DAMP = 0.0_dp
5036
      MASS = 0.0_dp
5037
      FORCE = 0.0_dp
5038
      
5039
      CALL ExecLocalAssembly( ProcPntr, CurrentModel, Solver, &
5040
          dt, Transient, MASS, DAMP, STIFF, FORCE, Element, &
5041
          Nrow, Ncol )
5042

5043
      IF ( Transient ) THEN
5044
        IF( Solver % TimeOrder == 1 ) THEN
5045
          CALL Default1stOrderTime( MASS,STIFF,FORCE)
5046
        ELSE IF( Solver % TimeOrder == 2) THEN
5047
          CALL Default2ndOrderTime( MASS,DAMP,STIFF,FORCE )
5048
        END IF
5049
      END IF
5050
      
5051
      IF ( Solver % NOFEigenValues > 0 ) THEN
5052
        IF( Solver % TimeOrder == 1 ) THEN
5053
          CALL DefaultUpdateMass(MASS)
5054
        ELSE IF( Solver % TimeOrder == 2) THEN
5055
          CALL DefaultUpdateDamp(DAMP)
5056
          CALL DefaultUpdateMass(MASS)
5057
        END IF
5058
      END IF
5059

5060
      
5061
      ! Assemble the matrix with offset
5062
      !-----------------------------------------------
5063
      IF( .TRUE. .OR. ColInd0 > 0 .OR. RowInd0 > 0 ) THEN
5064
        CALL GlueLocalSubMatrix( Amat, &
5065
            RowInd0,ColInd0,Nrow,Ncol,RowInds,ColInds,&
5066
            RowDofs,ColDofs,STIFF )
5067
        
5068
        ! Assemble the r.h.s with offset
5069
        !-----------------------------------------------
5070
        DO i=1,Nrow
5071
          DO j=1,RowDofs
5072
            Row = RowInd0 + RowDofs * (RowInds(i)-1) + j
5073
            ForceVector(Row) = ForceVector(Row) + &
5074
                FORCE(RowDofs*(i-1)+j)
5075
          END DO
5076
        END DO
5077
        
5078
        ! For some constraints assemble also the transpose
5079
        !--------------------------------------------------
5080
        IF( AssemblySymmetric ) THEN
5081
          CALL GlueLocalSubMatrix( Amat, &
5082
              ColInd0,RowInd0,Ncol,Nrow,ColInds,RowInds,&
5083
              ColDofs,RowDofs,TRANSPOSE(STIFF) )               
5084
        ELSE IF( AssemblyAntisymmetric ) THEN
5085
          CALL GlueLocalSubMatrix( Amat, &
5086
              ColInd0,RowInd0,Ncol,Nrow,ColInds,RowInds,&
5087
              ColDofs,RowDofs,-TRANSPOSE(STIFF) )               
5088
        END IF
5089
      ELSE              
5090
        CALL DefaultUpdateEquations( STIFF, FORCE )
5091
      END IF
5092
      
5093
    END DO
5094
    
5095
    IF(BulkMode) THEN
5096
      ! CALL Info( Caller, 'Bulk assembly done for blocks', Level=4 )
5097
      BulkMode = .FALSE.
5098
!      GOTO 100
5099
    ELSE 
5100
      ! CALL Info( Caller, 'Boundary assembly done for blocks', Level=4 )
5101
    END IF
5102

5103
  END SUBROUTINE BlockSystemAssembly
5104
!-----------------------------------------------------------------------------------
5105

5106

5107
  SUBROUTINE ExecSolverInSteps( Model, Solver, dt, TransientSimulation )
5108
!------------------------------------------------------------------------------
5109
    TYPE(Model_t)  :: Model
5110
    TYPE(Solver_t),POINTER :: Solver
5111
    LOGICAL :: TransientSimulation
5112
    REAL(KIND=dp) :: dt
5113
!------------------------------------------------------------------------------
5114
    INTEGER(KIND=AddrInt) :: SolverAddr
5115
    CHARACTER(LEN=MAX_NAME_LEN) :: ProcName
5116
    INTEGER :: iter, MaxIter
5117
    LOGICAL :: Found
5118
!------------------------------------------------------------------------------
5119
    INTEGER :: NColours, col
5120
    REAL(KIND=dp) :: Norm
5121

5122
    CALL Info('ExecSolverInSteps','Performing solution in steps',Level=6)
5123
    ProcName = ListGetString( Solver % Values,'Procedure', Found )
5124

5125
    MaxIter = ListGetInteger( Solver % Values,'Nonlinear System Max Iterations', Found ) 
5126
    IF( .NOT. Found ) MaxIter = 1
5127

5128
    DO iter = 1, MaxIter
5129
      CALL DefaultInitialize( Solver )
5130
      
5131
      IF( ASSOCIATED( Solver % ColourIndexList ) ) THEN
5132
        ncolours = Solver % ColourIndexList % n
5133
      ELSE
5134
        ncolours = 1 
5135
      END IF
5136
      Solver % CurrentColour = 0
5137

5138
      SolverAddr = Solver % PROCEDURE
5139
      DO col=1,ncolours
5140
        ! The > CurrentColour < is advanced by GetNOFActive() routine
5141
        ! to have similar interface as for non-steps solver
5142
        CALL ExecSolver( SolverAddr, Model, Solver, dt, TransientSimulation)
5143
      END DO
5144
      CALL DefaultFinishBulkAssembly( Solver )
5145

5146
      SolverAddr = GetProcAddr( TRIM(ProcName)//'_boundary', abort=.FALSE. )
5147
      IF( SolverAddr /= 0 ) THEN
5148
        CALL ExecSolver( SolverAddr, Model, Solver, dt, TransientSimulation)
5149
      END IF
5150

5151
      CALL DefaultFinishBoundaryAssembly( Solver )
5152
      CALL DefaultFinishAssembly( Solver )
5153
      CALL DefaultDirichletBCs( Solver )
5154

5155
      Norm = DefaultSolve( Solver )
5156
      
5157
      IF( Solver % Variable % NonlinConverged > 0 ) EXIT
5158
    END DO
5159

5160
    !SolverAddr = GetProcAddr( TRIM(ProcName)//'_post', abort=.FALSE. )
5161
    !IF( SolverAddr /= 0 ) THEN
5162
    !  CALL ExecSolver( SolverAddr, Model, Solver, dt, TransientSimulation)
5163
    !END IF
5164
    
5165

5166
  END SUBROUTINE ExecSolverInSteps
5167

5168
  
5169
!------------------------------------------------------------------------------
5170
!> This executes the original line of solvers (legacy solvers) where each solver 
5171
!> includes looping over elements and the convergence control. From generality
5172
!> point of view this misses some opportunities to have control of the nonlinear
5173
!> system. 
5174
!------------------------------------------------------------------------------
5175
  RECURSIVE SUBROUTINE SingleSolver( Model, Solver, dt, TransientSimulation )
5176
!------------------------------------------------------------------------------
5177
     TYPE(Model_t)  :: Model
5178
     TYPE(Solver_t),POINTER :: Solver
5179
     LOGICAL :: TransientSimulation
5180
     REAL(KIND=dp) :: dt
5181
!------------------------------------------------------------------------------
5182
     LOGICAL :: stat, Found, GB, MeActive, GotIt
5183
     INTEGER :: i, j, k, l, col, row, n, BDOFs, maxdim, dsize, size0
5184
     TYPE(Element_t), POINTER :: CurrentElement
5185
     TYPE(ValueList_t), POINTER :: SolverParams
5186
     INTEGER(KIND=AddrInt) :: SolverAddr
5187
     CHARACTER(:), ALLOCATABLE :: EquationName
5188

5189
     INTEGER, ALLOCATABLE :: memb(:)
5190
     TYPE(Matrix_t), POINTER :: M
5191
     INTEGER :: comm_active, group_active, group_world, ierr
5192

5193
     LOGICAL :: ApplyMortar, FoundMortar, SlaveNotParallel, Parallel, UseOrigMesh
5194
     TYPE(Matrix_t), POINTER :: CM, CM0, CM1, CMP
5195
     TYPE(Mesh_t), POINTER :: Mesh
5196
     LOGICAL :: DoBC, DoBulk
5197
     CHARACTER(*), PARAMETER :: Caller="SingleSolver"   
5198
     
5199
!------------------------------------------------------------------------------
5200
     MeActive = ASSOCIATED(Solver % Matrix)
5201
     IF ( MeActive ) MeActive = (Solver % Matrix % NumberOfRows > 0)
5202

5203
     Parallel = Solver % Parallel 
5204
     !------------------------------------------------------------------------------
5205

5206
     IF( Solver % Mesh % Changed .OR. Solver % NumberOfActiveElements <= 0 ) THEN
5207
       Solver % NumberOFActiveElements = 0
5208
       EquationName = ListGetString( Solver % Values, 'Equation', Found)
5209

5210
       IF ( Found ) THEN
5211

5212
         CALL SetActiveElementsTable( Model, Solver, MaxDim  ) 
5213
         CALL ListAddInteger( Solver % Values, 'Active Mesh Dimension', Maxdim )
5214
         
5215
         ! Calculate accumulated integration weights for bulk if requested          
5216
         DoBulk = ListGetLogical( Solver % Values,'Calculate Weights',Found )
5217
         ! Calculate weight for boundary 
5218
         DoBC = ListGetLogical( Solver % Values,'Calculate Boundary Weights',Found )
5219

5220
         ! In parallel we have to prepare the communicator already for the weights
5221
         IF(DoBulk .OR. DoBC ) THEN
5222
           IF ( Parallel .AND. MeActive ) THEN
5223
             IF ( ASSOCIATED(Solver % Mesh % ParallelInfo % GInterface) ) THEN
5224
               IF (.NOT. ASSOCIATED(Solver % Matrix % ParMatrix) ) &
5225
                   CALL ParallelInitMatrix(Solver, Solver % Matrix )               
5226
               ParEnv => Solver % Matrix % ParMatrix % ParEnv
5227
               ParEnv % ActiveComm = Solver % Matrix % Comm
5228
             END IF
5229
           END IF           
5230
         END IF
5231
         
5232
         IF(DoBulk) CALL CalculateNodalWeights(Solver,.FALSE.)
5233
         IF(DoBC) CALL CalculateNodalWeights(Solver,.TRUE.) 
5234
       END IF
5235
     END IF
5236
!------------------------------------------------------------------------------
5237
     UseOrigMesh = ListGetLogical(Solver % Values,'Use Original Coordinates',Found )
5238
     IF(UseOrigMesh ) THEN
5239
       Mesh => Solver % Mesh
5240
       IF(.NOT. ASSOCIATED(Mesh % NodesOrig)) THEN
5241
         CALL Fatal(Caller,'Cannot toggle between meshes: NodesOrig not associated!')
5242
       END IF
5243
       IF(.NOT. ASSOCIATED(Mesh % NodesMapped)) THEN
5244
         CALL Fatal(Caller,'Cannot toggle between meshes: NodesMapped not associated!')
5245
       END IF
5246
       Mesh % Nodes => Mesh % NodesOrig
5247
       CALL Info(Caller,'Using stored original coordinate in solver')
5248
     END IF
5249

5250
     SlaveNotParallel = ListGetLogical( Solver % Values, 'Slave not parallel',Found )
5251

5252
     IF ( Parallel .AND. .NOT. SlaveNotParallel ) THEN
5253
       ! Set the communicator and active info partitions.
5254

5255
BLOCK
5256
       LOGICAL :: ChangedActiveParts
5257

5258
       ChangedActiveParts = .FALSE.
5259

5260
       !block partitions containing ONLY halos
5261
       IF( MeActive ) THEN
5262
         IF ( ListGetLogical( Solver % Values, 'Skip Halo Only Partitions', Found) ) THEN
5263
           MeActive = .FALSE.
5264
           DO i=1,Solver % NumberOfActiveElements
5265
             IF(Solver % Mesh % Elements(Solver % ActiveElements(i)) % PartIndex==ParEnv % myPE) THEN
5266
               MeActive = .TRUE.; EXIT
5267
             END IF
5268
           END DO
5269
           IF(.NOT. MeActive) ChangedActiveParts = .TRUE.
5270
         END IF
5271
       END IF
5272

5273
       IF ( ASSOCIATED(Solver  % Matrix) ) THEN
5274
          IF ( ASSOCIATED(Solver  % Matrix % ParMatrix) ) THEN
5275
            ParEnv => Solver % Matrix % ParMatrix % ParEnv
5276
          END IF
5277
       END IF
5278

5279
       CALL ParallelActive( MeActive )
5280
       n = COUNT(ParEnv % Active)
5281

5282
       IF ( n>0 .AND. n<ParEnv % PEs ) THEN
5283
         IF ( ASSOCIATED(Solver % Matrix) ) THEN
5284
           IF ( Solver % Matrix % Comm /= ELMER_COMM_WORLD .AND. Solver % Matrix % Comm /= MPI_COMM_NULL ) &
5285
             CALL MPI_Comm_Free( Solver % Matrix % Comm, ierr )
5286
         END IF
5287

5288
         CALL MPI_Comm_group( ELMER_COMM_WORLD, group_world, ierr )
5289
         ALLOCATE(memb(n))
5290
         n = 0
5291
         DO i=1,ParEnv % PEs
5292
           IF ( ParEnv % Active(i) ) THEN
5293
             n=n+1
5294
             memb(n)=i-1
5295
           END IF
5296
         END DO
5297
         CALL MPI_Group_incl( group_world, n, memb, group_active, ierr)
5298
         DEALLOCATE(memb)
5299
         CALL MPI_Comm_create( ELMER_COMM_WORLD, group_active, &
5300
                 comm_active, ierr)
5301

5302
         M => Solver % Matrix
5303
         DO WHILE(ASSOCIATED(M))
5304
           M % Comm = comm_active
5305
           M => M % Parent
5306
         END DO
5307

5308
         IF( ANY( ParEnv % Active(MinOutputPE+1:MIN(MaxOutputPE+1,ParEnv % PEs)) ) ) THEN
5309
           ! If any of the active output partitions in active just use it.
5310
           ! Typically the 1st one. Others are passive. 
5311
           IF( ParEnv % MyPe >= MinOutputPE .AND. ParEnv % MyPe <= MaxOutputPE ) THEN 
5312
             OutputPE = ParEnv % MyPE
5313
           ELSE
5314
             OutputPE = -1
5315
           END IF
5316
         ELSE         
5317
           ! Otherwise find the 1st active partition and if found use it.
5318
           ! Otherwise use the 0:th partition. 
5319
           DO i=1,ParEnv % PEs
5320
             IF ( ParEnv % Active(i) ) EXIT
5321
           END DO
5322

5323
           OutputPE = -1
5324
           IF ( i-1 == ParEnv % MyPE ) THEN
5325
             OutputPE = i-1 
5326
           ELSE IF( i > ParEnv % PEs .AND. ParEnv % myPE == 0 ) THEN
5327
             OutputPE = 0
5328
           END IF
5329
         END IF
5330
       ELSE
5331
         M => Solver % Matrix
5332
         DO WHILE( ASSOCIATED(M) )
5333
           M % Comm = ELMER_COMM_WORLD
5334
           M => M % Parent
5335
         END DO
5336

5337
         IF(.NOT.ASSOCIATED(Solver % Matrix)) ParEnv % Active = .TRUE.
5338

5339
         ! Here set the default partitions active. 
5340
         IF( ParEnv % MyPe >= MinOutputPE .AND. &
5341
             ParEnv % MyPe <= MaxOutputPE ) THEN 
5342
           OutputPE = ParEnv % MyPE
5343
         ELSE
5344
           OutputPE = -1
5345
         END IF
5346
       END IF
5347

5348
       ! POTENTIAL INCOMPATIBILITY: don't execute solvers for non-active partitions
5349
       ! (usually this is done within the solver by:
5350
       ! IF (.NOT. ASSOCIATED(Solver % Matrix) ) RETURN
5351
       ! or some such .... )
5352

5353
       IF(.NOT. MeActive .AND. n>0) RETURN
5354
END BLOCK
5355
     END IF
5356

5357
     IF( ListGetLogical( Solver % Values,'CutFEM',Found ) ) GOTO 1
5358

5359
     IF ( ASSOCIATED(Solver % Matrix) ) THEN
5360
       IF ( Parallel .AND. MeActive ) THEN
5361
         IF ( ASSOCIATED(Solver % Mesh % ParallelInfo % GInterface) ) THEN
5362
           ParEnv % ActiveComm = Solver % Matrix % Comm
5363

5364
           IF (.NOT. ASSOCIATED(Solver % Matrix % ParMatrix) ) &
5365
             CALL ParallelInitMatrix(Solver, Solver % Matrix )
5366

5367
           ParEnv => Solver % Matrix % ParMatrix % ParEnv
5368
           ParEnv % ActiveComm = Solver % Matrix % Comm
5369

5370
#if 0
5371
           ! This one is mainly for debugging of parallel problems
5372
           BLOCK
5373
             TYPE(Variable_t), POINTER :: Gvar
5374
             
5375
             GVar => VariableGet( Solver % Mesh % Variables,&
5376
                 Solver % Variable % Name(Solver % Variable NameLen)//' Gdofs')      
5377
             IF ( ASSOCIATED(GVar) ) THEN               
5378
               DO i = 1, SIZE( GVar % Perm )
5379
                 j = GVar % Perm(i)
5380
                 IF( j == 0 ) CYCLE
5381
                 k = Solver % Matrix % ParallelInfo % GlobalDOFs(j)
5382
                 GVar % Values(j) = 1.0_dp * k
5383
               END DO
5384
             END IF
5385
           END BLOCK             
5386
#endif
5387
           
5388
         END IF
5389
       END IF
5390
     ELSE IF (.NOT.SlaveNotParallel) THEN
5391
       Parenv % ActiveComm = ELMER_COMM_WORLD
5392
     END IF
5393

5394

5395
     ! This is more featured version than the original one with just one flag.
5396
     ! This way different solvers can detect when their mesh has been updated. 
5397
1    Solver % MeshChanged = Solver % Mesh % Changed
5398
     IF( Solver % MeshTag /= Solver % Mesh % MeshTag ) THEN
5399
       Solver % MeshChanged = .TRUE.
5400
       Solver % MeshTag = Solver % Mesh % MeshTag
5401
     END IF       
5402
     
5403
     ! Linear constraints from mortar BCs:
5404
     ! -----------------------------------
5405
     IF(.NOT. ListGetLogical(Solver % Values,'Mortar BCs Fixed', Found ) ) THEN
5406
       CALL GenerateProjectors(Model,Solver,Nonlinear = .FALSE. )
5407
     END IF
5408
       
5409
     CALL Info(Caller, "Attempting to call solver: "//I2S(Solver % SolverId), level=8)
5410
     SolverParams => ListGetSolverParams(Solver)
5411
     EquationName = GetString(SolverParams, 'Equation', GotIt)
5412
     IF (GotIt) THEN
5413
        Message = 'Solver Equation string is: '//TRIM(EquationName)
5414
        CALL Info(Caller, Message, level=8)
5415
     END IF
5416

5417
     IF( Solver % SolverMode == SOLVER_MODE_STEPS ) THEN
5418
       CALL ExecSolverinSteps( Model, Solver, dt, TransientSimulation)
5419
     ELSE
5420
       SolverAddr = Solver % PROCEDURE
5421
       CALL ExecSolver( SolverAddr, Model, Solver, dt, TransientSimulation)
5422
     END IF
5423

5424

5425
     ! Special slot for post-processing solvers
5426
     ! This makes it convenient to separate the solution and postprocessing.
5427
     ! This solver must use the same structures as the primary solver.
5428
     ! If the postprocessing solver uses different element basis this is
5429
     ! not a good idea. 
5430
     !-----------------------------------------------------------------------
5431
     BLOCK 
5432
       CHARACTER(LEN=MAX_NAME_LEN) :: ProcName
5433
       LOGICAL :: PostActive
5434

5435
       PostActive = ListGetLogical( Solver % Values,'PostSolver Active',Found )
5436

5437
       IF( PostActive ) THEN
5438
         ProcName = ListGetString( Solver % Values,'Procedure', Found )
5439
         SolverAddr = GetProcAddr( TRIM(ProcName)//'_post', abort=.FALSE. )
5440
         IF( SolverAddr /= 0 ) THEN
5441
           CALL Info(Caller,'Calling solver for postprocessing',Level=10)
5442
           CALL ExecSolver( SolverAddr, Model, Solver, dt, TransientSimulation)
5443
         END IF
5444
       END IF
5445
     END BLOCK
5446

5447
     IF( ListGetLogical( Solver % Values,'Library Adaptivity', Found ) ) THEN
5448
#ifdef LIBRARY_ADAPTIVITY
5449
       ! Do adaptive meshing, whether to do this before or after "_post" is a matter  of taste i guess
5450
       BLOCK 
5451
         USE, INTRINSIC :: ISO_C_BINDING
5452

5453
         CHARACTER(LEN=MAX_NAME_LEN) :: ProcName
5454
         LOGICAL :: AdaptiveActive
5455
         TYPE(Variable_t), POINTER :: Var
5456
         INTEGER(KIND=AddrInt) :: IResidual, EResidual, BResidual
5457

5458
         INTERFACE
5459
           SUBROUTINE RefineMeshExt(Model,Solver,Quant,Perm,InsideResidual,EdgeResidual,BoundaryResidual)
5460
             USE Types
5461

5462
             TYPE( Model_t ) :: Model
5463
             TYPE(Solver_t), TARGET :: Solver
5464
             REAL(KIND=dp) :: Quant(:)
5465
             INTEGER :: Perm(:)
5466

5467
             INTERFACE
5468
                SUBROUTINE BoundaryResidual( Model,Edge,Mesh,Quant,Perm,Gnorm,Indicator )
5469
                   USE Types
5470
                   TYPE(Element_t), POINTER :: Edge
5471
                   TYPE(Model_t) :: Model
5472
                   TYPE(Mesh_t), POINTER :: Mesh
5473
                   REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm
5474
                   INTEGER :: Perm(:)
5475
                END SUBROUTINE BoundaryResidual
5476

5477
                SUBROUTINE EdgeResidual( Model,Edge,Mesh,Quant,Perm, Indicator)
5478
                   USE Types
5479
                   TYPE(Element_t), POINTER :: Edge
5480
                   TYPE(Model_t) :: Model
5481
                   TYPE(Mesh_t), POINTER :: Mesh
5482
                   REAL(KIND=dp) :: Quant(:), Indicator(2)
5483
                   INTEGER :: Perm(:)
5484
                END SUBROUTINE EdgeResidual
5485

5486

5487
                SUBROUTINE InsideResidual( Model,Element,Mesh,Quant,Perm,Fnorm, Indicator)
5488
                   USE Types
5489
                   TYPE(Element_t), POINTER :: Element
5490
                   TYPE(Model_t) :: Model
5491
                   TYPE(Mesh_t), POINTER :: Mesh
5492
                   REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm
5493
                   INTEGER :: Perm(:)
5494
                END SUBROUTINE InsideResidual
5495
             END INTERFACE
5496
           END SUBROUTINE RefineMeshExt
5497

5498
           SUBROUTINE BoundaryResidual( Model,Edge,Mesh,Quant,Perm,Gnorm,Indicator)
5499
             USE Types
5500
             TYPE(Element_t), POINTER :: Edge
5501
             TYPE(Model_t) :: Model
5502
             TYPE(Mesh_t), POINTER :: Mesh
5503
             REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm
5504
             INTEGER :: Perm(:)
5505
           END SUBROUTINE BoundaryResidual
5506

5507
           SUBROUTINE EdgeResidual( Model,Edge,Mesh,Quant,Perm,Indicator)
5508
             USE Types
5509
             TYPE(Element_t), POINTER :: Edge
5510
             TYPE(Model_t) :: Model
5511
             TYPE(Mesh_t), POINTER :: Mesh
5512
             REAL(KIND=dp) :: Quant(:), Indicator(2)
5513
             INTEGER :: Perm(:)
5514
           END SUBROUTINE EdgeResidual
5515

5516
           SUBROUTINE InsideResidual( Model,Element,Mesh,Quant,Perm,Fnorm,Indicator)
5517
             USE Types
5518
             TYPE(Element_t), POINTER :: Element
5519
             TYPE(Model_t) :: Model
5520
             TYPE(Mesh_t), POINTER :: Mesh
5521
             REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm
5522
             INTEGER :: Perm(:)
5523
           END SUBROUTINE InsideResidual
5524
         END INTERFACE
5525

5526
         PROCEDURE(InsideResidual), POINTER :: InsidePtr
5527
         PROCEDURE(EdgeResidual), POINTER :: EdgePtr
5528
         PROCEDURE(BoundaryResidual), POINTER :: BoundaryPtr
5529

5530
         TYPE(C_FUNPTR) :: IResFunC, EResFunC, BresFunC
5531

5532
         AdaptiveActive = ListGetLogical(Solver % Values, 'Adaptive Mesh Refinement', Found)
5533

5534
         IF (AdaptiveActive) THEN
5535
           ProcName = ListGetString( Solver % Values,'Procedure', Found )
5536

5537
           IResidual = GetProcAddr( TRIM(ProcName)//'_inside_residual', abort=.FALSE. )
5538
           EResidual   = GetProcAddr( TRIM(ProcName)//'_edge_residual', abort=.FALSE. )
5539
           BResidual   = GetProcAddr( TRIM(ProcName)//'_boundary_residual', abort=.FALSE. )
5540

5541
           IF( IResidual/=0 .AND. EResidual /= 0 .AND. BResidual /= 0 ) THEN
5542
             IResFunC = TRANSFER( Iresidual, IresFunC )
5543
             EResFunC = TRANSFER( Eresidual, EresFunC )
5544
             BResFunC = TRANSFER( Bresidual, BresFunC )
5545

5546
             CALL C_F_PROCPOINTER(IresFunC, InsidePtr)
5547
             CALL C_F_PROCPOINTER(EResFunC, EdgePtr )
5548
             CALL C_F_PROCPOINTER(BResFunC, BoundaryPtr )
5549

5550
             Var => Solver % Variable
5551
             CALL RefineMeshExt( Model, Solver, Var % Values, Var % Perm, InsidePtr, EdgePtr, BoundaryPtr )
5552
           END IF
5553
         END IF
5554
       END BLOCK
5555
#else
5556
       CALL Fatal(Caller,'Library version of adaptivity residuals not compiled with!')
5557
#endif
5558
     END IF            
5559

5560
     ! Compute all dependent fields, components and derivatives related to the primary solver.
5561
     !-----------------------------------------------------------------------   
5562
     CALL UpdateDependentObjects( Solver, .TRUE. ) 
5563

5564
     IF( UseOrigMesh ) THEN
5565
       CALL Info(Caller,'Reverting back to current coordinates',Level=12)
5566
       Mesh % Nodes => Mesh % NodesMapped
5567
     END IF
5568
!------------------------------------------------------------------------------
5569
   END SUBROUTINE SingleSolver
5570
!------------------------------------------------------------------------------
5571

5572

5573
!------------------------------------------------------------------------------
5574
!> Runs a solver as defined in the command file. There are several
5575
!> ways how to skip the execution. The are also three different ways
5576
!> how the matrices may be assembled and solved: standard (single), coupled and
5577
!> block. 
5578
!------------------------------------------------------------------------------
5579
  RECURSIVE SUBROUTINE SolverActivate( Model, Solver, dt, TransientSimulation )
5580
!------------------------------------------------------------------------------
5581
     TYPE(Model_t)  :: Model
5582
     TYPE(Solver_t),POINTER :: Solver
5583
     LOGICAL :: TransientSimulation
5584
     REAL(KIND=dp) :: dt
5585
!------------------------------------------------------------------------------
5586
     REAL(KIND=dp) :: OrigDT, DTScal
5587
     LOGICAL :: stat, Found, TimeDerivativeActive, Timing, IsPassiveBC, &
5588
         GotCoordTransform, NamespaceFound
5589
     INTEGER :: i, j, k, n, BDOFs, timestep, timei, timei0, PassiveBcId, Execi
5590
     INTEGER, POINTER :: ExecIntervals(:),ExecIntervalsOffset(:)
5591
     REAL(KIND=dp) :: tcond, t0, rt0, st, rst, ct
5592
     TYPE(Variable_t), POINTER :: TimeVar, IterV
5593
     TYPE(ValueList_t), POINTER :: Params
5594
     INTEGER, POINTER :: UpdateComponents(:)
5595

5596
     INTEGER :: ScanningLoops, scan, sOutputPE
5597
     LOGICAL :: GotLoops
5598
     TYPE(Variable_t), POINTER :: ScanVar, Var
5599
     CHARACTER(:), ALLOCATABLE :: str, CoordTransform
5600
     TYPE(Mesh_t), POINTER :: Mesh, pMesh
5601

5602
     SAVE TimeVar
5603
!------------------------------------------------------------------------------
5604
     sOutputPE = OutputPE
5605

5606

5607
     Mesh => Solver % Mesh 
5608

5609
     IF( ASSOCIATED( Mesh % Child ) .AND. .NOT. Mesh % OutputActive ) THEN
5610
       i = 0
5611
       pMesh => Mesh
5612
       DO WHILE( ASSOCIATED( pMesh % Child ) )
5613
         pMesh => pMesh % Child 
5614
         i = i+1 
5615
         IF(pMesh % OutputActive) EXIT 
5616
       END DO
5617
       IF( .NOT. ASSOCIATED(pMesh,Mesh) .AND. ASSOCIATED(pMesh) ) THEN
5618
         IF( .NOT. ListCheckPresent( Solver % Values,'Relative Mesh Level') ) THEN
5619
           CALL Info('SolverActivate','By some logic the mesh is switched here to child mesh!!!')
5620
           CALL Info('SolverActivate','Changing Solver '//I2S(Solver % SolverId)//&
5621
               ' mesh to be the '//TRIM(I2S(i))//'th Child mesh: '&
5622
               //TRIM(pMesh % Name),Level=7)
5623
           Solver % Mesh => pMesh
5624
         END IF
5625
       END IF
5626
     END IF
5627
     
5628
     CALL SetCurrentMesh( Model, Solver % Mesh )
5629

5630
     Model % Solver => Solver
5631
     Params => ListGetSolverParams(Solver)
5632

5633
     CoordTransform = ListGetString(Params,'Coordinate Transformation',&
5634
         GotCoordTransform )
5635
     IF( GotCoordTransform ) THEN
5636
       CALL CoordinateTransformation( Solver % Mesh, CoordTransform, &
5637
           Params, .FALSE. )
5638
     END IF
5639

5640
     ! Some properties might be rotated by "Property Rotate"
5641
     !--------------------------------------------------------------
5642
     CALL SetRotatedProperties(Model, Solver % Mesh)
5643

5644
     ! Some keywords may be normalized by area or volume
5645
     !--------------------------------------------------------------
5646
     CALL SetNormalizedKeywords(Model, Solver % Mesh)
5647

5648
!------------------------------------------------------------------------------
5649
! The solver may be skipped if the exec condition is negative. 
5650
! This allows for complicated conditions, and the earlier simple ones 
5651
! have therefore been replaced by this keyword.
5652
!------------------------------------------------------------------------------
5653
     IF( ListCheckPresent( Params,'Start Time') ) &
5654
       CALL Fatal('SolverActivate','Use > Exec Condition = Real < instead of > Start Time <')
5655

5656
     IF( ListCheckPresent( Params,'Stop Time') ) & 
5657
       CALL Fatal('SolverActivate','Use > Exec Condition = Real < instead of > Stop Time <')
5658

5659
     st = ListGetCReal( Params,'Exec Condition',Found) 
5660
     IF( Found .AND. st < 0.0 ) RETURN
5661

5662
!---------------------------------------------------------------------------------   
5663
! There may also be predefined discrete intervals for the execution of the solver.
5664
!---------------------------------------------------------------------------------   
5665
     execi = 1
5666
     ExecIntervals => ListGetIntegerArray( Params,'Exec Intervals', Found )
5667
     IF( .NOT. Found ) THEN
5668
       ExecIntervals =>  ListGetIntegerArray( Params,'Exec Interval', Found )
5669
     END IF
5670
     IF ( Found ) THEN
5671
       TimeVar => VariableGet( Model % Variables, 'Timestep Interval' )
5672
       timei = NINT(Timevar % Values(1))
5673

5674
       IF( ExecIntervals(timei) == 0 ) RETURN
5675

5676
       TimeVar => VariableGet( Model % Variables, 'Timestep' )
5677
       timestep = NINT(TimeVar % Values(1))
5678

5679
       ExecIntervalsOffset =>  ListGetIntegerArray( Params,&
5680
           'Exec Intervals Offset', Found )
5681
       IF( Found ) THEN
5682
         timei0 = ExecIntervalsOffset(timei)
5683
       ELSE
5684
         timei0 = 0
5685
       END IF
5686
       
5687
       execi = ExecIntervals(timei)
5688
       IF( MOD( timestep-1-timei0, execi) /= 0 ) RETURN              
5689
     END IF
5690

5691
!-------------------------------------------------------------------------------
5692
! Set solver parameters to avoid list operations during assembly
5693
!-------------------------------------------------------------------------------
5694
     Solver % DG = ListGetLogical(Params, 'Discontinuous Galerkin', Found)
5695
     Solver % GlobalBubbles = ListGetLogical(Params, 'Bubbles in Global System', Found)
5696
     IF(.NOT. Found) Solver % GlobalBubbles = .TRUE.
5697
     IF(GetString(Params, 'Linear System Direct Method', Found) == 'permon') THEN
5698
       Solver % DirectMethod = DIRECT_PERMON
5699
     END IF
5700
     
5701
     str = ListGetString( Params, 'Boundary Element Procedure', Found)
5702
     IF(Found) THEN
5703
       Solver % BoundaryElementProcedure = GetProcAddr( Str, abort=.FALSE., quiet=.TRUE. )
5704
     ELSE
5705
       Solver % BoundaryElementProcedure = 0
5706
     END IF
5707

5708
     str = ListGetString( Params, 'Bulk Element Procedure', Found)
5709
     IF(Found) THEN
5710
       Solver % BulkElementProcedure = GetProcAddr( Str, abort=.FALSE., quiet=.TRUE. )
5711
     ELSE
5712
       Solver % BulkElementProcedure = 0
5713
     END IF
5714

5715
!------------------------------------------------------------------------------
5716
! If solver timing is requested start the watches
5717
!------------------------------------------------------------------------------
5718
     Timing = ListCheckPrefix(Params,'Solver Timing')
5719
     IF( Timing ) THEN
5720
       t0 = CPUTime()
5721
       rt0 = RealTime()
5722
     END IF
5723

5724
     IF (.NOT. Solver % Mesh % AdaptiveDepth > 0) Solver % Mesh % OutputActive = .TRUE.
5725
     TimeDerivativeActive = TransientSimulation
5726

5727
!----------------------------------------------------------------------
5728
! This is to avoid resetting of certain info that could be interesting
5729
! when saving data i.e. using an auxiliary solver.
5730
!----------------------------------------------------------------------
5731
     IF(.NOT. ListGetLogical( Params,'Auxiliary Solver',Found)) THEN
5732
       DTScal = ListGetConstReal( Params, 'Timestep Scale', Found )
5733
       IF ( .NOT. Found ) DTScal = 1.0_dp
5734

5735
       IF( ListGetLogical( Params,'Timestep Over Intervals',Found) ) THEN
5736
         DTScal = 1.0_dp * Execi
5737
       END IF
5738

5739
       Solver % dt = DtScal * dt 
5740

5741
       IF ( TransientSimulation ) THEN
5742
         TimeDerivativeActive = &
5743
           ListGetLogical( Params, 'Time Derivative Active', Found )
5744

5745
         IF ( .NOT. Found ) THEN
5746
           TimeDerivativeActive = .TRUE.
5747
           tcond = ListGetCReal(Params,'Time Derivative Condition',Found)
5748
           IF ( Found ) TimeDerivativeActive = TimeDerivativeActive .AND. tcond>0
5749
         END IF
5750
       END IF
5751

5752
       str = ListGetString( Params, 'Namespace', NamespaceFound )
5753
       IF (NamespaceFound) CALL ListPushNamespace(TRIM(str))
5754
     END IF
5755

5756
 !------------------------------------------------------------------------------
5757
 ! Check for passive-active boundaries
5758
 !------------------------------------------------------------------------------
5759
     PassiveBcId = 0
5760
     IsPassiveBC = .FALSE.
5761
     DO j=1,Model % NumberOfBCs
5762
       IsPassiveBC = ListGetLogical( Model % BCs(j) % Values, &
5763
           'Passive Target',Found)
5764
       IF (IsPassiveBC) THEN
5765
         PassiveBcId = j
5766
         EXIT
5767
       END IF
5768
     END DO
5769
     IF ( IsPassiveBC ) THEN
5770
       CALL GetPassiveBoundary( Model, Model % Mesh, PassiveBcId )
5771
       Message = 'Passive element BC no. '//I2S(j)//' assigned to BC-ID no. '// &
5772
            I2S(PassiveBcId)
5773
       CALL Info('MainUtils',Message,Level=6)
5774
     END IF
5775

5776
     ScanningLoops = ListGetInteger( Params,'Scanning Loops',GotLoops)
5777
     IF( GotLoops ) THEN
5778
       ScanVar => VariableGet( Solver % Mesh % Variables,'scan')
5779
       IF(.NOT. ASSOCIATED( ScanVar ) ) THEN
5780
         CALL Fatal('SolverActivate','For scanning we should have scanning variable!')
5781
       END IF
5782
     ELSE
5783
       ScanningLoops = 1
5784
     END IF
5785

5786
     CALL SwapRefElemNodes( ANY(Solver % Def_Dofs(:,:,6)>0) )
5787

5788
     DO scan = 1, ScanningLoops        
5789
       !----------------------------------------------------------------------
5790
       ! This is to avoid resetting iteration that may be interesting and we
5791
       ! don't want to override it.
5792
       !----------------------------------------------------------------------
5793
       IF(.NOT. ListGetLogical( Params,'Auxiliary Solver',Found)) THEN
5794
         iterV => VariableGet( Solver % Mesh % Variables, 'nonlin iter' )
5795
         IF(ASSOCIATED(iterV)) iterV % Values(1) = 1
5796
       END IF
5797

5798
       IF( GotLoops ) THEN
5799
         ScanVar % Values(1) = scan
5800
       END IF
5801
       
5802
       !---------------------------------------------------------------------
5803
       ! Call the correct type of solver: standard (single), coupled or block
5804
       ! This is where everything really happens!
5805
       !---------------------------------------------------------------------
5806
       IF( Solver % SolverMode == SOLVER_MODE_COUPLED .OR. &
5807
           Solver % SolverMode == SOLVER_MODE_ASSEMBLY ) THEN
5808
         CALL CoupledSolver( Model, Solver, DTScal * dt, TimeDerivativeActive )
5809
       ELSE IF( Solver % SolverMode == SOLVER_MODE_BLOCK ) THEN
5810
         CALL BlockSolver( Model, Solver, DTScal * dt, TimeDerivativeActive )
5811
       ELSE 
5812
         CALL SingleSolver( Model, Solver, DTScal * dt, TimeDerivativeActive )
5813
       END IF
5814

5815
       IF( GotLoops ) THEN
5816
         IF( ListGetLogical( Params,'Save Scanning Modes',Found ) ) THEN
5817
           n = SIZE( Solver % Variable % Values )
5818
           IF ( .NOT. ASSOCIATED( Solver % Variable % EigenValues ) ) THEN
5819
             CALL Info('MainUtils','Creating modes over scanned fields',Level=8)
5820
             ALLOCATE( Solver % Variable % EigenValues(ScanningLoops) )
5821
             ALLOCATE( Solver % Variable % EigenVectors(ScanningLoops,n) )             
5822

5823
             IF( Solver % Variable % Dofs > 1 ) THEN
5824
               DO k=1,Solver % Variable % DOFs
5825
                 str = ComponentName( Solver % Variable % Name, k )
5826
                 Var => VariableGet( Solver % Mesh % Variables, str, .TRUE. )
5827
                 IF ( ASSOCIATED( Var ) ) THEN
5828
                   Var % EigenValues => Solver % Variable % EigenValues
5829
                   Var % EigenVectors =>  & 
5830
                       Solver % Variable % EigenVectors(:,k::Solver % Variable % DOFs )
5831
                 END IF
5832
               END DO
5833
             END IF
5834
           END IF
5835
           Solver % Variable % EigenValues(scan) = 1.0_dp * scan
5836
           Solver % Variable % EigenVectors(scan,:) = Solver % Variable % Values
5837
         END IF
5838
       END IF
5839
       
5840
       Solver % TimesVisited = Solver % TimesVisited + 1
5841
     END DO
5842
       
5843
     IF(.NOT. ListGetLogical( Params,'Auxiliary Solver',Found)) THEN
5844
       IF(NamespaceFound) CALL ListPopNamespace()
5845
     END IF
5846
     Solver % dt = dt
5847

5848
     IF( GotCoordTransform ) THEN
5849
       CALL BackCoordinateTransformation( Solver % Mesh )
5850
     END IF
5851

5852
!------------------------------------------------------------------------------
5853
! After solution register the timing, if requested
5854
!------------------------------------------------------------------------------
5855
     IF( Timing ) THEN
5856
       st  = CPUTime() - t0
5857
       rst = RealTime() - rt0
5858

5859
       str = ListGetString( Params,'Equation',Found)
5860
       WRITE(Message,'(a,f8.2,f8.2,a)') 'Solver time (CPU,REAL) for '&
5861
           //TRIM(str)//': ',st,rst,' (s)'
5862
       CALL Info('SolverActivate',Message,Level=4)    
5863
      
5864
       IF( ListGetLogical(Params,'Solver Timing',Found)) THEN
5865
         CALL ListAddConstReal(CurrentModel % Simulation,'res: solver cpu time '&
5866
             //TRIM(str),st)
5867
         CALL ListAddConstReal(CurrentModel % Simulation,'res: solver real time '&
5868
             //TRIM(str),rst)
5869
       END IF
5870
         
5871
       IF( ListGetLogical(Params,'Solver Timing Cumulative',Found)) THEN
5872
          ct = ListGetConstReal(CurrentModel % Simulation,'res: cum solver cpu time '&
5873
                //TRIM(str),Found)
5874
          st = st + ct
5875
          ct = ListGetConstReal(CurrentModel % Simulation,'res: cum solver real time '&
5876
                //TRIM(str),Found)
5877
          rst = rst + ct
5878
          CALL ListAddConstReal(CurrentModel % Simulation,'res: cum solver cpu time '&
5879
              //TRIM(str),st)
5880
          CALL ListAddConstReal(CurrentModel % Simulation,'res: cum solver real time '&
5881
              //TRIM(str),rst)
5882
        END IF 
5883
      END IF
5884

5885
      OutputPE = sOutputPE
5886
!------------------------------------------------------------------------------
5887
   END SUBROUTINE SolverActivate
5888
!------------------------------------------------------------------------------
5889

5890
!------------------------------------------------------------------------------
5891
!  The predictor-corrector scheme to change dt
5892
!------------------------------------------------------------------------------
5893
!------------------------------------------------------------------------------
5894
    SUBROUTINE PredictorCorrectorControl( Model, dt, RealTimestep )
5895
!------------------------------------------------------------------------------
5896
      
5897
      TYPE(Model_t), INTENT(IN) :: Model
5898
      REAL(KIND=dp), INTENT(INOUT) :: dt
5899
      INTEGER, OPTIONAL :: RealTimestep
5900
!------------------------------------------------------------------------------
5901
      TYPE(Solver_t), POINTER :: Solver
5902
      TYPE(ValueList_t), POINTER :: SolverParams
5903
      INTEGER :: PredCorrOrder, i, predcorrIndex = 0
5904
      REAL(KIND=dp) :: epsilon, beta1, beta2
5905
      LOGICAL :: Found, OutputFlag = .FALSE.
5906

5907
      REAL(KIND=dp) :: timeError, timeErrorMax, timeError2Norm, eta
5908

5909
      REAL(KIND=dp), SAVE:: dtOld, etaOld, zeta
5910

5911

5912
      DO i=1,Model % NumberOFSolvers
5913
        Solver => Model % Solvers(i)
5914
        !> Find the Solver for adaptive predictor, there should be only one solver as predictor
5915
        IF (Solver % SolverExecWhen == SOLVER_EXEC_PREDCORR) THEN
5916
          predcorrIndex = i
5917
        END IF 
5918
      END DO
5919

5920
      IF ( predcorrIndex == 0) THEN 
5921
        CALL Fatal('Predictor-Corrector Control','Predictor-Corrector Solver is not found!')
5922
      ELSE
5923
        ! Do Predictor-Corrector
5924
        Solver => Model % Solvers(predcorrIndex)
5925
        SolverParams => ListGetSolverParams(Solver)
5926

5927
        IF (RealTimestep == 1) THEN 
5928
        ! Do nothing on the first step
5929
          dtOld = dt
5930
          dt = 0.0_dp
5931
          zeta = 1.0_dp
5932
        ELSE IF (RealTimestep == 2) THEN
5933
        ! Use the initial time step, force to use first order time schemes 
5934
          dt = dtOld
5935
          zeta = 1.0_dp
5936
        ELSE IF (RealTimestep > 2) THEN
5937
        ! Use local error estimate and PI control 
5938

5939
          ! Read in the settings
5940
          CALL ReadPredCorrParams( Model, SolverParams, PredCorrOrder, epsilon, beta1, beta2 )
5941

5942
          ! Compute the error  |H-\tilde{H}|_inf
5943

5944
          timeErrorMax  =  MAXVAL( ABS(Solver % Variable % Values(:) - Solver % Variable % PrevValues(:,1)))
5945
          timeErrorMax = ParallelReduction(timeErrorMax,2)
5946

5947
          timeError = timeErrorMax
5948

5949
          ! 1st order error estimate for the first control step
5950
          IF (RealTimestep == 3 ) THEN 
5951
            PredCorrOrder = 1
5952
          END IF
5953

5954
          ! Estimate local truncation error use old zeta
5955
          CALL PredCorrErrorEstimate( eta, dtOld, PredCorrOrder, timeError, zeta )
5956
          IF (RealTimestep == 3 ) THEN 
5957
            etaOld =eta
5958
          END IF
5959
          ! PI controller
5960
          CALL TimeStepController( dt, dtOld, eta, etaOld, epsilon, beta1, beta2 )
5961
          ! Compute new zeta and for predictor time scheme
5962
          zeta = dt / dtOld
5963
          CALL ListAddConstReal(Solver % Values, 'Adams Zeta', zeta)
5964
          ! Save old eta
5965
          etaOld = eta
5966

5967

5968
          !> Save the time errors!     
5969
          OutputFlag = ListGetLogical(SolverParams, 'Predictor-Corrector Save Error', Found)   
5970
          IF ( OutputFlag ) THEN                 
5971
            OPEN (unit=135, file="ErrorPredictorCorrector.dat", POSITION='APPEND')
5972
            WRITE(135, *) dtOld, eta, timeError                                                
5973
            CLOSE(135)
5974
          END IF
5975

5976
          !> Output
5977
          WRITE (Message,*) "---------------- Predictor-Corrector Control ----------------------"
5978
          CALL Info('Predictor-Corrector', Message, Level=4)
5979
          WRITE (Message,*) "current dt=", dtOld, "next dt=",  dt
5980
          CALL Info('Predictor-Corrector', Message, Level=4)
5981
          WRITE (Message,*) "zeta=", zeta, "eta=",  eta, "terr=", timeError
5982
          CALL Info('Predictor-Corrector', Message, Level=6)
5983
          dtOld = dt
5984
        END IF 
5985
      END IF
5986
      
5987
!------------------------------------------------------------------------------
5988
    END SUBROUTINE PredictorCorrectorControl
5989
!------------------------------------------------------------------------------
5990

5991

5992
!------------------------------------------------------------------------------
5993
!  The adaptive controller for the predictor-corrector scheme
5994
!------------------------------------------------------------------------------
5995
!------------------------------------------------------------------------------
5996
    SUBROUTINE TimeStepController( dt, dtOld, eta, etaOld, epsilon, beta1, beta2 )
5997
!------------------------------------------------------------------------------
5998

5999
      REAL(KIND=dp), INTENT(INOUT):: dt  
6000
      REAL(KIND=dp), INTENT(IN):: dtOld, eta, etaOld, epsilon, beta1, beta2     
6001

6002
      REAL(KIND=dp) :: gfactor
6003

6004
      IF ((eta /= 0.0_dp) .AND. (etaOld /= 0.0_dp)) THEN 
6005
        gfactor = ((epsilon/eta)**beta1) * ((epsilon/etaOld)**beta2)
6006
        CALL TimeStepLimiter(dtOld, dt, gfactor)
6007
      ELSE 
6008
        dt = dtOld
6009
      END IF
6010

6011
!------------------------------------------------------------------------------
6012
    END  SUBROUTINE TimeStepController
6013
!------------------------------------------------------------------------------
6014

6015
!------------------------------------------------------------------------------
6016
!  Time step limiter for Adaptive TimeStepping
6017
!------------------------------------------------------------------------------
6018
!------------------------------------------------------------------------------
6019
    SUBROUTINE TimeStepLimiter(dtOld, dt, gfactor, k)
6020
!------------------------------------------------------------------------------
6021
      REAL(KIND=dp), INTENT(IN) :: dtOld, gfactor
6022
      REAL(KIND=dp), INTENT(OUT) :: dt
6023
      REAL(KIND=dp) :: xfactor
6024
      INTEGER, OPTIONAL :: k 
6025

6026
      IF( PRESENT(k) ) THEN
6027
        xfactor = 1.0 + k * ATAN((gfactor-1)/k)
6028
      ELSE
6029
        xfactor = 1.0 + 2.0 * ATAN((gfactor-1)*0.5)
6030
      END IF 
6031

6032
      dt = dtOld * xfactor
6033
!------------------------------------------------------------------------------
6034
    END SUBROUTINE TimeStepLimiter
6035
!------------------------------------------------------------------------------
6036

6037

6038
!------------------------------------------------------------------------------
6039
!  Local truncation error for AB-AM 1st and 2nd order methods
6040
!------------------------------------------------------------------------------
6041
!------------------------------------------------------------------------------
6042
    SUBROUTINE  PredCorrErrorEstimate(eta, dt, PredCorrOrder, timeError, zeta)
6043
!------------------------------------------------------------------------------
6044
      REAL(KIND=dp), INTENT(IN) :: dt, timeError, zeta
6045
      INTEGER, INTENT(IN) :: PredCorrOrder
6046
      REAL(KIND=dp), INTENT(OUT) :: eta
6047
    
6048
      IF (dt > 0.0_dp) THEN
6049
        IF (PredCorrOrder == 2) THEN
6050
          eta = timeError * zeta / dt / (zeta + 1.0_dp) / 3.0_dp
6051
        ELSE
6052
          eta = timeError / dt / 2.0_dp
6053
        END IF
6054
      ELSE 
6055
        CALL WARN('Predictor-Corrector','Time Step is 0 in Local error estimate!')        
6056
        eta =0.0_dp
6057
      END IF
6058
     
6059
!------------------------------------------------------------------------------
6060
    END SUBROUTINE PredCorrErrorEstimate
6061
!------------------------------------------------------------------------------
6062

6063

6064

6065
!------------------------------------------------------------------------------
6066
!  Set Default / Read in the settings for Predictor-Corrector scheme
6067
!------------------------------------------------------------------------------
6068
!------------------------------------------------------------------------------
6069
    SUBROUTINE ReadPredCorrParams(Model, SolverParams, outOrder, outEps, outB1, outB2)
6070

6071
      TYPE(Model_t), INTENT(IN) :: Model
6072
      TYPE(ValueList_t), POINTER, INTENT(INOUT):: SolverParams
6073
      INTEGER, OPTIONAL, INTENT(OUT) :: outOrder
6074
      REAL(KIND=dp), OPTIONAL, INTENT(OUT) :: outEps, outB1, outB2
6075

6076
      INTEGER :: PredCorrOrder
6077
      REAL(KIND=dp) :: epsilon, beta1, beta2
6078

6079
      LOGICAL :: Found
6080

6081

6082
      PredCorrOrder = ListGetInteger( SolverParams,  &
6083
                      'Predictor-Corrector Scheme Order', Found)
6084
      IF ( .NOT. Found ) THEN
6085
        PredCorrOrder = ListGetInteger( Model % Simulation,  &
6086
                        'Predictor-Corrector Scheme Order', Found)
6087
        IF ( .NOT. Found ) THEN
6088
          PredCorrOrder = 2
6089
        END IF
6090
        CALL ListAddInteger( SolverParams,  &
6091
                        'Predictor-Corrector Scheme Order', PredCorrOrder )
6092
      END IF
6093

6094
      epsilon = ListGetCReal( SolverParams, &
6095
                        'Predictor-Corrector Control Tolerance', Found )
6096
      IF ( .NOT. Found ) THEN
6097
        epsilon = ListGetCReal( Model % Simulation,  &
6098
                        'Predictor-Corrector Control Tolerance', Found )
6099
        IF ( .NOT. Found ) THEN
6100
          epsilon = 1.0e-6
6101
        END IF
6102
        CALL ListAddConstReal( SolverParams, &
6103
                        'Predictor-Corrector Control Tolerance', epsilon )
6104
      END IF
6105

6106
      beta1 = ListGetCReal( SolverParams, &
6107
                        'Predictor-Corrector Control Beta 1', Found )
6108
      IF ( .NOT. Found ) THEN
6109
        beta1 = ListGetCReal( Model % Simulation,  &
6110
                        'Predictor-Corrector Control Beta 1', Found )
6111
        IF ( .NOT. Found ) THEN
6112
          beta1 = 0.6_dp / ( PredCorrOrder + 1.0_dp )
6113
        END IF
6114
        CALL ListAddConstReal( SolverParams, &
6115
                        'Predictor-Corrector Control Beta 1', beta1 )
6116
      END IF
6117

6118
      beta2 = ListGetCReal( SolverParams, &
6119
                        'Predictor-Corrector Control Beta 2', Found )
6120
      IF ( .NOT. Found ) THEN
6121
        beta2 = ListGetCReal( Model % Simulation,  &
6122
                        'Predictor-Corrector Control Beta 2', Found )
6123
        IF ( .NOT. Found ) THEN
6124
          beta2 = -0.2_dp / ( PredCorrOrder + 1.0_dp )
6125
        END IF
6126
        CALL ListAddConstReal( SolverParams, &
6127
                        'Predictor-Corrector Control Beta 2', beta2 )
6128
      END IF
6129

6130
      ! Output if required
6131
      IF ( PRESENT( outOrder ) ) outOrder = PredCorrOrder
6132
      IF ( PRESENT( outEps ) ) outEps = epsilon
6133
      IF ( PRESENT( outB1 ) ) outB1 = beta1
6134
      IF ( PRESENT( outB2 ) ) outB2 = beta2
6135
!------------------------------------------------------------------------------
6136
    END  SUBROUTINE ReadPredCorrParams
6137
!------------------------------------------------------------------------------
6138

6139
END MODULE MainUtils
6140

6141
!> \} ElmerLib
6142

6143
Product

Resources

Company
