Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/do concurrent #74

Merged
merged 3 commits into from
Nov 20, 2024
Merged

Feature/do concurrent #74

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2cd9705
Replace tightly nested loops with do-concurrent
fluidnumerics-joe Nov 15, 2024
70bfa3c
Merge remote-tracking branch 'origin/main' into feature/do-concurrent
fluidnumerics-joe Nov 18, 2024
2496357
Merge remote-tracking branch 'origin/main' into feature/do-concurrent
fluidnumerics-joe Nov 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 24 additions & 35 deletions src/SELF_DGModel1D_t.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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 -- !
Expand All @@ -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')

Expand Down
153 changes: 67 additions & 86 deletions src/SELF_DGModel2D_t.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -336,16 +336,13 @@
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
Expand All @@ -364,21 +361,18 @@
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
Expand All @@ -392,16 +386,13 @@
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
Expand All @@ -416,44 +407,42 @@
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
Expand All @@ -469,48 +458,46 @@
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)

Check warning on line 479 in src/SELF_DGModel2D_t.f90

View check run for this annotation

Codecov / codecov/patch

src/SELF_DGModel2D_t.f90#L478-L479 

Added lines #L478 - L479 were not covered by tests

dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2)
dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2)

Check warning on line 481 in src/SELF_DGModel2D_t.f90

View check run for this annotation

Codecov / codecov/patch

src/SELF_DGModel2D_t.f90#L481 

Added line #L481 was not covered by tests

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)

Check warning on line 484 in src/SELF_DGModel2D_t.f90

View check run for this annotation

Codecov / codecov/patch

src/SELF_DGModel2D_t.f90#L483-L484 

Added lines #L483 - L484 were not covered by tests
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
Expand Down Expand Up @@ -575,20 +562,17 @@
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)

Expand All @@ -601,19 +585,16 @@
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 -- !

Expand Down
Loading
Loading