/FluidMechanics/Stokes/RoutineCheck/Dynamic/src/DynamicExample.f90

  1!> \file
  2!> \author Sebastian Krittian
  3!> \brief This is an example program to solve a dynamic Stokes equation using OpenCMISS calls.
  4!>
  5!> \section LICENSE
  6!>
  7!> Version: MPL 1.1/GPL 2.0/LGPL 2.1
  8!>
  9!> The contents of this file are subject to the Mozilla Public License
 10!> Version 1.1 (the "License"); you may not use this file except in
 11!> compliance with the License. You may obtain a copy of the License at
 12!> http://www.mozilla.org/MPL/
 13!>
 14!> Software distributed under the License is distributed on an "AS IS"
 15!> basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
 16!> License for the specific language governing rights and limitations
 17!> under the License.
 18!>
 19!> The Original Code is OpenCMISS
 20!>
 21!> The Initial Developer of the Original Code is University of Auckland,
 22!> Auckland, New Zealand and University of Oxford, Oxford, United
 23!> Kingdom. Portions created by the University of Auckland and University
 24!> of Oxford are Copyright (C) 2007 by the University of Auckland and
 25!> the University of Oxford. All Rights Reserved.
 26!>
 27!> Contributor(s):
 28!>
 29!> Alternatively, the contents of this file may be used under the terms of
 30!> either the GNU General Public License Version 2 or later (the "GPL"), or
 31!> the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
 32!> in which case the provisions of the GPL or the LGPL are applicable instead
 33!> of those above. If you wish to allow use of your version of this file only
 34!> under the terms of either the GPL or the LGPL, and not to allow others to
 35!> use your version of this file under the terms of the MPL, indicate your
 36!> decision by deleting the provisions above and replace them with the notice
 37!> and other provisions required by the GPL or the LGPL. If you do not delete
 38!> the provisions above, a recipient may use your version of this file under
 39!> the terms of any one of the MPL, the GPL or the LGPL.
 40!>
 41
 42!> \example FluidMechanics/Stokes/RoutineCheck/Dynamic/src/DynamicExample.f90
 43!! Example program to solve a dynamic Stokes equation using OpenCMISS calls.
 44!! \htmlinclude FluidMechanics/Stokes/RoutineCheck/Dynamic/history.html
 45!!
 46!<
 47
 48!> Main program
 49
 50PROGRAM STOKESDYNAMICEXAMPLE
 51
 52  !
 53  !================================================================================================================================
 54  !
 55
 56  !PROGRAM LIBRARIES
 57
 58  USE OPENCMISS
 59  USE FLUID_MECHANICS_IO_ROUTINES
 60  USE MPI
 61
 62#ifdef WIN32
 63  USE IFQWINCMISS
 64#endif
 65
 66  !
 67  !================================================================================================================================
 68  !
 69
 70  !PROGRAM VARIABLES AND TYPES
 71
 72  IMPLICIT NONE
 73
 74  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumber=1337
 75  TYPE(CMISSFieldType) :: EquationsSetField
 76
 77
 78  !Test program parameters
 79
 80  INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
 81  INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=2
 82  INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=3
 83  INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=4
 84  INTEGER(CMISSIntg), PARAMETER :: GeometricFieldUserNumber=5
 85  INTEGER(CMISSIntg), PARAMETER :: DependentFieldUserNumberStokes=6
 86  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberStokes=7
 87  INTEGER(CMISSIntg), PARAMETER :: EquationsSetUserNumberStokes=8
 88  INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=9
 89
 90  INTEGER(CMISSIntg), PARAMETER :: DomainUserNumber=2
 91  INTEGER(CMISSIntg), PARAMETER :: SolverStokesUserNumber=1
 92  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberStokesMu=1
 93  INTEGER(CMISSIntg), PARAMETER :: MaterialsFieldUserNumberStokesRho=2
 94
 95  !Program types
 96
 97  TYPE(EXPORT_CONTAINER):: CM
 98
 99  !Program variables
100
101  INTEGER(CMISSIntg) :: NUMBER_OF_DIMENSIONS
102  
103  INTEGER(CMISSIntg) :: BASIS_TYPE
104  INTEGER(CMISSIntg) :: BASIS_NUMBER_SPACE
105  INTEGER(CMISSIntg) :: BASIS_NUMBER_VELOCITY
106  INTEGER(CMISSIntg) :: BASIS_NUMBER_PRESSURE
107  INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_SPACE
108  INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_VELOCITY
109  INTEGER(CMISSIntg) :: BASIS_XI_GAUSS_PRESSURE
110  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_SPACE
111  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_VELOCITY
112  INTEGER(CMISSIntg) :: BASIS_XI_INTERPOLATION_PRESSURE
113  INTEGER(CMISSIntg) :: MESH_NUMBER_OF_COMPONENTS
114  INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_SPACE
115  INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_VELOCITY
116  INTEGER(CMISSIntg) :: MESH_COMPONENT_NUMBER_PRESSURE
117  INTEGER(CMISSIntg) :: NUMBER_OF_NODES_SPACE
118  INTEGER(CMISSIntg) :: NUMBER_OF_NODES_VELOCITY
119  INTEGER(CMISSIntg) :: NUMBER_OF_NODES_PRESSURE
120  INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_SPACE
121  INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_VELOCITY
122  INTEGER(CMISSIntg) :: NUMBER_OF_ELEMENT_NODES_PRESSURE
123  INTEGER(CMISSIntg) :: TOTAL_NUMBER_OF_NODES
124  INTEGER(CMISSIntg) :: TOTAL_NUMBER_OF_ELEMENTS
125  INTEGER(CMISSIntg) :: MAXIMUM_ITERATIONS
126  INTEGER(CMISSIntg) :: RESTART_VALUE
127!   INTEGER(CMISSIntg) :: MPI_IERROR
128  INTEGER(CMISSIntg) :: NUMBER_OF_FIXED_WALL_NODES_STOKES
129  INTEGER(CMISSIntg) :: NUMBER_OF_INLET_WALL_NODES_STOKES
130
131  INTEGER(CMISSIntg) :: EQUATIONS_STOKES_OUTPUT
132  INTEGER(CMISSIntg) :: COMPONENT_NUMBER
133  INTEGER(CMISSIntg) :: NODE_NUMBER
134  INTEGER(CMISSIntg) :: ELEMENT_NUMBER
135  INTEGER(CMISSIntg) :: NODE_COUNTER
136  INTEGER(CMISSIntg) :: CONDITION
137
138  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_STOKES_OUTPUT_FREQUENCY
139  INTEGER(CMISSIntg) :: DYNAMIC_SOLVER_STOKES_OUTPUT_TYPE
140  INTEGER(CMISSIntg) :: LINEAR_SOLVER_STOKES_OUTPUT_TYPE
141
142  INTEGER, ALLOCATABLE, DIMENSION(:):: FIXED_WALL_NODES_STOKES
143  INTEGER, ALLOCATABLE, DIMENSION(:):: INLET_WALL_NODES_STOKES
144
145
146  REAL(CMISSDP) :: INITIAL_FIELD_STOKES(3)
147  REAL(CMISSDP) :: BOUNDARY_CONDITIONS_STOKES(3)
148  REAL(CMISSDP) :: DIVERGENCE_TOLERANCE
149  REAL(CMISSDP) :: RELATIVE_TOLERANCE
150  REAL(CMISSDP) :: ABSOLUTE_TOLERANCE
151  REAL(CMISSDP) :: LINESEARCH_ALPHA
152  REAL(CMISSDP) :: VALUE
153  REAL(CMISSDP) :: MU_PARAM_STOKES
154  REAL(CMISSDP) :: RHO_PARAM_STOKES
155
156  REAL(CMISSDP) :: DYNAMIC_SOLVER_STOKES_START_TIME
157  REAL(CMISSDP) :: DYNAMIC_SOLVER_STOKES_STOP_TIME
158  REAL(CMISSDP) :: DYNAMIC_SOLVER_STOKES_THETA
159  REAL(CMISSDP) :: DYNAMIC_SOLVER_STOKES_TIME_INCREMENT
160
161  LOGICAL :: EXPORT_FIELD_IO
162  LOGICAL :: LINEAR_SOLVER_STOKES_DIRECT_FLAG
163  LOGICAL :: FIXED_WALL_NODES_STOKES_FLAG
164  LOGICAL :: INLET_WALL_NODES_STOKES_FLAG
165
166  !CMISS variables
167
168  !Regions
169  TYPE(CMISSRegionType) :: Region
170  TYPE(CMISSRegionType) :: WorldRegion
171  !Coordinate systems
172  TYPE(CMISSCoordinateSystemType) :: CoordinateSystem
173  TYPE(CMISSCoordinateSystemType) :: WorldCoordinateSystem
174  !Basis
175  TYPE(CMISSBasisType) :: BasisSpace
176  TYPE(CMISSBasisType) :: BasisVelocity
177  TYPE(CMISSBasisType) :: BasisPressure
178  !Nodes
179  TYPE(CMISSNodesType) :: Nodes
180  !Elements
181  TYPE(CMISSMeshElementsType) :: MeshElementsSpace
182  TYPE(CMISSMeshElementsType) :: MeshElementsVelocity
183  TYPE(CMISSMeshElementsType) :: MeshElementsPressure
184  !Meshes
185  TYPE(CMISSMeshType) :: Mesh
186  !Decompositions
187  TYPE(CMISSDecompositionType) :: Decomposition
188  !Fields
189  TYPE(CMISSFieldsType) :: Fields
190  !Field types
191  TYPE(CMISSFieldType) :: GeometricField
192  TYPE(CMISSFieldType) :: DependentFieldStokes
193  TYPE(CMISSFieldType) :: MaterialsFieldStokes
194  !Boundary conditions
195  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditionsStokes
196  !Equations sets
197  TYPE(CMISSEquationsSetType) :: EquationsSetStokes
198  !Equations
199  TYPE(CMISSEquationsType) :: EquationsStokes
200  !Problems
201  TYPE(CMISSProblemType) :: Problem
202  !Control loops
203  TYPE(CMISSControlLoopType) :: ControlLoop
204  !Solvers
205  TYPE(CMISSSolverType) :: DynamicSolverStokes
206  TYPE(CMISSSolverType) :: LinearSolverStokes
207  !Solver equations
208  TYPE(CMISSSolverEquationsType) :: SolverEquationsStokes
209
210#ifdef WIN32
211  !Quickwin type
212  LOGICAL :: QUICKWIN_STATUS=.FALSE.
213  TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
214#endif
215  
216  !Generic CMISS variables
217
218  INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber,BoundaryNodeDomain
219  INTEGER(CMISSIntg) :: EquationsSetIndex
220  INTEGER(CMISSIntg) :: Err
221  
222#ifdef WIN32
223  !Initialise QuickWin
224  QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
225  QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
226  QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
227  !Set the window parameters
228  QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
229  !If attempt fails set with system estimated values
230  IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
231#endif
232
233  !
234  !================================================================================================================================
235  !
236
237  !INITIALISE OPENCMISS
238
239  CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
240
241  CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
242
243  !
244  !================================================================================================================================
245  !
246
247  !CHECK COMPUTATIONAL NODE
248
249  !Get the computational nodes information
250  CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
251  CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
252
253  !
254  !================================================================================================================================
255  !
256
257  !PROBLEM CONTROL PANEL
258
259  !Import cmHeart mesh information
260  CALL FLUID_MECHANICS_IO_READ_CMHEART(CM,Err)  
261
262  BASIS_NUMBER_SPACE=CM%ID_M
263  BASIS_NUMBER_VELOCITY=CM%ID_V
264  BASIS_NUMBER_PRESSURE=CM%ID_P
265  NUMBER_OF_DIMENSIONS=CM%D
266  BASIS_TYPE=CM%IT_T
267  BASIS_XI_INTERPOLATION_SPACE=CM%IT_M
268  BASIS_XI_INTERPOLATION_VELOCITY=CM%IT_V
269  BASIS_XI_INTERPOLATION_PRESSURE=CM%IT_P
270  NUMBER_OF_NODES_SPACE=CM%N_M
271  NUMBER_OF_NODES_VELOCITY=CM%N_V
272  NUMBER_OF_NODES_PRESSURE=CM%N_P
273  TOTAL_NUMBER_OF_NODES=CM%N_T
274  TOTAL_NUMBER_OF_ELEMENTS=CM%E_T
275  NUMBER_OF_ELEMENT_NODES_SPACE=CM%EN_M
276  NUMBER_OF_ELEMENT_NODES_VELOCITY=CM%EN_V
277  NUMBER_OF_ELEMENT_NODES_PRESSURE=CM%EN_P
278  !Set initial values
279  INITIAL_FIELD_STOKES(1)=0.0_CMISSDP
280  INITIAL_FIELD_STOKES(2)=0.0_CMISSDP
281  INITIAL_FIELD_STOKES(3)=0.0_CMISSDP
282  !Set boundary conditions
283  FIXED_WALL_NODES_STOKES_FLAG=.TRUE.
284  INLET_WALL_NODES_STOKES_FLAG=.TRUE.
285  IF(FIXED_WALL_NODES_STOKES_FLAG) THEN
286    NUMBER_OF_FIXED_WALL_NODES_STOKES=80
287    ALLOCATE(FIXED_WALL_NODES_STOKES(NUMBER_OF_FIXED_WALL_NODES_STOKES))
288    FIXED_WALL_NODES_STOKES=(/1,2,3,4,5,7,9,10,11,12,13,14,17,20,24,28,29,30,31,32,33,34,35,37,39, & 
289    & 41,44,46,47,48,50,51,52,53,54,57,60,64,65,66,67,68,70,72,74,76,77,78,79,80,83,86, & 
290    & 89,90,91,92,93,94,95,97,99,101,102,103,104,105,106,107,108,111,114,115,116,117,118, & 
291    & 120,122,123,124,125/)
292  ENDIF
293  IF(INLET_WALL_NODES_STOKES_FLAG) THEN
294    NUMBER_OF_INLET_WALL_NODES_STOKES=9
295    ALLOCATE(INLET_WALL_NODES_STOKES(NUMBER_OF_INLET_WALL_NODES_STOKES))
296    INLET_WALL_NODES_STOKES=(/6,15,16,23,36,42,81,82,96/)
297    !Set initial boundary conditions
298    BOUNDARY_CONDITIONS_STOKES(1)=0.0_CMISSDP
299    BOUNDARY_CONDITIONS_STOKES(2)=1.0_CMISSDP
300    BOUNDARY_CONDITIONS_STOKES(3)=0.0_CMISSDP
301  ENDIF
302  !Set material parameters
303  MU_PARAM_STOKES=1.0_CMISSDP
304  RHO_PARAM_STOKES=1.0_CMISSDP
305  !Set interpolation parameters
306  BASIS_XI_GAUSS_SPACE=3
307  BASIS_XI_GAUSS_VELOCITY=3
308  BASIS_XI_GAUSS_PRESSURE=3
309  !Set output parameter
310  !(NoOutput/ProgressOutput/TimingOutput/SolverOutput/SolverMatrixOutput)
311  DYNAMIC_SOLVER_STOKES_OUTPUT_TYPE=CMISS_SOLVER_NO_OUTPUT
312  LINEAR_SOLVER_STOKES_OUTPUT_TYPE=CMISS_SOLVER_NO_OUTPUT
313  !(NoOutput/TimingOutput/MatrixOutput/ElementOutput)
314  EQUATIONS_STOKES_OUTPUT=CMISS_EQUATIONS_NO_OUTPUT
315  !Set time parameter
316  DYNAMIC_SOLVER_STOKES_START_TIME=0.0_CMISSDP
317  DYNAMIC_SOLVER_STOKES_STOP_TIME=5.0_CMISSDP 
318  DYNAMIC_SOLVER_STOKES_TIME_INCREMENT=1.0_CMISSDP
319  DYNAMIC_SOLVER_STOKES_THETA=1.0_CMISSDP
320  !Set result output parameter
321  DYNAMIC_SOLVER_STOKES_OUTPUT_FREQUENCY=1
322  !Set solver parameters
323  LINEAR_SOLVER_STOKES_DIRECT_FLAG=.FALSE.
324  RELATIVE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-05_CMISSDP
325  ABSOLUTE_TOLERANCE=1.0E-10_CMISSDP !default: 1.0E-10_CMISSDP
326  DIVERGENCE_TOLERANCE=1.0E20 !default: 1.0E5
327  MAXIMUM_ITERATIONS=100000 !default: 100000
328  RESTART_VALUE=3000 !default: 30
329  LINESEARCH_ALPHA=1.0
330
331  !
332  !================================================================================================================================
333  !
334
335  !COORDINATE SYSTEM
336
337  !Start the creation of a new RC coordinate system
338  CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
339  CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
340  !Set the coordinate system dimension
341  CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NUMBER_OF_DIMENSIONS,Err)
342  !Finish the creation of the coordinate system
343  CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
344
345  !
346  !================================================================================================================================
347  !
348
349  !REGION
350
351  !Start the creation of a new region
352  CALL CMISSRegion_Initialise(Region,Err)
353  CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
354  !Set the regions coordinate system as defined above
355  CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
356  !Finish the creation of the region
357  CALL CMISSRegion_CreateFinish(Region,Err)
358
359  !
360  !================================================================================================================================
361  !
362
363  !BASES
364
365  !Start the creation of new bases
366  MESH_NUMBER_OF_COMPONENTS=1
367  CALL CMISSBasis_Initialise(BasisSpace,Err)
368  CALL CMISSBasis_CreateStart(BASIS_NUMBER_SPACE,BasisSpace,Err)
369  !Set the basis type (Lagrange/Simplex)
370  CALL CMISSBasis_TypeSet(BasisSpace,BASIS_TYPE,Err)
371  !Set the basis xi number
372  CALL CMISSBasis_NumberOfXiSet(BasisSpace,NUMBER_OF_DIMENSIONS,Err)
373  !Set the basis xi interpolation and number of Gauss points
374  IF(NUMBER_OF_DIMENSIONS==2) THEN
375    CALL CMISSBasis_InterpolationXiSet(BasisSpace,(/BASIS_XI_INTERPOLATION_SPACE,BASIS_XI_INTERPOLATION_SPACE/),Err)
376    CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisSpace,(/BASIS_XI_GAUSS_SPACE,BASIS_XI_GAUSS_SPACE/),Err)
377  ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
378    CALL CMISSBasis_InterpolationXiSet(BasisSpace,(/BASIS_XI_INTERPOLATION_SPACE,BASIS_XI_INTERPOLATION_SPACE, & 
379      & BASIS_XI_INTERPOLATION_SPACE/),Err)                         
380    CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisSpace,(/BASIS_XI_GAUSS_SPACE,BASIS_XI_GAUSS_SPACE,BASIS_XI_GAUSS_SPACE/),Err)
381  ENDIF
382  !Finish the creation of the basis
383  CALL CMISSBasis_CreateFinish(BasisSpace,Err)
384  !Start the creation of another basis
385  IF(BASIS_XI_INTERPOLATION_VELOCITY==BASIS_XI_INTERPOLATION_SPACE) THEN
386    BasisVelocity=BasisSpace
387  ELSE
388    MESH_NUMBER_OF_COMPONENTS=MESH_NUMBER_OF_COMPONENTS+1
389    !Initialise a new velocity basis
390    CALL CMISSBasis_Initialise(BasisVelocity,Err)
391    !Start the creation of a basis
392    CALL CMISSBasis_CreateStart(BASIS_NUMBER_VELOCITY,BasisVelocity,Err)
393    !Set the basis type (Lagrange/Simplex)
394    CALL CMISSBasis_TypeSet(BasisVelocity,BASIS_TYPE,Err)
395    !Set the basis xi number
396    CALL CMISSBasis_NumberOfXiSet(BasisVelocity,NUMBER_OF_DIMENSIONS,Err)
397    !Set the basis xi interpolation and number of Gauss points
398    IF(NUMBER_OF_DIMENSIONS==2) THEN
399      CALL CMISSBasis_InterpolationXiSet(BasisVelocity,(/BASIS_XI_INTERPOLATION_VELOCITY,BASIS_XI_INTERPOLATION_VELOCITY/),Err)
400      CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisVelocity,(/BASIS_XI_GAUSS_VELOCITY,BASIS_XI_GAUSS_VELOCITY/),Err)
401    ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
402      CALL CMISSBasis_InterpolationXiSet(BasisVelocity,(/BASIS_XI_INTERPOLATION_VELOCITY,BASIS_XI_INTERPOLATION_VELOCITY, & 
403        & BASIS_XI_INTERPOLATION_VELOCITY/),Err)                         
404      CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisVelocity,(/BASIS_XI_GAUSS_VELOCITY,BASIS_XI_GAUSS_VELOCITY, & 
405        & BASIS_XI_GAUSS_VELOCITY/),Err)
406    ENDIF
407    !Finish the creation of the basis
408    CALL CMISSBasis_CreateFinish(BasisVelocity,Err)
409  ENDIF
410  !Start the creation of another basis
411  IF(BASIS_XI_INTERPOLATION_PRESSURE==BASIS_XI_INTERPOLATION_SPACE) THEN
412    BasisPressure=BasisSpace
413  ELSE IF(BASIS_XI_INTERPOLATION_PRESSURE==BASIS_XI_INTERPOLATION_VELOCITY) THEN
414    BasisPressure=BasisVelocity
415  ELSE
416    MESH_NUMBER_OF_COMPONENTS=MESH_NUMBER_OF_COMPONENTS+1
417    !Initialise a new pressure basis
418    CALL CMISSBasis_Initialise(BasisPressure,Err)
419    !Start the creation of a basis
420    CALL CMISSBasis_CreateStart(BASIS_NUMBER_PRESSURE,BasisPressure,Err)
421    !Set the basis type (Lagrange/Simplex)
422    CALL CMISSBasis_TypeSet(BasisPressure,BASIS_TYPE,Err)
423    !Set the basis xi number
424    CALL CMISSBasis_NumberOfXiSet(BasisPressure,NUMBER_OF_DIMENSIONS,Err)
425    !Set the basis xi interpolation and number of Gauss points
426    IF(NUMBER_OF_DIMENSIONS==2) THEN
427      CALL CMISSBasis_InterpolationXiSet(BasisPressure,(/BASIS_XI_INTERPOLATION_PRESSURE,BASIS_XI_INTERPOLATION_PRESSURE/),Err)
428      CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisPressure,(/BASIS_XI_GAUSS_PRESSURE,BASIS_XI_GAUSS_PRESSURE/),Err)
429    ELSE IF(NUMBER_OF_DIMENSIONS==3) THEN
430      CALL CMISSBasis_InterpolationXiSet(BasisPressure,(/BASIS_XI_INTERPOLATION_PRESSURE,BASIS_XI_INTERPOLATION_PRESSURE, & 
431        & BASIS_XI_INTERPOLATION_PRESSURE/),Err)                         
432      CALL CMISSBasis_QuadratureNumberOfGaussXiSet(BasisPressure,(/BASIS_XI_GAUSS_PRESSURE,BASIS_XI_GAUSS_PRESSURE, & 
433        & BASIS_XI_GAUSS_PRESSURE/),Err)
434    ENDIF
435    !Finish the creation of the basis
436    CALL CMISSBasis_CreateFinish(BasisPressure,Err)
437  ENDIF
438
439  !
440  !================================================================================================================================
441  !
442
443  !MESH
444
445  !Start the creation of mesh nodes
446  CALL CMISSNodes_Initialise(Nodes,Err)
447  CALL CMISSMesh_Initialise(Mesh,Err)
448  CALL CMISSNodes_CreateStart(Region,TOTAL_NUMBER_OF_NODES,Nodes,Err)
449  CALL CMISSNodes_CreateFinish(Nodes,Err)
450  !Start the creation of the mesh
451  CALL CMISSMesh_CreateStart(MeshUserNumber,Region,NUMBER_OF_DIMENSIONS,Mesh,Err)
452  !Set number of mesh elements
453  CALL CMISSMesh_NumberOfElementsSet(Mesh,TOTAL_NUMBER_OF_ELEMENTS,Err)
454  !Set number of mesh components
455  CALL CMISSMesh_NumberOfComponentsSet(Mesh,MESH_NUMBER_OF_COMPONENTS,Err)
456  !Specify spatial mesh component
457  CALL CMISSMeshElements_Initialise(MeshElementsSpace,Err)
458  CALL CMISSMeshElements_Initialise(MeshElementsVelocity,Err)
459  CALL CMISSMeshElements_Initialise(MeshElementsPressure,Err)
460  MESH_COMPONENT_NUMBER_SPACE=1
461  MESH_COMPONENT_NUMBER_VELOCITY=1
462  MESH_COMPONENT_NUMBER_PRESSURE=1
463  CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_SPACE,BasisSpace,MeshElementsSpace,Err)
464  DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
465    CALL CMISSMeshElements_NodesSet(MeshElementsSpace,ELEMENT_NUMBER,CM%M(ELEMENT_NUMBER,1:NUMBER_OF_ELEMENT_NODES_SPACE),Err)
466  ENDDO
467  CALL CMISSMeshElements_CreateFinish(MeshElementsSpace,Err)
468  !Specify velocity mesh component
469  IF(BASIS_XI_INTERPOLATION_VELOCITY==BASIS_XI_INTERPOLATION_SPACE) THEN
470    MeshElementsVelocity=MeshElementsSpace
471  ELSE
472    MESH_COMPONENT_NUMBER_VELOCITY=MESH_COMPONENT_NUMBER_SPACE+1
473    CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_VELOCITY,BasisVelocity,MeshElementsVelocity,Err)
474    DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
475      CALL CMISSMeshElements_NodesSet(MeshElementsVelocity,ELEMENT_NUMBER,CM%V(ELEMENT_NUMBER, & 
476        & 1:NUMBER_OF_ELEMENT_NODES_VELOCITY),Err)
477    ENDDO
478    CALL CMISSMeshElements_CreateFinish(MeshElementsVelocity,Err)
479  ENDIF
480  !Specify pressure mesh component
481  IF(BASIS_XI_INTERPOLATION_PRESSURE==BASIS_XI_INTERPOLATION_SPACE) THEN
482    MeshElementsPressure=MeshElementsSpace
483    MESH_COMPONENT_NUMBER_PRESSURE=MESH_COMPONENT_NUMBER_SPACE
484  ELSE IF(BASIS_XI_INTERPOLATION_PRESSURE==BASIS_XI_INTERPOLATION_VELOCITY) THEN
485    MeshElementsPressure=MeshElementsVelocity
486    MESH_COMPONENT_NUMBER_PRESSURE=MESH_COMPONENT_NUMBER_VELOCITY
487  ELSE
488    MESH_COMPONENT_NUMBER_PRESSURE=MESH_COMPONENT_NUMBER_VELOCITY+1
489    CALL CMISSMeshElements_CreateStart(Mesh,MESH_COMPONENT_NUMBER_PRESSURE,BasisPressure,MeshElementsPressure,Err)
490    DO ELEMENT_NUMBER=1,TOTAL_NUMBER_OF_ELEMENTS
491      CALL CMISSMeshElements_NodesSet(MeshElementsPressure,ELEMENT_NUMBER,CM%P(ELEMENT_NUMBER, & 
492        & 1:NUMBER_OF_ELEMENT_NODES_PRESSURE),Err)
493    ENDDO
494    CALL CMISSMeshElements_CreateFinish(MeshElementsPressure,Err)
495  ENDIF
496  !Finish the creation of the mesh
497  CALL CMISSMesh_CreateFinish(Mesh,Err)
498
499  !
500  !================================================================================================================================
501  !
502
503  !GEOMETRIC FIELD
504
505  !Create a decomposition
506  CALL CMISSDecomposition_Initialise(Decomposition,Err)
507  CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
508  !Set the decomposition to be a general decomposition with the specified number of domains
509  CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
510  CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfComputationalNodes,Err)
511  !Finish the decomposition
512  CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
513
514  !Start to create a default (geometric) field on the region
515  CALL CMISSField_Initialise(GeometricField,Err)
516  CALL CMISSField_CreateStart(GeometricFieldUserNumber,Region,GeometricField,Err)
517  !Set the field type
518  CALL CMISSField_TypeSet(GeometricField,CMISS_FIELD_GEOMETRIC_TYPE,Err)
519  !Set the decomposition to use
520  CALL CMISSField_MeshDecompositionSet(GeometricField,Decomposition,Err)
521  !Set the scaling to use
522  CALL CMISSField_ScalingTypeSet(GeometricField,CMISS_FIELD_NO_SCALING,Err)
523  !Set the mesh component to be used by the field components.
524  DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
525    CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, & 
526      & MESH_COMPONENT_NUMBER_SPACE,Err)
527  ENDDO
528  !Finish creating the field
529  CALL CMISSField_CreateFinish(GeometricField,Err)
530  !Update the geometric field parameters
531  DO NODE_NUMBER=1,NUMBER_OF_NODES_SPACE
532    DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
533      VALUE=CM%N(NODE_NUMBER,COMPONENT_NUMBER)
534      CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE_NUMBER,1,BoundaryNodeDomain,Err)
535      IF(BoundaryNodeDomain==ComputationalNodeNumber) THEN
536        CALL CMISSField_ParameterSetUpdateNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
537          & 1,CMISS_NO_GLOBAL_DERIV,NODE_NUMBER,COMPONENT_NUMBER,VALUE,Err)
538      ENDIF
539    ENDDO
540  ENDDO
541  CALL CMISSField_ParameterSetUpdateStart(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
542  CALL CMISSField_ParameterSetUpdateFinish(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,Err)
543
544
545
546  !
547  !================================================================================================================================
548  !
549
550  !EQUATIONS SETS
551
552  !Create the equations set for dynamic Stokes
553  CALL CMISSEquationsSet_Initialise(EquationsSetStokes,Err)
554    CALL CMISSField_Initialise(EquationsSetField,Err)
555CALL CMISSEquationsSet_CreateStart(EquationsSetUserNumberStokes,Region,GeometricField,CMISS_EQUATIONS_SET_FLUID_MECHANICS_CLASS, &
556    & CMISS_EQUATIONS_SET_STOKES_EQUATION_TYPE,CMISS_EQUATIONS_SET_TRANSIENT_STOKES_SUBTYPE,EquationsSetFieldUserNumber, &
557      & EquationsSetField, &
558    & EquationsSetStokes,Err)
559  !Set the equations set to be a dynamic Stokes problem
560  
561  !Finish creating the equations set
562  CALL CMISSEquationsSet_CreateFinish(EquationsSetStokes,Err)
563
564
565  !
566  !================================================================================================================================
567  !
568
569  !DEPENDENT FIELDS
570
571  !Create the equations set dependent field variables for dynamic Stokes
572  CALL CMISSField_Initialise(DependentFieldStokes,Err)
573  CALL CMISSEquationsSet_DependentCreateStart(EquationsSetStokes,DependentFieldUserNumberStokes, & 
574    & DependentFieldStokes,Err)
575  !Set the mesh component to be used by the field components.
576  DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
577    CALL CMISSField_ComponentMeshComponentSet(DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, & 
578      & MESH_COMPONENT_NUMBER_VELOCITY,Err)
579    CALL CMISSField_ComponentMeshComponentSet(DependentFieldStokes,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, & 
580      & MESH_COMPONENT_NUMBER_VELOCITY,Err)
581  ENDDO
582  COMPONENT_NUMBER=NUMBER_OF_DIMENSIONS+1
583    CALL CMISSField_ComponentMeshComponentSet(DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,COMPONENT_NUMBER, & 
584      & MESH_COMPONENT_NUMBER_PRESSURE,Err)
585    CALL CMISSField_ComponentMeshComponentSet(DependentFieldStokes,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,COMPONENT_NUMBER, & 
586      & MESH_COMPONENT_NUMBER_PRESSURE,Err)
587  !Finish the equations set dependent field variables
588  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSetStokes,Err)
589
590  !Initialise dependent field
591  DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
592    CALL CMISSField_ComponentValuesInitialise(DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
593      & COMPONENT_NUMBER,INITIAL_FIELD_STOKES(COMPONENT_NUMBER),Err)
594  ENDDO
595
596
597  !
598  !================================================================================================================================
599  !
600
601  !MATERIALS FIELDS
602
603  !Create the equations set materials field variables for dynamic Stokes
604  CALL CMISSField_Initialise(MaterialsFieldStokes,Err)
605  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSetStokes,MaterialsFieldUserNumberStokes, & 
606    & MaterialsFieldStokes,Err)
607  !Finish the equations set materials field variables
608  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSetStokes,Err)
609  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
610    & MaterialsFieldUserNumberStokesMu,MU_PARAM_STOKES,Err)
611  CALL CMISSField_ComponentValuesInitialise(MaterialsFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, & 
612    & MaterialsFieldUserNumberStokesRho,RHO_PARAM_STOKES,Err)
613
614  !
615  !================================================================================================================================
616  !
617
618  !EQUATIONS
619
620
621  !Create the equations set equations
622  CALL CMISSEquations_Initialise(EquationsStokes,Err)
623  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSetStokes,EquationsStokes,Err)
624  !Set the equations matrices sparsity type
625  CALL CMISSEquations_SparsityTypeSet(EquationsStokes,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
626  !Set the equations lumping type
627  CALL CMISSEquations_LumpingTypeSet(EquationsStokes,CMISS_EQUATIONS_UNLUMPED_MATRICES,Err)
628  !Set the equations set output
629  CALL CMISSEquations_OutputTypeSet(EquationsStokes,EQUATIONS_STOKES_OUTPUT,Err)
630  !Finish the equations set equations
631  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSetStokes,Err)
632
633
634  !
635  !================================================================================================================================
636  !
637
638  !PROBLEMS
639
640  !Start the creation of a problem.
641  CALL CMISSProblem_Initialise(Problem,Err)
642  CALL CMISSControlLoop_Initialise(ControlLoop,Err)
643  CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err)
644  !Set the problem to be a dynamic Stokes problem
645  CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_FLUID_MECHANICS_CLASS,CMISS_PROBLEM_STOKES_EQUATION_TYPE, &
646    & CMISS_PROBLEM_TRANSIENT_STOKES_SUBTYPE,Err)
647  !Finish the creation of a problem.
648  CALL CMISSProblem_CreateFinish(Problem,Err)
649  !Start the creation of the problem control loop
650  CALL CMISSProblem_ControlLoopCreateStart(Problem,Err)
651  !Get the control loop
652  CALL CMISSProblem_ControlLoopGet(Problem,CMISS_CONTROL_LOOP_NODE,ControlLoop,Err)
653  !Set the times
654  CALL CMISSControlLoop_TimesSet(ControlLoop,DYNAMIC_SOLVER_STOKES_START_TIME,DYNAMIC_SOLVER_STOKES_STOP_TIME, & 
655    & DYNAMIC_SOLVER_STOKES_TIME_INCREMENT,Err)
656  !Set the output timing
657  CALL CMISSControlLoop_TimeOutputSet(ControlLoop,DYNAMIC_SOLVER_STOKES_OUTPUT_FREQUENCY,Err)
658  !Finish creating the problem control loop
659  CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err)
660
661  !
662  !================================================================================================================================
663  !
664
665  !SOLVERS
666
667  !Start the creation of the problem solvers
668  CALL CMISSSolver_Initialise(DynamicSolverStokes,Err)
669  CALL CMISSSolver_Initialise(LinearSolverStokes,Err)
670  CALL CMISSProblem_SolversCreateStart(Problem,Err)
671  !Get the dynamic dymamic solver
672  CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,SolverStokesUserNumber,DynamicSolverStokes,Err)
673  !Set the output type
674  CALL CMISSSolver_OutputTypeSet(DynamicSolverStokes,DYNAMIC_SOLVER_STOKES_OUTPUT_TYPE,Err)
675  !Set theta
676  CALL CMISSSolver_DynamicThetaSet(DynamicSolverStokes,DYNAMIC_SOLVER_STOKES_THETA,Err)
677!   CALL CMISSSolverDynamicDynamicSet(DynamicSolverStokes,.TRUE.,Err)
678  !Get the dynamic linear solver
679  CALL CMISSSolver_DynamicLinearSolverGet(DynamicSolverStokes,LinearSolverStokes,Err)
680  !Set the output type
681  CALL CMISSSolver_OutputTypeSet(LinearSolverStokes,LINEAR_SOLVER_STOKES_OUTPUT_TYPE,Err)
682  !Set the solver settings
683  IF(LINEAR_SOLVER_STOKES_DIRECT_FLAG) THEN
684    CALL CMISSSolver_LinearTypeSet(LinearSolverStokes,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
685    CALL CMISSSolver_LibraryTypeSet(LinearSolverStokes,CMISS_SOLVER_MUMPS_LIBRARY,Err)
686  ELSE
687    CALL CMISSSolver_LinearTypeSet(LinearSolverStokes,CMISS_SOLVER_LINEAR_ITERATIVE_SOLVE_TYPE,Err)
688    CALL CMISSSolver_LinearIterativeMaximumIterationsSet(LinearSolverStokes,MAXIMUM_ITERATIONS,Err)
689    CALL CMISSSolver_LinearIterativeDivergenceToleranceSet(LinearSolverStokes,DIVERGENCE_TOLERANCE,Err)
690    CALL CMISSSolver_LinearIterativeRelativeToleranceSet(LinearSolverStokes,RELATIVE_TOLERANCE,Err)
691    CALL CMISSSolver_LinearIterativeAbsoluteToleranceSet(LinearSolverStokes,ABSOLUTE_TOLERANCE,Err)
692    CALL CMISSSolver_LinearIterativeGMRESRestartSet(LinearSolverStokes,RESTART_VALUE,Err)
693  ENDIF
694  !Finish the creation of the problem solver
695  CALL CMISSProblem_SolversCreateFinish(Problem,Err)
696
697  !
698  !================================================================================================================================
699  !
700
701  !SOLVER EQUATIONS
702
703  !Start the creation of the problem solver equations
704  CALL CMISSSolver_Initialise(DynamicSolverStokes,Err)
705  CALL CMISSSolverEquations_Initialise(SolverEquationsStokes,Err)
706
707  CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err)
708  !Get the dynamic solver equations
709  CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,SolverStokesUserNumber,DynamicSolverStokes,Err)
710  CALL CMISSSolver_SolverEquationsGet(DynamicSolverStokes,SolverEquationsStokes,Err)
711  !Set the solver equations sparsity
712  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquationsStokes,CMISS_SOLVER_SPARSE_MATRICES,Err)
713  !Add in the equations set
714  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquationsStokes,EquationsSetStokes,EquationsSetIndex,Err)
715  !Finish the creation of the problem solver equations
716  CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err)
717
718  !
719  !================================================================================================================================
720  !
721
722  !BOUNDARY CONDITIONS
723
724  !Start the creation of the equations set boundary conditions for Stokes
725  CALL CMISSBoundaryConditions_Initialise(BoundaryConditionsStokes,Err)
726  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquationsStokes,BoundaryConditionsStokes,Err)
727  !Set fixed wall nodes
728  IF(FIXED_WALL_NODES_STOKES_FLAG) THEN
729    DO NODE_COUNTER=1,NUMBER_OF_FIXED_WALL_NODES_STOKES
730      NODE_NUMBER=FIXED_WALL_NODES_STOKES(NODE_COUNTER)
731      CONDITION=CMISS_BOUNDARY_CONDITION_FIXED_WALL
732      CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE_NUMBER,1,BoundaryNodeDomain,Err)
733      IF(BoundaryNodeDomain==ComputationalNodeNumber) THEN
734        DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
735          VALUE=0.0_CMISSDP
736          CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsStokes,DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,1, &
737            & CMISS_NO_GLOBAL_DERIV, &
738            & NODE_NUMBER,COMPONENT_NUMBER,CONDITION,VALUE,Err)
739        ENDDO
740      ENDIF
741    ENDDO
742  ENDIF
743  !Set velocity boundary conditions
744  IF(INLET_WALL_NODES_STOKES_FLAG) THEN
745    DO NODE_COUNTER=1,NUMBER_OF_INLET_WALL_NODES_STOKES
746      NODE_NUMBER=INLET_WALL_NODES_STOKES(NODE_COUNTER)
747      CONDITION=CMISS_BOUNDARY_CONDITION_FIXED_INLET
748      CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE_NUMBER,1,BoundaryNodeDomain,Err)
749      IF(BoundaryNodeDomain==ComputationalNodeNumber) THEN
750        DO COMPONENT_NUMBER=1,NUMBER_OF_DIMENSIONS
751          VALUE=BOUNDARY_CONDITIONS_STOKES(COMPONENT_NUMBER)
752          CALL CMISSBoundaryConditions_SetNode(BoundaryConditionsStokes,DependentFieldStokes,CMISS_FIELD_U_VARIABLE_TYPE,1, &
753            & CMISS_NO_GLOBAL_DERIV, &
754            & NODE_NUMBER,COMPONENT_NUMBER,CONDITION,VALUE,Err)
755        ENDDO
756      ENDIF
757    ENDDO
758  ENDIF
759  !Finish the creation of the equations set boundary conditions
760  CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquationsStokes,Err)
761
762  !
763  !================================================================================================================================
764  !
765
766  !RUN SOLVERS
767
768  !Turn of PETSc error handling
769  !CALL PETSC_ERRORHANDLING_SET_ON(ERR,ERROR,*999)
770
771  !Solve the problem
772  WRITE(*,'(A)') "Solving problem..."
773  CALL CMISSProblem_Solve(Problem,Err)
774  WRITE(*,'(A)') "Problem solved!"
775
776  !
777  !================================================================================================================================
778  !
779
780  !OUTPUT
781
782  EXPORT_FIELD_IO=.TRUE.
783  IF(EXPORT_FIELD_IO) THEN
784    WRITE(*,'(A)') "Exporting fields..."
785    CALL CMISSFields_Initialise(Fields,Err)
786    CALL CMISSFields_Create(Region,Fields,Err)
787    CALL CMISSFields_NodesExport(Fields,"DynamicStokes","FORTRAN",Err)
788    CALL CMISSFields_ElementsExport(Fields,"DynamicStokes","FORTRAN",Err)
789    CALL CMISSFields_Finalise(Fields,Err)
790    WRITE(*,'(A)') "Field exported!"
791  ENDIF
792  
793  !Finialise CMISS
794!   CALL CMISSFinalise(Err)
795
796  WRITE(*,'(A)') "Program successfully completed."
797  
798  STOP
799
800END PROGRAM STOKESDYNAMICEXAMPLE