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 ODY solver with global dofs.
27! *
28! *  Authors: Peter Råback
29! *  Email:   Peter.Raback@csc.fi
30! *  Web:     http://www.csc.fi/elmer
31! *  Address: CSC - IT Center for Science Ltd.
32! *           Keilaranta 14
33! *           02101 Espoo, Finland
34! *
35! *  Original Date: 4.5.2015
36! *
37! *****************************************************************************/
38
39!> \ingroup Solvers
40!> \{
41
42
43!------------------------------------------------------------------------------
44!> Initialization for the primary solver
45!------------------------------------------------------------------------------
46SUBROUTINE OdeSolver_init( Model,Solver,dt,TransientSimulation )
47!------------------------------------------------------------------------------
48  USE DefUtils
49  IMPLICIT NONE
50!------------------------------------------------------------------------------
51  TYPE(Solver_t) :: Solver       !< Linear & nonlinear equation solver options
52  TYPE(Model_t) :: Model         !< All model information (mesh, materials, BCs, etc...)
53  REAL(KIND=dp) :: dt            !< Timestep size for time dependent simulations
54  LOGICAL :: TransientSimulation !< Steady state or transient simulation
55!------------------------------------------------------------------------------
56  TYPE(ValueList_t), POINTER :: Params
57
58  Params => GetSolverParams()
59
60  IF(.NOT. ListCheckPresent( Params,'Time Derivative Order') ) THEN
61    CALL ListAddInteger( Params, 'Time derivative order', 2 )
62  END IF
63  CALL ListAddLogical( Params, 'Use Global Mass Matrix', .TRUE. )
64  CALL ListAddLogical( Params, 'Ode Matrix', .TRUE. )
65  CALL ListAddLogical( Params, 'Variable Global', .TRUE. )
66
67!------------------------------------------------------------------------------
68END SUBROUTINE OdeSolver_init
69!------------------------------------------------------------------------------
70
71
72!------------------------------------------------------------------------------
73!> This is a simple solver for ordinary differential equations (ODE)
74!> utilizing many of the same subroutines as the default PDE solution.
75!> The idea is to easily allow the use of same time integration schemes
76!> etc. There is plenty of room for improvement and polishing.
77!------------------------------------------------------------------------------
78SUBROUTINE OdeSolver( Model,Solver,dt,TransientSimulation )
79!------------------------------------------------------------------------------
80  USE DefUtils
81  IMPLICIT NONE
82!------------------------------------------------------------------------------
83  TYPE(Solver_t) :: Solver       !< Linear & nonlinear equation solver options
84  TYPE(Model_t) :: Model         !< All model information (mesh, materials, BCs, etc...)
85  REAL(KIND=dp) :: dt            !< Timestep size for time dependent simulations
86  LOGICAL :: TransientSimulation !< Steady state or transient simulation
87!------------------------------------------------------------------------------
88! Local variables
89!------------------------------------------------------------------------------
90  LOGICAL :: Found
91  REAL(KIND=dp) :: Norm
92  INTEGER :: iter, MaxIter,TimeOrder
93  CHARACTER(LEN=MAX_NAME_LEN) :: str
94  TYPE(ValueList_t), POINTER :: Params
95!------------------------------------------------------------------------------
96
97  CALL Info('OdeSolver','Solving ordinary differential equation',Level=4)
98
99
100  Params => GetSolverParams()
101
102  MaxIter = GetInteger( Params,'Nonlinear System Max Iterations',Found )
103  IF(.NOT. Found) MaxIter = 1
104
105  TimeOrder = GetInteger( Params,'Time Derivative Order' )
106
107  DO iter = 1, MaxIter
108    IF( MaxIter > 1 ) THEN
109      CALL Info('OdeSolver','Nonlinear iteration: '//TRIM(I2S(iter)),Level=5)
110    END IF
111
112    CALL DefaultInitialize()
113
114    CALL OdeMatrix()
115
116    CALL DefaultFinishBulkAssembly()
117    CALL DefaultFinishBoundaryAssembly()
118
119    ! This sets the time-integration and is therefore imperative
120    CALL DefaultFinishAssembly()
121
122    ! This is probably not needed...
123    !CALL DefaultDirichletBCs()
124
125    Norm = DefaultSolve()
126
127    IF( Solver % Variable % NonlinConverged == 1 ) EXIT
128  END DO
129
130  CALL Info('OdeSolver','All done',Level=5)
131
132
133CONTAINS
134
135
136  ! Creates the local matrix equation for the ODY before time integration.
137  ! This optionally uses > Active Components < where the component entry
138  ! could be a natural place for parameters of the ODE system.
139  !-----------------------------------------------------------------------
140
141  SUBROUTINE OdeMatrix()
142
143    TYPE(Matrix_t), POINTER :: A
144    INTEGER :: i,j,Dofs
145    REAL(KIND=dp), POINTER CONTIG :: SaveValues(:)
146    REAL(KIND=dp) :: val
147    TYPE(ValueList_t), POINTER :: OdeList
148    INTEGER, POINTER :: ActiveComponents(:)
149
150
151    ActiveComponents => ListGetIntegerArray( Params, &
152        'Active Components', Found )
153    IF( Found ) THEN
154      IF( SIZE( ActiveComponents ) > 1 ) THEN
155        CALL Fatal('OdeSolver','Currently implemented only for one component!')
156      END IF
157      i = ActiveComponents(1)
158      IF( i > CurrentModel % NumberOfComponents ) THEN
159        CALL Fatal('OdeSolver','Active Component index out of range')
160      END IF
161      OdeList => CurrentModel % Components(i) % Values
162      CALL Info('OdeSolver','Using active component: '//TRIM(I2S(i)),Level=10)
163    ELSE
164      OdeList => Params
165      CALL Info('OdeSolver','Using solver section',Level=10)
166    END IF
167
168    A => Solver % Matrix
169    Dofs = Solver % Variable % Dofs
170
171    IF ( .NOT. ASSOCIATED( A % MassValues ) ) THEN
172      CALL Info('OdeSolver','Allocating mass matrix',Level=10)
173      ALLOCATE( A % MassValues(SIZE(A % Values)) )
174      A % MassValues = 0.0_dp
175    END IF
176
177    IF( TimeOrder >= 2 ) THEN
178      IF ( .NOT. ASSOCIATED( A % DampValues ) ) THEN
179        CALL Info('OdeSolver','Allocating damping matrix',Level=10)
180        ALLOCATE( A % DampValues(SIZE(A % Values)) )
181        A % DampValues = 0.0_dp
182      END IF
183    END IF
184
185    ! Set stiffness matrix values
186    DO i=1,Dofs
187      DO j=1,Dofs
188        str = 'Stiffness Matrix '//TRIM(I2S(i))//TRIM(I2S(j))
189        val = ListGetCReal( OdeList, str, Found )
190        CALL CRS_SetMatrixElement( A,i,j,val )
191      END DO
192    END DO
193
194    ! Set mass matrix values
195    SaveValues => A % Values
196    A % Values => A % MassValues
197    DO i=1,Dofs
198      DO j=1,Dofs
199
200        ! In Elmer the convention is to call the highest time-derivatie always mass
201        ! Even for 1st order PDEs. However, in this solver we like to call the coefficients
202        ! mass, damping and spring.
203        IF( TimeOrder >= 2 ) THEN
204          str = 'Mass Matrix '//TRIM(I2S(i))//TRIM(I2S(j))
205        ELSE
206          str = 'Damping Matrix '//TRIM(I2S(i))//TRIM(I2S(j))
207        END IF
208        val = ListGetCReal( OdeList, str, Found )
209        CALL CRS_SetMatrixElement( A,i,j,val )
210      END DO
211    END DO
212    A % Values => SaveValues
213
214    IF( TimeOrder >= 2 ) THEN
215      ! Set damp matrix values
216      SaveValues => A % Values
217      A % Values => A % DampValues
218      DO i=1,Dofs
219        DO j=1,Dofs
220          str = 'Damping Matrix '//TRIM(I2S(i))//TRIM(I2S(j))
221          val = ListGetCReal( OdeList, str, Found )
222          CALL CRS_SetMatrixElement( A,i,j,val )
223        END DO
224      END DO
225      A % Values => SaveValues
226    END IF
227
228    DO i=1,Dofs
229      str = 'Force '//TRIM(I2S(i))
230      val = ListGetCReal( OdeList, str, Found )
231      A % Rhs(i) = val
232    END DO
233
234    IF(.FALSE.) THEN
235      PRINT *,'Dofs:',Dofs
236      PRINT *,'A % Values',A % values
237      PRINT *,'A % DampValues',A % Dampvalues
238      IF( TimeOrder >= 2 ) THEN
239        PRINT *,'A % MassValues',A % Massvalues
240      END IF
241      PRINT *,'A % Rhs',A % Rhs
242    END IF
243
244  END SUBROUTINE OdeMatrix
245
246!------------------------------------------------------------------------------
247END SUBROUTINE OdeSolver
248!------------------------------------------------------------------------------
249