diff --git a/euler_2d_shockbubble_amrclaw/Makefile b/euler/shockbubble/2d_axisym/Makefile similarity index 97% rename from euler_2d_shockbubble_amrclaw/Makefile rename to euler/shockbubble/2d_axisym/Makefile index 507c36a4..57287a0c 100644 --- a/euler_2d_shockbubble_amrclaw/Makefile +++ b/euler/shockbubble/2d_axisym/Makefile @@ -56,6 +56,7 @@ SOURCES = \ $(AMRLIB)/advanc.f \ $(AMRLIB)/bound.f90 \ $(AMRLIB)/stepgrid.f \ + $(AMRLIB)/stepgrid_dimSplit.f \ $(AMRLIB)/cellave.f \ $(AMRLIB)/fss.f \ $(AMRLIB)/zeroin.f \ @@ -126,7 +127,10 @@ SOURCES = \ $(AMRLIB)/icall.f \ $(AMRLIB)/preicall.f \ $(AMRLIB)/step2.f90 \ + $(AMRLIB)/step2x.f90 \ + $(AMRLIB)/step2y.f90 \ $(AMRLIB)/flux2.f \ + $(AMRLIB)/flux2_dimSplit.f \ $(AMRLIB)/inlinelimiter.f \ $(AMRLIB)/cstore.f \ $(AMRLIB)/saveqc.f \ diff --git a/euler/shockbubble/2d_axisym/README.rst b/euler/shockbubble/2d_axisym/README.rst new file mode 100644 index 00000000..5867011d --- /dev/null +++ b/euler/shockbubble/2d_axisym/README.rst @@ -0,0 +1,28 @@ + +.. _apps_euler_shockbubble_2d_axisym: + +2D Axisymmetric shock-bubble intereaction +========================================== + + +Shock-bubble interaction solved with 2-dimensional axisymmetric Euler +equations using AMR. + +A circular (i.e. spherical) bubble of gas is hit by a shock wave. +The same gas is inside and outside the bubble (ideal gas with gamma = 1.4). +The gas inside the bubble has a lower density but the same pressure +initially as the gas outside. Parameters are specified in `setrun.py`. + +A passive tracer is also advected to show the motion of the gas originally inside +the bubble. + +Version history: +---------------- + +- This version works with Clawpack 5.3.0 +- 28 Dec 2014: + + - Updated `Makefile` to include dimensional splitting + option introduced in 5.3. + - Moved from `apps/euler_2d_shockbubble` to `apps/euler/shockbubble/2d_axisym` + diff --git a/euler_2d_shockbubble_amrclaw/fdisc.f b/euler/shockbubble/2d_axisym/fdisc.f similarity index 100% rename from euler_2d_shockbubble_amrclaw/fdisc.f rename to euler/shockbubble/2d_axisym/fdisc.f diff --git a/euler_2d_shockbubble_amrclaw/qinit.f b/euler/shockbubble/2d_axisym/qinit.f similarity index 100% rename from euler_2d_shockbubble_amrclaw/qinit.f rename to euler/shockbubble/2d_axisym/qinit.f diff --git a/euler_2d_shockbubble_amrclaw/setaux.f b/euler/shockbubble/2d_axisym/setaux.f similarity index 100% rename from euler_2d_shockbubble_amrclaw/setaux.f rename to euler/shockbubble/2d_axisym/setaux.f diff --git a/euler_2d_shockbubble_amrclaw/setplot.py b/euler/shockbubble/2d_axisym/setplot.py similarity index 99% rename from euler_2d_shockbubble_amrclaw/setplot.py rename to euler/shockbubble/2d_axisym/setplot.py index a66e02f6..1ed0f1cc 100644 --- a/euler_2d_shockbubble_amrclaw/setplot.py +++ b/euler/shockbubble/2d_axisym/setplot.py @@ -119,6 +119,7 @@ def aa(current_data): plotfigure = plotdata.new_plotfigure(name='u-velocity', figno=3) + plotfigure.show = False plotfigure.kwargs = {'figsize':(16,5)} # Set up for axes in this figure: diff --git a/euler_2d_shockbubble_amrclaw/setprob.f b/euler/shockbubble/2d_axisym/setprob.f similarity index 100% rename from euler_2d_shockbubble_amrclaw/setprob.f rename to euler/shockbubble/2d_axisym/setprob.f diff --git a/euler_2d_shockbubble_amrclaw/setrun.py b/euler/shockbubble/2d_axisym/setrun.py similarity index 99% rename from euler_2d_shockbubble_amrclaw/setrun.py rename to euler/shockbubble/2d_axisym/setrun.py index efdce8d1..fa8075f3 100644 --- a/euler_2d_shockbubble_amrclaw/setrun.py +++ b/euler/shockbubble/2d_axisym/setrun.py @@ -182,7 +182,7 @@ def setrun(claw_pkg='amrclaw'): # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters clawdata.order = 2 - # Use dimensional splitting? (not yet available for AMR) + # Use dimensional splitting? clawdata.dimensional_split = 'unsplit' # For unsplit method, transverse_waves can be diff --git a/euler_2d_shockbubble_amrclaw/src1d.f b/euler/shockbubble/2d_axisym/src1d.f similarity index 100% rename from euler_2d_shockbubble_amrclaw/src1d.f rename to euler/shockbubble/2d_axisym/src1d.f diff --git a/euler_2d_shockbubble_amrclaw/src2.f b/euler/shockbubble/2d_axisym/src2.f similarity index 100% rename from euler_2d_shockbubble_amrclaw/src2.f rename to euler/shockbubble/2d_axisym/src2.f diff --git a/euler/shockbubble/3d/Makefile b/euler/shockbubble/3d/Makefile new file mode 100644 index 00000000..c9eff0ef --- /dev/null +++ b/euler/shockbubble/3d/Makefile @@ -0,0 +1,144 @@ + +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = amrclaw # Clawpack package to use +EXE = xamr # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +OVERWRITE ?= True # False ==> make a copy of OUTDIR first +RESTART ?= False # Should = clawdata.restart in setrun + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# List of sources for this program: +# --------------------------------- + +AMRLIB = $(CLAW)/amrclaw/src/3d + +MODULES = \ + $(AMRLIB)/amr_module.f90 +#$(AMRLIB)/gauges_module.f90 \ +#$(AMRLIB)/regions_module.f90 \ + +SOURCES = \ + qinit.f \ + setprob.f \ + cellave.f \ + fdisc.f \ + $(CLAW)/riemann/src/rpn3_euler.f90 \ + $(CLAW)/riemann/src/rpt3_euler.f90 \ + $(CLAW)/riemann/src/rptt3_euler.f90 \ + $(AMRLIB)/amr3.f90 \ + $(AMRLIB)/setaux.f90 \ + $(AMRLIB)/bc3amr.f \ + $(AMRLIB)/opendatafile.f \ + $(AMRLIB)/b4step3.f90 \ + $(AMRLIB)/qad.f \ + $(AMRLIB)/src3.f90 \ + $(AMRLIB)/src1d.f90 \ + $(AMRLIB)/advanc.f \ + $(AMRLIB)/bound.f \ + $(AMRLIB)/stepgrid.f \ + $(AMRLIB)/stepgrid_dimSplit.f \ + $(AMRLIB)/auxcoarsen.f \ + $(AMRLIB)/fixcapaq.f \ + $(AMRLIB)/estdt.f \ + $(AMRLIB)/igetsp.f \ + $(AMRLIB)/reclam.f \ + $(AMRLIB)/birect.f \ + $(AMRLIB)/cleanup.f \ + $(AMRLIB)/colate.f \ + $(AMRLIB)/errest.f \ + $(AMRLIB)/flag2refine.f \ + $(AMRLIB)/allowflag.f \ + $(AMRLIB)/bufnst.f \ + $(AMRLIB)/spest.f \ + $(AMRLIB)/errf1.f \ + $(AMRLIB)/gfixup.f \ + $(AMRLIB)/filval.f \ + $(AMRLIB)/filpatch.f \ + $(AMRLIB)/prefilp.f \ + $(AMRLIB)/flglvl.f \ + $(AMRLIB)/fluxad.f \ + $(AMRLIB)/fluxsv.f \ + $(AMRLIB)/ginit.f \ + $(AMRLIB)/grdfit.f \ + $(AMRLIB)/intfil.f \ + $(AMRLIB)/moment.f \ + $(AMRLIB)/nestck.f \ + $(AMRLIB)/prepc.f \ + $(AMRLIB)/prepf.f \ + $(AMRLIB)/projec.f \ + $(AMRLIB)/signs.f \ + $(AMRLIB)/findcut.f \ + $(AMRLIB)/smartbis.f \ + $(AMRLIB)/putnod.f \ + $(AMRLIB)/putsp.f \ + $(AMRLIB)/regrid.f \ + $(AMRLIB)/setgrd.f \ + $(AMRLIB)/setuse.f \ + $(AMRLIB)/stst1.f \ + $(AMRLIB)/tick.f \ + $(AMRLIB)/trimbd.f \ + $(AMRLIB)/update.f \ + $(AMRLIB)/nodget.f \ + $(AMRLIB)/upbnd.f \ + $(AMRLIB)/basic.f \ + $(AMRLIB)/outval.f \ + $(AMRLIB)/copysol.f \ + $(AMRLIB)/outvar.f \ + $(AMRLIB)/outmsh.f \ + $(AMRLIB)/outtre.f \ + $(AMRLIB)/domain.f \ + $(AMRLIB)/setflags.f \ + $(AMRLIB)/shiftset.f \ + $(AMRLIB)/conck.f \ + $(AMRLIB)/domshrink.f \ + $(AMRLIB)/domprep.f \ + $(AMRLIB)/domup.f \ + $(AMRLIB)/domcopy.f \ + $(AMRLIB)/coarsen.f \ + $(AMRLIB)/intcopy.f \ + $(AMRLIB)/preintcopy.f \ + $(AMRLIB)/icall.f \ + $(AMRLIB)/preicall.f \ + $(AMRLIB)/step3.f \ + $(AMRLIB)/step3x.f \ + $(AMRLIB)/step3y.f \ + $(AMRLIB)/step3z.f \ + $(AMRLIB)/flux3.f \ + $(AMRLIB)/flux3_dimSplit.f \ + $(AMRLIB)/limiter.f \ + $(AMRLIB)/philim.f \ + $(AMRLIB)/cstore.f \ + $(AMRLIB)/saveqc.f \ + $(AMRLIB)/valout.f \ + $(AMRLIB)/check.f \ + $(AMRLIB)/restrt.f \ + $(AMRLIB)/quick_sort1.f \ + $(AMRLIB)/init_alloc.f90 \ + $(AMRLIB)/resize_alloc.f90 \ + $(AMRLIB)/restrt_alloc.f90 + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + diff --git a/euler/shockbubble/3d/README.rst b/euler/shockbubble/3d/README.rst new file mode 100644 index 00000000..fe2f73c6 --- /dev/null +++ b/euler/shockbubble/3d/README.rst @@ -0,0 +1,29 @@ + +.. _apps_euler_shockbubble_3d: + +3D shock-bubble intereaction +============================== + + +Shock-bubble interaction solved with 3-dimensional axisymmetric Euler +equations using AMR. + +A spherical bubble of gas is hit by a shock wave. +The same gas is inside and outside the bubble (ideal gas with gamma = 1.4). +The gas inside the bubble has a lower density but the same pressure +initially as the gas outside. Parameters are specified in `setrun.py`. + +A passive tracer is also advected to show the motion of the gas originally inside +the bubble. + +Version history: +---------------- + +- This version works with Clawpack 5.3.0 +- 28 Dec 2014: + + - Updated `Makefile` to include dimensional splitting + option introduced in 5.3. + - Moved from `amrclaw/examples/euler_3d_shockbubble` to + `apps/euler/shockbubble/3d` + diff --git a/euler/shockbubble/3d/addaxes.m b/euler/shockbubble/3d/addaxes.m new file mode 100644 index 00000000..0207e32c --- /dev/null +++ b/euler/shockbubble/3d/addaxes.m @@ -0,0 +1,7 @@ +% add axes to 3d plots if grid is turned off +hold on +plot3([0 1.2],[.5 .5],[0 0],'k') +plot3([1.2 1.2],[.5 .5],[0 .5],'k') +plot3([1.2 1.2],[0 .5],[0 0],'k') +hold off + diff --git a/euler/shockbubble/3d/afterframe.m b/euler/shockbubble/3d/afterframe.m new file mode 100644 index 00000000..dcd9da4f --- /dev/null +++ b/euler/shockbubble/3d/afterframe.m @@ -0,0 +1,28 @@ +if PlotType ~= 4 + axis([0 2 0 0.5 0 0.5]); + daspect([1 1 1]); +end + +if PlotType==1 + yrbcolormap + caxis([0,3]) +end + +if PlotType==3 + set(gca,'box','on'); +end +camlight; + +% showcubes; +% setcubecolor('k',1); +% setcubecolor('b',2); +% setcubecolor('r',3); + +if PlotType==5 + camlight left; + grid off; +end + +shg +colorbar; +clear afterframe diff --git a/euler/shockbubble/3d/cellave.f b/euler/shockbubble/3d/cellave.f new file mode 100644 index 00000000..205a334f --- /dev/null +++ b/euler/shockbubble/3d/cellave.f @@ -0,0 +1,238 @@ +c +c +c +c +c ================================================= + subroutine cellave(xlow,ylow,dx,dy,wl) +c ================================================= + implicit double precision (a-h,o-z) + external fss + logical fl(5),alll,allr + dimension x(10),y(10),xx(5),yy(5) + common/fsscorn/ xc0,yc0,xc1,yc1 +c +c # compute wl, fraction of cell that lies in left state. +c # For initial data with two states ql and qr separated by a +c # discontinuity. The curve along which the discontinuity lies is +c # specified by the function fdisc, which should return a value that +c # is negative on the side where ql lies and positive on the qr side. +c +c # xlow,ylow is the coordinate of the lower left corner of the cell. +c # dx, dy are grid spacing in x and y. +c + xx(1) = xlow + xx(2) = xlow + xx(3) = xlow+dx + xx(4) = xlow+dx + xx(5) = xx(1) + yy(1) = ylow + yy(2) = ylow+dy + yy(3) = ylow+dy + yy(4) = ylow + yy(5) = yy(1) + alll = .true. + allr = .true. +c + do 20 i=1,4 + fl(i) = fdisc(xx(i),yy(i)) .lt. 0.d0 + alll = alll .and. fl(i) + allr = allr .and. (.not. fl(i)) + 20 continue + fl(5) = fl(1) +c + if (alll) then + wl = 1.d0 + return + endif + if (allr) then + wl = 0.d0 + return + endif +c + iv = 0 + do 40 i=1,4 + if (fl(i)) then + iv = iv+1 + x(iv) = xx(i) + y(iv) = yy(i) + endif + if (fl(i).neqv.fl(i+1)) then + iv = iv+1 + xc0 = xx(i) + yc0 = yy(i) + xc1 = xx(i+1) + yc1 = yy(i+1) + ss = zeroin(0.d0, 1.d0, fss, 1d-8) +c write(27,*) 'xc,yc,ss:',xc0,yc0,xc1,yc1,ss + x(iv) = xx(i) + ss*(xx(i+1)-xx(i)) + y(iv) = yy(i) + ss*(yy(i+1)-yy(i)) + endif + 40 continue +c +c # compute area: +c + if (iv.eq.0) then + wl = 0.d0 + return + endif +c + x(iv+1) = x(1) + y(iv+1) = y(1) + area = 0.d0 + do 50 i=1,iv + area = area + .5d0*(y(i)+y(i+1))*(x(i+1)-x(i)) +c write(27,*) ' x,y:',x(i),y(i) + 50 continue +c + wl = area / (dx*dy) +c write(27,*) 'area,wl:',area,wl +c + return + end +c +c +c +c +c +c ================================================= + function fss(s) +c ================================================= + implicit double precision (a-h,o-z) + common/fsscorn/ xc0,yc0,xc1,yc1 +c +c # compute fdisc at distance s between corners (xc0,yc0) and (xc1,yc1) +c + x = xc0 + s*(xc1-xc0) + y = yc0 + s*(yc1-yc0) + fss = fdisc(x,y) + return + end +c +c +c +c ================================================= + function zeroin(ax,bx,f,tol) +c ================================================= + implicit double precision (a-h,o-z) + external f +c +c a zero of the function f(x) is computed in the interval ax,bx . +c (Standard routine from netlib) +c +c input.. +c +c ax left endpoint of initial interval +c bx right endpoint of initial interval +c f function subprogram which evaluates f(x) for any x in +c the interval ax,bx +c tol desired length of the interval of uncertainty of the +c final result ( .ge. 0.d0) +c +c +c output.. +c +c zeroin abcissa approximating a zero of f in the interval ax,bx +c +c +c it is assumed that f(ax) and f(bx) have opposite signs +c without a check. zeroin returns a zero x in the given interval +c ax,bx to within a tolerance 4*macheps*dabs(x) + tol, where macheps +c is the relative machine precision. +c this function subprogram is a slightly modified translation of +c the algol 60 procedure zero given in richard brent, algorithms for +c minimization without derivatives, prentice - hall, inc. (1973). +c +c +c +c compute eps, the relative machine precision +c + eps = 1.d0 + 10 eps = eps/2.d0 + tol1 = 1.d0 + eps + if (tol1 .gt. 1.d0) go to 10 +c +c initialization +c + a = ax + b = bx + fa = f(a) + fb = f(b) +c +c begin step +c + 20 c = a + fc = fa + d = b - a + e = d + 30 if (dabs(fc) .ge. dabs(fb)) go to 40 + a = b + b = c + c = a + fa = fb + fb = fc + fc = fa +c +c convergence test +c + 40 tol1 = 2.d0*eps*dabs(b) + 0.5*tol + xm = .5*(c - b) + if (dabs(xm) .le. tol1) go to 90 + if (fb .eq. 0.d0) go to 90 +c +c is bisection necessary +c + if (dabs(e) .lt. tol1) go to 70 + if (dabs(fa) .le. dabs(fb)) go to 70 +c +c is quadratic interpolation possible +c + if (a .ne. c) go to 50 +c +c linear interpolation +c + s = fb/fa + p = 2.d0*xm*s + q = 1.d0 - s + go to 60 +c +c inverse quadratic interpolation +c + 50 q = fa/fc + r = fb/fc + s = fb/fa + p = s*(2.d0*xm*q*(q - r) - (b - a)*(r - 1.d0)) + q = (q - 1.d0)*(r - 1.d0)*(s - 1.d0) +c +c adjust signs +c + 60 if (p .gt. 0.d0) q = -q + p = dabs(p) +c +c is interpolation acceptable +c + if ((2.d0*p) .ge. (3.d0*xm*q - dabs(tol1*q))) go to 70 + if (p .ge. dabs(0.5*e*q)) go to 70 + e = d + d = p/q + go to 80 +c +c bisection +c + 70 d = xm + e = d +c +c complete step +c + 80 a = b + fa = fb + if (dabs(d) .gt. tol1) b = b + d + if (dabs(d) .le. tol1) b = b + dsign(tol1, xm) + fb = f(b) + if ((fb*(fc/dabs(fc))) .gt. 0.d0) go to 20 + go to 30 +c +c done +c + 90 zeroin = b + return + end diff --git a/euler/shockbubble/3d/fdisc.f b/euler/shockbubble/3d/fdisc.f new file mode 100644 index 00000000..be189648 --- /dev/null +++ b/euler/shockbubble/3d/fdisc.f @@ -0,0 +1,22 @@ +c +c ================================================= + function fdisc(x,y) +c ================================================= + implicit double precision (a-h,o-z) + common/cdisc/ x0,y0, r0, radius +c +c # for computing cell averages for initial data that has a +c # discontinuity along some curve. fdisc should be negative to the +c # left of the curve and positive to the right +c # idisc specifies the nature of the discontinuity for two +c # particular cases (a straight line and circle) but this routine +c # can be modified for any other curve. +c + +c # circle of radius r0: + + + fdisc = (x-x0)**2 + (y-y0)**2 - radius**2 +c + return + end diff --git a/euler/shockbubble/3d/qinit.f b/euler/shockbubble/3d/qinit.f new file mode 100644 index 00000000..d59fc261 --- /dev/null +++ b/euler/shockbubble/3d/qinit.f @@ -0,0 +1,70 @@ +c +c ===================================================== + subroutine qinit(meqn,mbc,mx,my,mz,xlower,ylower,zlower, + & dx,dy,dz,q,maux,aux) +c ===================================================== +c +c # Set initial conditions for q. +c + implicit double precision (a-h,o-z) + dimension q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc) + dimension aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc,1-mbc:mz+mbc) + + + double precision xc(1-mbc:mx+mbc) + double precision xe(1-mbc:mx+mbc) + + double precision yc(1-mbc:my+mbc) + double precision ye(1-mbc:my+mbc) + + double precision zc(1-mbc:mz+mbc) + double precision ze(1-mbc:mz+mbc) + + common /cparam/ gamma + common /comic/ qin(5),qout(5) + common /cominf/ rinf,vinf,einf + common /cdisc/ x0,y0,r0,radius + + do i = 1-mbc,mx+mbc + xe(i) = xlower + (i - 1 )*dx + xc(i) = xlower + (i - 0.5)*dx + enddo + + do j = 1-mbc,my+mbc + ye(j) = ylower + (j - 1.0)*dy + yc(j) = ylower + (j - 0.5)*dy + enddo + + do k = 1-mbc,mz+mbc + ze(k) = zlower + (k - 1.0)*dz + zc(k) = zlower + (k - 0.5)*dz + enddo + + + do k = 1,mz + do j = 1,my + do i = 1,mx + if (zc(k) <= r0) then + radius = dsqrt(r0**2 - zc(k)**2) + call cellave(xe(i),ye(j),dx,dy,win) + else + win = 0.0 + endif + do m=1,meqn + q(m,i,j,k) = win*qin(m) + (1.d0-win)*qout(m) + enddo + + if (xe(i) < 0.2) then +c # behind shock: + q(1,i,j,k) = rinf + q(2,i,j,k) = rinf*vinf + q(3,i,j,k) = 0.d0 + q(4,i,j,k) = 0.d0 + q(5,i,j,k) = einf + endif + enddo + enddo + enddo + + return + end diff --git a/euler/shockbubble/3d/setplot.py b/euler/shockbubble/3d/setplot.py new file mode 100644 index 00000000..df25deaf --- /dev/null +++ b/euler/shockbubble/3d/setplot.py @@ -0,0 +1,47 @@ + +""" +Dummy setplot does nothing currently in 3d. +Still need to implement Python plotting in 3d. + +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + + +#-------------------------- +def setplot(plotdata): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of clawpack.visclaw.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + plotdata.clearfigures() # clear any old figures,axes,items data + + print "**** Python plotting tools not yet implemented in 3d" + print "**** No plots will be generated." + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via clawpack.visclaw.frametools.printframes: + + plotdata.printfigs = False # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = [] # list of frames to print + plotdata.print_fignos = [] # list of figures to print + plotdata.html = False # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.html_movie = 'JSAnimation' # new style, or "4.x" for old style + plotdata.latex = False # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + + return plotdata + diff --git a/euler/shockbubble/3d/setplot3.m b/euler/shockbubble/3d/setplot3.m new file mode 100644 index 00000000..b95e940f --- /dev/null +++ b/euler/shockbubble/3d/setplot3.m @@ -0,0 +1,59 @@ +% setplot3.m +% called in plotclaw3 before plotting to set various parameters + +OutputDir = '_output'; + +setopengl; +setviews % set viewpoints so that view(xSlice), for example, can be used. + +PlotType = 1; % type of plot to produce: + % 1 = pcolor on slices (with optional contours) + % 2 = contour lines in 3d on transparent slices + % 3 = Schlieren plot on slices + % 4 = scatter plot of q vs. r + +mq = 1; % which component of q to plot +UserVariable = 0; % set to 1 to specify a user-defined variable +UserVariableFile = ' '; % name of m-file mapping data to q +MappedGrid = 0; % set to 1 if mapc2p.m exists for nonuniform grid +MaxFrames = 1000; % max number of frames to loop over +MaxLevels = 6; % max number of AMR levels + +PlotData = [1 1 1 0 0 0]; % Data on refinement level k is plotted only + % if k'th component is nonzero +PlotGrid = [1 0 0 0 0 0]; % Plot grid lines on each level? +PlotGridEdges = [1 1 0 0 0 0]; % Plot edges of patches of each grid at + % this level on slices? +PlotCubeEdges = [0 0 0 0 0 0]; % Plot edges of cube of refinement patch at + % this level? + + +% ContourValues is a vector of contour lines that can be used with +% PlotType = 1,2,3. Empty ==> no contour lines drawn + ContourValues = []; + +% The next three parameters are vectors of x,y,z coordinates of 2d slices +% to be displayed for PlotType = 1,2,3. Empty ==> no slices in that direction. + xSliceCoords = []; + ySliceCoords = 0.5; + zSliceCoords = 0.; + +% For PlotType = 4 (Scatter plot) +% plot q(r) vs. r = sqrt((x-x0)^2 + (y-y0)^2 + (z-z0)^2); +% Use symbol PlotStyle(k) at refinement level k. + x0 = 0.5; + y0 = 0.5; + z0 = 0.5; + PlotStyle = setplotstyle('o','x','.','s','v','^'); + +% Note: Lengths of SurfTransparency and SurfColors must greater than or +% equal to the length of SurfValues.. + + IsosurfValues = [0.5]; % Plot surfaces at q = surfValue(i). + + IsosurfAlphas = [0.5]; % Transparency of each surface + % (0=clear; 1=opaque) + % NOTE: Your system must be able to + % use the OpenGL Renderer. + + IsosurfColors = strvcat('y'); % Colors for each surface. diff --git a/euler/shockbubble/3d/setprob.f b/euler/shockbubble/3d/setprob.f new file mode 100644 index 00000000..e7708e6e --- /dev/null +++ b/euler/shockbubble/3d/setprob.f @@ -0,0 +1,65 @@ +c ================== + subroutine setprob +c ================== + + implicit double precision (a-h,o-z) + character*12 fname + common /comic/ qin(5),qout(5) + common /cparam/ gamma + common /cdisc/ x0,y0,r0,radius + common /cominf/ rinf,vinf,einf +c +c + iunit = 7 + fname = 'setprob.data' +c # open the unit with new routine from Clawpack 4.4 to skip over +c # comment lines starting with #: + call opendatafile(iunit, fname) + +c # set idisc for cellave routines (see function fdisc) + idisc = 2 +c + read(7,*) gamma + gamma1 = gamma - 1.d0 + +c # read center and radius of bubble: + read(7,*) x0 + read(7,*) y0 + read(7,*) r0 + +c # density in bubble: + read(7,*) rhoin + +c # pressure behind shock: + read(7,*) pinf +c +c # density outside bubble and pressure ahead of shock are fixed: + rhoout = 1.d0 + pout = 1.d0 + pin = 1.d0 + + qin(1) = rhoin + qin(2) = 0.d0 + qin(3) = 0.d0 + qin(4) = 0.d0 + qin(5) = pin/gamma1 + + qout(1) = rhoout + qout(2) = 0.d0 + qout(3) = 0.d0 + qout(4) = 0.d0 + qout(5) = pout/gamma1 +c +c # Compute density and velocity behind shock from Hugoniot relations: +c # ------------------------------------------------------------------ + + rinf = ( gamma1 + (gamma+1)*pinf )/( (gamma+1) + gamma1*pinf ) + vinf = (1.0d0/sqrt(gamma)) * (pinf - 1.d0)/ + & sqrt( 0.5*((gamma+1)/gamma) * pinf + 0.5*gamma1/gamma ) + einf = 0.5*rinf*vinf*vinf + pinf/gamma1 + + write(6,601) pinf,rinf,vinf,einf + 601 format('pinf,rinf,vinf,einf:',/, 4e14.6) + + return + end diff --git a/euler/shockbubble/3d/setrun.py b/euler/shockbubble/3d/setrun.py new file mode 100644 index 00000000..03696145 --- /dev/null +++ b/euler/shockbubble/3d/setrun.py @@ -0,0 +1,345 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +import os +import numpy as np + +#------------------------------ +def setrun(claw_pkg='amrclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "amrclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + + assert claw_pkg.lower() == 'amrclaw', "Expected claw_pkg = 'amrclaw'" + + num_dim = 3 + rundata = data.ClawRunData(claw_pkg, num_dim) + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + probdata.add_param('gamma', 1.4, 'gamma') + probdata.add_param('x0', .5, 'x0') + probdata.add_param('y0', .5, 'y0') + probdata.add_param('r0', .2, 'r0') + probdata.add_param('rhoin', .1, 'rhoin') + probdata.add_param('pinf', 5.0, 'pinf') + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amrclaw.data for AMR) + #------------------------------------------------------------------ + + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + clawdata.lower[0] = 0.000000e+00 # xlower + clawdata.upper[0] = 2.000000e+00 # xupper + clawdata.lower[1] = 0.000000e+00 # ylower + clawdata.upper[1] = 0.500000e+00 # yupper + clawdata.lower[2] = 0.000000e+00 # zlower + clawdata.upper[2] = 0.500000e+00 # zupper + + # Number of grid cells: + clawdata.num_cells[0] = 40 # mx + clawdata.num_cells[1] = 10 # my + clawdata.num_cells[2] = 10 # mz + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 5 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 0 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 0 + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.000000 + + + # Restart from checkpoint file of a previous run? + # Note: If restarting, you must also change the Makefile to set: + # RESTART = True + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = 'fort.chk00006' # File to use for restart data + + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output ntimes frames at equally spaced times up to tfinal: + # Can specify num_output_times = 0 for no output + clawdata.num_output_times = 7 + clawdata.tfinal = .7 + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list or numpy array of output times: + # Include t0 if you want output at the initial time. + clawdata.output_times = [0., 1.0, 2.0, 3.0] + + elif clawdata.output_style == 3: + # Output every step_interval timesteps over total_steps timesteps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 10 + clawdata.output_t0 = True # output at initial (or restart) time? + + + clawdata.output_format = 'ascii' # 'ascii', 'binary', 'netcdf' + + clawdata.output_q_components = 'all' # could be list such as [True,True] + clawdata.output_aux_components = 'none' # could be list + clawdata.output_aux_onlyonce = True # output aux arrays only at t0 + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==True: variable time steps used based on cfl_desired, + # if dt_variable==False: fixed time steps dt = dt_initial always used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # (If dt_variable==0 then dt=dt_initial for all steps) + clawdata.dt_initial = 5e-3 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1.000000e+99 + + # Desired Courant number if variable dt used + clawdata.cfl_desired = 0.850000 + # max Courant number to allow without retaking step with a smaller dt: + clawdata.cfl_max = 1.000000 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 50000 + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? + clawdata.dimensional_split = 'godunov' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 22 + + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'vanleer' ==> van Leer + # 4 or 'mc' ==> MC limiter + clawdata.limiter = ['mc','mc','mc'] + + clawdata.use_fwaves = False # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 0 + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 or 'user' => user specified (must modify bcNamr.f to use this option) + # 1 or 'extrap' => extrapolation (non-reflecting outflow) + # 2 or 'periodic' => periodic (must specify this at both boundaries) + # 3 or 'wall' => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' # at xlower + clawdata.bc_upper[0] = 'extrap' # at xupper + + clawdata.bc_lower[1] = 'extrap' # at ylower + clawdata.bc_upper[1] = 'wall' # at yupper + + clawdata.bc_lower[2] = 'wall' # at zlower + clawdata.bc_upper[2] = 'extrap' # at zupper + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif clawdata.checkpt_style == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = [0.1,0.15] + + elif clawdata.checkpt_style == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + + # --------------- + # AMR parameters: + # --------------- + + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 3 + + # List of refinement ratios at each level (length at least amr_level_max-1) + amrdata.refinement_ratios_x = [4, 2] + amrdata.refinement_ratios_y = [4, 2] + amrdata.refinement_ratios_z = [4, 2] + amrdata.refinement_ratios_t = [4, 2] + + + # Specify type of each aux variable in clawdata.auxtype. + # This must be a list of length num_aux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + amrdata.aux_type = [] + + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = False # use Richardson? + amrdata.flag_richardson_tol = 0.001000e+00 # Richardson tolerance + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = True # use this? + amrdata.flag2refine_tol = 0.2 # tolerance used in this routine + # User can modify flag2refine to change the criterion for flagging. + # Default: check maximum absolute difference of first component of q + # between a cell and each of its neighbors. + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 2 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 3 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.7 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 1 + + + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + return rundata + + # end of function setrun + # ---------------------- + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() + diff --git a/euler_2d_shockbubble_amrclaw/README.txt b/euler_2d_shockbubble_amrclaw/README.txt deleted file mode 100644 index a7186eff..00000000 --- a/euler_2d_shockbubble_amrclaw/README.txt +++ /dev/null @@ -1,4 +0,0 @@ - -Shock-bubble interaction solved with 2-dimensional axi-symmetric Euler -equations and AMR. -