diff --git a/src/SELF_DGModel1D_t.f90 b/src/SELF_DGModel1D_t.f90 index 99a3177f3..886bc1689 100644 --- a/src/SELF_DGModel1D_t.f90 +++ b/src/SELF_DGModel1D_t.f90 @@ -452,23 +452,21 @@ subroutine BoundaryFlux_DGModel1D_t(this) real(prec) :: fout(1:this%solution%nvar) real(prec) :: dfdx(1:this%solution%nvar),nhat - do iel = 1,this%mesh%nelem - do iside = 1,2 - - ! set the normal velocity - if(iside == 1) then - nhat = -1.0_prec - else - nhat = 1.0_prec - endif - - fin = this%solution%boundary(iside,iel,1:this%solution%nvar) ! interior solution - fout = this%solution%extboundary(iside,iel,1:this%solution%nvar) ! exterior solution - dfdx = this%solutionGradient%avgboundary(iside,iel,1:this%solution%nvar) ! average solution gradient (with direction taken into account) - this%flux%boundarynormal(iside,iel,1:this%solution%nvar) = & - this%riemannflux1d(fin,fout,dfdx,nhat) + do concurrent(iside=1:2,iel=1:this%mesh%nElem) + + ! set the normal velocity + if(iside == 1) then + nhat = -1.0_prec + else + nhat = 1.0_prec + endif + + fin = this%solution%boundary(iside,iel,1:this%solution%nvar) ! interior solution + fout = this%solution%extboundary(iside,iel,1:this%solution%nvar) ! exterior solution + dfdx = this%solutionGradient%avgboundary(iside,iel,1:this%solution%nvar) ! average solution gradient (with direction taken into account) + this%flux%boundarynormal(iside,iel,1:this%solution%nvar) = & + this%riemannflux1d(fin,fout,dfdx,nhat) - enddo enddo endsubroutine BoundaryFlux_DGModel1D_t @@ -481,16 +479,14 @@ subroutine fluxmethod_DGModel1D_t(this) integer :: i real(prec) :: f(1:this%solution%nvar),dfdx(1:this%solution%nvar) - do iel = 1,this%mesh%nelem - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,iel=1:this%mesh%nElem) - f = this%solution%interior(i,iel,1:this%solution%nvar) - dfdx = this%solutionGradient%interior(i,iel,1:this%solution%nvar) + f = this%solution%interior(i,iel,1:this%solution%nvar) + dfdx = this%solutionGradient%interior(i,iel,1:this%solution%nvar) - this%flux%interior(i,iel,1:this%solution%nvar) = & - this%flux1d(f,dfdx) + this%flux%interior(i,iel,1:this%solution%nvar) = & + this%flux1d(f,dfdx) - enddo enddo endsubroutine fluxmethod_DGModel1D_t @@ -503,16 +499,14 @@ subroutine sourcemethod_DGModel1D_t(this) integer :: i real(prec) :: f(1:this%solution%nvar),dfdx(1:this%solution%nvar) - do iel = 1,this%mesh%nelem - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,iel=1:this%mesh%nElem) - f = this%solution%interior(i,iel,1:this%solution%nvar) - dfdx = this%solutionGradient%interior(i,iel,1:this%solution%nvar) + f = this%solution%interior(i,iel,1:this%solution%nvar) + dfdx = this%solutionGradient%interior(i,iel,1:this%solution%nvar) - this%source%interior(i,iel,1:this%solution%nvar) = & - this%source1d(f,dfdx) + this%source%interior(i,iel,1:this%solution%nvar) = & + this%source1d(f,dfdx) - enddo enddo endsubroutine sourcemethod_DGModel1D_t @@ -578,19 +572,16 @@ subroutine Write_DGModel1D_t(this,fileName) call Open_HDF5(pickupFile,H5F_ACC_TRUNC_F,fileId) ! Write the interpolant to the file - INFO("Writing interpolant data to file") call this%solution%interp%WriteHDF5(fileId) ! In this section, we write the solution and geometry on the control (quadrature) grid ! which can be used for model pickup runs or post-processing ! Write the model state to file - INFO("Writing control grid solution to file") call CreateGroup_HDF5(fileId,'/controlgrid') call this%solution%WriteHDF5(fileId,'/controlgrid/solution') ! Write the geometry to file - INFO("Writing control grid geometry to file") call CreateGroup_HDF5(fileId,'/controlgrid/geometry') call this%geometry%x%WriteHDF5(fileId,'/controlgrid/geometry/x') ! -- END : writing solution on control grid -- ! @@ -614,12 +605,10 @@ subroutine Write_DGModel1D_t(this,fileName) call this%solution%GridInterp(solution%interior) ! Write the model state to file - INFO("Writing target grid solution to file") call CreateGroup_HDF5(fileId,'/targetgrid') call solution%WriteHDF5(fileId,'/targetgrid/solution') ! Write the geometry to file - INFO("Writing target grid geometry to file") call CreateGroup_HDF5(fileId,'/targetgrid/geometry') call x%WriteHDF5(fileId,'/targetgrid/geometry/x') diff --git a/src/SELF_DGModel2D_t.f90 b/src/SELF_DGModel2D_t.f90 index 297f6ddaa..11225f0ae 100644 --- a/src/SELF_DGModel2D_t.f90 +++ b/src/SELF_DGModel2D_t.f90 @@ -336,16 +336,13 @@ subroutine fluxmethod_DGModel2D_t(this) integer :: j real(prec) :: s(1:this%nvar),dsdx(1:this%nvar,1:2) - do iel = 1,this%mesh%nelem - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iel=1:this%mesh%nElem) - s = this%solution%interior(i,j,iel,1:this%nvar) - dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) - this%flux%interior(i,j,iel,1:this%nvar,1:2) = this%flux2d(s,dsdx) + s = this%solution%interior(i,j,iel,1:this%nvar) + dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) + this%flux%interior(i,j,iel,1:this%nvar,1:2) = this%flux2d(s,dsdx) - enddo - enddo enddo endsubroutine fluxmethod_DGModel2D_t @@ -364,21 +361,18 @@ subroutine BoundaryFlux_DGModel2D_t(this) real(prec) :: dsdx(1:this%nvar,1:2) real(prec) :: nhat(1:2),nmag - do iEl = 1,this%solution%nElem - do j = 1,4 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:4, & + iel=1:this%mesh%nElem) - ! Get the boundary normals on cell edges from the mesh geometry - nhat = this%geometry%nHat%boundary(i,j,iEl,1,1:2) - sL = this%solution%boundary(i,j,iel,1:this%nvar) ! interior solution - sR = this%solution%extboundary(i,j,iel,1:this%nvar) ! exterior solution - dsdx = this%solutiongradient%avgboundary(i,j,iel,1:this%nvar,1:2) - nmag = this%geometry%nScale%boundary(i,j,iEl,1) + ! Get the boundary normals on cell edges from the mesh geometry + nhat = this%geometry%nHat%boundary(i,j,iEl,1,1:2) + sL = this%solution%boundary(i,j,iel,1:this%nvar) ! interior solution + sR = this%solution%extboundary(i,j,iel,1:this%nvar) ! exterior solution + dsdx = this%solutiongradient%avgboundary(i,j,iel,1:this%nvar,1:2) + nmag = this%geometry%nScale%boundary(i,j,iEl,1) - this%flux%boundaryNormal(i,j,iEl,1:this%nvar) = this%riemannflux2d(sL,sR,dsdx,nhat)*nmag + this%flux%boundaryNormal(i,j,iEl,1:this%nvar) = this%riemannflux2d(sL,sR,dsdx,nhat)*nmag - enddo - enddo enddo endsubroutine BoundaryFlux_DGModel2D_t @@ -392,16 +386,13 @@ subroutine sourcemethod_DGModel2D_t(this) integer :: j real(prec) :: s(1:this%nvar),dsdx(1:this%nvar,1:2) - do iel = 1,this%mesh%nelem - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iel=1:this%mesh%nElem) - s = this%solution%interior(i,j,iel,1:this%nvar) - dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) - this%source%interior(i,j,iel,1:this%nvar) = this%source2d(s,dsdx) + s = this%solution%interior(i,j,iel,1:this%nvar) + dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) + this%source%interior(i,j,iel,1:this%nvar) = this%source2d(s,dsdx) - enddo - enddo enddo endsubroutine sourcemethod_DGModel2D_t @@ -416,44 +407,42 @@ subroutine setboundarycondition_DGModel2D_t(this) integer :: i,iEl,j,e2,bcid real(prec) :: nhat(1:2),x(1:2) - do iEl = 1,this%solution%nElem ! Loop over all elements - do j = 1,4 ! Loop over all sides + do concurrent(j=1:4,iel=1:this%mesh%nElem) - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID + bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID + e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then + if(e2 == 0) then + if(bcid == SELF_BC_PRESCRIBED) then - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,iEl,1,1:2) + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + x = this%geometry%x%boundary(i,j,iEl,1,1:2) - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_Prescribed(x,this%t) - enddo + this%solution%extBoundary(i,j,iEl,1:this%nvar) = & + this%hbc2d_Prescribed(x,this%t) + enddo - elseif(bcid == SELF_BC_RADIATION) then + elseif(bcid == SELF_BC_RADIATION) then - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_Radiation(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) - enddo + this%solution%extBoundary(i,j,iEl,1:this%nvar) = & + this%hbc2d_Radiation(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then + elseif(bcid == SELF_BC_NONORMALFLOW) then - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_NoNormalFlow(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) - enddo + this%solution%extBoundary(i,j,iEl,1:this%nvar) = & + this%hbc2d_NoNormalFlow(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) + enddo - endif endif + endif - enddo enddo endsubroutine setboundarycondition_DGModel2D_t @@ -469,48 +458,46 @@ subroutine setgradientboundarycondition_DGModel2D_t(this) real(prec) :: dsdx(1:this%nvar,1:2) real(prec) :: nhat(1:2),x(1:2) - do iEl = 1,this%solution%nElem ! Loop over all elements - do j = 1,4 ! Loop over all sides + do concurrent(j=1:4,iel=1:this%mesh%nElem) - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID + bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID + e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then + if(e2 == 0) then + if(bcid == SELF_BC_PRESCRIBED) then - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,iEl,1,1:2) + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + x = this%geometry%x%boundary(i,j,iEl,1,1:2) - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_Prescribed(x,this%t) - enddo + this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & + this%pbc2d_Prescribed(x,this%t) + enddo - elseif(bcid == SELF_BC_RADIATION) then + elseif(bcid == SELF_BC_RADIATION) then - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) + dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_Radiation(dsdx,nhat) - enddo + this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & + this%pbc2d_Radiation(dsdx,nhat) + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then + elseif(bcid == SELF_BC_NONORMALFLOW) then - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) + dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_NoNormalFlow(dsdx,nhat) - enddo + this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & + this%pbc2d_NoNormalFlow(dsdx,nhat) + enddo - endif endif + endif - enddo enddo endsubroutine setgradientboundarycondition_DGModel2D_t @@ -575,20 +562,17 @@ subroutine Write_DGModel2D_t(this,fileName) call Open_HDF5(pickupFile,H5F_ACC_TRUNC_F,fileId,this%mesh%decomp%mpiComm) ! Write the interpolant to the file - print*,__FILE__//" : Writing interpolant data to file" call this%solution%interp%WriteHDF5(fileId) ! In this section, we write the solution and geometry on the control (quadrature) grid ! which can be used for model pickup runs or post-processing ! Write the model state to file - print*,__FILE__//" : Writing control grid solution to file" call CreateGroup_HDF5(fileId,'/controlgrid') print*," offset, nglobal_elem : ",this%mesh%decomp%offsetElem(this%mesh%decomp%rankId+1),this%mesh%decomp%nElem call this%solution%WriteHDF5(fileId,'/controlgrid/solution', & this%mesh%decomp%offsetElem(this%mesh%decomp%rankId+1),this%mesh%decomp%nElem) ! Write the geometry to file - print*,__FILE__//" : Writing control grid geometry to file" call this%geometry%x%WriteHDF5(fileId,'/controlgrid/geometry', & this%mesh%decomp%offsetElem(this%mesh%decomp%rankId+1),this%mesh%decomp%nElem) @@ -601,19 +585,16 @@ subroutine Write_DGModel2D_t(this,fileName) call Open_HDF5(pickupFile,H5F_ACC_TRUNC_F,fileId) ! Write the interpolant to the file - print*,__FILE__//" : Writing interpolant data to file" call this%solution%interp%WriteHDF5(fileId) ! In this section, we write the solution and geometry on the control (quadrature) grid ! which can be used for model pickup runs or post-processing ! Write the model state to file - print*,__FILE__//" : Writing control grid solution to file" call CreateGroup_HDF5(fileId,'/controlgrid') call this%solution%WriteHDF5(fileId,'/controlgrid/solution') ! Write the geometry to file - print*,__FILE__//" : Writing control grid geometry to file" call this%geometry%x%WriteHDF5(fileId,'/controlgrid/geometry') ! -- END : writing solution on control grid -- ! diff --git a/src/SELF_DGModel3D_t.f90 b/src/SELF_DGModel3D_t.f90 index e85b1872c..cefd73e3a 100644 --- a/src/SELF_DGModel3D_t.f90 +++ b/src/SELF_DGModel3D_t.f90 @@ -335,18 +335,13 @@ subroutine fluxmethod_DGModel3D_t(this) integer :: i,j,k real(prec) :: s(1:this%nvar),dsdx(1:this%nvar,1:3) - do iel = 1,this%mesh%nelem - do k = 1,this%solution%interp%N+1 - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + k=1:this%solution%N+1,iel=1:this%mesh%nElem) - s = this%solution%interior(i,j,k,iel,1:this%nvar) - dsdx = this%solutionGradient%interior(i,j,k,iel,1:this%nvar,1:3) - this%flux%interior(i,j,k,iel,1:this%nvar,1:3) = this%flux3d(s,dsdx) + s = this%solution%interior(i,j,k,iel,1:this%nvar) + dsdx = this%solutionGradient%interior(i,j,k,iel,1:this%nvar,1:3) + this%flux%interior(i,j,k,iel,1:this%nvar,1:3) = this%flux3d(s,dsdx) - enddo - enddo - enddo enddo endsubroutine fluxmethod_DGModel3D_t @@ -363,21 +358,17 @@ subroutine BoundaryFlux_DGModel3D_t(this) real(prec) :: dsdx(1:this%nvar,1:3) real(prec) :: nhat(1:3),nmag - do iEl = 1,this%solution%nElem - do k = 1,6 - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 - ! Get the boundary normals on cell edges from the mesh geometry - nhat = this%geometry%nHat%boundary(i,j,k,iEl,1,1:3) - sL = this%solution%boundary(i,j,k,iel,1:this%nvar) ! interior solution - sR = this%solution%extboundary(i,j,k,iel,1:this%nvar) ! exterior solution - dsdx = this%solutiongradient%avgboundary(i,j,k,iel,1:this%nvar,1:3) - nmag = this%geometry%nScale%boundary(i,j,k,iEl,1) - - this%flux%boundaryNormal(i,j,k,iEl,1:this%nvar) = this%riemannflux3d(sL,sR,dsdx,nhat)*nmag - enddo - enddo - enddo + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + k=1:6,iel=1:this%mesh%nElem) + ! Get the boundary normals on cell edges from the mesh geometry + nhat = this%geometry%nHat%boundary(i,j,k,iEl,1,1:3) + sL = this%solution%boundary(i,j,k,iel,1:this%nvar) ! interior solution + sR = this%solution%extboundary(i,j,k,iel,1:this%nvar) ! exterior solution + dsdx = this%solutiongradient%avgboundary(i,j,k,iel,1:this%nvar,1:3) + nmag = this%geometry%nScale%boundary(i,j,k,iEl,1) + + this%flux%boundaryNormal(i,j,k,iEl,1:this%nvar) = this%riemannflux3d(sL,sR,dsdx,nhat)*nmag + enddo endsubroutine BoundaryFlux_DGModel3D_t @@ -389,18 +380,13 @@ subroutine sourcemethod_DGModel3D_t(this) integer :: i,j,k,iel real(prec) :: s(1:this%nvar),dsdx(1:this%nvar,1:3) - do iel = 1,this%mesh%nelem - do k = 1,this%solution%interp%N+1 - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + k=1:this%solution%N+1,iel=1:this%mesh%nElem) - s = this%solution%interior(i,j,k,iel,1:this%nvar) - dsdx = this%solutionGradient%interior(i,j,k,iel,1:this%nvar,1:3) - this%source%interior(i,j,k,iel,1:this%nvar) = this%source3d(s,dsdx) + s = this%solution%interior(i,j,k,iel,1:this%nvar) + dsdx = this%solutionGradient%interior(i,j,k,iel,1:this%nvar,1:3) + this%source%interior(i,j,k,iel,1:this%nvar) = this%source3d(s,dsdx) - enddo - enddo - enddo enddo endsubroutine sourcemethod_DGModel3D_t @@ -415,50 +401,48 @@ subroutine setboundarycondition_DGModel3D_t(this) integer :: i,iEl,j,k,e2,bcid real(prec) :: nhat(1:3),x(1:3) - do iEl = 1,this%solution%nElem ! Loop over all elements - do k = 1,6 ! Loop over all sides + do concurrent(k=1:6,iel=1:this%mesh%nElem) - bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID + bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID + e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then + if(e2 == 0) then + if(bcid == SELF_BC_PRESCRIBED) then - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solution%interp%N+1 ! Loop over quadrature points + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + x = this%geometry%x%boundary(i,j,k,iEl,1,1:3) - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_Prescribed(x,this%t) - enddo + this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & + this%hbc3d_Prescribed(x,this%t) enddo + enddo - elseif(bcid == SELF_BC_RADIATION) then + elseif(bcid == SELF_BC_RADIATION) then - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solution%interp%N+1 ! Loop over quadrature points + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_Radiation(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) - enddo + this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & + this%hbc3d_Radiation(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) enddo + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then + elseif(bcid == SELF_BC_NONORMALFLOW) then - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solution%interp%N+1 ! Loop over quadrature points + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_NoNormalFlow(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) - enddo + this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & + this%hbc3d_NoNormalFlow(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) enddo + enddo - endif endif + endif - enddo enddo endsubroutine setboundarycondition_DGModel3D_t @@ -474,54 +458,52 @@ subroutine setgradientboundarycondition_DGModel3D_t(this) real(prec) :: dsdx(1:this%nvar,1:3) real(prec) :: nhat(1:3),x(1:3) - do iEl = 1,this%solution%nElem ! Loop over all elements - do k = 1,6 ! Loop over all sides + do concurrent(k=1:6,iel=1:this%mesh%nElem) - bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID + bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID + e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then + if(e2 == 0) then + if(bcid == SELF_BC_PRESCRIBED) then - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - x = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + x = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_Prescribed(x,this%t) - enddo + this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & + this%pbc3d_Prescribed(x,this%t) enddo + enddo - elseif(bcid == SELF_BC_RADIATION) then + elseif(bcid == SELF_BC_RADIATION) then - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) + dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_Radiation(dsdx,nhat) - enddo + this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & + this%pbc3d_Radiation(dsdx,nhat) enddo + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then + elseif(bcid == SELF_BC_NONORMALFLOW) then - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) + dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_NoNormalFlow(dsdx,nhat) - enddo + this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & + this%pbc3d_NoNormalFlow(dsdx,nhat) enddo + enddo - endif endif + endif - enddo enddo endsubroutine setgradientboundarycondition_DGModel3D_t @@ -585,19 +567,16 @@ subroutine Write_DGModel3D_t(this,fileName) call Open_HDF5(pickupFile,H5F_ACC_TRUNC_F,fileId,this%mesh%decomp%mpiComm) ! Write the interpolant to the file - print*,__FILE__//" : Writing interpolant data to file" call this%solution%interp%WriteHDF5(fileId) ! In this section, we write the solution and geometry on the control (quadrature) grid ! which can be used for model pickup runs or post-processing ! Write the model state to file - print*,__FILE__//" : Writing control grid solution to file" call CreateGroup_HDF5(fileId,'/controlgrid') call this%solution%WriteHDF5(fileId,'/controlgrid/solution', & this%mesh%decomp%offsetElem(this%mesh%decomp%rankId+1),this%mesh%decomp%nElem) ! Write the geometry to file - print*,__FILE__//" : Writing control grid geometry to file" call this%geometry%x%WriteHDF5(fileId,'/controlgrid/geometry', & this%mesh%decomp%offsetElem(this%mesh%decomp%rankId+1),this%mesh%decomp%nElem) @@ -610,19 +589,16 @@ subroutine Write_DGModel3D_t(this,fileName) call Open_HDF5(pickupFile,H5F_ACC_TRUNC_F,fileId) ! Write the interpolant to the file - print*,__FILE__//" : Writing interpolant data to file" call this%solution%interp%WriteHDF5(fileId) ! In this section, we write the solution and geometry on the control (quadrature) grid ! which can be used for model pickup runs or post-processing ! Write the model state to file - print*,__FILE__//" : Writing control grid solution to file" call CreateGroup_HDF5(fileId,'/controlgrid') call this%solution%WriteHDF5(fileId,'/controlgrid/solution') ! Write the geometry to file - print*,__FILE__//" : Writing control grid geometry to file" call this%geometry%x%WriteHDF5(fileId,'/controlgrid/geometry') ! -- END : writing solution on control grid -- ! diff --git a/src/gpu/SELF_DGModel1D.f90 b/src/gpu/SELF_DGModel1D.f90 index 66d456575..55d42de39 100644 --- a/src/gpu/SELF_DGModel1D.f90 +++ b/src/gpu/SELF_DGModel1D.f90 @@ -340,23 +340,21 @@ subroutine BoundaryFlux_DGModel1D(this) this%solutiongradient%avgboundary_gpu,sizeof(this%solutiongradient%avgboundary), & hipMemcpyDeviceToHost)) - do iel = 1,this%mesh%nelem - do iside = 1,2 - - ! set the normal velocity - if(iside == 1) then - nhat = -1.0_prec - else - nhat = 1.0_prec - endif - - fin = this%solution%boundary(iside,iel,1:this%solution%nvar) ! interior solution - fout = this%solution%extboundary(iside,iel,1:this%solution%nvar) ! exterior solution - dfdx = this%solutionGradient%avgboundary(iside,iel,1:this%solution%nvar) ! average solution gradient (with direction taken into account) - this%flux%boundarynormal(iside,iel,1:this%solution%nvar) = & - this%riemannflux1d(fin,fout,dfdx,nhat) + do concurrent(iside=1:2,iel=1:this%mesh%nElem) + + ! set the normal velocity + if(iside == 1) then + nhat = -1.0_prec + else + nhat = 1.0_prec + endif + + fin = this%solution%boundary(iside,iel,1:this%solution%nvar) ! interior solution + fout = this%solution%extboundary(iside,iel,1:this%solution%nvar) ! exterior solution + dfdx = this%solutionGradient%avgboundary(iside,iel,1:this%solution%nvar) ! average solution gradient (with direction taken into account) + this%flux%boundarynormal(iside,iel,1:this%solution%nvar) = & + this%riemannflux1d(fin,fout,dfdx,nhat) - enddo enddo call gpuCheck(hipMemcpy(this%flux%boundarynormal_gpu, & @@ -374,16 +372,14 @@ subroutine fluxmethod_DGModel1D(this) integer :: i real(prec) :: f(1:this%solution%nvar),dfdx(1:this%solution%nvar) - do iel = 1,this%mesh%nelem - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,iel=1:this%mesh%nElem) - f = this%solution%interior(i,iel,1:this%solution%nvar) - dfdx = this%solutionGradient%interior(i,iel,1:this%solution%nvar) + f = this%solution%interior(i,iel,1:this%solution%nvar) + dfdx = this%solutionGradient%interior(i,iel,1:this%solution%nvar) - this%flux%interior(i,iel,1:this%solution%nvar) = & - this%flux1d(f,dfdx) + this%flux%interior(i,iel,1:this%solution%nvar) = & + this%flux1d(f,dfdx) - enddo enddo call gpuCheck(hipMemcpy(this%flux%interior_gpu, & @@ -409,16 +405,14 @@ subroutine sourcemethod_DGModel1D(this) this%solutiongradient%interior_gpu,sizeof(this%solutiongradient%interior), & hipMemcpyDeviceToHost)) - do iel = 1,this%mesh%nelem - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,iel=1:this%mesh%nElem) - f = this%solution%interior(i,iel,1:this%solution%nvar) - dfdx = this%solutionGradient%interior(i,iel,1:this%solution%nvar) + f = this%solution%interior(i,iel,1:this%solution%nvar) + dfdx = this%solutionGradient%interior(i,iel,1:this%solution%nvar) - this%source%interior(i,iel,1:this%solution%nvar) = & - this%source1d(f,dfdx) + this%source%interior(i,iel,1:this%solution%nvar) = & + this%source1d(f,dfdx) - enddo enddo call gpuCheck(hipMemcpy(this%source%interior_gpu, & diff --git a/src/gpu/SELF_DGModel2D.f90 b/src/gpu/SELF_DGModel2D.f90 index 7d5da4330..30a651b37 100644 --- a/src/gpu/SELF_DGModel2D.f90 +++ b/src/gpu/SELF_DGModel2D.f90 @@ -194,16 +194,13 @@ subroutine fluxmethod_DGModel2D(this) integer :: j real(prec) :: s(1:this%nvar),dsdx(1:this%nvar,1:2) - do iel = 1,this%mesh%nelem - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iel=1:this%mesh%nElem) - s = this%solution%interior(i,j,iel,1:this%nvar) - dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) - this%flux%interior(i,j,iel,1:this%nvar,1:2) = this%flux2d(s,dsdx) + s = this%solution%interior(i,j,iel,1:this%nvar) + dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) + this%flux%interior(i,j,iel,1:this%nvar,1:2) = this%flux2d(s,dsdx) - enddo - enddo enddo call gpuCheck(hipMemcpy(this%flux%interior_gpu, & @@ -239,21 +236,18 @@ subroutine BoundaryFlux_DGModel2D(this) this%solutiongradient%avgboundary_gpu,sizeof(this%solutiongradient%avgboundary), & hipMemcpyDeviceToHost)) - do iEl = 1,this%solution%nElem - do j = 1,4 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:4, & + iel=1:this%mesh%nElem) - ! Get the boundary normals on cell edges from the mesh geometry - nhat = this%geometry%nHat%boundary(i,j,iEl,1,1:2) - sL = this%solution%boundary(i,j,iel,1:this%nvar) ! interior solution - sR = this%solution%extboundary(i,j,iel,1:this%nvar) ! exterior solution - dsdx = this%solutiongradient%avgboundary(i,j,iel,1:this%nvar,1:2) - nmag = this%geometry%nScale%boundary(i,j,iEl,1) + ! Get the boundary normals on cell edges from the mesh geometry + nhat = this%geometry%nHat%boundary(i,j,iEl,1,1:2) + sL = this%solution%boundary(i,j,iel,1:this%nvar) ! interior solution + sR = this%solution%extboundary(i,j,iel,1:this%nvar) ! exterior solution + dsdx = this%solutiongradient%avgboundary(i,j,iel,1:this%nvar,1:2) + nmag = this%geometry%nScale%boundary(i,j,iEl,1) - this%flux%boundaryNormal(i,j,iEl,1:this%nvar) = this%riemannflux2d(sL,sR,dsdx,nhat)*nmag + this%flux%boundaryNormal(i,j,iEl,1:this%nvar) = this%riemannflux2d(sL,sR,dsdx,nhat)*nmag - enddo - enddo enddo call gpuCheck(hipMemcpy(this%flux%boundarynormal_gpu, & @@ -280,16 +274,13 @@ subroutine sourcemethod_DGModel2D(this) this%solutiongradient%interior_gpu,sizeof(this%solutiongradient%interior), & hipMemcpyDeviceToHost)) - do iel = 1,this%mesh%nelem - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iel=1:this%mesh%nElem) - s = this%solution%interior(i,j,iel,1:this%nvar) - dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) - this%source%interior(i,j,iel,1:this%nvar) = this%source2d(s,dsdx) + s = this%solution%interior(i,j,iel,1:this%nvar) + dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) + this%source%interior(i,j,iel,1:this%nvar) = this%source2d(s,dsdx) - enddo - enddo enddo call gpuCheck(hipMemcpy(this%source%interior_gpu, & @@ -317,44 +308,42 @@ subroutine setboundarycondition_DGModel2D(this) this%solution%extboundary_gpu,sizeof(this%solution%extboundary), & hipMemcpyDeviceToHost)) - do iEl = 1,this%solution%nElem ! Loop over all elements - do j = 1,4 ! Loop over all sides + do concurrent(j=1:4,iel=1:this%mesh%nElem) - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID + bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID + e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then + if(e2 == 0) then + if(bcid == SELF_BC_PRESCRIBED) then - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,iEl,1,1:2) + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + x = this%geometry%x%boundary(i,j,iEl,1,1:2) - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_Prescribed(x,this%t) - enddo + this%solution%extBoundary(i,j,iEl,1:this%nvar) = & + this%hbc2d_Prescribed(x,this%t) + enddo - elseif(bcid == SELF_BC_RADIATION) then + elseif(bcid == SELF_BC_RADIATION) then - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_Radiation(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) - enddo + this%solution%extBoundary(i,j,iEl,1:this%nvar) = & + this%hbc2d_Radiation(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then + elseif(bcid == SELF_BC_NONORMALFLOW) then - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_NoNormalFlow(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) - enddo + this%solution%extBoundary(i,j,iEl,1:this%nvar) = & + this%hbc2d_NoNormalFlow(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) + enddo - endif endif + endif - enddo enddo call gpuCheck(hipMemcpy(this%solution%extBoundary_gpu, & @@ -383,48 +372,46 @@ subroutine setgradientboundarycondition_DGModel2D(this) this%solutiongradient%extboundary_gpu,sizeof(this%solutiongradient%extboundary), & hipMemcpyDeviceToHost)) - do iEl = 1,this%solution%nElem ! Loop over all elements - do j = 1,4 ! Loop over all sides + do concurrent(j=1:4,iel=1:this%mesh%nElem) - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID + bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID + e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then + if(e2 == 0) then + if(bcid == SELF_BC_PRESCRIBED) then - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,iEl,1,1:2) + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + x = this%geometry%x%boundary(i,j,iEl,1,1:2) - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_Prescribed(x,this%t) - enddo + this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & + this%pbc2d_Prescribed(x,this%t) + enddo - elseif(bcid == SELF_BC_RADIATION) then + elseif(bcid == SELF_BC_RADIATION) then - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) + dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_Radiation(dsdx,nhat) - enddo + this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & + this%pbc2d_Radiation(dsdx,nhat) + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then + elseif(bcid == SELF_BC_NONORMALFLOW) then - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) + dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_NoNormalFlow(dsdx,nhat) - enddo + this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & + this%pbc2d_NoNormalFlow(dsdx,nhat) + enddo - endif endif + endif - enddo enddo call gpuCheck(hipMemcpy(this%solutiongradient%extBoundary_gpu, & diff --git a/src/gpu/SELF_DGModel3D.f90 b/src/gpu/SELF_DGModel3D.f90 index 471c0181b..294149cad 100644 --- a/src/gpu/SELF_DGModel3D.f90 +++ b/src/gpu/SELF_DGModel3D.f90 @@ -199,18 +199,13 @@ subroutine fluxmethod_DGModel3D(this) integer :: i,j,k real(prec) :: s(1:this%nvar),dsdx(1:this%nvar,1:3) - do iel = 1,this%mesh%nelem - do k = 1,this%solution%interp%N+1 - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + k=1:this%solution%N+1,iel=1:this%mesh%nElem) - s = this%solution%interior(i,j,k,iel,1:this%nvar) - dsdx = this%solutionGradient%interior(i,j,k,iel,1:this%nvar,1:3) - this%flux%interior(i,j,k,iel,1:this%nvar,1:3) = this%flux3d(s,dsdx) + s = this%solution%interior(i,j,k,iel,1:this%nvar) + dsdx = this%solutionGradient%interior(i,j,k,iel,1:this%nvar,1:3) + this%flux%interior(i,j,k,iel,1:this%nvar,1:3) = this%flux3d(s,dsdx) - enddo - enddo - enddo enddo call gpuCheck(hipMemcpy(this%flux%interior_gpu, & @@ -243,21 +238,18 @@ subroutine BoundaryFlux_DGModel3D(this) call gpuCheck(hipMemcpy(c_loc(this%solutiongradient%avgboundary), & this%solutiongradient%avgboundary_gpu,sizeof(this%solutiongradient%avgboundary), & hipMemcpyDeviceToHost)) - do iEl = 1,this%solution%nElem - do k = 1,6 - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 - ! Get the boundary normals on cell edges from the mesh geometry - nhat = this%geometry%nHat%boundary(i,j,k,iEl,1,1:3) - sL = this%solution%boundary(i,j,k,iel,1:this%nvar) ! interior solution - sR = this%solution%extboundary(i,j,k,iel,1:this%nvar) ! exterior solution - dsdx = this%solutiongradient%avgboundary(i,j,k,iel,1:this%nvar,1:3) - nmag = this%geometry%nScale%boundary(i,j,k,iEl,1) - - this%flux%boundaryNormal(i,j,k,iEl,1:this%nvar) = this%riemannflux3d(sL,sR,dsdx,nhat)*nmag - enddo - enddo - enddo + + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + k=1:6,iel=1:this%mesh%nElem) + ! Get the boundary normals on cell edges from the mesh geometry + nhat = this%geometry%nHat%boundary(i,j,k,iEl,1,1:3) + sL = this%solution%boundary(i,j,k,iel,1:this%nvar) ! interior solution + sR = this%solution%extboundary(i,j,k,iel,1:this%nvar) ! exterior solution + dsdx = this%solutiongradient%avgboundary(i,j,k,iel,1:this%nvar,1:3) + nmag = this%geometry%nScale%boundary(i,j,k,iEl,1) + + this%flux%boundaryNormal(i,j,k,iEl,1:this%nvar) = this%riemannflux3d(sL,sR,dsdx,nhat)*nmag + enddo call gpuCheck(hipMemcpy(this%flux%boundarynormal_gpu, & @@ -282,18 +274,13 @@ subroutine sourcemethod_DGModel3D(this) this%solutiongradient%interior_gpu,sizeof(this%solutiongradient%interior), & hipMemcpyDeviceToHost)) - do iel = 1,this%mesh%nelem - do k = 1,this%solution%interp%N+1 - do j = 1,this%solution%interp%N+1 - do i = 1,this%solution%interp%N+1 + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + k=1:this%solution%N+1,iel=1:this%mesh%nElem) - s = this%solution%interior(i,j,k,iel,1:this%nvar) - dsdx = this%solutionGradient%interior(i,j,k,iel,1:this%nvar,1:3) - this%source%interior(i,j,k,iel,1:this%nvar) = this%source3d(s,dsdx) + s = this%solution%interior(i,j,k,iel,1:this%nvar) + dsdx = this%solutionGradient%interior(i,j,k,iel,1:this%nvar,1:3) + this%source%interior(i,j,k,iel,1:this%nvar) = this%source3d(s,dsdx) - enddo - enddo - enddo enddo call gpuCheck(hipMemcpy(this%source%interior_gpu, & @@ -321,50 +308,48 @@ subroutine setboundarycondition_DGModel3D(this) this%solution%extboundary_gpu,sizeof(this%solution%extboundary), & hipMemcpyDeviceToHost)) - do iEl = 1,this%solution%nElem ! Loop over all elements - do k = 1,6 ! Loop over all sides + do concurrent(k=1:6,iel=1:this%mesh%nElem) - bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID + bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID + e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then + if(e2 == 0) then + if(bcid == SELF_BC_PRESCRIBED) then - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solution%interp%N+1 ! Loop over quadrature points + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + x = this%geometry%x%boundary(i,j,k,iEl,1,1:3) - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_Prescribed(x,this%t) - enddo + this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & + this%hbc3d_Prescribed(x,this%t) enddo + enddo - elseif(bcid == SELF_BC_RADIATION) then + elseif(bcid == SELF_BC_RADIATION) then - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solution%interp%N+1 ! Loop over quadrature points + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_Radiation(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) - enddo + this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & + this%hbc3d_Radiation(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) enddo + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then + elseif(bcid == SELF_BC_NONORMALFLOW) then - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solution%interp%N+1 ! Loop over quadrature points + do i = 1,this%solution%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_NoNormalFlow(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) - enddo + this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & + this%hbc3d_NoNormalFlow(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) enddo + enddo - endif endif + endif - enddo enddo call gpuCheck(hipMemcpy(this%solution%extBoundary_gpu, & @@ -393,54 +378,52 @@ subroutine setgradientboundarycondition_DGModel3D(this) this%solutiongradient%extboundary_gpu,sizeof(this%solutiongradient%extboundary), & hipMemcpyDeviceToHost)) - do iEl = 1,this%solution%nElem ! Loop over all elements - do k = 1,6 ! Loop over all sides + do concurrent(k=1:6,iel=1:this%mesh%nElem) - bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID + bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID + e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then + if(e2 == 0) then + if(bcid == SELF_BC_PRESCRIBED) then - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - x = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + x = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_Prescribed(x,this%t) - enddo + this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & + this%pbc3d_Prescribed(x,this%t) enddo + enddo - elseif(bcid == SELF_BC_RADIATION) then + elseif(bcid == SELF_BC_RADIATION) then - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) + dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_Radiation(dsdx,nhat) - enddo + this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & + this%pbc3d_Radiation(dsdx,nhat) enddo + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then + elseif(bcid == SELF_BC_NONORMALFLOW) then - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points + nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) + dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_NoNormalFlow(dsdx,nhat) - enddo + this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & + this%pbc3d_NoNormalFlow(dsdx,nhat) enddo + enddo - endif endif + endif - enddo enddo call gpuCheck(hipMemcpy(this%solutiongradient%extBoundary_gpu, &