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! * EXPERIMENTAL 27! * 28! * Module for exporting parallel results in Xdmf/HDF5 file format 29! * 30! ****************************************************************************** 31! * 32! * Authors: Mikko Lyly 33! * 34! * Email: Juha.Ruokolainen@csc.fi 35! * Web: http://www.csc.fi/elmer 36! * Address: CSC - IT Center for Science Ltd. 37! * Keilaranta 14 38! * 02101 Espoo, Finland 39! * 40! * Original Date: 11 Feb 2011 41! * 42! *****************************************************************************/ 43! 44!------------------------------------------------------------------------------ 45SUBROUTINE XdmfWriter(Model, Solver, dt, TransientSimulation) 46!------------------------------------------------------------------------------ 47 USE HDF5 48 USE DefUtils 49 IMPLICIT NONE 50!------------------------------------------------------------------------------ 51 TYPE(Solver_t) :: Solver 52 TYPE(Model_t) :: Model 53 REAL(KIND=dp) :: dt 54 LOGICAL :: TransientSimulation 55!------------------------------------------------------------------------------ 56 INTEGER, PARAMETER :: MaxFields = 1000 57 INTEGER :: NofScalarFields, NofVectorFields 58 CHARACTER(LEN=MAX_NAME_LEN) :: ScalarFieldNames(MaxFields) 59 CHARACTER(LEN=MAX_NAME_LEN) :: VectorFieldNames(MaxFields) 60 CHARACTER(LEN=MAX_NAME_LEN) :: BaseFileName 61 TYPE(Mesh_t), POINTER :: Mesh 62 INTEGER :: PEs, MyPE, i, ierr 63 INTEGER, ALLOCATABLE :: NofNodes(:), NofElements(:), NofStorage(:), itmp(:) 64 INTEGER(HID_T) :: file_id, plist_id, fill_type_id 65 LOGICAL :: Found 66 CHARACTER(LEN=MAX_NAME_LEN) :: Str 67 REAL(KIND=dp) :: RealTime, StartTime, TotalTime 68 INTEGER:: Order203(3), Order510(10), Order820(20) 69 REAL(KIND=dp), ALLOCATABLE :: TimeValues(:), dtmp(:) 70 TYPE(Variable_t), POINTER :: Var 71 INTEGER :: Counter = 1 72 SAVE TimeValues, Counter 73!------------------------------------------------------------------------------ 74 StartTime = RealTime() 75 Mesh => GetMesh() 76 PEs = ParEnv % PEs 77 MyPE = ParEnv % MyPE + 1 78 79 ! Save simulation time history: 80 !------------------------------- 81 IF(Counter == 1) THEN 82 ALLOCATE(TimeValues(1)) 83 TimeValues(1) = GetTime() 84 ELSE 85 ALLOCATE(dtmp(SIZE(TimeValues))) 86 dtmp = TimeValues 87 DEALLOCATE(TimeValues) 88 ALLOCATE(TimeValues(SIZE(dtmp)+1)) 89 TimeValues(1:SIZE(dtmp)) = dtmp 90 TimeValues(SIZE(dtmp)+1) = GetTime() 91 DEALLOCATE(dtmp) 92 END IF 93 94 ! Determine the base file name and field variables: 95 !--------------------------------------------------- 96 BaseFileName = ListGetString(Solver % Values, 'base file name', Found) 97 IF(.NOT.Found) BaseFileName = 'results' 98 CALL INFO('XdmfWriter', 'Base file name: '//TRIM(BaseFileName)) 99 100 CALL FindScalarFields(NofScalarFields, ScalarFieldNames) 101 DO i = 1, NofScalarFields 102 CALL INFO('XdmfWriter', 'Scalar field: '//TRIM(ScalarFieldNames(i))) 103 END DO 104 105 CALL FindVectorFields(NofVectorFields, VectorFieldNames) 106 DO i = 1, NofVectorFields 107 CALL INFO('XdmfWriter', 'Vector field: '//TRIM(VectorFieldNames(i))) 108 END DO 109 110 ! Set up node permutation vectors for quadratic elements: 111 !--------------------------------------------------------- 112 Order203(:) = (/ 1,3,2 /) 113 Order510(:) = (/ 1,2,4,3,5,9,8,7,6,10/) 114 Order820(:) = (/ 1,2,3,4,5,6,7,8,9,10,11,12,17,18,19,20,13,14,15,16 /) 115 116 ! Determine Nof nodes, Nof elements and mixed xdmf storage size for all PEs: 117 !---------------------------------------------------------------------------- 118 ALLOCATE(itmp(PEs)) 119 120 ALLOCATE(NofNodes(PEs), NofElements(PEs), NofStorage(PEs)) 121 122 NofNodes = 0; itmp = 0; itmp(MyPE) = Mesh % NumberOfNodes 123 CALL MPI_ALLREDUCE(itmp, NofNodes, PEs, MPI_INTEGER, MPI_SUM, ELMER_COMM_WORLD, ierr) 124 125 NofElements = 0; itmp = 0; itmp(MyPE) = GetNofActive() 126 CALL MPI_ALLREDUCE(itmp, NofElements, PEs, MPI_INTEGER, MPI_SUM, ELMER_COMM_WORLD, ierr) 127 128 NofStorage = 0; itmp = 0; itmp(MyPE) = GetElementStorageSize() 129 CALL MPI_ALLREDUCE(itmp, NofStorage, PEs, MPI_INTEGER, MPI_SUM, ELMER_COMM_WORLD, ierr) 130 131 DEALLOCATE(itmp) 132 133 ! Initialize and create/open the hdf5 file collectively: 134 !-------------------------------------------------------- 135 CALL h5open_f(ierr) 136 CALL h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, ierr) 137 CALL h5pset_fapl_mpio_f(plist_id, ELMER_COMM_WORLD, MPI_INFO_NULL, ierr) 138 139 IF(Counter == 1) THEN 140 CALL h5fcreate_f(TRIM(BaseFileName)//'.h5', H5F_ACC_TRUNC_F, file_id, ierr, access_prp = plist_id) 141 ELSE 142 CALL h5fopen_f(TRIM(BaseFileName)//'.h5', H5F_ACC_RDWR_F, file_id, ierr, access_prp = plist_id) 143 END IF 144 145 CALL h5pclose_f(plist_id, ierr) 146 147 ! Determine output precision for reals: 148 !--------------------------------------- 149 fill_type_id = H5T_NATIVE_DOUBLE 150 IF(ListGetLogical(Solver % Values, 'single precision', Found)) fill_type_id = H5T_NATIVE_REAL 151 IF(.NOT.Found) fill_type_id = H5T_NATIVE_DOUBLE 152 153 ! Write nodes, elements and part numbers (first time only): 154 !----------------------------------------------------------- 155 IF(Counter == 1) THEN 156 CALL WriteNodes(file_id, PEs, MyPE, NofNodes, fill_type_id) 157 CALL WriteElements(file_id, PEs, MyPE, NofStorage) 158 CALL WriteParts(file_id, PEs, MyPE, NofNodes, fill_type_id) 159 END IF 160 161 ! Export field variables: 162 !------------------------- 163 DO i = 1, NofScalarFields 164 CALL WriteScalars(file_id, PEs, MyPE, NofNodes, ScalarFieldNames(i), fill_type_id) 165 END DO 166 167 DO i = 1, NofVectorFields 168 CALL WriteVectors(file_id, PEs, MyPE, NofNodes, VectorFieldNames(i), fill_type_id) 169 END DO 170 171 ! Rewrite the xdmf-file: 172 !------------------------ 173 IF(MyPE == 1) CALL WriteXdmfFile(PEs, NofNodes, NofElements, & 174 NofStorage, NofScalarFields, ScalarFieldNames, NofVectorFields, & 175 VectorFieldNames, BaseFileName, fill_type_id) 176 177 ! Finalize: 178 !----------- 179 CALL h5fclose_f(file_id, ierr) 180 DEALLOCATE(NofElements, NofNodes, NofStorage) 181 TotalTime = RealTime() - StartTime 182 WRITE(Str, *) TotalTime 183 CALL INFO('XdmfWriter', 'Total write time (REAL): '//TRIM(ADJUSTL(Str))) 184 Counter = Counter + 1 185 186CONTAINS 187 188!------------------------------------------------------------------------------ 189 SUBROUTINE FindScalarFields(NofScalarFields, ScalarFieldNames) 190!------------------------------------------------------------------------------ 191 IMPLICIT NONE 192 INTEGER, INTENT(OUT) :: NofScalarFields 193 CHARACTER(LEN=MAX_NAME_LEN), INTENT(OUT) :: ScalarFieldNames(:) 194 195 INTEGER :: id 196 LOGICAL :: Found 197 CHARACTER(LEN=MAX_NAME_LEN) :: LHS, RHS, Tmp1 198 TYPE(Variable_t), POINTER :: Variable 199 200 NofScalarFields = 0 201 202 DO id = 1, SIZE(ScalarFieldNames) 203 WRITE(Tmp1, *) id 204 205 WRITE(LHS, '(A)') 'scalar field '//TRIM(ADJUSTL(Tmp1)) 206 207 RHS = ListGetString(Solver % Values, TRIM(LHS), Found) 208 209 IF(.NOT.Found) CYCLE 210 211 Variable => VariableGet(Solver % Mesh % Variables, TRIM(RHS)) 212 213 IF(.NOT.ASSOCIATED(Variable)) THEN 214 CALL INFO('XdmfWriter', 'Bad scalar field: '//TRIM(RHS)) 215 CYCLE 216 END IF 217 218 NofScalarFields = NofScalarFields + 1 219 ScalarFieldNames(NofScalarFields) = TRIM(RHS) 220 END DO 221!------------------------------------------------------------------------------ 222 END SUBROUTINE FindScalarFields 223!------------------------------------------------------------------------------ 224 225!------------------------------------------------------------------------------ 226 SUBROUTINE FindVectorFields(NofVectorFields, VectorFieldNames) 227!------------------------------------------------------------------------------ 228 IMPLICIT NONE 229 INTEGER, INTENT(OUT) :: NofVectorFields 230 CHARACTER(LEN=MAX_NAME_LEN), INTENT(OUT) :: VectorFieldNames(:) 231 232 INTEGER :: id 233 LOGICAL :: Found 234 CHARACTER(LEN=MAX_NAME_LEN) :: LHS, RHS, Tmp1 235 TYPE(Variable_t), POINTER :: Variable 236 237 NofVectorFields = 0 238 239 DO id = 1, SIZE(VectorFieldNames) 240 WRITE(Tmp1, *) id 241 242 WRITE(LHS, '(A)') 'vector field '//TRIM(ADJUSTL(Tmp1)) 243 244 RHS = ListGetString(Solver % Values, TRIM(LHS), Found) 245 246 IF(.NOT.Found) CYCLE 247 248 Variable => VariableGet(Solver % Mesh % Variables, TRIM(RHS)) 249 250 IF(.NOT.ASSOCIATED(Variable)) THEN 251 ! Try componentwise: 252 !-------------------- 253 WRITE(Tmp1, '(A)') TRIM(RHS)//' 1' 254 Variable => VariableGet(Solver % Mesh % Variables, TRIM(Tmp1)) 255 IF(.NOT.ASSOCIATED(Variable)) THEN 256 CALL INFO('XdmfWriter', 'Bad vector field: '//TRIM(RHS)) 257 CYCLE 258 END IF 259 END IF 260 261 NofVectorFields = NofVectorFields + 1 262 VectorFieldNames(NofVectorFields) = TRIM(RHS) 263 END DO 264!------------------------------------------------------------------------------ 265 END SUBROUTINE FindVectorFields 266!------------------------------------------------------------------------------ 267 268!------------------------------------------------------------------------------ 269 INTEGER FUNCTION GetXdmfCode(ElementCode) RESULT(XdmfCode) 270!------------------------------------------------------------------------------ 271 IMPLICIT NONE 272 INTEGER :: ElementCode, XdmfCode 273 274 SELECT CASE(ElementCode) 275 CASE(202) 276 XdmfCode = 2 ! polyline 277 CASE(303) 278 XdmfCode = 4 ! linear triangle 279 CASE(404) 280 XdmfCode = 5 ! linear quadrilateral 281 CASE(504) 282 XdmfCode = 6 ! linear tetrahedron 283 CASE(510) 284 XdmfCode = 38 ! quadratic tetrahedron 285 CASE(808) 286 XdmfCode = 9 ! linear hexahedron 287 CASE(820) 288 XdmfCode = 48 ! quadratic hexahedron 289 CASE DEFAULT 290 XdmfCode = -1 ! not supported, yet 291 END SELECT 292 293!XDMF_NOTOPOLOGY 0x0 294!XDMF_POLYVERTEX 0x1 295!XDMF_POLYLINE 0x2 296!XDMF_POLYGON 0x3 297!XDMF_TRI 0x4 298!XDMF_QUAD 0x5 299!XDMF_TET 0x6 300!XDMF_PYRAMID 0x7 301!XDMF_WEDGE 0x8 302!XDMF_HEX 0x9 303!XDMF_EDGE_3 0x0022 304!XDMF_TRI_6 0x0024 305!XDMF_QUAD_8 0x0025 306!XDMF_TET_10 0x0026 307!XDMF_PYRAMID_13 0x0027 308!XDMF_WEDGE_15 0x0028 309!XDMF_WEDGE_18 0x0029 310!XDMF_HEX_20 0x0030 311!XDMF_HEX_24 0x0031 312!XDMF_HEX_27 0x0032 313!XDMF_MIXED 0x0070 314!XDMF_2DSMESH 0x0100 315!XDMF_2DRECTMESH 0x0101 316!XDMF_2DCORECTMESH 0x0102 317!XDMF_3DSMESH 0x1100 318!XDMF_3DRECTMESH 0x1101 319!XDMF_3DCORECTMESH 0x1102 320 321!------------------------------------------------------------------------------ 322 END FUNCTION GetXdmfCode 323!------------------------------------------------------------------------------ 324 325!------------------------------------------------------------------------------ 326 INTEGER FUNCTION GetElementStorageSize() RESULT(StorageSize) 327!------------------------------------------------------------------------------ 328 INTEGER :: i, StorageSize 329 TYPE(Element_t), POINTER :: Element 330 INTEGER :: XdmfCode 331 332 StorageSize = 0 333 DO i = 1, GetNofActive() 334 Element => GetActiveElement(i) 335 IF(.NOT.ASSOCIATED(Element)) CYCLE 336 XdmfCode = GetXdmfCode(Element % Type % ElementCode) 337 IF(XdmfCode < 0) CYCLE ! unknown: skip this element 338 StorageSize = StorageSize + 1 339 IF(XdmfCode == 2) StorageSize = StorageSize + 1 ! polyline 340 StorageSize = StorageSize + GetElementNofNodes() 341 END DO 342 343!------------------------------------------------------------------------------ 344 END FUNCTION GetElementStorageSize 345!------------------------------------------------------------------------------ 346 347!------------------------------------------------------------------------------ 348 SUBROUTINE WriteNodes(file_id, PEs, MyPE, NofNodes, fill_type_id) 349!------------------------------------------------------------------------------ 350 IMPLICIT NONE 351 INTEGER :: file_id, PEs, MyPE, NofNodes(:) 352 INTEGER(HID_T) :: fill_type_id 353 354 INTEGER :: i, ierr 355 INTEGER(HSIZE_T) :: dims(2) 356 INTEGER(HID_T) :: dset_id(PEs), filespace, memspace, plist_id 357 REAL(KIND=dp), ALLOCATABLE :: data(:,:) 358 CHARACTER(LEN=MAX_NAME_LEN) :: Str 359 360 ! Create the datasets collectively: 361 !----------------------------------- 362 DO i = 1, PEs 363 dims(1) = 3 364 dims(2) = NofNodes(i) 365 366 WRITE(Str, *) i 367 WRITE(Str, '(A)') 'nodes_'//TRIM(ADJUSTL(Str)) 368 369 CALL h5screate_simple_f(2, dims, filespace, ierr) 370 CALL h5dcreate_f(file_id, TRIM(ADJUSTL(Str)), fill_type_id, filespace, dset_id(i), ierr) 371 CALL h5sclose_f(filespace, ierr) 372 END DO 373 374 ! Write the data independently: 375 !------------------------------- 376 ALLOCATE(data(3, NofNodes(MyPE))) 377 378 dims(1) = SIZE(data, 1) 379 dims(2) = SIZE(data, 2) 380 381 DO i = 1, dims(2) 382 data(1, i) = Mesh % Nodes % x(i) 383 data(2, i) = Mesh % Nodes % y(i) 384 data(3, i) = Mesh % Nodes % z(i) 385 END DO 386 387 CALL h5screate_simple_f(2, dims, memspace, ierr) 388 CALL h5dget_space_f(dset_id(MyPE), filespace, ierr) 389 390 CALL h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr) 391 CALL h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_INDEPENDENT_F, ierr) 392 393 CALL h5dwrite_f(dset_id(MyPE), H5T_NATIVE_DOUBLE, data, dims, ierr, & 394 file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id) 395 396 ! Finalize: 397 !----------- 398 DEALLOCATE(data) 399 CALL h5pclose_f(plist_id, ierr) 400 CALL h5sclose_f(filespace, ierr) 401 CALL h5sclose_f(memspace, ierr) 402 403 DO i = 1, PEs 404 CALL h5dclose_f(dset_id(i), ierr) 405 END DO 406 407!------------------------------------------------------------------------------ 408 END SUBROUTINE WriteNodes 409!------------------------------------------------------------------------------ 410 411!------------------------------------------------------------------------------ 412 SUBROUTINE WriteElements(file_id, PEs, MyPE, NofStorage) 413!------------------------------------------------------------------------------ 414 IMPLICIT NONE 415 INTEGER :: file_id, PEs, MyPE, NofStorage(:) 416 417 INTEGER :: i, j, k, ierr, XdmfCode 418 INTEGER(HSIZE_T) :: dims(2) 419 INTEGER(HID_T) :: dset_id(PEs), filespace, memspace, plist_id 420 INTEGER, ALLOCATABLE :: data(:,:) 421 CHARACTER(LEN=MAX_NAME_LEN) :: Str 422 TYPE(Element_t), POINTER :: Element 423 424 ! Create the datasets collectively: 425 !----------------------------------- 426 DO i = 1, PEs 427 dims(1) = 1 428 dims(2) = NofStorage(i) 429 430 WRITE(Str, *) i 431 WRITE(Str, '(A)') 'elements_'//TRIM(ADJUSTL(Str)) 432 433 CALL h5screate_simple_f(2, dims, filespace, ierr) 434 CALL h5dcreate_f(file_id, TRIM(ADJUSTL(Str)), H5T_NATIVE_INTEGER, filespace, dset_id(i), ierr) 435 CALL h5sclose_f(filespace, ierr) 436 END DO 437 438 ! Write the data independently: 439 !------------------------------- 440 ALLOCATE(data(1, NofStorage(MyPE))) 441 442 dims(1) = SIZE(data, 1) 443 dims(2) = SIZE(data, 2) 444 445 j = 0 446 DO i = 1, GetNofActive() 447 Element => GetActiveElement(i) 448 IF(.NOT.ASSOCIATED(Element)) CYCLE 449 XdmfCode = GetXdmfCode(Element % Type % ElementCode) 450 IF(XdmfCode < 0) CYCLE ! unknown: skip this element 451 452 j = j + 1 453 data(1, j) = XdmfCode 454 455 IF(XdmfCode == 2) THEN 456 j = j + 1 ! polyline: nof nodes 457 data(1, j) = GetElementNofNodes() 458 END IF 459 460 DO k = 1, GetElementNofNodes() 461 j = j + 1 462 463 ! Permuted C-style numbering 464 SELECT CASE(Element % Type % ElementCode) 465 CASE(203) 466 data(1, j) = Element % NodeIndexes(Order203(k)) - 1 467 CASE(510) 468 data(1, j) = Element % NodeIndexes(Order510(k)) - 1 469 CASE(820) 470 data(1, j) = Element % NodeIndexes(Order820(k)) - 1 471 CASE DEFAULT 472 data(1, j) = Element % NodeIndexes(k) - 1 473 END SELECT 474 475 END DO 476 END DO 477 478 IF(j /= NofStorage(MyPE)) CALL Fatal('XdmfWriter', 'Bad element numbering') 479 480 CALL h5screate_simple_f(2, dims, memspace, ierr) 481 CALL h5dget_space_f(dset_id(MyPE), filespace, ierr) 482 483 CALL h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr) 484 CALL h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_INDEPENDENT_F, ierr) 485 486 CALL h5dwrite_f(dset_id(MyPE), H5T_NATIVE_INTEGER, data, dims, ierr, & 487 file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id) 488 489 ! Finalize: 490 !----------- 491 DEALLOCATE(data) 492 CALL h5pclose_f(plist_id, ierr) 493 CALL h5sclose_f(filespace, ierr) 494 CALL h5sclose_f(memspace, ierr) 495 496 DO i = 1, PEs 497 CALL h5dclose_f(dset_id(i), ierr) 498 END DO 499 500!------------------------------------------------------------------------------ 501 END SUBROUTINE WriteElements 502!------------------------------------------------------------------------------ 503 504!------------------------------------------------------------------------------ 505 SUBROUTINE WriteParts(file_id, PEs, MyPE, NofNodes, fill_type_id) 506!------------------------------------------------------------------------------ 507 IMPLICIT NONE 508 INTEGER :: file_id, PEs, MyPE, NofNodes(:) 509 INTEGER(HID_T) :: fill_type_id 510 511 INTEGER :: i, ierr 512 INTEGER(HSIZE_T) :: dims(2) 513 INTEGER(HID_T) :: dset_id(PEs), filespace, memspace, plist_id 514 REAL(KIND=dp), ALLOCATABLE :: data(:,:) 515 CHARACTER(LEN=MAX_NAME_LEN) :: Str 516 517 ! Create the datasets collectively: 518 !----------------------------------- 519 DO i = 1, PEs 520 dims(1) = 1 521 dims(2) = NofNodes(i) 522 523 WRITE(Str, *) i 524 WRITE(Str, '(A)') 'part_number_'//TRIM(ADJUSTL(Str)) 525 526 CALL h5screate_simple_f(2, dims, filespace, ierr) 527 CALL h5dcreate_f(file_id, TRIM(ADJUSTL(Str)), fill_type_id, filespace, dset_id(i), ierr) 528 CALL h5sclose_f(filespace, ierr) 529 END DO 530 531 ! Write the data independently: 532 !------------------------------- 533 ALLOCATE(data(1, NofNodes(MyPE))) 534 535 dims(1) = SIZE(data, 1) 536 dims(2) = SIZE(data, 2) 537 538 data = MyPE 539 540 CALL h5screate_simple_f(2, dims, memspace, ierr) 541 CALL h5dget_space_f(dset_id(MyPE), filespace, ierr) 542 543 CALL h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr) 544 CALL h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_INDEPENDENT_F, ierr) 545 546 CALL h5dwrite_f(dset_id(MyPE), H5T_NATIVE_DOUBLE, data, dims, ierr, & 547 file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id) 548 549 ! Finalize: 550 !----------- 551 DEALLOCATE(data) 552 CALL h5pclose_f(plist_id, ierr) 553 CALL h5sclose_f(filespace, ierr) 554 CALL h5sclose_f(memspace, ierr) 555 556 DO i = 1, PEs 557 CALL h5dclose_f(dset_id(i), ierr) 558 END DO 559 560!------------------------------------------------------------------------------ 561 END SUBROUTINE WriteParts 562!------------------------------------------------------------------------------ 563 564!------------------------------------------------------------------------------ 565 SUBROUTINE WriteScalars(file_id, PEs, MyPE, NofNodes, ScalarFieldName, fill_type_id) 566!------------------------------------------------------------------------------ 567 IMPLICIT NONE 568 INTEGER :: file_id, PEs, MyPE, NofNodes(:) 569 CHARACTER(LEN=MAX_NAME_LEN) :: ScalarFieldName 570 INTEGER(HID_T) :: fill_type_id 571 572 INTEGER :: i, j, ierr 573 INTEGER(HSIZE_T) :: dims(2) 574 INTEGER(HID_T) :: dset_id(PEs), filespace, memspace, plist_id 575 REAL(KIND=dp), ALLOCATABLE :: data(:,:) 576 CHARACTER(LEN=MAX_NAME_LEN) :: Str, Tmp1, Tmp2 577 TYPE(Variable_t), POINTER :: Var 578 579 ! Create the datasets collectively: 580 !----------------------------------- 581 DO i = 1, PEs 582 dims(1) = 1 583 dims(2) = NofNodes(i) 584 585 WRITE(Tmp1, *) i 586 WRITE(Tmp2, *) Counter 587 588 WRITE(Str, '(A)') TRIM(ScalarFieldName)//'_'//TRIM(ADJUSTL(Tmp2))//'_'//TRIM(ADJUSTL(Tmp1)) 589 590 CALL h5screate_simple_f(2, dims, filespace, ierr) 591 CALL h5dcreate_f(file_id, TRIM(ADJUSTL(Str)), fill_type_id, filespace, dset_id(i), ierr) 592 CALL h5sclose_f(filespace, ierr) 593 END DO 594 595 ! Write the data independently: 596 !------------------------------- 597 ALLOCATE(data(1, NofNodes(MyPE))) 598 599 dims(1) = SIZE(data, 1) 600 dims(2) = SIZE(data, 2) 601 602 Var => VariableGet(Solver % Mesh % Variables, TRIM(ScalarFieldName)) 603 IF(.NOT.ASSOCIATED(Var)) CALL INFO('XdmfWriter', 'Scalar not found') 604 605 DO i = 1, dims(2) 606 j = Var % Perm(i) 607 data(1, i) = Var % Values(j) 608 END DO 609 610 CALL h5screate_simple_f(2, dims, memspace, ierr) 611 CALL h5dget_space_f(dset_id(MyPE), filespace, ierr) 612 613 CALL h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr) 614 CALL h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_INDEPENDENT_F, ierr) 615 616 CALL h5dwrite_f(dset_id(MyPE), H5T_NATIVE_DOUBLE, data, dims, ierr, & 617 file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id) 618 619 ! Finalize: 620 !----------- 621 DEALLOCATE(data) 622 CALL h5pclose_f(plist_id, ierr) 623 CALL h5sclose_f(filespace, ierr) 624 CALL h5sclose_f(memspace, ierr) 625 626 DO i = 1, PEs 627 CALL h5dclose_f(dset_id(i), ierr) 628 END DO 629 630!------------------------------------------------------------------------------ 631 END SUBROUTINE WriteScalars 632!------------------------------------------------------------------------------ 633 634!------------------------------------------------------------------------------ 635 SUBROUTINE WriteVectors(file_id, PEs, MyPE, NofNodes, VectorFieldName, fill_type_id) 636!------------------------------------------------------------------------------ 637 IMPLICIT NONE 638 INTEGER :: file_id, PEs, MyPE, NofNodes(:) 639 CHARACTER(LEN=MAX_NAME_LEN) :: VectorFieldName 640 INTEGER(HID_T) :: fill_type_id 641 642 INTEGER :: i, j, k, ierr, dofs 643 INTEGER(HSIZE_T) :: dims(2) 644 INTEGER(HID_T) :: dset_id(PEs), filespace, memspace, plist_id 645 REAL(KIND=dp), ALLOCATABLE :: data(:,:) 646 CHARACTER(LEN=MAX_NAME_LEN) :: Str, Tmp1, Tmp2 647 TYPE(Variable_t), POINTER :: Var 648 649 ! Create the datasets collectively: 650 !----------------------------------- 651 DO i = 1, PEs 652 dims(1) = 3 653 dims(2) = NofNodes(i) 654 655 WRITE(Tmp1, *) i 656 WRITE(Tmp2, *) Counter 657 658 WRITE(Str, '(A)') TRIM(VectorFieldName)//'_'//TRIM(ADJUSTL(Tmp2))//'_'//TRIM(ADJUSTL(Tmp1)) 659 660 CALL h5screate_simple_f(2, dims, filespace, ierr) 661 CALL h5dcreate_f(file_id, TRIM(ADJUSTL(Str)), fill_type_id, filespace, dset_id(i), ierr) 662 CALL h5sclose_f(filespace, ierr) 663 END DO 664 665 ! Write the data independently: 666 !------------------------------- 667 ALLOCATE(data(3, NofNodes(MyPE))) 668 669 dims(1) = SIZE(data, 1) 670 dims(2) = SIZE(data, 2) 671 672 data = 0.0d0 673 674 Var => VariableGet(Solver % Mesh % Variables, TRIM(VectorFieldName)) 675 676 IF(ASSOCIATED(Var)) THEN 677 ! Storage type: multiple DOFs per node: 678 !--------------------------------------- 679 dofs = Var % DOFs 680 DO i = 1, dims(2) 681 j = Var % Perm(i) 682 DO k = 1, MIN(3, dofs) 683 data(k, i) = Var % Values(dofs*(j - 1) + k) 684 END DO 685 END DO 686 ELSE 687 ! Storage type: multiple scalars per node: 688 !----------------------------------------- 689 DO k = 1, 3 690 WRITE(Tmp1, *) k 691 WRITE(Tmp2, '(A)') TRIM(VectorFieldName)//' '//TRIM(ADJUSTL(Tmp1)) 692 Var => VariableGet(Solver % Mesh % Variables, TRIM(Tmp2)) 693 IF(.NOT.ASSOCIATED(Var)) CYCLE 694 dofs = Var % DOFs 695 IF(dofs /= 1) CYCLE 696 DO i = 1, dims(2) 697 j = Var % Perm(i) 698 data(k, i) = Var % Values(j) 699 END DO 700 END DO 701 END IF 702 703 CALL h5screate_simple_f(2, dims, memspace, ierr) 704 CALL h5dget_space_f(dset_id(MyPE), filespace, ierr) 705 706 CALL h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr) 707 CALL h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_INDEPENDENT_F, ierr) 708 709 CALL h5dwrite_f(dset_id(MyPE), H5T_NATIVE_DOUBLE, data, dims, ierr, & 710 file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id) 711 712 ! Finalize: 713 !----------- 714 DEALLOCATE(data) 715 CALL h5pclose_f(plist_id, ierr) 716 CALL h5sclose_f(filespace, ierr) 717 CALL h5sclose_f(memspace, ierr) 718 719 DO i = 1, PEs 720 CALL h5dclose_f(dset_id(i), ierr) 721 END DO 722 723!------------------------------------------------------------------------------ 724 END SUBROUTINE WriteVectors 725!------------------------------------------------------------------------------ 726 727!------------------------------------------------------------------------------ 728 LOGICAL FUNCTION Dmp(fid, indent, str) RESULT(ok) 729!------------------------------------------------------------------------------ 730 IMPLICIT NONE 731 INTEGER :: fid, indent 732 CHARACTER(LEN=*) :: str 733 LOGICAL :: ok 734 735 INTEGER :: i 736 CHARACTER(LEN=MAX_NAME_LEN) :: line 737 738 WRITE(line, '(A)') REPEAT(' ', indent)//TRIM(ADJUSTL(str))//CHAR(10) 739 740 WRITE(fid) TRIM(line) 741 742 ok = .TRUE. ! TODO: Return false if WRITE fails 743!------------------------------------------------------------------------------ 744 END FUNCTION Dmp 745!------------------------------------------------------------------------------ 746 747!------------------------------------------------------------------------------ 748 SUBROUTINE WriteXdmfFile(PEs, NofNodes, NofElements, & 749 NofStorage, NofScalarFields, ScalarFieldNames, & 750 NofVectorFields, VectorFieldNames, BaseFileName, & 751 fill_type_id) 752!------------------------------------------------------------------------------ 753 IMPLICIT NONE 754 INTEGER :: PEs, NofScalarFields, NofVectorFields 755 INTEGER :: NofElements(:), NofNodes(:), NofStorage(:) 756 CHARACTER(LEN=MAX_NAME_LEN) :: ScalarFieldNames(:), VectorFieldNames(:) 757 CHARACTER(LEN=MAX_NAME_LEN) :: BaseFileName 758 INTEGER(HID_T) :: fill_type_id 759 760 CHARACTER(LEN=MAX_NAME_LEN) :: Tmp1, Tmp2, Tmp3, Tmp4, Tmp5, Tmp6 761 CHARACTER(LEN=MAX_NAME_LEN) :: FileName 762 CHARACTER(LEN=MAX_NAME_LEN) :: H5FileName 763 INTEGER :: i, j, k 764 LOGICAL :: ok 765 766 ! Initialize: 767 !------------- 768 FileName = TRIM(ADJUSTL(BaseFileName))//'.xmf' 769 H5FileName = TRIM(ADJUSTL(BaseFileName))//'.h5' 770 771 OPEN(UNIT=10, FILE=TRIM(FileName), FORM='unformatted', ACCESS='stream', STATUS='unknown') 772 773 ok = Dmp(10, 0, '<?xml version="1.0" ?>') 774 ok = Dmp(10, 0, '<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>') 775 ok = Dmp(10, 0, '') 776 ok = Dmp(10, 0, '<Xdmf>') 777 ok = Dmp(10, 2, '<Domain>') 778 779 ok = Dmp(10, 4, '<Grid Name="mesh" GridType="Collection" CollectionType="Temporal">') 780 781 DO j = 1, Counter 782 WRITE(Tmp1, *) j; WRITE(Tmp1, '(A)') TRIM(ADJUSTL(Tmp1)) 783 WRITE(Tmp2, *) TimeValues(j); WRITE(Tmp2, '(A)') TRIM(ADJUSTL(Tmp2)) 784 785 ok = Dmp(10, 6, '<Grid Name="mesh_'//TRIM(Tmp1)//'" GridType="Collection" CollectionType="Spatial">') 786 ok = Dmp(10, 8, '<Time Value="'//TRIM(Tmp2)//'" />') 787 788 DO i = 1, PEs 789 WRITE(Tmp1, *) i; WRITE(Tmp1, '(A)') TRIM(ADJUSTL(Tmp1)) 790 WRITE(Tmp2, *) NofElements(i); WRITE(Tmp2, '(A)') TRIM(ADJUSTL(Tmp2)) 791 WRITE(Tmp3, *) NofStorage(i); WRITE(Tmp3, '(A)') TRIM(ADJUSTL(Tmp3)) 792 WRITE(Tmp4, *) NofNodes(i); WRITE(Tmp4, '(A)') TRIM(ADJUSTL(Tmp4)) 793 WRITE(Tmp5, *) j; WRITE(Tmp5, '(A)') TRIM(ADJUSTL(Tmp5)) 794 WRITE(Tmp6, *) TimeValues(j); WRITE(Tmp6, '(A)') TRIM(ADJUSTL(Tmp6)) 795 796 ! Init part: 797 !------------ 798 ok = Dmp(10, 8, '<Grid Name="mesh_'//TRIM(Tmp5)//'_'//TRIM(Tmp1)//'">') 799 ok = Dmp(10, 10, '<Time Value="'//TRIM(Tmp6)//'" />') 800 801 ! Write elements: 802 !----------------- 803 ok = Dmp(10, 10, '<Topology Type="Mixed" NumberOfElements="'//TRIM(Tmp2)//'">') 804 ok = Dmp(10, 12, '<DataItem Format="HDF" DataType="Int" Dimensions="'//TRIM(Tmp3)//'">') 805 ok = Dmp(10, 14, TRIM(H5FileName)//':/elements_'//TRIM(Tmp1)) 806 ok = Dmp(10, 12, '</DataItem>') 807 ok = Dmp(10, 10, '</Topology>') 808 809 ! Write nodes: 810 !-------------- 811 ok = Dmp(10, 10, '<Geometry Type="XYZ">') 812 IF(fill_type_id == H5T_NATIVE_DOUBLE) THEN 813 ok = Dmp(10, 12, '<DataItem Format="HDF" DataType="Float" Precision="8" Dimensions="'//TRIM(Tmp4)//' 3">') 814 ELSE 815 ok = Dmp(10, 12, '<DataItem Format="HDF" DataType="Float" Precision="4" Dimensions="'//TRIM(Tmp4)//' 3">') 816 END IF 817 ok = Dmp(10, 14, TRIM(H5FileName)//':/nodes_'//TRIM(Tmp1)) 818 ok = Dmp(10, 12, '</DataItem>') 819 ok = Dmp(10, 10, '</Geometry>') 820 821 ! Write part number: 822 !-------------------- 823 ok = Dmp(10, 10, ' <Attribute Name="part_number" AttributeType="Scalar" Center="Node">') 824 IF(fill_type_id == H5T_NATIVE_DOUBLE) THEN 825 ok = Dmp(10, 12, ' <DataItem Format="HDF" DataType="Float" Precision="8" Dimensions="'//TRIM(Tmp4)//' 1">') 826 ELSE 827 ok = Dmp(10, 12, ' <DataItem Format="HDF" DataType="Float" Precision="4" Dimensions="'//TRIM(Tmp4)//' 1">') 828 END IF 829 ok = Dmp(10, 14, TRIM(H5FileName)//':/part_number_'//TRIM(Tmp1)) 830 ok = Dmp(10, 12, '</DataItem>') 831 ok = Dmp(10, 10, '</Attribute>') 832 833 ! Write scalar fields: 834 !---------------------- 835 DO k = 1, NofScalarFields 836 ok = Dmp(10, 10, ' <Attribute Name="'//TRIM(ScalarFieldNames(k))//'" AttributeType="Scalar" Center="Node">') 837 IF(fill_type_id == H5T_NATIVE_DOUBLE) THEN 838 ok = Dmp(10, 12, ' <DataItem Format="HDF" DataType="Float" Precision="8" Dimensions="'//TRIM(Tmp4)//' 1">') 839 ELSE 840 ok = Dmp(10, 12, ' <DataItem Format="HDF" DataType="Float" Precision="4" Dimensions="'//TRIM(Tmp4)//' 1">') 841 END IF 842 ok = Dmp(10, 14, TRIM(H5FileName)//':/'//TRIM(ScalarFieldNames(k))//'_'//TRIM(Tmp5)//'_'//TRIM(Tmp1)) 843 ok = Dmp(10, 12, '</DataItem>') 844 ok = Dmp(10, 10, '</Attribute>') 845 END DO 846 847 ! Write vector fields: 848 !---------------------- 849 DO k = 1, NofVectorFields 850 ok = Dmp(10, 10, ' <Attribute Name="'//TRIM(VectorFieldNames(k))//'" AttributeType="Vector" Center="Node">') 851 IF(fill_type_id == H5T_NATIVE_DOUBLE) THEN 852 ok = Dmp(10, 12, ' <DataItem Format="HDF" DataType="Float" Precision="8" Dimensions="'//TRIM(Tmp4)//' 3">') 853 ELSE 854 ok = Dmp(10, 12, ' <DataItem Format="HDF" DataType="Float" Precision="4" Dimensions="'//TRIM(Tmp4)//' 3">') 855 END IF 856 ok = Dmp(10, 14, TRIM(H5FileName)//':/'//TRIM(VectorFieldNames(k))//'_'//TRIM(Tmp5)//'_'//TRIM(Tmp1)) 857 ok = Dmp(10, 12, '</DataItem>') 858 ok = Dmp(10, 10, '</Attribute>') 859 END DO 860 861 ! Finalize part: 862 !---------------- 863 ok = Dmp(10, 8, '</Grid>') ! part 864 END DO 865 866 ok = Dmp(10, 6, '</Grid>') ! spatial collection 867 868 END DO 869 870 ok = Dmp(10, 4, '</Grid>') ! temporal collection 871 ok = Dmp(10, 2, '</Domain>') 872 ok = Dmp(10, 0, '</Xdmf>') 873 874 CLOSE(10) 875 876!------------------------------------------------------------------------------ 877 END SUBROUTINE WriteXdmfFile 878!------------------------------------------------------------------------------ 879 880!------------------------------------------------------------------------------ 881END SUBROUTINE XdmfWriter 882!------------------------------------------------------------------------------ 883