← Back to branch summary

~ubuntu-branches/ubuntu/maverick/elmerfem/maverick

« back to all changes in this revision

Viewing changes to fem/src/FluxSolver.src

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2009-04-08 16:58:46 UTC
  • Revision ID: james.westby@ubuntu.com-20090408165846-qydtxka1cvrjw2ix
Tags: upstream-5.5.0.svn.4096.dfsg
Import upstream version 5.5.0.svn.4096.dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
fem/src/FluxSolver.src
 
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 program is free software; you can redistribute it and/or
 
8
! *  modify it under the terms of the GNU General Public License
 
9
! *  as published by the Free Software Foundation; either version 2
 
10
! *  of the License, or (at your option) any later version.
 
11
! * 
 
12
! *  This program 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
 
15
! *  GNU General Public License for more details.
 
16
! *
 
17
! *  You should have received a copy of the GNU General Public License
 
18
! *  along with this program (in file fem/GPL-2); if not, write to the 
 
19
! *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 
 
20
! *  Boston, MA 02110-1301, USA.
 
21
! *
 
22
! *****************************************************************************/
 
23
!
 
24
!/******************************************************************************
 
25
! *
 
26
! *  Module for computing fluxes of Poisson type of equations
 
27
! *
 
28
! ******************************************************************************
 
29
! *
 
30
! *  Authors: Peter R�back, Juha Ruokolainen
 
31
! *  Email:   Peter.Raback@csc.fi
 
32
! *  Web:     http://www.csc.fi/elmer
 
33
! *  Address: CSC - IT Center for Science Ltd.
 
34
! *           Keilaranta 14
 
35
! *           02101 Espoo, Finland 
 
36
! *
 
37
! *  Original Date: 20.06.2007
 
38
! *
 
39
! *****************************************************************************/
 
40
 
 
41
!------------------------------------------------------------------------------
 
42
  SUBROUTINE FluxSolver_Init( Model,Solver,dt,Transient )
 
43
!------------------------------------------------------------------------------
 
44
    USE DefUtils
 
45
 
 
46
    TYPE(Model_t)  :: Model
 
47
    TYPE(Solver_t) :: Solver
 
48
    REAL(KIND=dp) :: DT
 
49
    LOGICAL :: Transient
 
50
!------------------------------------------------------------------------------
 
51
    TYPE(ValueList_t), POINTER :: SolverParams
 
52
    INTEGER :: dim
 
53
!------------------------------------------------------------------------------
 
54
    SolverParams => GetSolverParams()
 
55
    dim = CoordinateSystemDimension()
 
56
 
 
57
    IF ( .NOT. ListCheckPresent( SolverParams,'Variable') ) THEN
 
58
      CALL ListAddInteger( SolverParams, 'Variable DOFs', 1 )
 
59
      IF(ListCheckPresent(SolverParams,'Flux Component')) THEN
 
60
        CALL ListAddString( SolverParams,'Variable','Flux' )
 
61
      ELSE
 
62
        CALL ListAddString( SolverParams, 'Variable','-nooutput flux_temp' )
 
63
        CALL ListAddString(SolverParams, 'Flux Result Variable','F')
 
64
        IF(dim == 2) THEN
 
65
          CALL ListAddString(SolverParams, 'Exported Variable 1','F[Flux:2]')
 
66
        ELSE IF(dim == 3) THEN
 
67
          CALL ListAddString(SolverParams, 'Exported Variable 1','F[Flux:3]')
 
68
        ELSE
 
69
          CALL Fatal('VortictySolver_init','Flux computation makes sense only in 2D and 3D')
 
70
        END IF
 
71
      END IF
 
72
    END IF
 
73
    CALL ListAddInteger( SolverParams, 'Time derivative order', 0 )
 
74
 
 
75
    ! Add linear system defaults: cg+diagonal
 
76
    IF(.NOT. ListCheckPresent(SolverParams,'Linear System Solver')) &
 
77
      CALL ListAddString(SolverParams,'Linear System Solver','Iterative')
 
78
    IF(.NOT. ListCheckPresent(SolverParams,'Linear System Iterative Method')) &
 
79
      CALL ListAddString(SolverParams,'Linear System Iterative Method','cg')
 
80
    IF(.NOT. ListCheckPresent(SolverParams,'Linear System Preconditioning')) &
 
81
      CALL ListAddString(SolverParams,'Linear System Preconditioning','diagonal')
 
82
    IF(.NOT. ListCheckPresent(SolverParams,'Linear System Max Iterations')) &
 
83
      CALL ListAddInteger(SolverParams,'Linear System Max Iterations',500)
 
84
    IF(.NOT. ListCheckPresent(SolverParams,'Linear System Residual Output')) &
 
85
      CALL ListAddInteger(SolverParams,'Linear System Residual Output',10)
 
86
    IF(.NOT. ListCheckPresent(SolverParams,'Linear System Convergence Tolerance')) &
 
87
      CALL ListAddConstReal(SolverParams,'Linear System Convergence Tolerance',1.0e-10_dp)
 
88
 
 
89
!------------------------------------------------------------------------------
 
90
  END SUBROUTINE FluxSolver_Init
 
91
!------------------------------------------------------------------------------
 
92
 
 
93
 
 
94
 
 
95
 
 
96
!------------------------------------------------------------------------------
 
97
SUBROUTINE FluxSolver( Model,Solver,dt,Transient )
 
98
!------------------------------------------------------------------------------
 
99
 
 
100
  USE CoordinateSystems
 
101
  USE DefUtils
 
102
 
 
103
  IMPLICIT NONE
 
104
!------------------------------------------------------------------------------
 
105
!******************************************************************************
 
106
!
 
107
!  Solve stress equations for one timestep
 
108
!
 
109
!  ARGUMENTS:
 
110
!
 
111
!  TYPE(Model_t) :: Model,  
 
112
!     INPUT: All model information (mesh,materials,BCs,etc...)
 
113
!
 
114
!  TYPE(Solver_t) :: Solver
 
115
!     INPUT: Linear equation solver options
 
116
!
 
117
!  REAL(KIND=dp) :: dt,
 
118
!     INPUT: Timestep size for time dependent simulations (NOTE: Not used
 
119
!            currently)
 
120
!
 
121
!******************************************************************************
 
122
 
 
123
  TYPE(Model_t)  :: Model
 
124
  TYPE(Solver_t), TARGET :: Solver
 
125
  LOGICAL ::  Transient
 
126
  REAL(KIND=dp) :: dt
 
127
!------------------------------------------------------------------------------
 
128
!    Local variables
 
129
!------------------------------------------------------------------------------
 
130
  TYPE(ValueList_t),POINTER :: SolverParams
 
131
  CHARACTER(LEN=MAX_NAME_LEN) :: VarName, CondName
 
132
  INTEGER :: i,j,dim,DOFs,ActiveDir
 
133
  LOGICAL :: ConstantBulkMatrix, ConstantBulkMatrixInUse, CSymmetry
 
134
  LOGICAL :: GotIt, Visited = .FALSE., DirMask(3)
 
135
  REAL(KIND=dp) :: Unorm, Totnorm
 
136
  REAL(KIND=dp), POINTER :: ForceVector(:,:), SaveRHS(:)
 
137
  REAL(KIND=dp) :: at0,at1,at2,CPUTime,RealTime
 
138
  TYPE(Variable_t), POINTER :: FluxSol
 
139
  
 
140
  SAVE Visited
 
141
 
 
142
 
 
143
  CALL Info( 'FluxSolver', '-------------------------------------',Level=4 )
 
144
  CALL Info( 'FluxSolver','Computing the flux',Level=4 )
 
145
  CALL Info( 'FluxSolver', '-------------------------------------',Level=4 )
 
146
 
 
147
  dim = CoordinateSystemDimension()
 
148
!------------------------------------------------------------------------------
 
149
!    Get variables needed for solution
 
150
!------------------------------------------------------------------------------
 
151
  IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN
 
152
  IF ( COUNT( Solver % Variable % Perm > 0 ) <= 0 ) RETURN
 
153
  
 
154
  SolverParams => GetSolverParams()
 
155
 
 
156
!-------------------------------------------------------------------------------
 
157
! If only one component is used use the scalar equation, otherwise use an
 
158
! auxiliary variable to store all the dimensions
 
159
!-------------------------------------------------------------------------------
 
160
  ActiveDir = GetInteger( SolverParams, 'Flux Component', GotIt) 
 
161
  IF(GotIt) THEN
 
162
    Dofs = 1
 
163
  ELSE
 
164
    VarName = GetString(SolverParams,'Flux Result Variable',GotIt )
 
165
    IF(.NOT. gotIt) VarName = 'Flux'
 
166
    FluxSol => VariableGet( Solver % Mesh % Variables, VarName )
 
167
    IF(ASSOCIATED(FluxSol)) THEN
 
168
      Dofs = FluxSol % DOFs
 
169
      IF(Dofs /= DIM) THEN
 
170
        CALL Fatal('FluxSolver','The flux should have DOFs equal to DIM')
 
171
      END IF
 
172
    ELSE
 
173
      CALL Fatal('FluxSolver','Flux Result Variable is missing: '//TRIM(VarName))      
 
174
    END IF
 
175
    DirMask = .TRUE.
 
176
  END IF
 
177
  
 
178
  CSymmetry = CurrentCoordinateSystem() == AxisSymmetric .OR. &
 
179
      CurrentCoordinateSystem() == CylindricSymmetric
 
180
  
 
181
  VarName = GetString(SolverParams,'Flux Variable',GotIt )
 
182
  IF(.NOT. gotIt) VarName = TRIM('Temperature')
 
183
 
 
184
  CondName = ListGetString(SolverParams,'Flux Coefficient',GotIt )
 
185
  IF(.NOT. gotIt) CondName = TRIM('Heat Conductivity')
 
186
  
 
187
  at0 = RealTime()
 
188
  
 
189
  ConstantBulkMatrix = GetLogical( SolverParams, 'Constant Bulk Matrix', GotIt )
 
190
  ConstantBulkMatrixInUse = ConstantBulkMatrix .AND. &
 
191
      ASSOCIATED(Solver % Matrix % BulkValues)
 
192
  
 
193
  IF ( ConstantBulkMatrixInUse ) THEN
 
194
    Solver % Matrix % Values = Solver % Matrix % BulkValues        
 
195
    Solver % Matrix % rhs = 0.0_dp
 
196
  ELSE
 
197
    CALL DefaultInitialize()
 
198
  END IF
 
199
  
 
200
  IF(Dofs > 1) THEN
 
201
    ALLOCATE(ForceVector(Dofs-1,SIZE(Solver % Matrix % RHS)))  
 
202
    ForceVector = 0.0_dp
 
203
    SaveRHS => Solver % Matrix % RHS
 
204
  END IF
 
205
  
 
206
  CALL BulkAssembly()
 
207
  CALL DefaultFinishAssembly()
 
208
  
 
209
  at1 = RealTime()
 
210
  WRITE(Message,* ) 'Assembly Time: ',at1-at0
 
211
  CALL Info( 'FluxSolver', Message, Level=5 )
 
212
!        
 
213
!------------------------------------------------------------------------------     
 
214
 
 
215
  IF(Dofs > 1) THEN
 
216
    TotNorm = 0.0_dp
 
217
    DO i=1,Dofs
 
218
      IF(i==1) THEN
 
219
        Solver % Matrix % RHS => SaveRHS
 
220
      ELSE
 
221
        Solver % Matrix % RHS => ForceVector(i-1,:)
 
222
      END IF
 
223
      UNorm = DefaultSolve()
 
224
      TotNorm = TotNorm + Unorm ** 2
 
225
      DO j=1,Solver % Matrix % NumberOfRows
 
226
        FluxSol % Values(DOFs*(j-1)+i) = Solver % Variable % Values(j)
 
227
      END DO
 
228
    END DO
 
229
    DEALLOCATE( ForceVector )  
 
230
    Solver % Matrix % RHS => SaveRHS
 
231
    TotNorm = SQRT(TotNorm)
 
232
    Solver % Variable % Norm = Totnorm
 
233
  ELSE     
 
234
    TotNorm = DefaultSolve()
 
235
  END IF
 
236
!------------------------------------------------------------------------------     
 
237
 
 
238
  at2 = RealTime()
 
239
  WRITE(Message,* ) 'Solution Time: ',at2-at1
 
240
  CALL Info( 'FluxSolver', Message, Level=5 )
 
241
  
 
242
  WRITE( Message, * ) 'Result Norm: ',TotNorm
 
243
  CALL Info( 'FluxSolver', Message, Level=4 )
 
244
  
 
245
CONTAINS
 
246
 
 
247
 
 
248
!------------------------------------------------------------------------------
 
249
  SUBROUTINE BulkAssembly()
 
250
!------------------------------------------------------------------------------
 
251
       
 
252
    INTEGER :: elem,t,i,j,p,q,n,nd, Rank
 
253
    REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), FORCE(:,:)
 
254
    TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
 
255
    TYPE(Nodes_t) :: Nodes
 
256
    TYPE(Element_t), POINTER :: Element
 
257
    REAL(KIND=dp) :: weight,grad(3),C(3,3),coeff,detJ
 
258
    REAL(KIND=dp), ALLOCATABLE :: Basis(:), dBasisdx(:,:)
 
259
    REAL(KIND=dp), ALLOCATABLE :: LocalPotential(:)
 
260
    LOGICAL :: Found
 
261
    TYPE(ValueList_t), POINTER :: Material
 
262
    REAL(KIND=dp), POINTER :: Conductivity(:,:,:)=>NULL()
 
263
    
 
264
    SAVE Conductivity, Nodes
 
265
    
 
266
    n = MAX( Solver % Mesh % MaxElementDOFs, Solver % Mesh % MaxElementNodes )
 
267
    ALLOCATE( STIFF(n,n), FORCE(dim,n) )
 
268
    ALLOCATE( LocalPotential(n), Basis(n), dBasisdx(n,3) )
 
269
 
 
270
    DO elem = 1,Solver % NumberOFActiveElements
 
271
         
 
272
      ! Element information
 
273
      ! ---------------------
 
274
      Element => GetActiveElement(elem)
 
275
      CALL GetElementNodes( Nodes )
 
276
      nd = GetElementNOFDOFs()
 
277
      n  = GetElementNOFNodes()
 
278
      
 
279
      CALL GetRealArray( GetMaterial(), Conductivity, CondName, Found )
 
280
      Rank = 0
 
281
      IF ( Found ) Rank = GetTensorRank(Conductivity)
 
282
      CALL GetScalarLocalSolution( LocalPotential, VarName )
 
283
      
 
284
      ! Integrate local stresses:
 
285
      ! -------------------------
 
286
      IntegStuff = GaussPoints( Element )
 
287
      STIFF  = 0.0_dp
 
288
      FORCE  = 0.0_dp
 
289
 
 
290
      C = 0.0_dp
 
291
      DO i=1,dim
 
292
        C(i,i) = 1.0_dp
 
293
      END DO
 
294
      
 
295
      DO t=1,IntegStuff % n
 
296
        Found = ElementInfo( Element, Nodes, IntegStuff % u(t), &
 
297
            IntegStuff % v(t), IntegStuff % w(t), detJ, Basis, dBasisdx )
 
298
        
 
299
        Weight = IntegStuff % s(t) * detJ
 
300
        IF ( CSymmetry ) Weight = Weight * SUM( Basis(1:n) * Nodes % x(1:n) )
 
301
        
 
302
        IF ( .NOT. ConstantBulkMatrixInUse ) THEN
 
303
          DO p=1,nd
 
304
            DO q=1,nd
 
305
              STIFF(p,q) = STIFF(p,q) + Weight * Basis(q) * Basis(p)
 
306
            END DO
 
307
          END DO
 
308
        END IF
 
309
        
 
310
        Grad(1:dim) = MATMUL( LocalPotential(1:nd), dBasisdx(1:nd,1:dim) )
 
311
        
 
312
        SELECT CASE(Rank)
 
313
        CASE(0)
 
314
        CASE(1)
 
315
          DO i=1,dim
 
316
            C(i,i) = SUM( Basis(1:n) * Conductivity(1,1,1:n) )
 
317
          END DO
 
318
        CASE(2)
 
319
          DO i=1,dim
 
320
            C(i,i) = SUM( Basis(1:n) * Conductivity(i,1,1:n) )
 
321
          END DO
 
322
        CASE DEFAULT
 
323
          DO i=1,dim
 
324
            DO j=1,dim
 
325
              C(i,j) = SUM( Basis(1:n) * Conductivity(i,j,1:n) )
 
326
            END DO
 
327
          END DO
 
328
        END SELECT
 
329
        
 
330
        IF(Dofs == 1) THEN
 
331
          Coeff = Weight * SUM( C(ActiveDir,1:dim) * Grad(1:dim) )
 
332
          FORCE(1,1:nd) = FORCE(1,1:nd) - Coeff*Basis(1:nd)
 
333
        ELSE            
 
334
          DO i=1,dim    
 
335
            Coeff = Weight * SUM( C(i,1:dim) * Grad(1:dim) )
 
336
            FORCE(i,1:nd) = FORCE(i,1:nd) - Coeff*Basis(1:nd)
 
337
          END DO
 
338
        END IF
 
339
      END DO
 
340
 
 
341
!------------------------------------------------------------------------------
 
342
!      Update global matrices from local matrices 
 
343
!------------------------------------------------------------------------------
 
344
      IF ( .NOT. ConstantBulkMatrixInUse ) THEN
 
345
        IF(Dofs > 1) Solver % Matrix % RHS => SaveRHS
 
346
        CALL DefaultUpdateEquations( STIFF, FORCE(1,1:nd), &
 
347
            BulkUpdate=ConstantBulkMatrix )
 
348
      ELSE
 
349
        CALL DefaultUpdateForce( FORCE(1,1:nd) )       
 
350
      END IF
 
351
 
 
352
      DO i=2,Dofs
 
353
        Solver % Matrix % RHS => ForceVector(i-1,:)
 
354
        CALL DefaultUpdateForce( FORCE(i,1:nd) )
 
355
      END DO
 
356
    END DO
 
357
    
 
358
    DEALLOCATE( LocalPotential, STIFF, FORCE, Basis, dBasisdx )
 
359
!------------------------------------------------------------------------------
 
360
  END SUBROUTINE BulkAssembly
 
361
!------------------------------------------------------------------------------
 
362
 
 
363
 
 
364
!------------------------------------------------------------------------------
 
365
  FUNCTION GetTensorRank( Tensor ) RESULT ( Rank )
 
366
!------------------------------------------------------------------------------
 
367
    REAL(KIND=dp), POINTER :: Tensor(:,:,:)
 
368
    INTEGER :: Rank
 
369
    
 
370
    IF ( SIZE(Tensor,1) == 1 ) THEN
 
371
      Rank = 1
 
372
    ELSE IF ( SIZE(Tensor,2) == 1 ) THEN
 
373
      Rank = 2
 
374
    ELSE
 
375
      Rank = 3
 
376
    END IF
 
377
!-----------------------------------------------------------------------------
 
378
  END FUNCTION GetTensorRank
 
379
!------------------------------------------------------------------------------
 
380
 
 
381
 
 
382
!------------------------------------------------------------------------------
 
383
END SUBROUTINE FluxSolver
 
384
!------------------------------------------------------------------------------
 
385