Merge pull request #424 from sanromd/dev-tfluct

tfluct as input to solver, following rp

examples/euler_1d/Makefile

@@ -2,17 +2,21 @@ F2PY = f2py
 
 SHARPCLAW = ../../src/pyclaw/sharpclaw
 
-ONE_D_SHARPCLAW_SOURCES = evec.f90 tfluct.f90 $(SHARPCLAW)/flux1.f90 $(SHARPCLAW)/ClawParams.f90 $(SHARPCLAW)/workspace.f90 $(SHARPCLAW)/weno.f90 $(SHARPCLAW)/reconstruct.f90
+ONE_D_SHARPCLAW_SOURCES = evec.f90 $(SHARPCLAW)/flux1.f90 $(SHARPCLAW)/ClawParams.f90 $(SHARPCLAW)/workspace.f90 $(SHARPCLAW)/weno.f90 $(SHARPCLAW)/reconstruct.f90
 
 all: 
-	make sharpclaw1.so
+    make sharpclaw1.so
+    make euler_tfluct.so
 
 sharpclaw1.so: $(ONE_D_SHARPCLAW_SOURCES) 
-	${F2PY} -m sharpclaw1 -c $(ONE_D_SHARPCLAW_SOURCES) 
+    ${F2PY} -m sharpclaw1 -c $(ONE_D_SHARPCLAW_SOURCES) 
+
+euler_tfluct.so: euler_tfluct.f90
+    f2py -m $(basename $(notdir $@)) -c $^
 
 clean:
-	rm -f *.o *.so *.pyc *.log
+    rm -f *.o *.so *.pyc *.log
 
 clobber: clean
-	rm -rf _output/
-	rm -rf _plots/
+    rm -rf _output/
+    rm -rf _plots/

examples/euler_1d/tfluct.f90 → examples/euler_1d/euler_tfluct.f90

@@ -1,7 +1,6 @@
 ! ============================================================================
-subroutine tfluct(ixyz,maxnx,num_eqn,num_waves,num_ghost,mx,ql,qr,auxl,auxr,adq)
+subroutine tfluct1(maxnx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,adq)
 ! ============================================================================
-!   
 !   "Internal" Riemann solver for the euler equations in 1D.
 !   The riemann problem is solved by assuming a discontinuity at the
 !   center of the i'th cell.
@@ -11,8 +10,8 @@ subroutine tfluct(ixyz,maxnx,num_eqn,num_waves,num_ghost,mx,ql,qr,auxl,auxr,adq)
 !             auxl contains the auxiliary vector at the left edge of each cell
 !             auxr contains the state vector at the right edge of each cell
 !             maxnx is the number of physical points (without ghost cells)
-!             num_ghost is the number of ghost cells
-!             num_eqn is the number of equations
+!             mbc is the number of ghost cells
+!             meqn is the number of equations
 !             ixyz is the dimension index
 !             mx is the size of the patch for the dimension corresponding
 !               to the value of ixyz
@@ -28,16 +27,21 @@ subroutine tfluct(ixyz,maxnx,num_eqn,num_waves,num_ghost,mx,ql,qr,auxl,auxr,adq)
 !             |_  ( gamma*q3 - 0.5(gamma-1)q2^2/q1 ) * q2/g1  _|
 
 
-    implicit double precision (a-h,o-z)
-    integer,          intent(in) :: maxnx, num_eqn, num_waves, num_ghost, mx, ixyz
-    double precision, intent(in) :: ql(num_eqn,1-num_ghost:maxnx+num_ghost)
-    double precision, intent(in) :: qr(num_eqn,1-num_ghost:maxnx+num_ghost)
-    double precision, intent(out) :: adq(num_eqn,1-num_ghost:maxnx+num_ghost)
+    implicit none
+    integer,          intent(in)  :: maxnx, mx, meqn, mwaves, maux, mbc
+    double precision, intent(in)  :: auxl(maux,1-mbc:maxnx+mbc)
+    double precision, intent(in)  :: auxr(maux,1-mbc:maxnx+mbc)
+    double precision, intent(in)  ::   ql(meqn,1-mbc:maxnx+mbc)
+    double precision, intent(in)  ::   qr(meqn,1-mbc:maxnx+mbc)
+
+    double precision, intent(out) :: adq(meqn,1-mbc:maxnx+mbc)
+
+    integer :: i
+    double precision :: gamma, gamma1
 
 !   local storage
 !   ---------------
     common /cparam/  gamma1
-    double precision :: gamma
 
     gamma = gamma1 + 1.d0
 
@@ -50,5 +54,4 @@ subroutine tfluct(ixyz,maxnx,num_eqn,num_waves,num_ghost,mx,ql,qr,auxl,auxr,adq)
     enddo
 
     return
-end subroutine tfluct
-
+end subroutine tfluct1

examples/euler_1d/setup.py

@@ -24,7 +24,8 @@ def configuration(parent_package='',top_path=None):
                                                             'reconstruct.f90','flux1.f90']]
 
     config.add_extension('sharpclaw1',
-                         ['evec.f90','tfluct.f90'] + sharpclaw_srcs)
+                         ['evec.f90'] + sharpclaw_srcs)
+    config.add_extension('euler_tfluct','euler_tfluct.f90')
 
     return config
 

examples/euler_1d/shocksine.py

@@ -38,7 +38,7 @@
 c = np.array([0., .3772689153313680, .7545378306627360, .7289856616121880, .6992261359316680])
 
 def setup(use_petsc=False,iplot=False,htmlplot=False,outdir='./_output',solver_type='sharpclaw',
-        kernel_language='Fortran',use_char_decomp=False):
+        kernel_language='Fortran',use_char_decomp=False,tfluct_solver=True):
 
     if use_petsc:
         import clawpack.petclaw as pyclaw
@@ -60,11 +60,25 @@ def setup(use_petsc=False,iplot=False,htmlplot=False,outdir='./_output',solver_t
             try:
                 import sharpclaw1               # Import custom Fortran code
                 solver.fmod = sharpclaw1
-                solver.tfluct_solver = True     # Use total fluctuation solver for efficiency
-                solver.lim_type = 2             # WENO reconstruction 
-                solver.char_decomp = 2          # characteristic-wise reconstruction
+                solver.tfluct_solver = tfluct_solver     # Use total fluctuation solver for efficiency
+                if solver.tfluct_solver:
+                    try:
+                        import euler_tfluct
+                        solver.tfluct = euler_tfluct
+                    except ImportError:
+                        import logging
+                        logger = logging.getLogger()
+                        logger.error('Unable to load tfluct solver, did you run make?')
+                        print 'Unable to load tfluct solver, did you run make?'
+                        raise
             except ImportError:
+                import logging
+                logger = logging.getLogger()
+                logger.error('Unable to load sharpclaw1 solver, did you run make?')
+                print 'Unable to load sharpclaw1 solver, did you run make?'
                 pass
+            solver.lim_type = 2             # WENO reconstruction
+            solver.char_decomp = 2          # characteristic-wise reconstruction
     else:
         solver = pyclaw.ClawSolver1D(rs)
 

examples/euler_1d/woodward_colella_blast.py

@@ -41,12 +41,13 @@
         import logging
         logger = logging.getLogger()
         logger.warn('unable to compile Fortran modules; some SharpClaw options will not be available for this example')
+        print 'unable to compile Fortran modules; some SharpClaw options will not be available for this example'
         raise
 
 gamma = 1.4 # Ratio of specific heats
 gamma1 = gamma - 1.
 
-def setup(use_petsc=False,outdir='./_output',solver_type='sharpclaw',kernel_language='Fortran'):
+def setup(use_petsc=False,outdir='./_output',solver_type='sharpclaw',kernel_language='Fortran',tfluct_solver=True):
 
     if use_petsc:
         import clawpack.petclaw as pyclaw
@@ -63,12 +64,22 @@ def setup(use_petsc=False,outdir='./_output',solver_type='sharpclaw',kernel_lang
         solver.time_integrator = 'SSP33'
         solver.cfl_max = 0.65
         solver.cfl_desired = 0.6
+        solver.tfluct_solver = tfluct_solver
+        if solver.tfluct_solver:
+            try:
+                import euler_tfluct
+                solver.tfluct = euler_tfluct
+            except ImportError:
+                import logging
+                logger = logging.getLogger()
+                logger.error('Unable to load tfluct solver, did you run make?')
+                print 'Unable to load tfluct solver, did you run make?'
+                raise
+        solver.lim_type = 1
+        solver.char_decomp = 2
         try:
             import sharpclaw1
             solver.fmod = sharpclaw1
-            solver.tfluct_solver = True
-            solver.lim_type = 1     # TVD reconstruction 
-            solver.char_decomp = 2  # characteristic-wise reconstructiong
         except ImportError:
             pass
     elif solver_type=='classic':
@@ -80,13 +91,13 @@ def setup(use_petsc=False,outdir='./_output',solver_type='sharpclaw',kernel_lang
     solver.bc_lower[0]=pyclaw.BC.wall
     solver.bc_upper[0]=pyclaw.BC.wall
 
-    mx=800;
+    mx = 800;
     x = pyclaw.Dimension('x',0.0,1.0,mx)
     domain = pyclaw.Domain([x])
     state = pyclaw.State(domain,solver.num_eqn)
 
-    state.problem_data['gamma']= gamma
-    state.problem_data['gamma1']= gamma1
+    state.problem_data['gamma'] = gamma
+    state.problem_data['gamma1'] = gamma1
     if kernel_language =='Python':
         state.problem_data['efix'] = False
 

src/pyclaw/sharpclaw/flux1.f90

@@ -1,5 +1,5 @@
 ! ===================================================================
-subroutine flux1(q1d,dq1d,aux,dt,cfl,t,ixyz,num_aux,num_eqn,mx,num_ghost,maxnx,rp)
+subroutine flux1(q1d,dq1d,aux,dt,cfl,t,ixyz,num_aux,num_eqn,mx,num_ghost,maxnx,rp,tfluct)
 ! ===================================================================
 !
 !     # Evaluate (delta t) * dq(t)/dt
@@ -36,7 +36,7 @@ subroutine flux1(q1d,dq1d,aux,dt,cfl,t,ixyz,num_aux,num_eqn,mx,num_ghost,maxnx,r
     implicit none
 
     ! Input (dummy) variables
-    external :: rp
+    external :: rp, tfluct
     integer, intent(in) :: num_aux, num_eqn, num_ghost, maxnx, mx, ixyz
     double precision, intent(in) :: q1d(num_eqn,1-num_ghost:mx+num_ghost)
     double precision, intent(inout) :: dq1d(num_eqn,1-num_ghost:maxnx+num_ghost)
@@ -149,8 +149,13 @@ subroutine flux1(q1d,dq1d,aux,dt,cfl,t,ixyz,num_aux,num_eqn,mx,num_ghost,maxnx,r
         ! and right state qr(i), and returns a total fluctuation in amdq2
         ! NOTE that here amdq2 is really a total fluctuation (should be
         ! called adq); we do it this way just to avoid declaring more storage
-        call tfluct(ixyz,maxnx,num_eqn,num_waves,num_ghost,mx,ql,qr, &
-                     aux,aux,amdq2)
+        if (num_dim.eq.1) then
+            call tfluct(maxnx,num_eqn,num_waves,num_aux,num_ghost,mx,&
+                        ql,qr,aux,aux,amdq2)
+        else
+            call tfluct(ixyz,maxnx,num_eqn,num_waves,num_aux,num_ghost,mx,&
+                        ql,qr,aux,aux,amdq2)
+        endif
 
         ! Modify q using fluctuations:
         ! Note this may not correspond to a conservative flux-differencing

src/pyclaw/sharpclaw/flux2.f90

@@ -1,5 +1,5 @@
 ! ==========================================================
-subroutine flux2(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,rpn2)
+subroutine flux2(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,rpn2,tfluct2)
 ! ==========================================================
 !    Evaluate (delta t) *dq/dt
 !    On entry, q should give the initial data for this step
@@ -10,7 +10,7 @@ subroutine flux2(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,rpn2)
     implicit none
 
     ! Input (dummy) variables
-    external :: rpn2
+    external :: rpn2, tfluct2
     integer, intent(in) :: num_aux, num_eqn, num_ghost, maxnx, mx, my
     double precision, intent(in) :: dt,t
     double precision, target, intent(in) :: q(num_eqn, 1-num_ghost:mx+num_ghost, 1-num_ghost:my+num_ghost)
@@ -25,6 +25,7 @@ subroutine flux2(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,rpn2)
 ! Dummy interface just so f2py doesn't complain:
 !f2py real(DP) x
 !f2py x=rpn2(x)
+!f2py x=tfluct2(x)
 
     ! Local variables
     integer :: i,j,m
@@ -45,7 +46,7 @@ subroutine flux2(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,rpn2)
         aux1d => aux(:,:,j)
 
         ! compute modification dq1d along this slice:
-        call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,1,num_aux,num_eqn,mx,num_ghost,maxnx,rpn2)
+        call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,1,num_aux,num_eqn,mx,num_ghost,maxnx,rpn2,tfluct2)
         cfl = dmax1(cfl,cfl1d)
 
         forall(i=1:mx, m=1:num_eqn)
@@ -62,7 +63,7 @@ subroutine flux2(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,rpn2)
         q1d => q(:,i,:)
         aux1d => aux(:,i,:)
 
-        call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,2,num_aux,num_eqn,my,num_ghost,maxnx,rpn2)
+        call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,2,num_aux,num_eqn,my,num_ghost,maxnx,rpn2,tfluct2)
         cfl = dmax1(cfl,cfl1d)
 
         forall(j=1:my, m=1:num_eqn)

src/pyclaw/sharpclaw/flux3.f90

@@ -1,5 +1,5 @@
 ! ==========================================================
-subroutine flux3(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,mz,rpn3)
+subroutine flux3(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,mz,rpn3,tfluct3)
 ! ==========================================================
 
     ! Evaluate (delta t) *dq/dt
@@ -13,7 +13,7 @@ subroutine flux3(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,mz,rpn3
     use ClawParams
     implicit none
 
-    external :: rpn3
+    external :: rpn3, tfluct3
     integer, intent(in) :: num_aux,num_eqn,num_ghost,maxnx,mx,my,mz
     double precision, intent(in) :: dt,t
     double precision, target, intent(in) :: q(num_eqn,1-num_ghost:mx+num_ghost,1-num_ghost:my+num_ghost,1-num_ghost:mz+num_ghost)
@@ -28,6 +28,7 @@ subroutine flux3(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,mz,rpn3
 ! Dummy interface just so f2py doesn't complain:
 !f2py real(DP) x
 !f2py x=rpn3(x)
+!f2py x=tfluct3(x)
 
     !local variables
     integer :: i,j,k,m
@@ -48,7 +49,7 @@ subroutine flux3(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,mz,rpn3
             aux1d => aux(:,:,j,k)
 
             ! compute modification dq1d along this slice:
-            call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,1,num_aux,num_eqn,mx,num_ghost,maxnx,rpn3)
+            call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,1,num_aux,num_eqn,mx,num_ghost,maxnx,rpn3,tfluct3)
             cfl = dmax1(cfl,cfl1d)
 
             forall(i=1:mx, m=1:num_eqn)
@@ -67,7 +68,7 @@ subroutine flux3(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,mz,rpn3
             q1d => q(:,i,:,k)
             aux1d => aux(:,i,:,k)
 
-            call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,2,num_aux,num_eqn,my,num_ghost,maxnx,rpn3)
+            call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,2,num_aux,num_eqn,my,num_ghost,maxnx,rpn3,tfluct3)
             cfl = dmax1(cfl,cfl1d)
 
             forall(j=1:my,m=1:num_eqn)
@@ -85,7 +86,7 @@ subroutine flux3(q,dq,aux,dt,cfl,t,num_aux,num_eqn,num_ghost,maxnx,mx,my,mz,rpn3
             q1d => q(:,i,j,:)
             aux1d => aux(:,i,j,:)
 
-            call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,3,num_aux,num_eqn,mz,num_ghost,maxnx,rpn3)
+            call flux1(q1d,dq1d,aux1d,dt,cfl1d,t,3,num_aux,num_eqn,mz,num_ghost,maxnx,rpn3,tfluct3)
             cfl = dmax1(cfl,cfl1d)
 
             forall(k=1:mz,m=1:num_eqn)

src/pyclaw/sharpclaw/setup.py

@@ -6,16 +6,16 @@ def configuration(parent_package='',top_path=None):
 
     config.add_extension('sharpclaw1',
                          ['ClawParams.f90','weno.f90','reconstruct.f90',
-                          'tfluct.f90','evec.f90','workspace.f90','flux1.f90'])
+                          'evec.f90','workspace.f90','flux1.f90'])
 
     config.add_extension('sharpclaw2',
                          ['ClawParams.f90','weno.f90','reconstruct.f90',
-                          'tfluct.f90','evec.f90','workspace.f90','flux2.f90',
+                          'evec.f90','workspace.f90','flux2.f90',
                           'flux1.f90'])
 
     config.add_extension('sharpclaw3',
                          ['ClawParams.f90','weno.f90','reconstruct.f90',
-                          'tfluct.f90','evec.f90','workspace.f90','flux3.f90',
+                          'evec.f90','workspace.f90','flux3.f90',
                           'flux1.f90'])
 
     return config