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