diff --git a/.github/workflows/CI_parallel.yml b/.github/workflows/CI_parallel.yml index ec47eb797..d79d98cfe 100644 --- a/.github/workflows/CI_parallel.yml +++ b/.github/workflows/CI_parallel.yml @@ -782,7 +782,7 @@ jobs: if: '!cancelled()' - name: Run Cylinder_pAdaptationRL - working-directory: ./Solver/test/NavierStokes/Cylinder_pAdaptationRL/SETUP + working-directory: ./Solver/test/NavierStokes/Cylinder_pAdaptationRL run: | source /opt/intel/oneapi/setvars.sh || true mpiexec -n 8 ./horses3d.ns Cylinder_pAdaptationRL.control diff --git a/Solver/src/libs/mesh/HexMesh.f90 b/Solver/src/libs/mesh/HexMesh.f90 index cbae4a2ab..933faef98 100644 --- a/Solver/src/libs/mesh/HexMesh.f90 +++ b/Solver/src/libs/mesh/HexMesh.f90 @@ -94,6 +94,7 @@ MODULE HexMeshClass procedure :: ExportBoundaryMesh => HexMesh_ExportBoundaryMesh procedure :: SaveSolution => HexMesh_SaveSolution procedure :: pAdapt => HexMesh_pAdapt + procedure :: pAdapt_MPI => HexMesh_pAdapt_MPI #if defined(NAVIERSTOKES) procedure :: SaveStatistics => HexMesh_SaveStatistics procedure :: ResetStatistics => HexMesh_ResetStatistics @@ -4116,7 +4117,7 @@ end subroutine HexMesh_ConvertPhaseFieldToDensity ! -------------------------------------------------------- ! Adapts a mesh to new polynomial orders NNew ! -------------------------------------------------------- -subroutine HexMesh_pAdapt (self, NNew, controlVariables) +subroutine HexMesh_pAdapt_MPI (self, NNew, controlVariables) implicit none !-arguments----------------------------------------- class(HexMesh), target , intent(inout) :: self @@ -4234,186 +4235,187 @@ subroutine HexMesh_pAdapt (self, NNew, controlVariables) nullify (e) nullify (f) -end subroutine HexMesh_pAdapt +end subroutine HexMesh_pAdapt_MPI ! !/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ------------------------------------------------------------------------ ! Adapts a mesh to new polynomial orders NNew. -! Optimized version of HexMesh_Adapt, but this one does not work with MPI +! Optimized version of HexMesh_Adapt_MPI, but this one does not work with MPI +! This subroutine is required for pAdaptationClassTE ! ------------------------------------------------------------------------ -! subroutine HexMesh_pAdapt_optimized (self, NNew, controlVariables) -! implicit none -! !-arguments----------------------------------------- -! class(HexMesh), target , intent(inout) :: self -! integer , intent(in) :: NNew(NDIM,self % no_of_elements) -! type(FTValueDictionary) , intent(in) :: controlVariables -! !-local-variables----------------------------------- -! integer :: eID, fID, zoneID -! logical :: saveGradients, FaceComputeQdot -! logical :: analyticalJac ! Do we need analytical Jacobian storage? -! type(IntegerDataLinkedList_t) :: elementList -! type(IntegerDataLinkedList_t) :: facesList -! type(IntegerDataLinkedList_t) :: zoneList -! integer , allocatable :: zoneArray(:) -! integer , allocatable :: facesArray(:) -! integer , allocatable :: elementArray(:) -! type(Zone_t) , pointer :: zone -! type(Element) , pointer :: e -! type(Face) , pointer :: f -! #if (!defined(NAVIERSTOKES)) -! logical, parameter :: computeGradients = .true. -! #endif -! !--------------------------------------------------- - -! ! ************************************** -! ! Check if resulting mesh is anisotropic -! ! ************************************** -! if ( maxval(NNew) /= minval(NNew) ) self % anisotropic = .TRUE. - -! self % NDOF = 0 -! do eID=1, self % no_of_elements -! self % NDOF = self % NDOF + product( NNew(:,eID) + 1 ) -! end do - -! ! ******************** -! ! Some initializations -! ! ******************** -! saveGradients = controlVariables % logicalValueForKey("save gradients with solution") -! FaceComputeQdot = controlVariables % containsKey("acoustic analogy") - -! facesList = IntegerDataLinkedList_t(.FALSE.) -! elementList = IntegerDataLinkedList_t(.FALSE.) -! zoneList = IntegerDataLinkedList_t(.FALSE.) -! analyticalJac = self % storage % anJacobian - -! ! ********************************************* -! ! Adapt individual elements (geometry excluded) -! ! ********************************************* -! !$omp parallel do schedule(runtime) private(fID, e) -! do eID=1, self % no_of_elements -! e => self % elements(eID) ! Associate fails(!) here -! if ( all( e % Nxyz == NNew(:,eID)) ) then -! cycle -! else -! call e % pAdapt ( NNew(:,eID), self % nodeType, saveGradients, self % storage % prevSol_num ) -! !$omp critical -! self % Nx(eID) = NNew(1,eID) -! self % Ny(eID) = NNew(2,eID) -! self % Nz(eID) = NNew(3,eID) -! call elementList % add (eID) -! do fID=1, 6 -! call facesList % add (e % faceIDs(fID)) -! if (self % faces(e % faceIDs(fID)) % FaceType /= HMESH_BOUNDARY) then -! call elementList % add ( mpi_partition % global2localeID (e % Connection(fID) % globID) ) -! end if -! end do -! !$omp end critical - -! end if -! !~ end associate -! end do -! !$omp end parallel do - -! call facesList % ExportToArray(facesArray, .TRUE.) - -! ! ************************* -! ! Adapt corresponding faces -! ! ************************* - -! ! Destruct faces storage -! ! ---------------------- -! !$omp parallel do schedule(runtime) -! do fID=1, size(facesArray) -! call self % faces( facesArray(fID) ) % storage % destruct -! end do -! !$omp end parallel do - -! ! Set connectivities -! ! ------------------ -! call self % SetConnectivitiesAndLinkFaces (self % nodeType, facesArray) - -! ! Construct faces storage -! ! ----------------------- -! !$omp parallel do private(f) schedule(runtime) -! do fID=1, size(facesArray) -! f => self % faces( facesArray(fID) ) ! associate fails here in intel compilers -! call f % storage(1) % Construct(NDIM, f % Nf, f % NelLeft , computeGradients, analyticalJac, FaceComputeQdot) -! call f % storage(2) % Construct(NDIM, f % Nf, f % NelRight, computeGradients, analyticalJac, FaceComputeQdot) -! end do -! !$omp end parallel do - -! ! ******************** -! ! Reconstruct geometry -! ! ******************** - -! !* 1. Adapted elements -! !* 2. Surrounding faces of adapted elements -! !* 3. Neighbor elements of adapted elements whose intermediate face's geometry was adapted -! !* 4. Faces and elements that share a boundary with a reconstructed face (3D non-conforming representations) - - -! if (self % anisotropic .and. (.not. self % meshIs2D) ) then - -! ! Check if any of the faces belongs to a boundary -! ! ----------------------------------------------- -! do fID=1, size(facesArray) -! associate (f => self % faces( facesArray(fID) ) ) -! if ( f % FaceType == HMESH_BOUNDARY ) then -! call zoneList % add (f % zone) -! end if -! end associate -! end do - -! ! Add the corresponding faces and elements -! ! ---------------------------------------- -! call zoneList % ExportToArray (zoneArray) - -! do zoneID=1, size(zoneArray) -! zone => self % zones( zoneArray(zoneID) ) ! Compiler bug(?): If zone was implemented as associate, gfortran would not compile -! do fID=1, zone % no_of_faces -! call facesList % add ( zone % faces(fID) ) -! call elementList % add ( self % faces(zone % faces(fID)) % elementIDs(1) ) -! end do -! end do -! deallocate (zoneArray ) -! end if - -! deallocate ( facesArray ) - -! call facesList % ExportToArray(facesArray , .TRUE.) -! call elementList % ExportToArray(elementArray, .TRUE.) - -! ! Destruct old -! ! ------------ -! do eID=1, size (elementArray) -! call self % elements (elementArray(eID)) % geom % destruct -! end do -! do fID=1, size (facesArray) -! call self % faces (facesArray(fID)) % geom % destruct -! end do - -! call self % ConstructGeometry(facesArray, elementArray) - -! #if defined(NAVIERSTOKES) -! call self % ComputeWallDistances(facesArray, elementArray) -! #endif - -! ! ********* -! ! Finish up -! ! ********* -! call self % PrepareForIO - -! call facesList % destruct -! call elementList % destruct -! call zoneList % destruct -! nullify (zone) -! nullify (e) -! nullify (f) -! deallocate (facesArray ) -! deallocate (elementArray) - -! end subroutine HexMesh_pAdapt_optimized + subroutine HexMesh_pAdapt (self, NNew, controlVariables) + implicit none + !-arguments----------------------------------------- + class(HexMesh), target , intent(inout) :: self + integer , intent(in) :: NNew(NDIM,self % no_of_elements) + type(FTValueDictionary) , intent(in) :: controlVariables + !-local-variables----------------------------------- + integer :: eID, fID, zoneID + logical :: saveGradients, FaceComputeQdot + logical :: analyticalJac ! Do we need analytical Jacobian storage? + type(IntegerDataLinkedList_t) :: elementList + type(IntegerDataLinkedList_t) :: facesList + type(IntegerDataLinkedList_t) :: zoneList + integer , allocatable :: zoneArray(:) + integer , allocatable :: facesArray(:) + integer , allocatable :: elementArray(:) + type(Zone_t) , pointer :: zone + type(Element) , pointer :: e + type(Face) , pointer :: f +#if (!defined(NAVIERSTOKES)) + logical, parameter :: computeGradients = .true. +#endif + !--------------------------------------------------- + +! ************************************** +! Check if resulting mesh is anisotropic +! ************************************** + if ( maxval(NNew) /= minval(NNew) ) self % anisotropic = .TRUE. + + self % NDOF = 0 + do eID=1, self % no_of_elements + self % NDOF = self % NDOF + product( NNew(:,eID) + 1 ) + end do + +! ******************** +! Some initializations +! ******************** + saveGradients = controlVariables % logicalValueForKey("save gradients with solution") + FaceComputeQdot = controlVariables % containsKey("acoustic analogy") + + facesList = IntegerDataLinkedList_t(.FALSE.) + elementList = IntegerDataLinkedList_t(.FALSE.) + zoneList = IntegerDataLinkedList_t(.FALSE.) + analyticalJac = self % storage % anJacobian + +! ********************************************* +! Adapt individual elements (geometry excluded) +! ********************************************* +!$omp parallel do schedule(runtime) private(fID, e) + do eID=1, self % no_of_elements + e => self % elements(eID) ! Associate fails(!) here + if ( all( e % Nxyz == NNew(:,eID)) ) then + cycle + else + call e % pAdapt ( NNew(:,eID), self % nodeType, saveGradients, self % storage % prevSol_num ) +!$omp critical + self % Nx(eID) = NNew(1,eID) + self % Ny(eID) = NNew(2,eID) + self % Nz(eID) = NNew(3,eID) + call elementList % add (eID) + do fID=1, 6 + call facesList % add (e % faceIDs(fID)) + if (self % faces(e % faceIDs(fID)) % FaceType /= HMESH_BOUNDARY) then + call elementList % add ( mpi_partition % global2localeID (e % Connection(fID) % globID) ) + end if + end do +!$omp end critical + + end if +!~ end associate + end do +!$omp end parallel do + + call facesList % ExportToArray(facesArray, .TRUE.) + +! ************************* +! Adapt corresponding faces +! ************************* + +! Destruct faces storage +! ---------------------- +!$omp parallel do schedule(runtime) + do fID=1, size(facesArray) + call self % faces( facesArray(fID) ) % storage % destruct + end do +!$omp end parallel do + +! Set connectivities +! ------------------ + call self % SetConnectivitiesAndLinkFaces (self % nodeType, facesArray) + +! Construct faces storage +! ----------------------- +!$omp parallel do private(f) schedule(runtime) + do fID=1, size(facesArray) + f => self % faces( facesArray(fID) ) ! associate fails here in intel compilers + call f % storage(1) % Construct(NDIM, f % Nf, f % NelLeft , computeGradients, analyticalJac, FaceComputeQdot) + call f % storage(2) % Construct(NDIM, f % Nf, f % NelRight, computeGradients, analyticalJac, FaceComputeQdot) + end do +!$omp end parallel do + +! ******************** +! Reconstruct geometry +! ******************** + + !* 1. Adapted elements + !* 2. Surrounding faces of adapted elements + !* 3. Neighbor elements of adapted elements whose intermediate face's geometry was adapted + !* 4. Faces and elements that share a boundary with a reconstructed face (3D non-conforming representations) + + + if (self % anisotropic .and. (.not. self % meshIs2D) ) then + +! Check if any of the faces belongs to a boundary +! ----------------------------------------------- + do fID=1, size(facesArray) + associate (f => self % faces( facesArray(fID) ) ) + if ( f % FaceType == HMESH_BOUNDARY ) then + call zoneList % add (f % zone) + end if + end associate + end do + +! Add the corresponding faces and elements +! ---------------------------------------- + call zoneList % ExportToArray (zoneArray) + + do zoneID=1, size(zoneArray) + zone => self % zones( zoneArray(zoneID) ) ! Compiler bug(?): If zone was implemented as associate, gfortran would not compile + do fID=1, zone % no_of_faces + call facesList % add ( zone % faces(fID) ) + call elementList % add ( self % faces(zone % faces(fID)) % elementIDs(1) ) + end do + end do + deallocate (zoneArray ) + end if + + deallocate ( facesArray ) + + call facesList % ExportToArray(facesArray , .TRUE.) + call elementList % ExportToArray(elementArray, .TRUE.) + +! Destruct old +! ------------ + do eID=1, size (elementArray) + call self % elements (elementArray(eID)) % geom % destruct + end do + do fID=1, size (facesArray) + call self % faces (facesArray(fID)) % geom % destruct + end do + + call self % ConstructGeometry(facesArray, elementArray) + +#if defined(NAVIERSTOKES) + call self % ComputeWallDistances(facesArray, elementArray) +#endif + +! ********* +! Finish up +! ********* + call self % PrepareForIO + + call facesList % destruct + call elementList % destruct + call zoneList % destruct + nullify (zone) + nullify (e) + nullify (f) + deallocate (facesArray ) + deallocate (elementArray) + + end subroutine HexMesh_pAdapt ! !/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! diff --git a/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 b/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 index 2c08364e0..0c2256594 100644 --- a/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 +++ b/Solver/src/libs/timeintegrator/pAdaptationClassRL.f90 @@ -423,7 +423,7 @@ subroutine pAdaptation_pAdapt(this, sem, itera, t, computeTimeDerivative, Comput ! ---------------------------------- ! call Stopwatch % Start("pAdapt: Adaptation") - call sem % mesh % pAdapt (NNew, controlVariables) + call sem % mesh % pAdapt_MPI (NNew, controlVariables) call Stopwatch % Pause("pAdapt: Adaptation") ! Reconstruct probes