From 8964cc5a1ddd70db3b927102728c01a4a6192bab Mon Sep 17 00:00:00 2001
From: oawile <oawile@7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9>
Date: Wed, 10 Mar 2010 05:25:08 +0000
Subject: [PATCH] - added hamjac_reinit_rusoo routines and updated module file,
but the DS needs to be still updated (-> #63)
git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@566 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9
---
src/ppm_hamjac_reinit_russo_3d.f | 332 ++++++++++++++++++++++++
src/ppm_hamjac_reinit_russo_step_3d.f | 360 ++++++++++++++++++++++++++
src/ppm_module_hamjac_reinit.f | 14 +
3 files changed, 706 insertions(+)
create mode 100644 src/ppm_hamjac_reinit_russo_3d.f
create mode 100644 src/ppm_hamjac_reinit_russo_step_3d.f
diff --git a/src/ppm_hamjac_reinit_russo_3d.f b/src/ppm_hamjac_reinit_russo_3d.f
new file mode 100644
index 0000000..9fbd1a5
--- /dev/null
+++ b/src/ppm_hamjac_reinit_russo_3d.f
@@ -0,0 +1,332 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_hamjac_reinit_3d
+ !-------------------------------------------------------------------------
+ !
+ ! Purpose : Solve Hamilton-Jacobi for Gowas reinit
+ !
+ ! Input :
+ !
+ ! Input/Output :
+ !
+ ! Output :
+ !
+ ! Remarks :
+ !
+ !
+ ! References :
+ !
+ ! Revisions :
+ !-------------------------------------------------------------------------
+ ! $Log: ppm_hamjac_reinit_3d.f,v $
+ ! Revision 1.1.1.1 2006/07/25 15:18:19 menahel
+ ! initial import
+ !
+ ! Revision 1.1 2005/07/25 00:34:02 ivos
+ ! Initial check-in.
+ !
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! Institute of Computational Science
+ ! ETH Zentrum, Hirschengraben 84
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+
+#if __MODE == __SCA
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_hamjac_reinit_russo_3ds (phi, phigrad, tol, maxstep, &
+ & topo_id, mesh_id, ghostsize, info, indx)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_hamjac_reinit_russo_3dd (phi, phigrad, tol, maxstep, &
+ & topo_id, mesh_id, ghostsize, info, indx)
+#endif
+#elif __MODE == __VEC
+#error VECTOR NOT IMPLEMENTED
+#endif
+
+ USE ppm_module_data
+ USE ppm_module_data_mesh
+ USE ppm_module_error
+ USE ppm_module_write
+ USE ppm_module_substart
+ USE ppm_module_alloc
+ USE ppm_module_substop
+ USE ppm_module_map
+ IMPLICIT NONE
+
+#if __KIND == __SINGLE_PRECISION
+ INTEGER, PARAMETER :: MK = ppm_kind_single
+#elif __KIND == __DOUBLE_PRECISION
+ INTEGER, PARAMETER :: MK = ppm_kind_double
+#endif
+
+#ifdef __MPI
+ INCLUDE 'mpif.h'
+#endif
+ !-----------------------------------------------------
+ ! Arguments
+ !-----------------------------------------------------
+ REAL(MK), DIMENSION(:,:,:,: ), POINTER :: phi
+ REAL(MK), DIMENSION(:,:,:,:,:), POINTER :: phigrad
+ INTEGER, INTENT(in) :: topo_id, mesh_id
+ INTEGER, DIMENSION(3), INTENT(in) :: ghostsize
+ INTEGER, INTENT(inout) :: info
+ INTEGER, INTENT(in) :: maxstep
+ REAL(mk), INTENT(in) :: tol
+ INTEGER,INTENT(IN),OPTIONAL :: indx
+ !-----------------------------------------------------
+ ! Aliases
+ !-----------------------------------------------------
+ INTEGER, DIMENSION(:), POINTER :: isublist
+ REAL(mk), DIMENSION(:,:,:,: ), POINTER :: phi0,phirhs,phiprev
+ INTEGER :: nsublist
+ INTEGER, DIMENSION(:,:), POINTER :: ndata
+ INTEGER :: topoid,meshid
+ REAL(mk), DIMENSION(:,:), POINTER :: min_phys, max_phys
+ INTEGER, DIMENSION(6) :: orgbcdef
+ INTEGER :: s2didx, mpi_prec
+ !-----------------------------------------------------
+ ! standard stuff
+ !-----------------------------------------------------
+ INTEGER :: isub,isubl,i,j,k,maptype,istep,iopt
+ INTEGER :: ldl(4), ldu(4), ndata_max(3)
+ REAL(mk) :: len_phys(3), dx(3)
+ REAL(mk) :: t0, res, gres,tau
+ CHARACTER(len=256) :: msg
+
+ CALL substart('ppm_hamjac_reinit_russo_3d',t0,info)
+
+ IF(PRESENT(indx)) THEN
+ s2didx=indx
+ ELSE
+ s2didx= -1
+ END IF
+
+#ifdef __MPI
+ IF (ppm_kind.EQ.ppm_kind_single) THEN
+ MPI_PREC = MPI_REAL
+ ELSE
+ MPI_PREC = MPI_DOUBLE_PRECISION
+ ENDIF
+#endif
+ !-----------------------------------------------------
+ ! Get the mesh data
+ !-----------------------------------------------------
+ topoid = ppm_internal_topoid(topo_id)
+ meshid = ppm_meshid(topoid)%internal(mesh_id)
+ nsublist = ppm_nsublist(topoid)
+ ndata => ppm_cart_mesh(meshid,topoid)%nnodes
+
+ isublist => ppm_isublist(:,topoid)
+#if __KIND == __SINGLE_PRECISION
+ min_phys => ppm_min_physs
+ max_phys => ppm_max_physs
+#elif __KIND == __DOUBLE_PRECISION
+ min_phys => ppm_min_physd
+ max_phys => ppm_max_physd
+#endif
+
+ len_phys(1) = max_phys(1,topoid) - min_phys(1,topoid)
+ len_phys(2) = max_phys(2,topoid) - min_phys(2,topoid)
+ len_phys(3) = max_phys(3,topoid) - min_phys(3,topoid)
+
+ dx(1) = len_phys(1)/REAL(ppm_cart_mesh(meshid,topoid)%nm(1)-1,mk)
+ dx(2) = len_phys(2)/REAL(ppm_cart_mesh(meshid,topoid)%nm(2)-1,mk)
+ dx(3) = len_phys(3)/REAL(ppm_cart_mesh(meshid,topoid)%nm(3)-1,mk)
+
+ ! timestep
+ tau = 0.25_mk*MINVAL(dx)
+
+ !-----------------------------------------------------
+ ! allocate temporary storage
+ !-----------------------------------------------------
+ ldl(1:3) = 1 - ghostsize(1:3); ldl(4) = 1
+ ndata_max(1) = MAXVAL(ndata(1,1:nsublist))
+ ndata_max(2) = MAXVAL(ndata(2,1:nsublist))
+ ndata_max(3) = MAXVAL(ndata(3,1:nsublist))
+ ldu(1) = ndata_max(1) + ghostsize(1)
+ ldu(2) = ndata_max(2) + ghostsize(2)
+ ldu(3) = ndata_max(3) + ghostsize(3)
+ ldu(4) = nsublist
+ iopt = ppm_param_alloc_fit
+ CALL ppm_alloc(phi0,ldl,ldu,iopt,info)
+ CALL ppm_alloc(phirhs,ldl,ldu,iopt,info)
+ CALL ppm_alloc(phiprev,ldl,ldu,iopt,info)
+ IF(info.NE.0) THEN
+ info = ppm_error_fatal
+ CALL ppm_error(ppm_err_alloc,'ppm_hamjac_reinit_russo_3d', &
+ & 'temp storage for hamjac',__LINE__,info)
+ GOTO 9999
+ END IF
+
+ ! fill phi0 with initial condition
+ DO isub=1,nsublist
+ isubl = isublist(isub)
+ DO k=1,ndata(3,isubl)
+ DO j=1,ndata(2,isubl)
+ DO i=1,ndata(1,isubl)
+ phi0(i,j,k,isub) = phi(i,j,k,isub)
+ phiprev(i,j,k,isub) = phi(i,j,k,isub)
+ phigrad(1:3,i,j,k,isub) = 0.0_mk
+ END DO
+ END DO
+ END DO
+ END DO
+
+ maptype = ppm_param_map_init
+ CALL ppm_map_field_ghost(phi0,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_ghost_get
+ CALL ppm_map_field_ghost(phi0,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_push
+ CALL ppm_map_field_ghost(phi0,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_send
+ CALL ppm_map_field_ghost(phi0,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_pop
+ CALL ppm_map_field_ghost(phi0,topo_id,mesh_id,ghostsize,maptype,info)
+
+ !-----------------------------------------------------
+ ! start time integration loop: using TVD RK3
+ !-----------------------------------------------------
+
+ DO istep=1,maxstep
+ !--- map the ghosts
+ maptype = ppm_param_map_init
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_ghost_get
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_push
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_send
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_pop
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+
+ ! TVDRK3 substep A
+
+ CALL ppm_hamjac_reinit_russo_step(phi,phi0,phirhs,phigrad,res,s2didx,topo_id,mesh_id&
+ &, ghostsize,info)
+
+ DO isub=1,nsublist
+ isubl = isublist(isub)
+ DO k=1,ndata(3,isubl); DO j=1,ndata(2,isubl);DO i=1,ndata(1,isubl)
+ phi(i,j,k,isub) = phi(i,j,k,isub) + tau*phirhs(i,j,k,isub)
+ END DO; END DO; END DO
+ END DO
+
+ !--- map the ghosts
+ maptype = ppm_param_map_init
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_ghost_get
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_push
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_send
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_pop
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+
+ ! TVDRK3 substep B
+
+ CALL ppm_hamjac_reinit_russo_step(phi,phi0,phirhs,phigrad,res,s2didx,topo_id,mesh_id&
+ &, ghostsize,info)
+
+ DO isub=1,nsublist
+ isubl = isublist(isub)
+ DO k=1,ndata(3,isubl); DO j=1,ndata(2,isubl);DO i=1,ndata(1,isubl)
+ phi(i,j,k,isub) = phi(i,j,k,isub) + tau*phirhs(i,j,k,isub)
+ phi(i,j,k,isub) = 0.25_mk*phiprev(i,j,k,isub) + 0.75_mk*phi(i,j,k,isub)
+ END DO; END DO; END DO
+ END DO
+
+ !--- map the ghosts
+ maptype = ppm_param_map_init
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_ghost_get
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_push
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_send
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_pop
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+
+ ! TVDRK3 substep C
+
+ CALL ppm_hamjac_reinit_russo_step(phi,phi0,phirhs,phigrad,res,s2didx,topo_id,mesh_id&
+ &, ghostsize,info)
+
+ DO isub=1,nsublist
+ isubl = isublist(isub)
+ DO k=1,ndata(3,isubl); DO j=1,ndata(2,isubl);DO i=1,ndata(1,isubl)
+ phi(i,j,k,isub) = phi(i,j,k,isub) + tau*phirhs(i,j,k,isub)
+ phi(i,j,k,isub) = (2.0_mk*phiprev(i,j,k,isub) + phi(i,j,k,isub))/3.0_mk
+ phiprev(i,j,k,isub) = phi(i,j,k,isub)
+ END DO; END DO; END DO
+ END DO
+#ifdef __MPI
+ CALL MPI_AllReduce(res,gres,1,mpi_prec,MPI_MIN,0,ppm_comm,info)
+ res = gres
+#endif
+ IF(res.LT.tol) GOTO 666
+
+ IF(ppm_rank.EQ.0.AND.MOD(istep,5).EQ.0) THEN
+ WRITE(msg,*) 'iteration #',istep,' res=',res
+ CALL ppm_write(ppm_Rank,'ppm_hamjac',msg,info)
+ END IF
+
+ !IF(res.LT.tol) GOTO 666 ! does not work in parallel with ghosting, because all
+ ! processors need to be in loop for ghosting to work
+ END DO
+
+ !info = ppm_error_warning
+ !CALL ppm_error(ppm_err_converge,'ppm_hamjac_reinit_russo_3d', &
+ ! & 'failed to reach target residual',__LINE__,info)
+
+666 CONTINUE
+ !--- map the ghosts for last time
+ maptype = ppm_param_map_init
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_ghost_get
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_push
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_send
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+ maptype = ppm_param_map_pop
+ CALL ppm_map_field_ghost(phi,topo_id,mesh_id,ghostsize,maptype,info)
+
+ IF(ppm_rank.EQ.0) THEN
+ WRITE(msg,*) 'ended after iteration #',istep,' res=',res
+ CALL ppm_write(ppm_Rank,'ppm_hamjac',msg,info)
+ END IF
+
+ iopt = ppm_param_dealloc
+ CALL ppm_alloc(phi0,ldl,ldu,iopt,info)
+ CALL ppm_alloc(phirhs,ldl,ldu,iopt,info)
+ CALL ppm_alloc(phiprev,ldl,ldu,iopt,info)
+ IF(info.NE.0) THEN
+ info = ppm_error_error
+ CALL ppm_error(ppm_err_dealloc,'ppm_hamjac_reinit_russo_3d', &
+ & 'temp storage for hamjac not freed',__LINE__,info)
+ GOTO 9999
+ END IF
+
+
+9999 CONTINUE
+
+ CALL substop('ppm_hamjac_reinit_russo_3d',t0,info)
+
+
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_hamjac_reinit_russo_3ds
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_hamjac_reinit_russo_3dd
+#endif
+
+
+
+
+
+
+
+
diff --git a/src/ppm_hamjac_reinit_russo_step_3d.f b/src/ppm_hamjac_reinit_russo_step_3d.f
new file mode 100644
index 0000000..2c0cbbd
--- /dev/null
+++ b/src/ppm_hamjac_reinit_russo_step_3d.f
@@ -0,0 +1,360 @@
+ !-------------------------------------------------------------------------
+ ! Subroutine : ppm_hamjac_reinit_step_3d
+ !-------------------------------------------------------------------------
+ !
+ ! Purpose : Solve Hamilton-Jacobi for Gowas reinit
+ !
+ ! Input :
+ !
+ ! Input/Output :
+ !
+ ! Output :
+ !
+ ! Remarks :
+ !
+ !
+ ! References :
+ !
+ ! Revisions :
+ !-------------------------------------------------------------------------
+ ! $Log: ppm_hamjac_reinit_step_3d.f,v $
+ ! Revision 1.1.1.1 2006/07/25 15:18:19 menahel
+ ! initial import
+ !
+ ! Revision 1.1 2005/07/25 00:34:04 ivos
+ ! Initial check-in.
+ !
+ !-------------------------------------------------------------------------
+ ! Parallel Particle Mesh Library (PPM)
+ ! Institute of Computational Science
+ ! ETH Zentrum, Hirschengraben 84
+ ! CH-8092 Zurich, Switzerland
+ !-------------------------------------------------------------------------
+
+
+#if __MODE == __SCA
+#if __KIND == __SINGLE_PRECISION
+ SUBROUTINE ppm_hamjac_reinit_russo_step_3ds (phi, phi0, phirhs,phigrad,res, s2didx,&
+ & topo_id, mesh_id, ghostsize, info)
+#elif __KIND == __DOUBLE_PRECISION
+ SUBROUTINE ppm_hamjac_reinit_russo_step_3dd (phi, phi0, phirhs,phigrad,res, s2didx,&
+ & topo_id, mesh_id, ghostsize, info)
+#endif
+#elif __MODE == __VEC
+#error VECTOR NOT IMPLEMENTED
+#endif
+
+ USE ppm_module_data
+ USE ppm_module_data_mesh
+ USE ppm_module_error
+ USE ppm_module_substart
+ USE ppm_module_substop
+
+ IMPLICIT NONE
+
+#if __KIND == __SINGLE_PRECISION
+ INTEGER, PARAMETER :: MK = ppm_kind_single
+#elif __KIND == __DOUBLE_PRECISION
+ INTEGER, PARAMETER :: MK = ppm_kind_double
+#endif
+
+ !-----------------------------------------------------
+ ! Arguments
+ !-----------------------------------------------------
+ REAL(MK), DIMENSION(:,:,:,: ), POINTER :: phi, phi0,phirhs
+ REAL(MK), DIMENSION(:,:,:,:,:), POINTER :: phigrad
+ INTEGER, INTENT(in) :: topo_id, mesh_id
+ INTEGER, DIMENSION(3), INTENT(in) :: ghostsize
+ INTEGER, INTENT(inout) :: info
+ REAL(mk),INTENT(out) :: res
+ INTEGER, INTENT(in) :: s2didx
+
+ !-----------------------------------------------------
+ ! Aliases
+ !-----------------------------------------------------
+ INTEGER, DIMENSION(:), POINTER :: isublist
+ INTEGER :: nsublist
+ INTEGER, DIMENSION(:,:), POINTER :: ndata
+ INTEGER :: topoid, meshid
+ REAL(mk), DIMENSION(:,:), POINTER :: min_phys, max_phys
+
+ !-----------------------------------------------------
+ ! standard stuff
+ !-----------------------------------------------------
+ INTEGER :: isub,isubl,i,j,k,kmin,kmax
+ REAL(mk) :: len_phys(3)
+ REAL(mk) :: delta_phis(3,3), delta(3), D
+ !-----------------------------------------------------
+ ! WENO stuff
+ !-----------------------------------------------------
+ REAL(mk) :: oneg(3), opos(3), closeeps, wenoeps, pbs
+ REAL(mk) :: laps(-1:1,3), rpos(3), rneg(3), dx(3), dxi(3)
+ REAL(mk) :: phip(3), phin(3), phimid(3), rms, dphi_dt
+ INTEGER :: ilap
+ INTEGER, PARAMETER, DIMENSION(3,3) :: offs &
+ & = RESHAPE((/2,1,0,1,0,-1,0,-1,-2/),(/3,3/))
+ REAL(mk) :: t0
+
+
+ CALL substart('ppm_hamjac_reinit_russo_step_3d',t0,info)
+
+ !-----------------------------------------------------
+ ! Get the mesh data
+ !-----------------------------------------------------
+ topoid = ppm_internal_topoid(topo_id)
+ meshid = ppm_meshid(topoid)%internal(mesh_id)
+ nsublist = ppm_nsublist(topoid)
+ ndata => ppm_cart_mesh(meshid,topoid)%nnodes
+ ! COMMENT Thu May 26 19:39:51 PDT 2005: experimental
+ isublist => ppm_isublist(:,topoid)
+#if __KIND == __SINGLE_PRECISION
+ min_phys => ppm_min_physs
+ max_phys => ppm_max_physs
+#elif __KIND == __DOUBLE_PRECISION
+ min_phys => ppm_min_physd
+ max_phys => ppm_max_physd
+#endif
+
+ len_phys(1) = max_phys(1,topoid) - min_phys(1,topoid)
+ len_phys(2) = max_phys(2,topoid) - min_phys(2,topoid)
+ len_phys(3) = max_phys(3,topoid) - min_phys(3,topoid)
+
+ dx(1) = len_phys(1)/REAL(ppm_cart_mesh(meshid,topoid)%nm(1)-1,mk)
+ dx(2) = len_phys(2)/REAL(ppm_cart_mesh(meshid,topoid)%nm(2)-1,mk)
+ dx(3) = len_phys(3)/REAL(ppm_cart_mesh(meshid,topoid)%nm(3)-1,mk)
+ dxi(1) = 1.0_mk/dx(1)
+ dxi(2) = 1.0_mk/dx(2)
+ dxi(3) = 1.0_mk/dx(3)
+
+ closeeps = 1.0e-6_mk
+ wenoeps = 1.0e-6_mk
+
+ rms = -HUGE(rms)
+
+ !-----------------------------------------------------
+ ! compute right hand side
+ !-----------------------------------------------------
+
+ DO isub=1,nsublist
+ isubl = isublist(isub)
+ IF(s2didx.LT.0) THEN
+ kmin = 1 + ghostsize(3)
+ kmax = ndata(3,isubl) - ghostsize(3)
+ ELSE
+ kmin = s2didx
+ kmax = s2didx
+ END IF
+ DO k=kmin,kmax
+ DO j=1,ndata(2,isubl)
+ DO i=1,ndata(1,isubl)
+
+ IF( phi0(i,j,k,isub)*phi0(i+1,j,k,isub).LT.0.0_mk.OR.&
+ & phi0(i,j,k,isub)*phi0(i-1,j,k,isub).LT.0.0_mk.OR.&
+ & phi0(i,j,k,isub)*phi0(i,j+1,k,isub).LT.0.0_mk.OR.&
+ & phi0(i,j,k,isub)*phi0(i,j-1,k,isub).LT.0.0_mk.OR.&
+ & phi0(i,j,k,isub)*phi0(i,j,k+1,isub).LT.0.0_mk.OR.&
+ & phi0(i,j,k,isub)*phi0(i,j,k-1,isub).LT.0.0_mk) THEN
+ ! near interface: apply russo smereka (2000) technique to
+ ! prevent interface from moving
+
+ delta_phis(1,1) = 0.5_mk*(phi0(i+1,j,k,isub)-phi0(i-1,j,k,isub))**2
+ delta_phis(1,2) = 0.5_mk*(phi0(i,j+1,k,isub)-phi0(i,j-1,k,isub))**2
+ delta_phis(1,3) = 0.5_mk*(phi0(i,j,k+1,isub)-phi0(i,j,k-1,isub))**2
+
+ delta_phis(2,1) = (phi0(i+1,j,k,isub)-phi0(i,j,k,isub))**2
+ delta_phis(2,2) = (phi0(i,j+1,k,isub)-phi0(i,j,k,isub))**2
+ delta_phis(2,3) = (phi0(i,j,k+1,isub)-phi0(i,j,k,isub))**2
+
+ delta_phis(3,1) = (phi0(i,j,k,isub)-phi0(i-1,j,k,isub))**2
+ delta_phis(3,2) = (phi0(i,j,k,isub)-phi0(i,j-1,k,isub))**2
+ delta_phis(3,3) = (phi0(i,j,k,isub)-phi0(i,j,k-1,isub))**2
+
+ delta(1) = MAX(MAX(delta_phis(1,1),delta_phis(2,1)), &
+ & MAX(delta_phis(3,1),closeeps))
+ delta(2) = MAX(MAX(delta_phis(1,2),delta_phis(2,2)), &
+ & MAX(delta_phis(3,2),closeeps))
+ delta(3) = MAX(MAX(delta_phis(1,3),delta_phis(2,3)), &
+ & MAX(delta_phis(3,3),closeeps))
+
+ !here assume dx(1)==dx(2) need to do the dx(1)~=dx(2) case
+ D = dx(1)*phi0(i,j,k,isub)/(SQRT(delta(1)+delta(2)+delta(3)))
+
+ IF(phi0(i,j,k,isub).GT.0.0_mk) THEN
+ phirhs(i,j,k,isub) = -dxi(1)*(ABS(phi(i,j,k,isub))-D)
+
+ ELSE IF(phi0(i,j,k,isub).LT.0.0_mk) THEN
+ phirhs(i,j,k,isub) = -dxi(1)*(-1.0_mk*ABS(phi(i,j,k,isub))-D)
+ ELSE
+ phirhs(i,j,k,isub) = 0.0_mk
+ END IF
+ ! gradient not yet correctly implemented
+ phigrad(1,i,j,k,isub) = dxi(1)*D
+ phigrad(2,i,j,k,isub) = dxi(2)*D
+ phigrad(3,i,j,k,isub) = dxi(3)*D
+
+ ELSE
+ ! away from interface: normal flux calculation
+
+ phimid(1) = phi(i+1,j,k,isub)-phi(i-1,j,k,isub)
+ phimid(2) = phi(i,j+1,k,isub)-phi(i,j-1,k,isub)
+ phimid(3) = phi(i,j,k+1,isub)-phi(i,j,k-1,isub)
+
+ laps(-1,1) = phi(i+2,j,k,isub) &
+ & -2.0_mk * phi(i+1,j,k,isub) &
+ & + phi(i,j,k,isub)
+ laps(-1,2) = phi(i,j+2,k,isub) &
+ & -2.0_mk * phi(i,j+1,k,isub) &
+ & + phi(i,j,k,isub)
+ laps(-1,3) = phi(i,j,k+2,isub) &
+ & -2.0_mk * phi(i,j,k+1,isub) &
+ & + phi(i,j,k,isub)
+
+ laps(0,1) = phi(i+1,j,k,isub) &
+ & -2.0_mk * phi(i,j,k,isub) &
+ & + phi(i-1,j,k,isub)
+ laps(0,2) = phi(i,j+1,k,isub) &
+ & -2.0_mk * phi(i,j,k,isub) &
+ & + phi(i,j-1,k,isub)
+ laps(0,3) = phi(i,j,k+1,isub) &
+ & -2.0_mk * phi(i,j,k,isub) &
+ & + phi(i,j,k-1,isub)
+
+ laps(1,1) = phi(i,j,k,isub) &
+ & -2.0_mk * phi(i-1,j,k,isub) &
+ & + phi(i-2,j,k,isub)
+ laps(1,2) = phi(i,j,k,isub) &
+ & -2.0_mk * phi(i,j-1,k,isub) &
+ & + phi(i,j-2,k,isub)
+ laps(1,3) = phi(i,j,k,isub) &
+ & -2.0_mk * phi(i,j,k-1,isub) &
+ & + phi(i,j,k-2,isub)
+
+
+ rpos(1) = (wenoeps + laps( 1,1)**2)/(wenoeps + laps(0,1)**2)
+ rneg(1) = (wenoeps + laps(-1,1)**2)/(wenoeps + laps(0,1)**2)
+ rpos(2) = (wenoeps + laps( 1,2)**2)/(wenoeps + laps(0,2)**2)
+ rneg(2) = (wenoeps + laps(-1,2)**2)/(wenoeps + laps(0,2)**2)
+ rpos(3) = (wenoeps + laps( 1,3)**2)/(wenoeps + laps(0,3)**2)
+ rneg(3) = (wenoeps + laps(-1,3)**2)/(wenoeps + laps(0,3)**2)
+
+ opos(1) = 1.0_mk/(1.0_mk+2.0_mk*rpos(1)**2)
+ opos(2) = 1.0_mk/(1.0_mk+2.0_mk*rpos(2)**2)
+ opos(3) = 1.0_mk/(1.0_mk+2.0_mk*rpos(3)**2)
+ oneg(1) = 1.0_mk/(1.0_mk+2.0_mk*rneg(1)**2)
+ oneg(2) = 1.0_mk/(1.0_mk+2.0_mk*rneg(2)**2)
+ oneg(3) = 1.0_mk/(1.0_mk+2.0_mk*rneg(3)**2)
+
+ phip(1) = 0.5_mk*(phimid(1) - &
+ & opos(1)*( &
+ & phi(i+2,j,k,isub) - &
+ & 3.0_mk*(phi(i+1,j,k,isub) - phi(i ,j,k,isub)) - &
+ & phi(i-1,j,k,isub)))*dxi(1)
+ phip(2) = 0.5_mk*(phimid(2) - &
+ & opos(2)*( &
+ & phi(i,j+2,k,isub) - &
+ & 3.0_mk*(phi(i,j+1,k,isub) - phi(i ,j,k,isub)) - &
+ & phi(i,j-1,k,isub)))*dxi(2)
+ phip(3) = 0.5_mk*(phimid(3) - &
+ & opos(3)*( &
+ & phi(i,j,k+2,isub) - &
+ & 3.0_mk*(phi(i,j,k+1,isub) - phi(i ,j,k,isub)) - &
+ & phi(i,j,k-1,isub)))*dxi(3)
+ phin(1) = 0.5_mk*(phimid(1) - &
+ & oneg(1)*( &
+ & phi(i+1,j,k,isub) - &
+ & 3.0_mk*(phi(i ,j,k,isub) - phi(i-1,j,k,isub)) - &
+ & phi(i-2,j,k,isub)))*dxi(1)
+ phin(2) = 0.5_mk*(phimid(2) - &
+ & oneg(2)*( &
+ & phi(i,j+1,k,isub) - &
+ & 3.0_mk*(phi(i,j ,k,isub) - phi(i,j-1,k,isub)) - &
+ & phi(i,j-2,k,isub)))*dxi(2)
+ phin(3) = 0.5_mk*(phimid(3) - &
+ & oneg(3)*( &
+ & phi(i,j,k+1,isub) - &
+ & 3.0_mk*(phi(i,j,k ,isub) - phi(i,j,k-1,isub)) - &
+ & phi(i,j,k-2,isub)))*dxi(3)
+
+ !--- collect
+ IF(phi0(i,j,k,isub).GT.0.0_mk) THEN
+ pbs = SQRT( &
+ & MAX(-MIN(phip(1),0.0_mk),MAX(phin(1),0.0_mk))**2+&
+ & MAX(-MIN(phip(2),0.0_mk),MAX(phin(2),0.0_mk))**2+&
+ & MAX(-MIN(phip(3),0.0_mk),MAX(phin(3),0.0_mk))**2)&
+ & - 1.0_mk
+
+ phirhs(i,j,k,isub) = -1.0_mk*pbs
+
+
+ ELSE IF(phi0(i,j,k,isub).LT.0.0_mk) THEN
+ pbs = SQRT( &
+ & MAX(MAX(phip(1),0.0_mk),-MIN(phin(1),0.0_mk))**2+&
+ & MAX(MAX(phip(2),0.0_mk),-MIN(phin(2),0.0_mk))**2+&
+ & MAX(MAX(phip(3),0.0_mk),-MIN(phin(3),0.0_mk))**2)&
+ & - 1.0_mk
+
+ phirhs(i,j,k,isub) = +1.0_mk*pbs
+
+ ELSE
+ pbs = 0.0_mk
+ phirhs(i,j,k,isub) = 0.0_mk
+ END IF
+ ! gradient not yet correctly implemented
+ phigrad(1,i,j,k,isub) = MAX(-MIN(phip(1),0.0_mk),MAX(phin(1),0.0_mk))
+ phigrad(2,i,j,k,isub) = MAX(-MIN(phip(2),0.0_mk),MAX(phin(2),0.0_mk))
+ phigrad(3,i,j,k,isub) = MAX(-MIN(phip(3),0.0_mk),MAX(phin(3),0.0_mk))
+
+
+ END IF
+ ! simple euler for now
+ !tphi(i,j,k,isub) = phi(i,j,k,isub) + reinit_tau * flux
+
+ rms = MAX(rms,ABS(phirhs(i,j,k,isub)))
+
+ END DO
+
+ END DO
+ END DO
+
+ END DO
+
+
+ IF(s2didx.GT.0) THEN
+ DO isub=1,nsublist
+ isubl = isublist(isub)
+ DO k=1,ndata(3,isubl)
+ DO j=1,ndata(2,isubl)
+ DO i=1,ndata(1,isubl)
+ phirhs(i,j,k,isub) = phirhs(i,j,s2didx,isub)
+ phigrad(1,i,j,k,isub) = phigrad(1,i,j,s2didx,isub)
+ phigrad(2,i,j,k,isub) = phigrad(2,i,j,s2didx,isub)
+ phigrad(3,i,j,k,isub) = phigrad(3,i,j,s2didx,isub)
+ END DO
+ END DO
+ END DO
+ END DO
+ END IF
+ res = rms
+
+ CALL substop('ppm_hamjac_reinit_russo_step_3d',t0,info)
+
+#if __KIND == __SINGLE_PRECISION
+ END SUBROUTINE ppm_hamjac_reinit_russo_step_3ds
+#elif __KIND == __DOUBLE_PRECISION
+ END SUBROUTINE ppm_hamjac_reinit_russo_step_3dd
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/ppm_module_hamjac_reinit.f b/src/ppm_module_hamjac_reinit.f
index fca1c34..9e983df 100644
--- a/src/ppm_module_hamjac_reinit.f
+++ b/src/ppm_module_hamjac_reinit.f
@@ -65,6 +65,16 @@
MODULE PROCEDURE ppm_hamjac_reinit_ref_3dd
END INTERFACE
+ INTERFACE ppm_hamjac_reinit_russo
+ MODULE PROCEDURE ppm_hamjac_reinit_russo_3ds
+ MODULE PROCEDURE ppm_hamjac_reinit_russo_3dd
+ END INTERFACE
+
+ INTERFACE ppm_hamjac_reinit_russo_step
+ MODULE PROCEDURE ppm_hamjac_reinit_russo_step_3ds
+ MODULE PROCEDURE ppm_hamjac_reinit_russo_step_3dd
+ END INTERFACE
+
CONTAINS
#define __DIME __3D
#define __MODE __SCA
@@ -72,11 +82,13 @@
! 3D SCA SINGLE
#include "ppm_hamjac_reinit_step_3d.f"
#include "ppm_hamjac_reinit_step_ref_3d.f"
+#include "ppm_hamjac_reinit_russo_step_3d.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
! 3D SCA SINGLE
#include "ppm_hamjac_reinit_step_3d.f"
#include "ppm_hamjac_reinit_step_ref_3d.f"
+#include "ppm_hamjac_reinit_russo_step_3d.f"
#undef __KIND
#undef __MODE
#undef __DIME
@@ -114,11 +126,13 @@
! 3D SCA SINGLE
#include "ppm_hamjac_reinit_3d.f"
#include "ppm_hamjac_reinit_ref_3d.f"
+#include "ppm_hamjac_reinit_russo_3d.f"
#undef __KIND
#define __KIND __DOUBLE_PRECISION
! 3D SCA DOUBLE
#include "ppm_hamjac_reinit_3d.f"
#include "ppm_hamjac_reinit_ref_3d.f"
+#include "ppm_hamjac_reinit_russo_3d.f"
#undef __KIND
#undef __MODE
#undef __DIME
--
GitLab