Skip to content

Commit 2747b57

Browse files
p-costap-costa
committed
Symmetrize Poisson/Helmholtz problem to exploit matrix symmetries.
1 parent 6f657b4 commit 2747b57

File tree

3 files changed

+68
-56
lines changed

3 files changed

+68
-56
lines changed

‎src/fillps.f90‎

Lines changed: 17 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -30,9 +30,23 @@ subroutine fillps(lo,hi,dxf,dyf,dzf,dt,up,vp,wp,p)
3030
do j=lo(2),hi(2)
3131
do i=lo(1),hi(1)
3232
p(i,j,k) = ( &
33-
(wp(i,j,k)-wp(i,j,k-1))/dzf(k) + &
34-
(vp(i,j,k)-vp(i,j-1,k))/dyf(j) + &
35-
(up(i,j,k)-up(i-1,j,k))/dxf(i) &
33+
#if defined(_FFT_Z)
34+
(wp(i,j,k)-wp(i,j,k-1))/dzf(k)*dxf(i)*dyf(j) + &
35+
(vp(i,j,k)-vp(i,j-1,k))*dxf(i) + &
36+
(up(i,j,k)-up(i-1,j,k))*dyf(j) &
37+
#elif defined(_FFT_Y)
38+
(wp(i,j,k)-wp(i,j,k-1))*dxf(i) + &
39+
(vp(i,j,k)-vp(i,j-1,k))/dyf(j)*dxf(i)*dzf(k) + &
40+
(up(i,j,k)-up(i-1,j,k))*dzf(k) &
41+
#elif defined(_FFT_X)
42+
(wp(i,j,k)-wp(i,j,k-1))*dyf(j) + &
43+
(vp(i,j,k)-vp(i,j-1,k))*dzf(k) + &
44+
(up(i,j,k)-up(i-1,j,k))/dxf(i)*dyf(j)*dzf(k) &
45+
#else
46+
(wp(i,j,k)-wp(i,j,k-1))*dxf(i)*dyf(j) + &
47+
(vp(i,j,k)-vp(i,j-1,k))*dxf(i)*dzf(k) + &
48+
(up(i,j,k)-up(i-1,j,k))*dyf(j)*dzf(k) &
49+
#endif
3650
)/dt
3751
end do
3852
end do

‎src/main.f90‎

Lines changed: 9 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -95,13 +95,9 @@ program snac
9595
integer , dimension( 3) :: hiu,hiv,hiw
9696
#endif
9797
type(hypre_solver) :: psolver
98-
logical :: is_symm_matrix_p
9998
#ifdef _IMPDIFF
10099
type(hypre_solver) :: usolver,vsolver,wsolver
101100
real(rp) :: alphai,alphaoi
102-
logical :: is_symm_matrix_u, &
103-
is_symm_matrix_v, &
104-
is_symm_matrix_w
105101
#endif
106102
!
107103
real(rp) :: dt,dtmax,time,dtrk,divtot,divmax
@@ -112,7 +108,6 @@ program snac
112108
dxc_g,dxf_g,xc_g,xf_g, &
113109
dyc_g,dyf_g,yc_g,yf_g, &
114110
dzc_g,dzf_g,zc_g,zf_g
115-
logical :: is_uniform_grid
116111
logical, dimension(3) :: is_centered
117112
!
118113
#if defined(_FFT_X) || defined(_FFT_Y) || defined(_FFT_Z)
@@ -501,21 +496,6 @@ program snac
501496
call MPI_ALLREDUCE(MPI_IN_PLACE,dzf_g(hi_g(3)+1),1,MPI_REAL_RP,MPI_MAX,comm_block)
502497
#endif
503498
#endif
504-
is_uniform_grid = all(dzf(:) == dzf(lo(3))) .and. &
505-
all(dyf(:) == dyf(lo(2))) .and. &
506-
all(dxf(:) == dxf(lo(1))) .and. &
507-
#ifdef _FFT_X
508-
dyf(lo(2)) == dzf(lo(3))
509-
#elif _FFT_Y
510-
dxf(lo(1)) == dzf(lo(3))
511-
#elif _FFT_Z
512-
dxf(lo(1)) == dyf(lo(2))
513-
#else
514-
dxf(lo(1)) == dyf(lo(2)) .and. &
515-
dxf(lo(1)) == dzf(lo(3)) .and. &
516-
dyf(lo(2)) == dzf(lo(3))
517-
#endif
518-
call mpi_allreduce(MPI_IN_PLACE,is_uniform_grid,1,MPI_LOGICAL,MPI_LAND,MPI_COMM_WORLD)
519499
!
520500
! initialization of the flow fields
521501
!
@@ -713,7 +693,6 @@ program snac
713693
allocate(comms_fft(hi_a(idir)-lo_a(idir)+1))
714694
allocate(lambda_p_a(hi_a(idir)-lo_a(idir)+1))
715695
is_bound_a(:,:) = is_bound(:,:)
716-
is_symm_matrix_p = is_uniform_grid
717696
#ifndef _FFT_USE_SLABS
718697
comms_fft(:) = MPI_COMM_WORLD
719698
lambda_p_a(:) = lambda_p
@@ -728,17 +707,17 @@ program snac
728707
#endif
729708
#ifndef _FFT_USE_SLICED_PENCILS
730709
call init_n_2d_matrices(cbcpre(:,il:iu:iskip),bcpre(:,il:iu:iskip),dl(:,il:iu:iskip), &
731-
is_symm_matrix_p,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
710+
.true.,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
732711
lo_a(idir),hi_a(idir),lo_a(il:iu:iskip),hi_a(il:iu:iskip),periods(il:iu:iskip), &
733712
dl1_1,dl1_2,dl2_1,dl2_2,alpha,alpha_bc,lambda_p_a,comms_fft,psolver_fft)
734713
#else
735-
call init_n_3d_matrices(idir,nslices,cbcpre,bcpre,dl,is_symm_matrix_p,is_bound,is_centered,lo,periods, &
714+
call init_n_3d_matrices(idir,nslices,cbcpre,bcpre,dl,.true.,is_bound,is_centered,lo,periods, &
736715
lo_sp,hi_sp,dxc,dxf,dyc,dyf,dzc,dzf,alpha,alpha_bc,lambda_p_a,psolver_fft)
737716
#endif
738717
call create_n_solvers(npsolvers,hypre_maxiter,hypre_tol,hypre_solver_i,psolver_fft)
739718
call setup_n_solvers(npsolvers,psolver_fft)
740719
#else
741-
call init_matrix_3d(cbcpre,bcpre,dl,is_symm_matrix_p,is_bound,is_centered,lo,hi,periods, &
720+
call init_matrix_3d(cbcpre,bcpre,dl,.true.,is_bound,is_centered,lo,hi,periods, &
742721
dxc,dxf,dyc,dyf,dzc,dzf,alpha,alpha_bc,psolver)
743722
call create_solver(hypre_maxiter,hypre_tol,hypre_solver_i,psolver)
744723
call setup_solver(psolver)
@@ -750,7 +729,6 @@ program snac
750729
hiu(:) = hi(:)
751730
if(is_bound(1,1)) hiu(:) = hiu(:)-[1,0,0]
752731
is_centered(:) = [.false.,.true.,.true.]
753-
is_symm_matrix_u = is_uniform_grid
754732
call init_bc_rhs(cbcvel(:,:,1),bcvel(:,:,1),dlu,is_bound,is_centered,lo,hiu,periods, &
755733
dxf,dxc,dyc,dyf,dzc,dzf,rhsu%x,rhsu%y,rhsu%z,bcu%x,bcu%y,bcu%z)
756734
#if defined(_FFT_X) || defined(_FFT_Y) || defined(_FFT_Z)
@@ -768,11 +746,11 @@ program snac
768746
if(periods(1) == 0) hiu_a(:) = hiu_a(:)-[1,0,0]
769747
#endif
770748
call init_n_2d_matrices(cbcvel(:,il:iu:iskip,1),bcvel(:,il:iu:iskip,1),dlu(:,il:iu:iskip), &
771-
is_symm_matrix_u,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
749+
.true.,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
772750
lo_a(idir),hiu_a(idir),lo_a(il:iu:iskip),hiu_a(il:iu:iskip),periods(il:iu:iskip), &
773751
dlu1_1,dlu1_2,dlu2_1,dlu2_2,alpha,alpha_bc,lambda_u_a,comms_fft,usolver_fft)
774752
#else
775-
call init_matrix_3d(cbcvel(:,:,1),bcvel(:,:,1),dlu,is_symm_matrix_u,is_bound,is_centered,lo,hiu,periods, &
753+
call init_matrix_3d(cbcvel(:,:,1),bcvel(:,:,1),dlu,.true.,is_bound,is_centered,lo,hiu,periods, &
776754
dxf,dxc,dyc,dyf,dzc,dzf,alpha,alpha_bc,usolver)
777755
#endif
778756
dlv = reshape([dxc_g(lo_g(1)-1),dxc_g(hi_g(1)), &
@@ -796,14 +774,13 @@ program snac
796774
lambda_v_a(:) = lambda_v(lo_s(idir)-lo(idir)+1:hi_s(idir)-lo(idir)+1)
797775
hiv_a(:) = hi_s(:)
798776
if(periods(2) == 0) hiv_a(:) = hiv_a(:)-[0,1,0]
799-
is_symm_matrix_v = is_uniform_grid
800777
#endif
801778
call init_n_2d_matrices(cbcvel(:,il:iu:iskip,2),bcvel(:,il:iu:iskip,2),dlv(:,il:iu:iskip), &
802-
is_symm_matrix_v,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
779+
.true.,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
803780
lo_a(idir),hiv_a(idir),lo_a(il:iu:iskip),hiv_a(il:iu:iskip),periods(il:iu:iskip), &
804781
dlv1_1,dlv1_2,dlv2_1,dlv2_2,alpha,alpha_bc,lambda_v_a,comms_fft,vsolver_fft)
805782
#else
806-
call init_matrix_3d(cbcvel(:,:,2),bcvel(:,:,2),dlv,is_symm_matrix_v,is_bound,is_centered,lo,hiv,periods, &
783+
call init_matrix_3d(cbcvel(:,:,2),bcvel(:,:,2),dlv,.true.,is_bound,is_centered,lo,hiv,periods, &
807784
dxc,dxf,dyf,dyc,dzc,dzf,alpha,alpha_bc,vsolver)
808785
#endif
809786
dlw = reshape([dxc_g(lo_g(1)-1),dxc_g(hi_g(1)), &
@@ -827,14 +804,13 @@ program snac
827804
lambda_w_a(:) = lambda_w(lo_s(idir)-lo(idir)+1:hi_s(idir)-lo(idir)+1)
828805
hiw_a(:) = hi_s(:)
829806
if(periods(3) == 0) hiw_a(:) = hiw_a(:)-[0,0,1]
830-
is_symm_matrix_w = is_uniform_grid
831807
#endif
832808
call init_n_2d_matrices(cbcvel(:,il:iu:iskip,3),bcvel(:,il:iu:iskip,3),dlw(:,il:iu:iskip), &
833-
is_symm_matrix_w,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
809+
.true.,is_bound_a(:,il:iu:iskip),is_centered(il:iu:iskip), &
834810
lo_a(idir),hiw_a(idir),lo_a(il:iu:iskip),hiw_a(il:iu:iskip),periods(il:iu:iskip), &
835811
dlw1_1,dlw1_2,dlw2_1,dlw2_2,alpha,alpha_bc,lambda_w_a,comms_fft,wsolver_fft)
836812
#else
837-
call init_matrix_3d(cbcvel(:,:,3),bcvel(:,:,3),dlw,is_symm_matrix_w,is_bound,is_centered,lo,hiw,periods, &
813+
call init_matrix_3d(cbcvel(:,:,3),bcvel(:,:,3),dlw,.true.,is_bound,is_centered,lo,hiw,periods, &
838814
dxc,dxf,dyc,dyf,dzf,dzc,alpha,alpha_bc,wsolver)
839815
#endif
840816
#endif

‎src/solver.f90‎

Lines changed: 42 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -241,29 +241,43 @@ subroutine init_matrix_3d(cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo,hi,p
241241
do j=lo(2),hi(2)
242242
do i=lo(1),hi(1)
243243
q = q + 1
244-
cxm = 1._rp/(dx1(i-1+qqq(1))*dx2(i))
245-
cxp = 1._rp/(dx1(i +qqq(1))*dx2(i))
246-
cym = 1._rp/(dy1(j-1+qqq(2))*dy2(j))
247-
cyp = 1._rp/(dy1(j +qqq(2))*dy2(j))
248-
czm = 1._rp/(dz1(k-1+qqq(3))*dz2(k))
249-
czp = 1._rp/(dz1(k +qqq(3))*dz2(k))
244+
cc = 0.
250245
#ifdef _FFT_X
251246
cxm = 0._rp
252247
cxp = 0._rp
248+
cym = 1._rp/(dy1(j-1+qqq(2)))*dz2(k)
249+
cyp = 1._rp/(dy1(j +qqq(2)))*dz2(k)
250+
czm = 1._rp/(dz1(k-1+qqq(3)))*dy2(j)
251+
czp = 1._rp/(dz1(k +qqq(3)))*dy2(j)
253252
qq = i - lo(1) + 1
253+
cc = lambda(qq)*dy2(j)*dz2(k)
254254
#elif _FFT_Y
255+
cxm = 1._rp/(dx1(i-1+qqq(1)))*dz2(k)
256+
cxp = 1._rp/(dx1(i +qqq(1)))*dz2(k)
255257
cym = 0._rp
256258
cyp = 0._rp
259+
czm = 1._rp/(dz1(k-1+qqq(3)))*dx2(i)
260+
czp = 1._rp/(dz1(k +qqq(3)))*dx2(i)
257261
qq = j - lo(2) + 1
262+
cc = lambda(qq)*dx2(i)*dz2(k)
258263
#elif _FFT_Z
264+
cxm = 1._rp/(dx1(i-1+qqq(1)))*dy2(j)
265+
cxp = 1._rp/(dx1(i +qqq(1)))*dy2(j)
266+
cym = 1._rp/(dy1(j-1+qqq(2)))*dx2(i)
267+
cyp = 1._rp/(dy1(j +qqq(2)))*dx2(i)
259268
czm = 0._rp
260269
czp = 0._rp
261270
qq = k - lo(3) + 1
271+
cc = lambda(qq)*dx2(i)*dy2(j)
272+
#else
273+
cxm = 1._rp/(dx1(i-1+qqq(1)))*dy2(j)*dz2(k)
274+
cxp = 1._rp/(dx1(i +qqq(1)))*dy2(j)*dz2(k)
275+
cym = 1._rp/(dy1(j-1+qqq(2)))*dx2(i)*dz2(k)
276+
cyp = 1._rp/(dy1(j +qqq(2)))*dx2(i)*dz2(k)
277+
czm = 1._rp/(dz1(k-1+qqq(3)))*dx2(i)*dy2(j)
278+
czp = 1._rp/(dz1(k +qqq(3)))*dx2(i)*dy2(j)
262279
#endif
263-
cc = -(cxm+cxp+cym+cyp+czm+czp) + alpha
264-
#if defined(_FFT_X) || defined(_FFT_Y) || defined(_FFT_Z)
265-
cc = cc + lambda(qq)
266-
#endif
280+
cc = cc - (cxm+cxp+cym+cyp+czm+czp) + alpha
267281
if(periods(1) == 0) then
268282
if(is_bound(0,1).and.i == lo(1)) then
269283
cc = cc + sgn(0,1)*cxm + alpha_bc(0,1)
@@ -358,14 +372,14 @@ subroutine create_solver(maxiter,maxerror,stype,asolver)
358372
call HYPRE_StructSMGCreate(asolver%comm_hypre,solver,ierr)
359373
call HYPRE_StructSMGSetMaxIter(solver,maxiter,ierr)
360374
call HYPRE_StructSMGSetTol(solver,maxerror,ierr)
361-
call hypre_structSMGsetLogging(solver,1,ierr)
362-
call HYPRE_StructSMGSetPrintLevel(solver,1,ierr)
375+
! call hypre_structSMGsetLogging(solver,1,ierr)
376+
! call HYPRE_StructSMGSetPrintLevel(solver,1,ierr)
363377
else if ( stype == HYPRESolverPFMG ) then
364378
call HYPRE_StructPFMGCreate(asolver%comm_hypre,solver,ierr)
365379
call HYPRE_StructPFMGSetMaxIter(solver,maxiter,ierr)
366380
call HYPRE_StructPFMGSetTol(solver,maxerror,ierr)
367-
call HYPRE_structPFMGsetLogging(solver,1,ierr)
368-
call HYPRE_StructPFMGSetPrintLevel(solver,1,ierr)
381+
! call HYPRE_structPFMGsetLogging(solver,1,ierr)
382+
! call HYPRE_StructPFMGSetPrintLevel(solver,1,ierr)
369383
call HYPRE_StructPFMGSetRelChange(solver,1,ierr)
370384
! Relaxiation Method: 2 is the fastest if symm matrix
371385
! 0: Jacobi
@@ -711,10 +725,10 @@ subroutine init_matrix_2d(cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo,hi,p
711725
q = 0
712726
do i2=lo(2),hi(2)
713727
do i1=lo(1),hi(1)
714-
c1m = 1._rp/(dl1_1(i1-1+qqq(1))*dl1_2(i1))
715-
c1p = 1._rp/(dl1_1(i1 +qqq(1))*dl1_2(i1))
716-
c2m = 1._rp/(dl2_1(i2-1+qqq(2))*dl2_2(i2))
717-
c2p = 1._rp/(dl2_1(i2 +qqq(2))*dl2_2(i2))
728+
c1m = 1._rp/(dl1_1(i1-1+qqq(1)))*dl2_2(i2)
729+
c1p = 1._rp/(dl1_1(i1 +qqq(1)))*dl2_2(i2)
730+
c2m = 1._rp/(dl2_1(i2-1+qqq(2)))*dl1_2(i1)
731+
c2p = 1._rp/(dl2_1(i2 +qqq(2)))*dl1_2(i1)
718732
cc = -(c1m+c1p+c2m+c2p) + alpha
719733
if(periods(1) == 0) then
720734
if(is_bound(0,1).and.i1 == lo(1)) then
@@ -928,16 +942,24 @@ subroutine init_n_2d_matrices(cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo_
928942
type(hypre_solver), intent(inout), dimension(:) :: asolver
929943
logical , intent(in ), optional, dimension(0:1,2) :: is_bound_outflow
930944
type(hypre_solver) :: asolver_aux
931-
integer :: i_out,q
945+
integer :: i_out,i1,i2,q,qq
932946
logical, dimension(0:1,2) :: is_bound_outflow_aux
947+
real(rp), dimension(product(hi(1:2)-lo(1:2)+1)) :: matvalues
933948
is_bound_outflow_aux(:,:) = .false.
934949
if( present(is_bound_outflow) ) is_bound_outflow_aux(:,:) = is_bound_outflow(:,:)
935950
do i_out=lo_out,hi_out
936951
q = i_out-lo_out+1
937952
asolver_aux = asolver(q)
938953
call init_matrix_2d(cbc,bc,dl,is_uniform_grid,is_bound,is_centered,lo,hi,periods, &
939954
dl1_1,dl1_2,dl2_1,dl2_2,alpha,alpha_bc,comm(q),asolver_aux,is_bound_outflow_aux)
940-
call add_constant_to_diagonal([lo(1),lo(2),1],[hi(1),hi(2),1],lambda(q),asolver_aux%mat)
955+
qq = 0
956+
do i2=lo(2),hi(2)
957+
do i1=lo(1),hi(1)
958+
qq = qq + 1
959+
matvalues(qq) = lambda(q)*dl1_2(i1)*dl2_2(i2)
960+
end do
961+
end do
962+
call HYPRE_StructMatrixAddToBoxValues(asolver_aux%mat,[lo(1),lo(2),1],[hi(1),hi(2),1],1,[0],matvalues,ierr)
941963
asolver(q) = asolver_aux
942964
end do
943965
end subroutine init_n_2d_matrices

0 commit comments

Comments
 (0)