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