diff --git a/examples/CosmoVolume/run.sh b/examples/CosmoVolume/run.sh
index 412456c3cb1ae9b869afa52b1046747e32c2eefe..e128aa8f88cc30e13c420bda2088c9958c5bfc31 100755
--- a/examples/CosmoVolume/run.sh
+++ b/examples/CosmoVolume/run.sh
@@ -7,4 +7,4 @@ then
./getIC.sh
fi
-../swift -s -t 16 cosmoVolume.yml
+../swift -s -t 16 cosmoVolume.yml 2>&1 | tee output.log
diff --git a/examples/EAGLE_12/run.sh b/examples/EAGLE_12/run.sh
index 58c92adc0a1afadf0d15f81613d749f6e736f211..21b965050ab1c786fd08c0dfd28b1612ac09e9e5 100755
--- a/examples/EAGLE_12/run.sh
+++ b/examples/EAGLE_12/run.sh
@@ -7,5 +7,5 @@ then
./getIC.sh
fi
-../swift -s -t 16 eagle_12.yml
+../swift -s -t 16 eagle_12.yml 2>&1 | tee output.log
diff --git a/examples/EAGLE_25/run.sh b/examples/EAGLE_25/run.sh
index f232cf98f772ab9900bd22b635090e3c103749b5..21d0c7baf4482bdc6df83c552253a02d47432ba2 100755
--- a/examples/EAGLE_25/run.sh
+++ b/examples/EAGLE_25/run.sh
@@ -7,5 +7,5 @@ then
./getIC.sh
fi
-../swift -s -t 16 eagle_25.yml
+../swift -s -t 16 eagle_25.yml 2>&1 | tee output.log
diff --git a/examples/EAGLE_50/run.sh b/examples/EAGLE_50/run.sh
index 87d522be6bd95fa503d437af861a68ce8bf85b4c..253e7e6eb11a04e1731c605e9ea067ba2bd13221 100755
--- a/examples/EAGLE_50/run.sh
+++ b/examples/EAGLE_50/run.sh
@@ -7,5 +7,5 @@ then
./getIC.sh
fi
-../swift -s -t 16 eagle_50.yml
+../swift -s -t 16 eagle_50.yml 2>&1 | tee output.log
diff --git a/examples/ExternalPointMass/run.sh b/examples/ExternalPointMass/run.sh
index 9bbe03738b472f48715ae1875dfd462ef577b3d9..9f90ca395a5c8cf83e67928b3fdbd4d8529ac254 100755
--- a/examples/ExternalPointMass/run.sh
+++ b/examples/ExternalPointMass/run.sh
@@ -7,4 +7,4 @@ then
python makeIC.py 10000
fi
-../swift -g -t 2 externalPointMass.yml
+../swift -g -t 2 externalPointMass.yml 2>&1 | tee output.log
diff --git a/examples/Glass/glass_50000.hdf5 b/examples/Glass/glass_50000.hdf5
deleted file mode 100644
index ad209b851f04ee26e527fe925b71322a2740bb55..0000000000000000000000000000000000000000
Binary files a/examples/Glass/glass_50000.hdf5 and /dev/null differ
diff --git a/examples/GreshoVortex/getGlass.sh b/examples/GreshoVortex_2D/getGlass.sh
similarity index 100%
rename from examples/GreshoVortex/getGlass.sh
rename to examples/GreshoVortex_2D/getGlass.sh
diff --git a/examples/GreshoVortex/gresho.yml b/examples/GreshoVortex_2D/gresho.yml
similarity index 96%
rename from examples/GreshoVortex/gresho.yml
rename to examples/GreshoVortex_2D/gresho.yml
index 066b433a948db232c9665fd7d007753a19e058fa..b4de5bde517556cacb94d996f5a11cbe05188bf9 100644
--- a/examples/GreshoVortex/gresho.yml
+++ b/examples/GreshoVortex_2D/gresho.yml
@@ -27,7 +27,7 @@ Statistics:
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
- max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
+ max_smoothing_length: 0.02 # Maximal smoothing length allowed (in internal units).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
diff --git a/examples/GreshoVortex/makeIC.py b/examples/GreshoVortex_2D/makeIC.py
similarity index 94%
rename from examples/GreshoVortex/makeIC.py
rename to examples/GreshoVortex_2D/makeIC.py
index 2f5bebc00ce0f86d3f4f3cccd030cfff5f90d51d..96fa7f098d5eb26c42a984e2b5ec94bafc710dc3 100644
--- a/examples/GreshoVortex/makeIC.py
+++ b/examples/GreshoVortex_2D/makeIC.py
@@ -1,7 +1,6 @@
###############################################################################
# This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
@@ -19,9 +18,7 @@
##############################################################################
import h5py
-import random
from numpy import *
-import sys
# Generates a swift IC file for the Gresho-Chan vortex in a periodic box
@@ -51,7 +48,6 @@ for i in range(numPart):
x = coords[i,0]
y = coords[i,1]
- z = coords[i,2]
r2 = (x - boxSize / 2)**2 + (y - boxSize / 2)**2
r = sqrt(r2)
@@ -83,7 +79,7 @@ fileOutput = h5py.File(fileOutputName, 'w')
# Header
grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
+grp.attrs["BoxSize"] = [boxSize, boxSize, 0.2]
grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
diff --git a/examples/GreshoVortex/plotSolution.py b/examples/GreshoVortex_2D/plotSolution.py
similarity index 100%
rename from examples/GreshoVortex/plotSolution.py
rename to examples/GreshoVortex_2D/plotSolution.py
diff --git a/examples/GreshoVortex/run.sh b/examples/GreshoVortex_2D/run.sh
similarity index 89%
rename from examples/GreshoVortex/run.sh
rename to examples/GreshoVortex_2D/run.sh
index 81c016f4775a1e5d693e3c2f78c7d1e5dc532701..6d537bcc96c68385434f685bd551a2d423f469e0 100755
--- a/examples/GreshoVortex/run.sh
+++ b/examples/GreshoVortex_2D/run.sh
@@ -13,7 +13,7 @@ then
fi
# Run SWIFT
-../swift -s -t 1 gresho.yml
+../swift -s -t 1 gresho.yml 2>&1 | tee output.log
# Plot the solution
python plotSolution.py 11
diff --git a/examples/IsothermalPotential/GravityOnly/legend.pro b/examples/IsothermalPotential/GravityOnly/legend.pro
deleted file mode 100755
index e510bdfcd7dcabed796cb090c6fa7d54536608b6..0000000000000000000000000000000000000000
--- a/examples/IsothermalPotential/GravityOnly/legend.pro
+++ /dev/null
@@ -1,459 +0,0 @@
-;+
-; NAME:
-; LEGEND
-; PURPOSE:
-; Create an annotation legend for a plot.
-; EXPLANATION:
-; This procedure makes a legend for a plot. The legend can contain
-; a mixture of symbols, linestyles, Hershey characters (vectorfont),
-; and filled polygons (usersym). A test procedure, legendtest.pro,
-; shows legend's capabilities. Placement of the legend is controlled
-; with keywords like /right, /top, and /center or by using a position
-; keyword for exact placement (position=[x,y]) or via mouse (/position).
-; CALLING SEQUENCE:
-; LEGEND [,items][,keyword options]
-; EXAMPLES:
-; The call:
-; legend,['Plus sign','Asterisk','Period'],psym=[1,2,3]
-; produces:
-; -----------------
-; | |
-; | + Plus sign |
-; | * Asterisk |
-; | . Period |
-; | |
-; -----------------
-; Each symbol is drawn with a plots command, so they look OK.
-; Other examples are given in optional output keywords.
-;
-; lines = indgen(6) ; for line styles
-; items = 'linestyle '+strtrim(lines,2) ; annotations
-; legend,items,linestyle=lines ; vertical legend---upper left
-; items = ['Plus sign','Asterisk','Period']
-; sym = [1,2,3]
-; legend,items,psym=sym ; ditto except using symbols
-; legend,items,psym=sym,/horizontal ; horizontal format
-; legend,items,psym=sym,box=0 ; sans border
-; legend,items,psym=sym,delimiter='=' ; embed '=' betw psym & text
-; legend,items,psym=sym,margin=2 ; 2-character margin
-; legend,items,psym=sym,position=[x,y] ; upper left in data coords
-; legend,items,psym=sym,pos=[x,y],/norm ; upper left in normal coords
-; legend,items,psym=sym,pos=[x,y],/device ; upper left in device coords
-; legend,items,psym=sym,/position ; interactive position
-; legend,items,psym=sym,/right ; at upper right
-; legend,items,psym=sym,/bottom ; at lower left
-; legend,items,psym=sym,/center ; approximately near center
-; legend,items,psym=sym,number=2 ; plot two symbols, not one
-; legend,items,/fill,psym=[8,8,8],colors=[10,20,30]; 3 filled squares
-; INPUTS:
-; items = text for the items in the legend, a string array.
-; For example, items = ['diamond','asterisk','square'].
-; You can omit items if you don't want any text labels.
-; OPTIONAL INPUT KEYWORDS:
-;
-; linestyle = array of linestyle numbers If linestyle[i] < 0, then omit
-; ith symbol or line to allow a multi-line entry. If
-; linestyle = -99 then text will be left-justified.
-; psym = array of plot symbol numbers. If psym[i] is negative, then a
-; line connects pts for ith item. If psym[i] = 8, then the
-; procedure usersym is called with vertices define in the
-; keyword usersym. If psym[i] = 88, then use the previously
-; defined user symbol
-; vectorfont = vector-drawn characters for the sym/line column, e.g.,
-; ['!9B!3','!9C!3','!9D!3'] produces an open square, a checkmark,
-; and a partial derivative, which might have accompanying items
-; ['BOX','CHECK','PARTIAL DERIVATIVE'].
-; There is no check that !p.font is set properly, e.g., -1 for
-; X and 0 for PostScript. This can produce an error, e.g., use
-; !20 with PostScript and !p.font=0, but allows use of Hershey
-; *AND* PostScript fonts together.
-; N. B.: Choose any of linestyle, psym, and/or vectorfont. If none is
-; present, only the text is output. If more than one
-; is present, all need the same number of elements, and normal
-; plot behaviour occurs.
-; By default, if psym is positive, you get one point so there is
-; no connecting line. If vectorfont[i] = '',
-; then plots is called to make a symbol or a line, but if
-; vectorfont[i] is a non-null string, then xyouts is called.
-; /help = flag to print header
-; /horizontal = flag to make the legend horizontal
-; /vertical = flag to make the legend vertical (D=vertical)
-; box = flag to include/omit box around the legend (D=include)
-; clear = flag to clear the box area before drawing the legend
-; delimiter = embedded character(s) between symbol and text (D=none)
-; colors = array of colors for plot symbols/lines (D=!P.color)
-; textcolors = array of colors for text (D=!P.color)
-; margin = margin around text measured in characters and lines
-; spacing = line spacing (D=bit more than character height)
-; pspacing = psym spacing (D=3 characters) (when number of symbols is
-; greater than 1)
-; charsize = just like !p.charsize for plot labels
-; charthick = just like !p.charthick for plot labels
-; thick = array of line thickness numbers (D = !P.thick), if used, then
-; linestyle must also be specified
-; position = data coordinates of the /top (D) /left (D) of the legend
-; normal = use normal coordinates for position, not data
-; device = use device coordinates for position, not data
-; number = number of plot symbols to plot or length of line (D=1)
-; usersym = 2-D array of vertices, cf. usersym in IDL manual.
-; (/USERSYM =square, default is to use existing USERSYM definition)
-; /fill = flag to fill the usersym
-; /left_legend = flag to place legend snug against left side of plot
-; window (D)
-; /right_legend = flag to place legend snug against right side of plot
-; window. If /right,pos=[x,y], then x is position of RHS and
-; text runs right-to-left.
-; /top_legend = flag to place legend snug against top of plot window (D)
-; /bottom = flag to place legend snug against bottom of plot window
-; /top,pos=[x,y] and /bottom,pos=[x,y] produce same positions.
-;
-; If LINESTYLE, PSYM, VECTORFONT, THICK, COLORS, or TEXTCOLORS are
-; supplied as scalars, then the scalar value is set for every line or
-; symbol in the legend.
-; Outputs:
-; legend to current plot device
-; OPTIONAL OUTPUT KEYWORDS:
-; corners = 4-element array, like !p.position, of the normalized
-; coords for the box (even if box=0): [llx,lly,urx,ury].
-; Useful for multi-column or multi-line legends, for example,
-; to make a 2-column legend, you might do the following:
-; c1_items = ['diamond','asterisk','square']
-; c1_psym = [4,2,6]
-; c2_items = ['solid','dashed','dotted']
-; c2_line = [0,2,1]
-; legend,c1_items,psym=c1_psym,corners=c1,box=0
-; legend,c2_items,line=c2_line,corners=c2,box=0,pos=[c1[2],c1[3]]
-; c = [c1[0]<c2[0],c1[1]<c2[1],c1[2]>c2[2],c1[3]>c2[3]]
-; plots,[c[0],c[0],c[2],c[2],c[0]],[c[1],c[3],c[3],c[1],c[1]],/norm
-; Useful also to place the legend. Here's an automatic way to place
-; the legend in the lower right corner. The difficulty is that the
-; legend's width is unknown until it is plotted. In this example,
-; the legend is plotted twice: the first time in the upper left, the
-; second time in the lower right.
-; legend,['1','22','333','4444'],linestyle=indgen(4),corners=corners
-; ; BOGUS LEGEND---FIRST TIME TO REPORT CORNERS
-; xydims = [corners[2]-corners[0],corners[3]-corners[1]]
-; ; SAVE WIDTH AND HEIGHT
-; chdim=[!d.x_ch_size/float(!d.x_size),!d.y_ch_size/float(!d.y_size)]
-; ; DIMENSIONS OF ONE CHARACTER IN NORMALIZED COORDS
-; pos = [!x.window[1]-chdim[0]-xydims[0] $
-; ,!y.window[0]+chdim[1]+xydims[1]]
-; ; CALCULATE POSITION FOR LOWER RIGHT
-; plot,findgen(10) ; SIMPLE PLOT; YOU DO WHATEVER YOU WANT HERE.
-; legend,['1','22','333','4444'],linestyle=indgen(4),pos=pos
-; ; REDO THE LEGEND IN LOWER RIGHT CORNER
-; You can modify the pos calculation to place the legend where you
-; want. For example to place it in the upper right:
-; pos = [!x.window[1]-chdim[0]-xydims[0],!y.window[1]-xydims[1]]
-; Common blocks:
-; none
-; Procedure:
-; If keyword help is set, call doc_library to print header.
-; See notes in the code. Much of the code deals with placement of the
-; legend. The main problem with placement is not being
-; able to sense the length of a string before it is output. Some crude
-; approximations are used for centering.
-; Restrictions:
-; Here are some things that aren't implemented.
-; - An orientation keyword would allow lines at angles in the legend.
-; - An array of usersyms would be nice---simple change.
-; - An order option to interchange symbols and text might be nice.
-; - Somebody might like double boxes, e.g., with box = 2.
-; - Another feature might be a continuous bar with ticks and text.
-; - There are no guards to avoid writing outside the plot area.
-; - There is no provision for multi-line text, e.g., '1st line!c2nd line'
-; Sensing !c would be easy, but !c isn't implemented for PostScript.
-; A better way might be to simply output the 2nd line as another item
-; but without any accompanying symbol or linestyle. A flag to omit
-; the symbol and linestyle is linestyle[i] = -1.
-; - There is no ability to make a title line containing any of titles
-; for the legend, for the symbols, or for the text.
-; Side Effects:
-; Modification history:
-; write, 24-25 Aug 92, F K Knight (knight@ll.mit.edu)
-; allow omission of items or omission of both psym and linestyle, add
-; corners keyword to facilitate multi-column legends, improve place-
-; ment of symbols and text, add guards for unequal size, 26 Aug 92, FKK
-; add linestyle(i)=-1 to suppress a single symbol/line, 27 Aug 92, FKK
-; add keyword vectorfont to allow characters in the sym/line column,
-; 28 Aug 92, FKK
-; add /top, /bottom, /left, /right keywords for automatic placement at
-; the four corners of the plot window. The /right keyword forces
-; right-to-left printing of menu. 18 Jun 93, FKK
-; change default position to data coords and add normal, data, and
-; device keywords, 17 Jan 94, FKK
-; add /center keyword for positioning, but it is not precise because
-; text string lengths cannot be known in advance, 17 Jan 94, FKK
-; add interactive positioning with /position keyword, 17 Jan 94, FKK
-; allow a legend with just text, no plotting symbols. This helps in
-; simply describing a plot or writing assumptions done, 4 Feb 94, FKK
-; added thick, symsize, and clear keyword Feb 96, W. Landsman HSTX
-; David Seed, HR Wallingford, d.seed@hrwallingford.co.uk
-; allow scalar specification of keywords, Mar 96, W. Landsman HSTX
-; added charthick keyword, June 96, W. Landsman HSTX
-; Made keyword names left,right,top,bottom,center longer,
-; Aug 16, 2000, Kim Tolbert
-; Added ability to have regular text lines in addition to plot legend
-; lines in legend. If linestyle is -99 that item is left-justified.
-; Previously, only option for no sym/line was linestyle=-1, but then text
-; was lined up after sym/line column. 10 Oct 2000, Kim Tolbert
-; Make default value of thick = !P.thick W. Landsman Jan. 2001
-; Don't overwrite existing USERSYM definition W. Landsman Mar. 2002
-;-
-pro legend, items, BOTTOM_LEGEND=bottom, BOX = box, CENTER_LEGEND=center, $
- CHARTHICK=charthick, CHARSIZE = charsize, CLEAR = clear, COLORS = colorsi, $
- CORNERS = corners, DATA=data, DELIMITER=delimiter, DEVICE=device, $
- FILL=fill, HELP = help, HORIZONTAL=horizontal,LEFT_LEGEND=left, $
- LINESTYLE=linestylei, MARGIN=margin, NORMAL=normal, NUMBER=number, $
- POSITION=position,PSPACING=pspacing, PSYM=psymi, RIGHT_LEGEND=right, $
- SPACING=spacing, SYMSIZE=symsize, TEXTCOLORS=textcolorsi, THICK=thicki, $
- TOP_LEGEND=top, USERSYM=usersym, VECTORFONT=vectorfonti, VERTICAL=vertical
-;
-; =====>> HELP
-;
-on_error,2
-if keyword_set(help) then begin & doc_library,'legend' & return & endif
-;
-; =====>> SET DEFAULTS FOR SYMBOLS, LINESTYLES, AND ITEMS.
-;
- ni = n_elements(items)
- np = n_elements(psymi)
- nl = n_elements(linestylei)
- nth = n_elements(thicki)
- nv = n_elements(vectorfonti)
- nlpv = max([np,nl,nv])
- n = max([ni,np,nl,nv]) ; NUMBER OF ENTRIES
-strn = strtrim(n,2) ; FOR ERROR MESSAGES
-if n eq 0 then message,'No inputs! For help, type legend,/help.'
-if ni eq 0 then begin
- items = replicate('',n) ; DEFAULT BLANK ARRAY
-endif else begin
- if size(items,/TNAME) NE 'STRING' then message, $
- 'First parameter must be a string array. For help, type legend,/help.'
- if ni ne n then message,'Must have number of items equal to '+strn
-endelse
-symline = (np ne 0) or (nl ne 0) ; FLAG TO PLOT SYM/LINE
- if (np ne 0) and (np ne n) and (np NE 1) then message, $
- 'Must have 0, 1 or '+strn+' elements in PSYM array.'
- if (nl ne 0) and (nl ne n) and (nl NE 1) then message, $
- 'Must have 0, 1 or '+strn+' elements in LINESTYLE array.'
- if (nth ne 0) and (nth ne n) and (nth NE 1) then message, $
- 'Must have 0, 1 or '+strn+' elements in THICK array.'
-
- case nl of
- 0: linestyle = intarr(n) ;Default = solid
- 1: linestyle = intarr(n) + linestylei
- else: linestyle = linestylei
- endcase
-
- case nth of
- 0: thick = replicate(!p.thick,n) ;Default = !P.THICK
- 1: thick = intarr(n) + thicki
- else: thick = thicki
- endcase
-
- case np of ;Get symbols
- 0: psym = intarr(n) ;Default = solid
- 1: psym = intarr(n) + psymi
- else: psym = psymi
- endcase
-
- case nv of
- 0: vectorfont = replicate('',n)
- 1: vectorfont = replicate(vectorfonti,n)
- else: vectorfont = vectorfonti
- endcase
-;
-; =====>> CHOOSE VERTICAL OR HORIZONTAL ORIENTATION.
-;
-if n_elements(horizontal) eq 0 then begin ; D=VERTICAL
- if n_elements(vertical) eq 0 then vertical = 1
-endif else begin
- if n_elements(vertical) eq 0 then vertical = not horizontal
-endelse
-;
-; =====>> SET DEFAULTS FOR OTHER OPTIONS.
-;
-if n_elements(box) eq 0 then box = 1
-if n_elements(clear) eq 0 then clear = 0
-
-if n_elements(margin) eq 0 then margin = 0.5
-if n_elements(delimiter) eq 0 then delimiter = ''
-if n_elements(charsize) eq 0 then charsize = !p.charsize
-if n_elements(charthick) eq 0 then charthick = !p.charthick
-if charsize eq 0 then charsize = 1
-if (n_elements (symsize) eq 0) then symsize= charsize + intarr(n)
-if n_elements(number) eq 0 then number = 1
- case N_elements(colorsi) of
- 0: colors = replicate(!P.color,n) ;Default is !P.COLOR
- 1: colors = replicate(colorsi,n)
- else: colors = colorsi
- endcase
-
- case N_elements(textcolorsi) of
- 0: textcolors = replicate(!P.color,n) ;Default is !P.COLOR
- 1: textcolors = replicate(textcolorsi,n)
- else: textcolors = textcolorsi
- endcase
- fill = keyword_set(fill)
-if n_elements(usersym) eq 1 then usersym = 2*[[0,0],[0,1],[1,1],[1,0],[0,0]]-1
-;
-; =====>> INITIALIZE SPACING
-;
-if n_elements(spacing) eq 0 then spacing = 1.2
-if n_elements(pspacing) eq 0 then pspacing = 3
-xspacing = !d.x_ch_size/float(!d.x_size) * (spacing > charsize)
-yspacing = !d.y_ch_size/float(!d.y_size) * (spacing > charsize)
-ltor = 1 ; flag for left-to-right
-if n_elements(left) eq 1 then ltor = left eq 1
-if n_elements(right) eq 1 then ltor = right ne 1
-ttob = 1 ; flag for top-to-bottom
-if n_elements(top) eq 1 then ttob = top eq 1
-if n_elements(bottom) eq 1 then ttob = bottom ne 1
-xalign = ltor ne 1 ; x alignment: 1 or 0
-yalign = -0.5*ttob + 1 ; y alignment: 0.5 or 1
-xsign = 2*ltor - 1 ; xspacing direction: 1 or -1
-ysign = 2*ttob - 1 ; yspacing direction: 1 or -1
-if not ttob then yspacing = -yspacing
-if not ltor then xspacing = -xspacing
-;
-; =====>> INITIALIZE POSITIONS: FIRST CALCULATE X OFFSET FOR TEXT
-;
-xt = 0
-if nlpv gt 0 then begin ; SKIP IF TEXT ITEMS ONLY.
-if vertical then begin ; CALC OFFSET FOR TEXT START
- for i = 0,n-1 do begin
- if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
- if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
- if psym[i] eq 0 then expand = 1 else expand = 2
- thisxt = (expand*pspacing*(num-1)*xspacing)
- if ltor then xt = thisxt > xt else xt = thisxt < xt
- endfor
-endif ; NOW xt IS AN X OFFSET TO ALIGN ALL TEXT ENTRIES.
-endif
-;
-; =====>> INITIALIZE POSITIONS: SECOND LOCATE BORDER
-;
-if !x.window[0] eq !x.window[1] then begin
- plot,/nodata,xstyle=4,ystyle=4,[0],/noerase
-endif
-; next line takes care of weirdness with small windows
-pos = [min(!x.window),min(!y.window),max(!x.window),max(!y.window)]
-case n_elements(position) of
- 0: begin
- if ltor then px = pos[0] else px = pos[2]
- if ttob then py = pos[3] else py = pos[1]
- if keyword_set(center) then begin
- if not keyword_set(right) and not keyword_set(left) then $
- px = (pos[0] + pos[2])/2. - xt
- if not keyword_set(top) and not keyword_set(bottom) then $
- py = (pos[1] + pos[3])/2. + n*yspacing
- endif
- position = [px,py] + [xspacing,-yspacing]
- end
- 1: begin ; interactive
- message,/inform,'Place mouse at upper left corner and click any mouse button.'
- cursor,x,y,/normal
- position = [x,y]
- end
- 2: begin ; convert upper left corner to normal coordinates
- if keyword_set(data) then $
- position = convert_coord(position,/to_norm) $
- else if keyword_set(device) then $
- position = convert_coord(position,/to_norm,/device) $
- else if not keyword_set(normal) then $
- position = convert_coord(position,/to_norm)
- end
- else: message,'Position keyword can have 0, 1, or 2 elements only. Try legend,/help.'
-endcase
-
-yoff = 0.25*yspacing*ysign ; VERT. OFFSET FOR SYM/LINE.
-
-x0 = position[0] + (margin)*xspacing ; INITIAL X & Y POSITIONS
-y0 = position[1] - margin*yspacing + yalign*yspacing ; WELL, THIS WORKS!
-;
-; =====>> OUTPUT TEXT FOR LEGEND, ITEM BY ITEM.
-; =====>> FOR EACH ITEM, PLACE SYM/LINE, THEN DELIMITER,
-; =====>> THEN TEXT---UPDATING X & Y POSITIONS EACH TIME.
-; =====>> THERE ARE A NUMBER OF EXCEPTIONS DONE WITH IF STATEMENTS.
-;
-for iclr = 0,clear do begin
- y = y0 ; STARTING X & Y POSITIONS
- x = x0
- if ltor then xend = 0 else xend = 1 ; SAVED WIDTH FOR DRAWING BOX
-
- if ttob then ii = [0,n-1,1] else ii = [n-1,0,-1]
- for i = ii[0],ii[1],ii[2] do begin
- if vertical then x = x0 else y = y0 ; RESET EITHER X OR Y
- x = x + xspacing ; UPDATE X & Y POSITIONS
- y = y - yspacing
- if nlpv eq 0 then goto,TEXT_ONLY ; FLAG FOR TEXT ONLY
- if (psym[i] eq 0) and (vectorfont[i] eq '') then num = (number + 1) > 3 else num = number
- if psym[i] lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE
- if psym[i] eq 0 then expand = 1 else expand = 2
- xp = x + expand*pspacing*indgen(num)*xspacing
- if (psym[i] gt 0) and (num eq 1) and vertical then xp = x + xt/2.
- yp = y + intarr(num)
- if vectorfont[i] eq '' then yp = yp + yoff
- if psym[i] eq 0 then begin
- xp = [min(xp),max(xp)] ; TO EXPOSE LINESTYLES
- yp = [min(yp),max(yp)] ; DITTO
- endif
- if (psym[i] eq 8) and (N_elements(usersym) GT 1) then $
- usersym,usersym,fill=fill,color=colors[i]
-;; extra by djseed .. psym=88 means use the already defined usersymbol
- if psym[i] eq 88 then psym[i] =8
- if vectorfont[i] ne '' then begin
-; if (num eq 1) and vertical then xp = x + xt/2 ; IF 1, CENTERED.
- xyouts,xp,yp,vectorfont[i],width=width,color=colors[i] $
- ,size=charsize,align=xalign,charthick = charthick,/norm
- xt = xt > width
- xp = xp + width/2.
- endif else begin
- if symline and (linestyle[i] ge 0) then plots,xp,yp,color=colors[i] $
- ,/normal,linestyle=linestyle[i],psym=psym[i],symsize=symsize[i], $
- thick=thick[i]
- endelse
-
- if vertical then x = x + xt else if ltor then x = max(xp) else x = min(xp)
- if symline then x = x + xspacing
- TEXT_ONLY:
- if vertical and (vectorfont[i] eq '') and symline and (linestyle[i] eq -99) then x=x0 + xspacing
- xyouts,x,y,delimiter,width=width,/norm,color=textcolors[i], $
- size=charsize,align=xalign,charthick = charthick
- x = x + width*xsign
- if width ne 0 then x = x + 0.5*xspacing
- xyouts,x,y,items[i],width=width,/norm,color=textcolors[i],size=charsize, $
- align=xalign,charthick=charthick
- x = x + width*xsign
- if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE
- xfinal = (x + xspacing*margin)
- if ltor then xend = xfinal > xend else xend = xfinal < xend ; UPDATE END X
- endfor
-
- if (iclr lt clear ) then begin
-; =====>> CLEAR AREA
- x = position[0]
- y = position[1]
- if vertical then bottom = n else bottom = 1
- ywidth = - (2*margin+bottom-0.5)*yspacing
- corners = [x,y+ywidth,xend,y]
- polyfill,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm,color=-1
-; plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],thick=2
- endif else begin
-
-;
-; =====>> OUTPUT BORDER
-;
- x = position[0]
- y = position[1]
- if vertical then bottom = n else bottom = 1
- ywidth = - (2*margin+bottom-0.5)*yspacing
- corners = [x,y+ywidth,xend,y]
- if box then plots,[x,xend,xend,x,x],y + [0,0,ywidth,ywidth,0],/norm
- return
- endelse
-endfor
-
-end
-
diff --git a/examples/IsothermalPotential/GravityOnly/run.sh b/examples/IsothermalPotential/GravityOnly/run.sh
index 329d04d30a6ccba7036f3802ff184c854ed00761..f6adfcececf4923485c0deabd97e9af9a6f64b05 100755
--- a/examples/IsothermalPotential/GravityOnly/run.sh
+++ b/examples/IsothermalPotential/GravityOnly/run.sh
@@ -7,4 +7,4 @@ then
python makeIC.py 1000 1
fi
-../../swift -g -t 2 isothermal.yml
+../../swift -g -t 2 isothermal.yml 2>&1 | tee output.log
diff --git a/examples/KelvinHelmholtz_2D/kelvinHelmholtz.yml b/examples/KelvinHelmholtz_2D/kelvinHelmholtz.yml
new file mode 100644
index 0000000000000000000000000000000000000000..38dd16880a209b885f7ad9c30c024988f4d8228f
--- /dev/null
+++ b/examples/KelvinHelmholtz_2D/kelvinHelmholtz.yml
@@ -0,0 +1,35 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 1.5 # The end time of the simulation (in internal units).
+ dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: kelvinHelmholtz # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.25 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ max_smoothing_length: 0.01 # Maximal smoothing length allowed (in internal units).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./kelvinHelmholtz.hdf5 # The file to read
diff --git a/examples/KelvinHelmholtz_2D/makeIC.py b/examples/KelvinHelmholtz_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..5c8632dea52ef301c453cfbf21c35923f12e2d5a
--- /dev/null
+++ b/examples/KelvinHelmholtz_2D/makeIC.py
@@ -0,0 +1,153 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+from numpy import *
+import sys
+
+# Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
+
+# Parameters
+L2 = 128 # Particles along one edge in the low-density region
+gamma = 5./3. # Gas adiabatic index
+P1 = 2.5 # Central region pressure
+P2 = 2.5 # Outskirts pressure
+v1 = 0.5 # Central region velocity
+v2 = -0.5 # Outskirts vlocity
+rho1 = 2 # Central density
+rho2 = 1 # Outskirts density
+omega0 = 0.1
+sigma = 0.05 / sqrt(2)
+fileOutputName = "kelvinHelmholtz.hdf5"
+#---------------------------------------------------
+
+# Start by generating grids of particles at the two densities
+numPart2 = L2 * L2
+L1 = int(sqrt(numPart2 / rho2 * rho1))
+numPart1 = L1 * L1
+
+#print "N2 =", numPart2, "N1 =", numPart1
+#print "L2 =", L2, "L1 = ", L1
+#print "rho2 =", rho2, "rho1 =", (float(L1*L1)) / (float(L2*L2))
+
+coords1 = zeros((numPart1, 3))
+coords2 = zeros((numPart2, 3))
+h1 = ones(numPart1) * 1.2348 / L1
+h2 = ones(numPart2) * 1.2348 / L2
+m1 = zeros(numPart1)
+m2 = zeros(numPart2)
+u1 = zeros(numPart1)
+u2 = zeros(numPart2)
+vel1 = zeros((numPart1, 3))
+vel2 = zeros((numPart2, 3))
+
+# Particles in the central region
+for i in range(L1):
+ for j in range(L1):
+
+ index = i * L1 + j
+
+ x = i / float(L1) + 1. / (2. * L1)
+ y = j / float(L1) + 1. / (2. * L1)
+
+ coords1[index, 0] = x
+ coords1[index, 1] = y
+ u1[index] = P1 / (rho1 * (gamma-1.))
+ vel1[index, 0] = v1
+
+# Particles in the outskirts
+for i in range(L2):
+ for j in range(L2):
+
+ index = i * L2 + j
+
+ x = i / float(L2) + 1. / (2. * L2)
+ y = j / float(L2) + 1. / (2. * L2)
+
+ coords2[index, 0] = x
+ coords2[index, 1] = y
+ u2[index] = P2 / (rho2 * (gamma-1.))
+ vel2[index, 0] = v2
+
+
+# Now concatenate arrays
+where1 = abs(coords1[:,1]-0.5) < 0.25
+where2 = abs(coords2[:,1]-0.5) > 0.25
+
+coords = append(coords1[where1, :], coords2[where2, :], axis=0)
+
+#print L2*(L2/2), L1*(L1/2)
+#print shape(coords), shape(coords1[where1,:]), shape(coords2[where2,:])
+#print shape(coords), shape(logical_not(coords1[where1,:])), shape(logical_not(coords2[where2,:]))
+
+vel = append(vel1[where1, :], vel2[where2, :], axis=0)
+h = append(h1[where1], h2[where2], axis=0)
+m = append(m1[where1], m2[where2], axis=0)
+u = append(u1[where1], u2[where2], axis=0)
+numPart = size(h)
+ids = linspace(1, numPart, numPart)
+m[:] = (0.5 * rho1 + 0.5 * rho2) / float(numPart)
+
+# Velocity perturbation
+vel[:,1] = omega0 * sin(4*pi*coords[:,0]) * (exp(-(coords[:,1]-0.25)**2 / (2 * sigma**2)) + exp(-(coords[:,1]-0.75)**2 / (2 * sigma**2)))
+
+#File
+fileOutput = h5py.File(fileOutputName, 'w')
+
+# Header
+grp = fileOutput.create_group("/Header")
+grp.attrs["BoxSize"] = [1., 1., 0.1]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFileOutputsPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+
+#Runtime parameters
+grp = fileOutput.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#Units
+grp = fileOutput.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+grp = fileOutput.create_group("/PartType0")
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
+ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = vel
+ds = grp.create_dataset('Masses', (numPart, 1), 'f')
+ds[()] = m.reshape((numPart,1))
+ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+ds[()] = h.reshape((numPart,1))
+ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
+ds[()] = u.reshape((numPart,1))
+ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
+ds[()] = ids.reshape((numPart,1))
+
+fileOutput.close()
+
+
diff --git a/examples/KelvinHelmholtz_2D/plotSolution.py b/examples/KelvinHelmholtz_2D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..9191f3ac7ec75c61d5fdab5d347c86222f787fab
--- /dev/null
+++ b/examples/KelvinHelmholtz_2D/plotSolution.py
@@ -0,0 +1,159 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the Gresho-Chan vortex and plots the SPH answer
+
+# Parameters
+gas_gamma = 5./3. # Gas adiabatic index
+P1 = 2.5 # Central region pressure
+P2 = 2.5 # Outskirts pressure
+v1 = 0.5 # Central region velocity
+v2 = -0.5 # Outskirts vlocity
+rho1 = 2 # Central density
+rho2 = 1 # Outskirts density
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+# Read the simulation data
+sim = h5py.File("kelvinHelmholtz_%03d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+pos = sim["/PartType0/Coordinates"][:,:]
+x = pos[:,0] - boxSize / 2
+y = pos[:,1] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:,:]
+v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
+rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+
+# Plot the interesting quantities
+figure()
+
+
+# Azimuthal velocity profile -----------------------------
+subplot(231)
+scatter(pos[:,0], pos[:,1], c=vel[:,0], cmap="PuBu", edgecolors='face', s=4, vmin=-1., vmax=1.)
+text(0.97, 0.97, "${\\rm{Velocity~along}}~x$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 1)
+
+# Radial density profile --------------------------------
+subplot(232)
+scatter(pos[:,0], pos[:,1], c=rho, cmap="PuBu", edgecolors='face', s=4, vmin=0.8, vmax=2.2)
+text(0.97, 0.97, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 1)
+
+# Radial pressure profile --------------------------------
+subplot(233)
+scatter(pos[:,0], pos[:,1], c=P, cmap="PuBu", edgecolors='face', s=4, vmin=1, vmax=4)
+text(0.97, 0.97, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 1)
+
+# Internal energy profile --------------------------------
+subplot(234)
+scatter(pos[:,0], pos[:,1], c=u, cmap="PuBu", edgecolors='face', s=4, vmin=1.5, vmax=5.)
+text(0.97, 0.97, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 1)
+
+# Radial entropy profile --------------------------------
+subplot(235)
+scatter(pos[:,0], pos[:,1], c=S, cmap="PuBu", edgecolors='face', s=4, vmin=0.5, vmax=3.)
+text(0.97, 0.97, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 1)
+
+# Image --------------------------------------------------
+#subplot(234)
+#scatter(pos[:,0], pos[:,1], c=v_norm, cmap="PuBu", edgecolors='face', s=4, vmin=0, vmax=1)
+#text(0.95, 0.95, "$|v|$", ha="right", va="top")
+#xlim(0,1)
+#ylim(0,1)
+#xlabel("$x$", labelpad=0)
+#ylabel("$y$", labelpad=0)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Kelvin-Helmholtz instability at $t=%.2f$"%(time), fontsize=10)
+text(-0.49, 0.8, "Centre:~~~ $(P, \\rho, v) = (%.3f, %.3f, %.3f)$"%(P1, rho1, v1), fontsize=10)
+text(-0.49, 0.7, "Outskirts: $(P, \\rho, v) = (%.3f, %.3f, %.3f)$"%(P2, rho2, v2), fontsize=10)
+plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.4, scheme, fontsize=10)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+savefig("KelvinHelmholtz.png", dpi=200)
diff --git a/examples/KelvinHelmholtz_2D/run.sh b/examples/KelvinHelmholtz_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..3d83bcc6dd7c226885ed83c3188f3177c4807154
--- /dev/null
+++ b/examples/KelvinHelmholtz_2D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e kelvinHelmholtz.hdf5 ]
+then
+ echo "Generating initial conditions for the Kelvin-Helmholtz example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 kelvinHelmholtz.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 6
diff --git a/examples/MultiTypes/run.sh b/examples/MultiTypes/run.sh
index 9782a27fe18a7fac94cc2a12353f81abb9efb06b..57465ce0ba6dde3988359df990f2a93323dbc617 100755
--- a/examples/MultiTypes/run.sh
+++ b/examples/MultiTypes/run.sh
@@ -7,4 +7,4 @@ then
python makeIC.py 50 60
fi
-../swift -s -g -t 16 multiTypes.yml
+../swift -s -g -t 16 multiTypes.yml 2>&1 | tee output.log
diff --git a/examples/PerturbedBox_3D/run.sh b/examples/PerturbedBox_3D/run.sh
index c15011ca526efb7250e93dc8cb002e49647b82b6..e20bff52d18322ce377fb589900fd9e13eefe64d 100755
--- a/examples/PerturbedBox_3D/run.sh
+++ b/examples/PerturbedBox_3D/run.sh
@@ -7,4 +7,4 @@ then
python makeIC.py 50
fi
-../swift -s -t 16 perturbedBox.yml
+../swift -s -t 16 perturbedBox.yml 2>&1 | tee output.log
diff --git a/examples/SedovBlast_3D/makeIC_fcc.py b/examples/SedovBlast_1D/makeIC.py
similarity index 51%
rename from examples/SedovBlast_3D/makeIC_fcc.py
rename to examples/SedovBlast_1D/makeIC.py
index 0d3a017a9b7f3b30b61e723e3d1646d7797b40a4..4bdf69eee99d98956d5e657be3f963d0cf9ea15b 100644
--- a/examples/SedovBlast_3D/makeIC_fcc.py
+++ b/examples/SedovBlast_1D/makeIC.py
@@ -1,7 +1,6 @@
###############################################################################
# This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
@@ -19,68 +18,42 @@
##############################################################################
import h5py
-import random
from numpy import *
# Generates a swift IC file for the Sedov blast test in a periodic cubic box
# Parameters
-periodic= 1 # 1 For periodic box
-boxSize = 5.
-L = 64 # Number of particles boxes along one axis
-rho = 1. # Density
-P = 1.e-5 # Pressure
-E0= 1.e2 # Energy of the explosion
-pert = 0.025
-gamma = 5./3. # Gas adiabatic index
-eta = 1.2349 # 48 ngbs with cubic spline kernel
+numPart = 1000
+gamma = 5./3. # Gas adiabatic index
+rho0 = 1. # Background density
+P0 = 1.e-6 # Background pressure
+E0= 1. # Energy of the explosion
+N_inject = 3 # Number of particles in which to inject energy
fileName = "sedov.hdf5"
-
#---------------------------------------------------
-numPart = 4*(L**3)
-mass = boxSize**3 * rho / numPart
-internalEnergy = P / ((gamma - 1.)*rho)
-off = array( [ [ 0.0 , 0.0 , 0.0 ] , [ 0.0 , 0.5 , 0.5 ] , [ 0.5 , 0.0 , 0.5 ] , [ 0.5 , 0.5 , 0.0 ] ] );
-hbox = boxSize / L
+coords = zeros((numPart, 3))
+h = zeros(numPart)
+vol = 1.
-# if L%2 == 0:
-# print "Number of particles along each dimension must be odd."
-# exit()
+for i in range(numPart):
+ coords[i,0] = i * vol/numPart + vol/(2.*numPart)
+ h[i] = 1.2348 * vol / numPart
-#Generate particles
-coords = zeros((numPart, 3))
-v = zeros((numPart, 3))
-m = zeros(numPart)
-h = zeros(numPart)
-u = zeros(numPart)
-ids = zeros(numPart, dtype='L')
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+r = zeros(numPart)
-for i in range(L):
- for j in range(L):
- for k in range(L):
- x = (i + 0.25) * hbox
- y = (j + 0.25) * hbox
- z = (k + 0.25) * hbox
- for ell in range(4):
- index = 4*(i*L*L + j*L + k) + ell
- coords[index,0] = x + off[ell,0] * hbox
- coords[index,1] = y + off[ell,1] * hbox
- coords[index,2] = z + off[ell,2] * hbox
- v[index,0] = 0.
- v[index,1] = 0.
- v[index,2] = 0.
- m[index] = mass
- h[index] = eta * hbox
- u[index] = internalEnergy
- ids[index] = index + 1
- if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox:
- u[index] = u[index] + E0 / (28. * mass)
- print "Particle " , index , " set to detonate."
- coords[index,0] += random.random() * pert * hbox
- coords[index,1] += random.random() * pert * hbox
- coords[index,2] += random.random() * pert * hbox
+r = abs(coords[:,0] - 0.5)
+m[:] = rho0 * vol / numPart
+u[:] = P0 / (rho0 * (gamma - 1))
+# Make the central particle detonate
+index = argsort(r)
+u[index[0:N_inject]] = E0 / (N_inject * m[0])
#--------------------------------------------------
@@ -89,7 +62,7 @@ file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
+grp.attrs["BoxSize"] = [1., 1., 1.]
grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
@@ -100,7 +73,7 @@ grp.attrs["Flag_Entropy_ICs"] = 0
#Runtime parameters
grp = file.create_group("/RuntimePars")
-grp.attrs["PeriodicBoundariesOn"] = periodic
+grp.attrs["PeriodicBoundariesOn"] = 1
#Units
grp = file.create_group("/Units")
diff --git a/examples/SedovBlast_1D/plotSolution.py b/examples/SedovBlast_1D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..a62775b012edda3217558031c266ed6e9b48f423
--- /dev/null
+++ b/examples/SedovBlast_1D/plotSolution.py
@@ -0,0 +1,279 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the 2D Sedov blast wave.
+# The script works for a given initial box and dumped energy and computes the solution at a later time t.
+
+# Parameters
+rho_0 = 1. # Background Density
+P_0 = 1.e-6 # Background Pressure
+E_0 = 1. # Energy of the explosion
+gas_gamma = 5./3. # Gas polytropic index
+
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("sedov_%03d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+pos = sim["/PartType0/Coordinates"][:,:]
+x = pos[:,0] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:,:]
+r = abs(x)
+v_r = x * vel[:,0] / r
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+
+# Now, work our the solution....
+
+from scipy.special import gamma as Gamma
+from numpy import *
+
+def calc_a(g,nu=3):
+ """
+ exponents of the polynomials of the sedov solution
+ g - the polytropic gamma
+ nu - the dimension
+ """
+ a = [0]*8
+
+ a[0] = 2.0 / (nu + 2)
+ a[2] = (1-g) / (2*(g-1) + nu)
+ a[3] = nu / (2*(g-1) + nu)
+ a[5] = 2 / (g-2)
+ a[6] = g / (2*(g-1) + nu)
+
+ a[1] = (((nu+2)*g)/(2.0+nu*(g-1.0)) ) * ( (2.0*nu*(2.0-g))/(g*(nu+2.0)**2) - a[2])
+ a[4] = a[1]*(nu+2) / (2-g)
+ a[7] = (2 + nu*(g-1))*a[1]/(nu*(2-g))
+ return a
+
+def calc_beta(v, g, nu=3):
+ """
+ beta values for the sedov solution (coefficients of the polynomials of the similarity variables)
+ v - the similarity variable
+ g - the polytropic gamma
+ nu- the dimension
+ """
+
+ beta = (nu+2) * (g+1) * array((0.25, (g/(g-1))*0.5,
+ -(2 + nu*(g-1))/2.0 / ((nu+2)*(g+1) -2*(2 + nu*(g-1))),
+ -0.5/(g-1)), dtype=float64)
+
+ beta = outer(beta, v)
+
+ beta += (g+1) * array((0.0, -1.0/(g-1),
+ (nu+2) / ((nu+2)*(g+1) -2.0*(2 + nu*(g-1))),
+ 1.0/(g-1)), dtype=float64).reshape((4,1))
+
+ return beta
+
+
+def sedov(t, E0, rho0, g, n=1000, nu=3):
+ """
+ solve the sedov problem
+ t - the time
+ E0 - the initial energy
+ rho0 - the initial density
+ n - number of points (10000)
+ nu - the dimension
+ g - the polytropic gas gamma
+ """
+ # the similarity variable
+ v_min = 2.0 / ((nu + 2) * g)
+ v_max = 4.0 / ((nu + 2) * (g + 1))
+
+ v = v_min + arange(n) * (v_max - v_min) / (n - 1.0)
+
+ a = calc_a(g, nu)
+ beta = calc_beta(v, g=g, nu=nu)
+ lbeta = log(beta)
+
+ r = exp(-a[0] * lbeta[0] - a[2] * lbeta[1] - a[1] * lbeta[2])
+ rho = ((g + 1.0) / (g - 1.0)) * exp(a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2])
+ p = exp(nu * a[0] * lbeta[0] + (a[5] + 1) * lbeta[3] + (a[4] - 2 * a[1]) * lbeta[2])
+ u = beta[0] * r * 4.0 / ((g + 1) * (nu + 2))
+ p *= 8.0 / ((g + 1) * (nu + 2) * (nu + 2))
+
+ # we have to take extra care at v=v_min, since this can be a special point.
+ # It is not a singularity, however, the gradients of our variables (wrt v) are.
+ # r -> 0, u -> 0, rho -> 0, p-> constant
+
+ u[0] = 0.0; rho[0] = 0.0; r[0] = 0.0; p[0] = p[1]
+
+ # volume of an n-sphere
+ vol = (pi ** (nu / 2.0) / Gamma(nu / 2.0 + 1)) * power(r, nu)
+
+ # note we choose to evaluate the integral in this way because the
+ # volumes of the first few elements (i.e near v=vmin) are shrinking
+ # very slowly, so we dramatically improve the error convergence by
+ # finding the volumes exactly. This is most important for the
+ # pressure integral, as this is on the order of the volume.
+
+ # (dimensionless) energy of the model solution
+ de = rho * u * u * 0.5 + p / (g - 1)
+ # integrate (trapezium rule)
+ q = inner(de[1:] + de[:-1], diff(vol)) * 0.5
+
+ # the factor to convert to this particular problem
+ fac = (q * (t ** nu) * rho0 / E0) ** (-1.0 / (nu + 2))
+
+ # shock speed
+ shock_speed = fac * (2.0 / (nu + 2))
+ rho_s = ((g + 1) / (g - 1)) * rho0
+ r_s = shock_speed * t * (nu + 2) / 2.0
+ p_s = (2.0 * rho0 * shock_speed * shock_speed) / (g + 1)
+ u_s = (2.0 * shock_speed) / (g + 1)
+
+ r *= fac * t
+ u *= fac
+ p *= fac * fac * rho0
+ rho *= rho0
+ return r, p, rho, u, r_s, p_s, rho_s, u_s, shock_speed
+
+
+# The main properties of the solution
+r_s, P_s, rho_s, v_s, r_shock, _, _, _, _ = sedov(time, E_0, rho_0, gas_gamma, 1000, 1)
+
+# Append points for after the shock
+r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock*1.5])
+rho_s = np.insert(rho_s, np.size(rho_s), [rho_0, rho_0])
+P_s = np.insert(P_s, np.size(P_s), [P_0, P_0])
+v_s = np.insert(v_s, np.size(v_s), [0, 0])
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+
+
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(r, v_r, '.', color='r', ms=2.)
+plot(r_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Radial~velocity}}~v_r$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-0.2, 3.8)
+
+# Density profile --------------------------------
+subplot(232)
+plot(r, rho, '.', color='r', ms=2.)
+plot(r_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=2)
+xlim(0, 1.3 * r_shock)
+ylim(-0.2, 5.2)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(r, P, '.', color='r', ms=2.)
+plot(r_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-1, 12.5)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(r, u, '.', color='r', ms=2.)
+plot(r_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-2, 22)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(r, S, '.', color='r', ms=2.)
+plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-5, 50)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Sedov blast with $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Background $\\rho_0=%.2f$"%(rho_0), fontsize=10)
+text(-0.49, 0.7, "Energy injected $E_0=%.2f$"%(E_0), fontsize=10)
+plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.4, scheme, fontsize=10)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+
+savefig("Sedov.png", dpi=200)
+
+
+
+
diff --git a/examples/SedovBlast_1D/run.sh b/examples/SedovBlast_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..4b9a84f069673bd6def3b96faec71b9d4fdd0dda
--- /dev/null
+++ b/examples/SedovBlast_1D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e sedov.hdf5 ]
+then
+ echo "Generating initial conditions for the Sedov blast example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 sedov.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 5
diff --git a/examples/SedovBlast_1D/sedov.yml b/examples/SedovBlast_1D/sedov.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6f519835d26ff5aa851ffb8999e650815c522cd3
--- /dev/null
+++ b/examples/SedovBlast_1D/sedov.yml
@@ -0,0 +1,36 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 5e-2 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sedov # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./sedov.hdf5 # The file to read
+
diff --git a/examples/SedovBlast_2D/makeIC.py b/examples/SedovBlast_2D/makeIC.py
index a30e727b4dd42c1f058827959cf12e3b4f152181..05233576f63f90b8d448aaa75fa6bfe7fce1f0e8 100644
--- a/examples/SedovBlast_2D/makeIC.py
+++ b/examples/SedovBlast_2D/makeIC.py
@@ -18,7 +18,6 @@
##############################################################################
import h5py
-import random
from numpy import *
# Generates a swift IC file for the Sedov blast test in a periodic cubic box
diff --git a/examples/SedovBlast_2D/run.sh b/examples/SedovBlast_2D/run.sh
index b6f6f73db3ff281492d2abc21eb5ded67e528516..a32c8f0d6f3116d5486fe1bd086bf8df49d06020 100755
--- a/examples/SedovBlast_2D/run.sh
+++ b/examples/SedovBlast_2D/run.sh
@@ -13,7 +13,7 @@ then
fi
# Run SWIFT
-../swift -s -t 1 sedov.yml
+../swift -s -t 1 sedov.yml 2>&1 | tee output.log
# Plot the solution
python plotSolution.py 5
diff --git a/examples/SedovBlast_3D/getGlass.sh b/examples/SedovBlast_3D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..d5c5f590ac37c9c9431d626a2ea61b0c12c1513c
--- /dev/null
+++ b/examples/SedovBlast_3D/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
diff --git a/examples/SedovBlast_3D/makeIC.py b/examples/SedovBlast_3D/makeIC.py
index e64942e8e92ee6fe67142f841f566019b1a668be..3c1e36a74b53ece8b886e2fcbe5d9178d9deefbc 100644
--- a/examples/SedovBlast_3D/makeIC.py
+++ b/examples/SedovBlast_3D/makeIC.py
@@ -1,7 +1,6 @@
###############################################################################
# This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
@@ -19,65 +18,42 @@
##############################################################################
import h5py
-import random
from numpy import *
# Generates a swift IC file for the Sedov blast test in a periodic cubic box
# Parameters
-periodic= 1 # 1 For periodic box
-boxSize = 10.
-L = 101 # Number of particles along one axis
-rho = 1. # Density
-P = 1.e-5 # Pressure
-E0= 1.e2 # Energy of the explosion
-pert = 0.1
-gamma = 5./3. # Gas adiabatic index
-eta = 1.2349 # 48 ngbs with cubic spline kernel
+gamma = 5./3. # Gas adiabatic index
+rho0 = 1. # Background density
+P0 = 1.e-6 # Background pressure
+E0= 1. # Energy of the explosion
+N_inject = 15 # Number of particles in which to inject energy
fileName = "sedov.hdf5"
-
#---------------------------------------------------
-numPart = L**3
-mass = boxSize**3 * rho / numPart
-internalEnergy = P / ((gamma - 1.)*rho)
-
-if L%2 == 0:
- print "Number of particles along each dimension must be odd."
- exit()
-
-#Generate particles
-coords = zeros((numPart, 3))
-v = zeros((numPart, 3))
-m = zeros(numPart)
-h = zeros(numPart)
-u = zeros(numPart)
-ids = zeros(numPart, dtype='L')
-
-for i in range(L):
- for j in range(L):
- for k in range(L):
- index = i*L*L + j*L + k
- x = i * boxSize / L + boxSize / (2*L)
- y = j * boxSize / L + boxSize / (2*L)
- z = k * boxSize / L + boxSize / (2*L)
- coords[index,0] = x
- coords[index,1] = y
- coords[index,2] = z
- v[index,0] = 0.
- v[index,1] = 0.
- v[index,2] = 0.
- m[index] = mass
- h[index] = eta * boxSize / L
- u[index] = internalEnergy
- ids[index] = index + 1
- if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L:
- u[index] = u[index] + E0 / (33. * mass)
- print "Particle " , index , " set to detonate."
- coords[index,0] = x + random.random() * pert * boxSize/(2.*L)
- coords[index,1] = y + random.random() * pert * boxSize/(2.*L)
- coords[index,2] = z + random.random() * pert * boxSize/(2.*L)
+glass = h5py.File("glassCube_64.hdf5", "r")
+
+# Read particle positions and h from the glass
+pos = glass["/PartType0/Coordinates"][:,:]
+h = glass["/PartType0/SmoothingLength"][:] * 0.3
+
+numPart = size(h)
+vol = 1.
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+r = zeros(numPart)
+
+r = sqrt((pos[:,0] - 0.5)**2 + (pos[:,1] - 0.5)**2 + (pos[:,2] - 0.5)**2)
+m[:] = rho0 * vol / numPart
+u[:] = P0 / (rho0 * (gamma - 1))
+# Make the central particles detonate
+index = argsort(r)
+u[index[0:N_inject]] = E0 / (N_inject * m[0])
#--------------------------------------------------
@@ -86,7 +62,7 @@ file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
+grp.attrs["BoxSize"] = [1., 1., 1.]
grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
@@ -97,7 +73,7 @@ grp.attrs["Flag_Entropy_ICs"] = 0
#Runtime parameters
grp = file.create_group("/RuntimePars")
-grp.attrs["PeriodicBoundariesOn"] = periodic
+grp.attrs["PeriodicBoundariesOn"] = 1
#Units
grp = file.create_group("/Units")
@@ -109,7 +85,7 @@ grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group
grp = file.create_group("/PartType0")
-grp.create_dataset('Coordinates', data=coords, dtype='d')
+grp.create_dataset('Coordinates', data=pos, dtype='d')
grp.create_dataset('Velocities', data=v, dtype='f')
grp.create_dataset('Masses', data=m, dtype='f')
grp.create_dataset('SmoothingLength', data=h, dtype='f')
diff --git a/examples/SedovBlast_3D/plotSolution.py b/examples/SedovBlast_3D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..f86ce17206ae1d15ff846fb14c61bbb6926e03bf
--- /dev/null
+++ b/examples/SedovBlast_3D/plotSolution.py
@@ -0,0 +1,281 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the 2D Sedov blast wave.
+# The script works for a given initial box and dumped energy and computes the solution at a later time t.
+
+# Parameters
+rho_0 = 1. # Background Density
+P_0 = 1.e-6 # Background Pressure
+E_0 = 1. # Energy of the explosion
+gas_gamma = 5./3. # Gas polytropic index
+
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("sedov_%03d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+pos = sim["/PartType0/Coordinates"][:,:]
+x = pos[:,0] - boxSize / 2
+y = pos[:,1] - boxSize / 2
+z = pos[:,2] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:,:]
+r = sqrt(x**2 + y**2 + z**2)
+v_r = (x * vel[:,0] + y * vel[:,1] + z * vel[:,2]) / r
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+
+# Now, work our the solution....
+
+from scipy.special import gamma as Gamma
+from numpy import *
+
+def calc_a(g,nu=3):
+ """
+ exponents of the polynomials of the sedov solution
+ g - the polytropic gamma
+ nu - the dimension
+ """
+ a = [0]*8
+
+ a[0] = 2.0 / (nu + 2)
+ a[2] = (1-g) / (2*(g-1) + nu)
+ a[3] = nu / (2*(g-1) + nu)
+ a[5] = 2 / (g-2)
+ a[6] = g / (2*(g-1) + nu)
+
+ a[1] = (((nu+2)*g)/(2.0+nu*(g-1.0)) ) * ( (2.0*nu*(2.0-g))/(g*(nu+2.0)**2) - a[2])
+ a[4] = a[1]*(nu+2) / (2-g)
+ a[7] = (2 + nu*(g-1))*a[1]/(nu*(2-g))
+ return a
+
+def calc_beta(v, g, nu=3):
+ """
+ beta values for the sedov solution (coefficients of the polynomials of the similarity variables)
+ v - the similarity variable
+ g - the polytropic gamma
+ nu- the dimension
+ """
+
+ beta = (nu+2) * (g+1) * array((0.25, (g/(g-1))*0.5,
+ -(2 + nu*(g-1))/2.0 / ((nu+2)*(g+1) -2*(2 + nu*(g-1))),
+ -0.5/(g-1)), dtype=float64)
+
+ beta = outer(beta, v)
+
+ beta += (g+1) * array((0.0, -1.0/(g-1),
+ (nu+2) / ((nu+2)*(g+1) -2.0*(2 + nu*(g-1))),
+ 1.0/(g-1)), dtype=float64).reshape((4,1))
+
+ return beta
+
+
+def sedov(t, E0, rho0, g, n=1000, nu=3):
+ """
+ solve the sedov problem
+ t - the time
+ E0 - the initial energy
+ rho0 - the initial density
+ n - number of points (10000)
+ nu - the dimension
+ g - the polytropic gas gamma
+ """
+ # the similarity variable
+ v_min = 2.0 / ((nu + 2) * g)
+ v_max = 4.0 / ((nu + 2) * (g + 1))
+
+ v = v_min + arange(n) * (v_max - v_min) / (n - 1.0)
+
+ a = calc_a(g, nu)
+ beta = calc_beta(v, g=g, nu=nu)
+ lbeta = log(beta)
+
+ r = exp(-a[0] * lbeta[0] - a[2] * lbeta[1] - a[1] * lbeta[2])
+ rho = ((g + 1.0) / (g - 1.0)) * exp(a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2])
+ p = exp(nu * a[0] * lbeta[0] + (a[5] + 1) * lbeta[3] + (a[4] - 2 * a[1]) * lbeta[2])
+ u = beta[0] * r * 4.0 / ((g + 1) * (nu + 2))
+ p *= 8.0 / ((g + 1) * (nu + 2) * (nu + 2))
+
+ # we have to take extra care at v=v_min, since this can be a special point.
+ # It is not a singularity, however, the gradients of our variables (wrt v) are.
+ # r -> 0, u -> 0, rho -> 0, p-> constant
+
+ u[0] = 0.0; rho[0] = 0.0; r[0] = 0.0; p[0] = p[1]
+
+ # volume of an n-sphere
+ vol = (pi ** (nu / 2.0) / Gamma(nu / 2.0 + 1)) * power(r, nu)
+
+ # note we choose to evaluate the integral in this way because the
+ # volumes of the first few elements (i.e near v=vmin) are shrinking
+ # very slowly, so we dramatically improve the error convergence by
+ # finding the volumes exactly. This is most important for the
+ # pressure integral, as this is on the order of the volume.
+
+ # (dimensionless) energy of the model solution
+ de = rho * u * u * 0.5 + p / (g - 1)
+ # integrate (trapezium rule)
+ q = inner(de[1:] + de[:-1], diff(vol)) * 0.5
+
+ # the factor to convert to this particular problem
+ fac = (q * (t ** nu) * rho0 / E0) ** (-1.0 / (nu + 2))
+
+ # shock speed
+ shock_speed = fac * (2.0 / (nu + 2))
+ rho_s = ((g + 1) / (g - 1)) * rho0
+ r_s = shock_speed * t * (nu + 2) / 2.0
+ p_s = (2.0 * rho0 * shock_speed * shock_speed) / (g + 1)
+ u_s = (2.0 * shock_speed) / (g + 1)
+
+ r *= fac * t
+ u *= fac
+ p *= fac * fac * rho0
+ rho *= rho0
+ return r, p, rho, u, r_s, p_s, rho_s, u_s, shock_speed
+
+
+# The main properties of the solution
+r_s, P_s, rho_s, v_s, r_shock, _, _, _, _ = sedov(time, E_0, rho_0, gas_gamma, 1000, 3)
+
+# Append points for after the shock
+r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock*1.5])
+rho_s = np.insert(rho_s, np.size(rho_s), [rho_0, rho_0])
+P_s = np.insert(P_s, np.size(P_s), [P_0, P_0])
+v_s = np.insert(v_s, np.size(v_s), [0, 0])
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+
+
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(r, v_r, '.', color='r', ms=1.)
+plot(r_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Radial~velocity}}~v_r$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-0.2, 3.8)
+
+# Density profile --------------------------------
+subplot(232)
+plot(r, rho, '.', color='r', ms=1.)
+plot(r_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=2)
+xlim(0, 1.3 * r_shock)
+ylim(-0.2, 5.2)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(r, P, '.', color='r', ms=1.)
+plot(r_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-1, 12.5)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(r, u, '.', color='r', ms=1.)
+plot(r_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-2, 22)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(r, S, '.', color='r', ms=1.)
+plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-5, 50)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Sedov blast with $\\gamma=%.3f$ in 2D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Background $\\rho_0=%.2f$"%(rho_0), fontsize=10)
+text(-0.49, 0.7, "Energy injected $E_0=%.2f$"%(E_0), fontsize=10)
+plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.4, scheme, fontsize=10)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+
+savefig("Sedov.png", dpi=200)
+
+
+
+
diff --git a/examples/SedovBlast_3D/profile.py b/examples/SedovBlast_3D/profile.py
deleted file mode 100644
index 1b1dfb389b2c9f6d0e5af02a17a07068f03bfd92..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_3D/profile.py
+++ /dev/null
@@ -1,46 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- #
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- #
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- # GNU General Public License for more details.
- #
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program. If not, see <http://www.gnu.org/licenses/>.
- #
- ##############################################################################
-
-import numpy as np
-import h5py
-import matplotlib
-matplotlib.use('Agg')
-import pylab as pl
-import glob
-import sedov
-
-# plot the radial density profile for all snapshots in the current directory,
-# as well as the analytical solution of Sedov
-
-for f in sorted(glob.glob("output*.hdf5")):
- file = h5py.File(f, "r")
- t = file["/Header"].attrs["Time"]
- coords = np.array(file["/PartType0/Coordinates"])
- rho = np.array(file["/PartType0/Density"])
-
- radius = np.sqrt( (coords[:,0]-5.)**2 + (coords[:,1]-5.)**2 + \
- (coords[:,2]-5.)**2 )
-
- r_theory, rho_theory = sedov.get_analytical_solution(100., 5./3., 3, t,
- max(radius))
-
- pl.plot(r_theory, rho_theory, "r-")
- pl.plot(radius, rho, "k.")
- pl.savefig("{name}.png".format(name = f[:-5]))
- pl.close()
diff --git a/examples/SedovBlast_3D/rdf.py b/examples/SedovBlast_3D/rdf.py
deleted file mode 100644
index 7f932cc814dc36e14e1bef52e33cf5ed1f527dfd..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_3D/rdf.py
+++ /dev/null
@@ -1,121 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- #
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- #
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- # GNU General Public License for more details.
- #
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program. If not, see <http://www.gnu.org/licenses/>.
- #
- ##############################################################################
-
-import h5py
-import random
-import sys
-import math
-from numpy import *
-
-# Reads the HDF5 output of SWIFT and generates a radial density profile
-# of the different physical values.
-
-# Input values?
-if len(sys.argv) < 3 :
- print "Usage: " , sys.argv[0] , " <filename> <nr. bins>"
- exit()
-
-# Get the input arguments
-fileName = sys.argv[1];
-nr_bins = int( sys.argv[2] );
-
-
-# Open the file
-fileName = sys.argv[1];
-file = h5py.File( fileName , 'r' )
-
-# Get the space dimensions.
-grp = file[ "/Header" ]
-boxsize = grp.attrs[ 'BoxSize' ]
-if len( boxsize ) > 0:
- boxsize = boxsize[0]
-
-# Get the particle data
-grp = file.get( '/PartType0' )
-ds = grp.get( 'Coordinates' )
-coords = ds[...]
-ds = grp.get( 'Velocities' )
-v = ds[...]
-ds = grp.get( 'Masses' )
-m = ds[...]
-ds = grp.get( 'SmoothingLength' )
-h = ds[...]
-ds = grp.get( 'InternalEnergy' )
-u = ds[...]
-ds = grp.get( 'ParticleIDs' )
-ids = ds[...]
-ds = grp.get( 'Density' )
-rho = ds[...]
-
-# Set the origin to be the middle of the box
-origin = [ boxsize / 2 , boxsize / 2 , boxsize / 2 ]
-
-# Get the maximum radius
-r_max = boxsize / 2
-
-# Init the bins
-nr_parts = coords.shape[0]
-bins_v = zeros( nr_bins )
-bins_m = zeros( nr_bins )
-bins_h = zeros( nr_bins )
-bins_u = zeros( nr_bins )
-bins_rho = zeros( nr_bins )
-bins_count = zeros( nr_bins )
-bins_P = zeros( nr_bins )
-
-# Loop over the particles and fill the bins.
-for i in range( nr_parts ):
-
- # Get the box index.
- r = coords[i] - origin
- r = sqrt( r[0]*r[0] + r[1]*r[1] + r[2]*r[2] )
- ind = floor( r / r_max * nr_bins )
-
- # Update the bins
- if ( ind < nr_bins ):
- bins_count[ind] += 1
- bins_v[ind] += sqrt( v[i,0]*v[i,0] + v[i,1]*v[i,1] + v[i,2]*v[i,2] )
- bins_m[ind] += m[i]
- bins_h[ind] += h[i]
- bins_u[ind] += u[i]
- bins_rho[ind] += rho[i]
- bins_P[ind] += (2.0/3)*u[i]*rho[i]
-
-# Loop over the bins and dump them
-print "# bucket left right count v m h u rho"
-for i in range( nr_bins ):
-
- # Normalize by the bin volume.
- r_left = r_max * i / nr_bins
- r_right = r_max * (i+1) / nr_bins
- vol = 4/3*math.pi*(r_right*r_right*r_right - r_left*r_left*r_left)
- ivol = 1.0 / vol
-
- print "%i %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e" % \
- ( i , r_left , r_right , \
- bins_count[i] * ivol , \
- bins_v[i] / bins_count[i] , \
- bins_m[i] * ivol , \
- bins_h[i] / bins_count[i] , \
- bins_u[i] / bins_count[i] , \
- bins_rho[i] / bins_count[i] ,
- bins_P[i] / bins_count[i] )
-
-
diff --git a/examples/SedovBlast_3D/run.sh b/examples/SedovBlast_3D/run.sh
index f71830eb6a66a9ea84e93fd1bed1261b1cb42b7b..00d5e5b91c31e64f824a3d2a28c8e1a126684a74 100755
--- a/examples/SedovBlast_3D/run.sh
+++ b/examples/SedovBlast_3D/run.sh
@@ -1,10 +1,19 @@
#!/bin/bash
-# Generate the initial conditions if they are not present.
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassCube_64.hdf5 ]
+then
+ echo "Fetching initial glass file for the Sedov blast example..."
+ ./getGlass.sh
+fi
if [ ! -e sedov.hdf5 ]
then
- echo "Generating initial conditions for the SedovBlast example..."
- python makeIC_fcc.py
+ echo "Generating initial conditions for the Sedov blast example..."
+ python makeIC.py
fi
-../swift -s -t 16 sedov.yml
+# Run SWIFT
+../swift -s -t 4 sedov.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 5
diff --git a/examples/SedovBlast_3D/sedov.py b/examples/SedovBlast_3D/sedov.py
deleted file mode 100755
index ca908fe65cd15ddd2b0e4fb41a402492432a0690..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_3D/sedov.py
+++ /dev/null
@@ -1,212 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- #
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- #
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- # GNU General Public License for more details.
- #
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program. If not, see <http://www.gnu.org/licenses/>.
- #
- ##############################################################################
-
-import numpy as np
-import scipy.integrate as integrate
-import scipy.optimize as optimize
-import os
-
-# Calculate the analytical solution of the Sedov-Taylor shock wave for a given
-# number of dimensions, gamma and time. We assume dimensionless units and a
-# setup with constant density 1, pressure and velocity 0. An energy 1 is
-# inserted in the center at t 0.
-#
-# The solution is a self-similar shock wave, which was described in detail by
-# Sedov (1959). We follow his notations and solution method.
-#
-# The position of the shock at time t is given by
-# r2 = (E/rho1)^(1/(2+nu)) * t^(2/(2+nu))
-# the energy E is related to the inserted energy E0 by E0 = alpha*E, with alpha
-# a constant which has to be calculated by imposing energy conservation.
-#
-# The density for a given radius at a certain time is determined by introducing
-# the dimensionless position coordinate lambda = r/r2. The density profile as
-# a function of lambda is constant in time and given by
-# rho = rho1 * R(V, gamma, nu)
-# and
-# V = V(lambda, gamma, nu)
-#
-# The function V(lambda, gamma, nu) is found by solving a differential equation
-# described in detail in Sedov (1959) chapter 4, section 5. Alpha is calculated
-# from the integrals in section 11 of the same chapter.
-#
-# Numerically, the complete solution requires the use of 3 quadratures and 1
-# root solver, which are implemented using the GNU Scientific Library (GSL).
-# Since some quadratures call functions that themselves contain quadratures,
-# the problem is highly non-linear and complex and takes a while to converge.
-# Therefore, we tabulate the alpha values and profile values the first time
-# a given set of gamma and nu values is requested and reuse these tabulated
-# values.
-#
-# Reference:
-# Sedov (1959): Sedov, L., Similitude et dimensions en mecanique (7th ed.;
-# Moscou: Editions Mir) - french translation of the original
-# book from 1959.
-
-# dimensionless variable z = gamma*P/R as a function of V (and gamma and nu)
-# R is a dimensionless density, while P is a dimensionless pressure
-# z is hence a sort of dimensionless sound speed
-# The expression below corresponds to eq. 11.9 in Sedov (1959), chapter 4
-def _z(V, gamma, nu):
- if V == 2./(nu+2.)/gamma:
- return 0.
- else:
- return (gamma-1.)*V*V*(V-2./(2.+nu))/2./(2./(2.+nu)/gamma-V)
-
-# differential equation that needs to be solved to obtain lambda for a given V
-# corresponds to eq. 5.11 in Sedov (1959), chapter 4 (omega = 0)
-def _dlnlambda_dV(V, gamma, nu):
- nom = _z(V, gamma, nu) - (V-2./(nu+2.))*(V-2./(nu+2.))
- denom = V*(V-1.)*(V-2./(nu+2.))+nu*(2./(nu+2.)/gamma-V)*_z(V, gamma, nu)
- return nom/denom
-
-# dimensionless variable lambda = r/r2 as a function of V (and gamma and nu)
-# found by solving differential equation 5.11 in Sedov (1959), chapter 4
-# (omega = 0)
-def _lambda(V, gamma, nu):
- if V == 2./(nu+2.)/gamma:
- return 0.
- else:
- V0 = 4./(nu+2.)/(gamma+1.)
- integral, err = integrate.quad(_dlnlambda_dV, V, V0, (gamma, nu),
- limit = 8000)
- return np.exp(-integral)
-
-# dimensionless variable R = rho/rho1 as a function of V (and gamma and nu)
-# found by inverting eq. 5.12 in Sedov (1959), chapter 4 (omega = 0)
-# the integration constant C is found by inserting the R, V and z values
-# at the shock wave, where lambda is 1. These correspond to eq. 11.8 in Sedov
-# (1959), chapter 4.
-def _R(V, gamma, nu):
- if V == 2./(nu+2.)/gamma:
- return 0.
- else:
- C = 8.*gamma*(gamma-1.)/(nu+2.)/(nu+2.)/(gamma+1.)/(gamma+1.) \
- *((gamma-1.)/(gamma+1.))**(gamma-2.) \
- *(4./(nu+2.)/(gamma+1.)-2./(nu+2.))
- lambda1 = _lambda(V, gamma, nu)
- lambda5 = lambda1**(nu+2)
- return (_z(V, gamma, nu)*(V-2./(nu+2.))*lambda5/C)**(1./(gamma-2.))
-
-# function of which we need to find the zero point to invert lambda(V)
-def _lambda_min_lambda(V, lambdax, gamma, nu):
- return _lambda(V, gamma, nu) - lambdax
-
-# dimensionless variable V = v*t/r as a function of lambda (and gamma and nu)
-# found by inverting the function lamdba(V) which is found by solving
-# differential equation 5.11 in Sedov (1959), chapter 4 (omega = 0)
-# the invertion is done by searching the zero point of the function
-# lambda_min_lambda defined above
-def _V_inv(lambdax, gamma, nu):
- if lambdax == 0.:
- return 2./(2.+nu)/gamma;
- else:
- return optimize.brentq(_lambda_min_lambda, 2./(nu+2.)/gamma,
- 4./(nu+2.)/(gamma+1.), (lambdax, gamma, nu))
-
-# integrand of the first integral in eq. 11.24 in Sedov (1959), chapter 4
-def _integrandum1(lambdax, gamma, nu):
- V = _V_inv(lambdax, gamma, nu)
- if nu == 2:
- return _R(V, gamma, nu)*V**2*lambdax**3
- else:
- return _R(V, gamma, nu)*V**2*lambdax**4
-
-# integrand of the second integral in eq. 11.24 in Sedov (1959), chapter 4
-def _integrandum2(lambdax, gamma, nu):
- V = _V_inv(lambdax, gamma, nu)
- if V == 2./(nu+2.)/gamma:
- P = 0.
- else:
- P = _z(V, gamma, nu)*_R(V, gamma, nu)/gamma
- if nu == 2:
- return P*lambdax**3
- else:
- return P*lambdax**4
-
-# calculate alpha = E0/E
-# this corresponds to eq. 11.24 in Sedov (1959), chapter 4
-def get_alpha(gamma, nu):
- integral1, err1 = integrate.quad(_integrandum1, 0., 1., (gamma, nu))
- integral2, err2 = integrate.quad(_integrandum2, 0., 1., (gamma, nu))
-
- if nu == 2:
- return np.pi*integral1+2.*np.pi/(gamma-1.)*integral2
- else:
- return 2.*np.pi*integral1+4.*np.pi/(gamma-1.)*integral2
-
-# get the analytical solution for the Sedov-Taylor blastwave given an input
-# energy E, adiabatic index gamma, and number of dimensions nu, at time t and
-# with a maximal outer region radius maxr
-def get_analytical_solution(E, gamma, nu, t, maxr = 1.):
- # we check for the existance of a datafile with precomputed alpha and
- # profile values
- # if it does not exist, we calculate it here and write it out
- # calculation of alpha and the profile takes a while...
- lvec = np.zeros(1000)
- Rvec = np.zeros(1000)
- fname = "sedov_profile_gamma_{gamma}_nu_{nu}.dat".format(gamma = gamma,
- nu = nu)
- if os.path.exists(fname):
- file = open(fname, "r")
- lines = file.readlines()
- alpha = float(lines[0])
- for i in range(len(lines)-1):
- data = lines[i+1].split()
- lvec[i] = float(data[0])
- Rvec[i] = float(data[1])
- else:
- alpha = get_alpha(gamma, nu)
- for i in range(1000):
- lvec[i] = (i+1)*0.001
- V = _V_inv(lvec[i], gamma, nu)
- Rvec[i] = _R(V, gamma, nu)
- file = open(fname, "w")
- file.write("{alpha}\n".format(alpha = alpha))
- for i in range(1000):
- file.write("{l}\t{R}\n".format(l = lvec[i], R = Rvec[i]))
-
- xvec = np.zeros(1002)
- rhovec = np.zeros(1002)
- if nu == 2:
- r2 = (E/alpha)**0.25*np.sqrt(t)
- else:
- r2 = (E/alpha)**0.2*t**0.4
-
- for i in range(1000):
- xvec[i] = lvec[i]*r2
- rhovec[i] = Rvec[i]
- xvec[1000] = 1.001*r2
- xvec[1001] = maxr
- rhovec[1000] = 1.
- rhovec[1001] = 1.
-
- return xvec, rhovec
-
-def main():
- E = 1.
- gamma = 1.66667
- nu = 3
- t = 0.001
- x, rho = get_analytical_solution(E, gamma, nu, t)
- for i in range(len(x)):
- print x[i], rho[i]
-
-if __name__ == "__main__":
- main()
diff --git a/examples/SedovBlast_3D/sedov.yml b/examples/SedovBlast_3D/sedov.yml
index eb95fd85d5599145b2ff037256dcde72fc306209..6f519835d26ff5aa851ffb8999e650815c522cd3 100644
--- a/examples/SedovBlast_3D/sedov.yml
+++ b/examples/SedovBlast_3D/sedov.yml
@@ -9,15 +9,15 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1. # The end time of the simulation (in internal units).
+ time_end: 5e-2 # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: sedov # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
+ delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
@@ -27,7 +27,7 @@ Statistics:
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
- max_smoothing_length: 1. # Maximal smoothing length allowed (in internal units).
+ max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
diff --git a/examples/SedovBlast_3D/solution.py b/examples/SedovBlast_3D/solution.py
deleted file mode 100644
index 9335e22bdd76585e8d53d3dba9f9a435197f92fc..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_3D/solution.py
+++ /dev/null
@@ -1,114 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- #
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- #
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- # GNU General Public License for more details.
- #
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program. If not, see <http://www.gnu.org/licenses/>.
- #
- ##############################################################################
-
-import random
-from numpy import *
-
-# Computes the analytical solution of the 3D Sedov blast wave.
-# The script works for a given initial box and dumped energy and computes the solution at a later time t.
-# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
-# entropic function on N radial points between r=0 and r=R_max.
-# Follows the solution in Landau & Lifschitz
-
-# Parameters
-rho_0 = 1. # Background Density
-P_0 = 1.e-5 # Background Pressure
-E_0 = 1.e2 # Energy of the explosion
-gamma = 5./3. # Gas polytropic index
-
-t = 0.275 # Time of the solution
-
-N = 1000 # Number of radial points
-R_max = 3. # Maximal radius
-
-# ---------------------------------------------------------------
-# Don't touch anything after this.
-# ---------------------------------------------------------------
-
-if gamma == 5./3.:
- alpha = 0.49
-else:
- print "Unknown value for alpha"
- exit
-
-# Position and velocity of the shock
-r_shock = (E_0 / (alpha * rho_0))**(1./5.) * t**(2./5.)
-v_shock = (1./5.) * (1./alpha)**(1./5.) * ((E_0 * t**2. / rho_0)**(-4./5.)) * E_0 * (t / rho_0)
-
-# Prepare arrays
-delta_r = R_max / N
-r_s = arange(0, R_max, delta_r)
-rho_s = ones(N) * rho_0
-P_s = ones(N) * P_0
-u_s = ones(N)
-v_s = zeros(N)
-
-# State on the shock
-rho_shock = rho_0 * (gamma+1.)/(gamma-1.)
-P_shock = 2./(gamma+1.) * rho_shock * v_shock**2
-
-# Integer position of the shock
-i_shock = min(floor(r_shock /delta_r), N)
-
-# Dimensionless velocity and its spatial derivative
-v_bar0 = (gamma+1.)/(2.*gamma)
-deltaV_bar = (1.0 - v_bar0) / (i_shock - 1)
-
-def rho_dimensionless(v_bar):
- power1 = (2./(gamma-2.))
- power2 = -(12.-7.*gamma+13.*gamma**2)/(2.-3.*gamma-11.*gamma**2+6.*gamma**3)
- power3 = 3./(1.+2.*gamma)
- term1 = ((1.+ gamma - 2.*v_bar)/(gamma-1.))**power1
- term2 = ((5.+5.*gamma+2.*v_bar-6.*gamma*v_bar)/(7.-gamma))**power2
- term3 = ((2.*gamma*v_bar - gamma -1.)/(gamma-1.))**power3
- return term1 * term2 * term3
-
-def P_dimensionless(v_bar):
- return (gamma+1. - 2.*v_bar)/(2.*gamma*v_bar - gamma - 1.)*v_bar**2 * rho_dimensionless(v_bar)
-
-def r_dimensionless(v_bar):
- power1 = (-12.+7.*gamma-13.*gamma**2)/(5.*(-1.+gamma+6.*gamma**2))
- power2 = (gamma - 1.)/(1.+2.*gamma)
- term1 = ((5.+5.*gamma+2.*v_bar-6.*gamma*v_bar)/(7.-gamma))**power1
- term2 = ((2.*gamma*v_bar-gamma-1.)/(gamma-1.))**power2
- return v_bar**(-2./5.)*term1*term2
-
-# Generate solution
-for i in range(1,int(i_shock)):
- v_bar = v_bar0 + (i-1)*deltaV_bar
- r_s[i] = r_dimensionless(v_bar) * r_shock
- rho_s[i] = rho_dimensionless(v_bar) * rho_shock
- P_s[i] = P_shock * (r_s[i]/r_shock)**2 * P_dimensionless(v_bar)
- v_s[i] = (4./(5.*(gamma+1.)) * r_s[i] / t * v_bar)
-
-u_s = P_s / ((gamma - 1.)*rho_s) # internal energy
-s_s = P_s / rho_s**gamma # entropic function
-rho_s[0] = 0.
-P_s[0] = P_s[1] # dirty...
-
-
-savetxt("rho.dat", column_stack((r_s, rho_s)))
-savetxt("P.dat", column_stack((r_s, P_s)))
-savetxt("v.dat", column_stack((r_s, v_s)))
-savetxt("u.dat", column_stack((r_s, u_s)))
-savetxt("s.dat", column_stack((r_s, s_s)))
-
-
-
diff --git a/examples/SodShock_1D/plotSolution.py b/examples/SodShock_1D/plotSolution.py
index 57f66fe29b9be86a3c4de0d90eafe615d0cb2dbb..0a7720f4a6cf26e5a8acda1101bd438850d8d553 100644
--- a/examples/SodShock_1D/plotSolution.py
+++ b/examples/SodShock_1D/plotSolution.py
@@ -22,8 +22,6 @@
# Generates the analytical solution for the Sod shock test case
# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
-# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
-# entropic function on N points between x_min and x_max.
# This follows the solution given in (Toro, 2009)
diff --git a/examples/SodShock_1D/run.sh b/examples/SodShock_1D/run.sh
index e3ac218c56caa81d0e7a6817e03a8db20bb575d5..4be4254baa4a87b105a5f3c1bfbf9059348a1e9e 100755
--- a/examples/SodShock_1D/run.sh
+++ b/examples/SodShock_1D/run.sh
@@ -8,7 +8,7 @@ then
fi
# Run SWIFT
-../swift -s -t 1 sodShock.yml
+../swift -s -t 1 sodShock.yml 2>&1 | tee output.log
# Plot the result
-python plotSolution.py 1
+python plotSolution.py 1
diff --git a/examples/SodShock_2D/plotSolution.py b/examples/SodShock_2D/plotSolution.py
index e74651c567bb711b3190662cf78d421a66134775..99ba7e9a6e9ae4b6d50688a1428f07e9a08b3b85 100644
--- a/examples/SodShock_2D/plotSolution.py
+++ b/examples/SodShock_2D/plotSolution.py
@@ -22,8 +22,6 @@
# Generates the analytical solution for the Sod shock test case
# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
-# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
-# entropic function on N points between x_min and x_max.
# This follows the solution given in (Toro, 2009)
diff --git a/examples/SodShock_2D/run.sh b/examples/SodShock_2D/run.sh
index 3d1f4beb3ac4c4becfda24ed39390efec8074da7..9e6bbfdf1c0a7c206ce6966fdca7b20a28047dd8 100755
--- a/examples/SodShock_2D/run.sh
+++ b/examples/SodShock_2D/run.sh
@@ -13,6 +13,6 @@ then
fi
# Run SWIFT
-../swift -s -t 1 sodShock.yml
+../swift -s -t 1 sodShock.yml 2>&1 | tee output.log
python plotSolution.py 1
diff --git a/examples/SodShock_3D/getGlass.sh b/examples/SodShock_3D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..f61b61d4e6c51b44576fd7cdd6242cb9f0133039
--- /dev/null
+++ b/examples/SodShock_3D/getGlass.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/SodShock_3D/glass_001.hdf5 b/examples/SodShock_3D/glass_001.hdf5
deleted file mode 100644
index a371826c2ef4c1d53ad5f50a6bc7eb590017220e..0000000000000000000000000000000000000000
Binary files a/examples/SodShock_3D/glass_001.hdf5 and /dev/null differ
diff --git a/examples/SodShock_3D/glass_002.hdf5 b/examples/SodShock_3D/glass_002.hdf5
deleted file mode 100644
index dffb8d343157a9ae8318e9572fc752eecd8955fb..0000000000000000000000000000000000000000
Binary files a/examples/SodShock_3D/glass_002.hdf5 and /dev/null differ
diff --git a/examples/SodShock_3D/makeIC.py b/examples/SodShock_3D/makeIC.py
index 8ae19050c11c0712579b44646c8870d7574d113b..84283732afc497825417546be8bc25e183ecb1cb 100644
--- a/examples/SodShock_3D/makeIC.py
+++ b/examples/SodShock_3D/makeIC.py
@@ -1,7 +1,6 @@
###############################################################################
# This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
@@ -19,87 +18,77 @@
##############################################################################
import h5py
-import random
from numpy import *
-# Generates a swift IC file for the Sod Shock in a periodic box
+# Generates a swift IC file for the 3D Sod Shock in a periodic box
# Parameters
-periodic= 1 # 1 For periodic box
-factor = 8
-boxSize = [ 1.0 , 1.0/factor , 1.0/factor ]
-L = 100 # Number of particles along one axis
-P1 = 1. # Pressure left state
-P2 = 0.1795 # Pressure right state
-gamma = 5./3. # Gas adiabatic index
+gamma = 5./3. # Gas adiabatic index
+x_min = -1.
+x_max = 1.
+rho_L = 1. # Density left state
+rho_R = 0.125 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
fileName = "sodShock.hdf5"
-vol = boxSize[0] * boxSize[1] * boxSize[2]
#---------------------------------------------------
-
-#Read in high density glass
-# glass1 = h5py.File("../Glass/glass_200000.hdf5")
-glass1 = h5py.File("glass_001.hdf5")
-pos1 = glass1["/PartType0/Coordinates"][:,:]
-pos1 = pos1 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-glass_h1 = glass1["/PartType0/SmoothingLength"][:] / factor
-
-#Read in high density glass
-# glass2 = h5py.File("../Glass/glass_50000.hdf5")
-glass2 = h5py.File("glass_002.hdf5")
-pos2 = glass2["/PartType0/Coordinates"][:,:]
-pos2 = pos2 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-glass_h2 = glass2["/PartType0/SmoothingLength"][:] / factor
-
-#Generate high density region
-rho1 = 1.
-coord1 = append(pos1, pos1 + [0.125, 0, 0], 0)
-coord1 = append(coord1, coord1 + [0.25, 0, 0], 0)
-# coord1 = append(pos1, pos1 + [0, 0.5, 0], 0)
-# coord1 = append(coord1, pos1 + [0, 0, 0.5], 0)
-# coord1 = append(coord1, pos1 + [0, 0.5, 0.5], 0)
-N1 = size(coord1)/3
-v1 = zeros((N1, 3))
-u1 = ones(N1) * P1 / ((gamma - 1.) * rho1)
-m1 = ones(N1) * vol * 0.5 * rho1 / N1
-h1 = append(glass_h1, glass_h1, 0)
-h1 = append(h1, h1, 0)
-
-#Generate low density region
-rho2 = 0.25
-coord2 = append(pos2, pos2 + [0.125, 0, 0], 0)
-coord2 = append(coord2, coord2 + [0.25, 0, 0], 0)
-# coord2 = append(pos2, pos2 + [0, 0.5, 0], 0)
-# coord2 = append(coord2, pos2 + [0, 0, 0.5], 0)
-# coord2 = append(coord2, pos2 + [0, 0.5, 0.5], 0)
-N2 = size(coord2)/3
-v2 = zeros((N2, 3))
-u2 = ones(N2) * P2 / ((gamma - 1.) * rho2)
-m2 = ones(N2) * vol * 0.5 * rho2 / N2
-h2 = append(glass_h2, glass_h2, 0)
-h2 = append(h2, h2, 0)
-
-#Merge arrays
-numPart = N1 + N2
-coords = append(coord1, coord2+[0.5, 0., 0.], 0)
-v = append(v1, v2,0)
-h = append(h1, h2,0)
-u = append(u1, u2,0)
-m = append(m1, m2,0)
-ids = zeros(numPart, dtype='L')
-for i in range(1, numPart+1):
- ids[i-1] = i
-
-#Final operation since we come from Gadget-2 cubic spline ICs
-h /= 1.825752
+boxSize = (x_max - x_min)
+
+glass_L = h5py.File("glassCube_64.hdf5", "r")
+glass_R = h5py.File("glassCube_32.hdf5", "r")
+
+pos_L = glass_L["/PartType0/Coordinates"][:,:] * 0.5
+pos_R = glass_R["/PartType0/Coordinates"][:,:] * 0.5
+h_L = glass_L["/PartType0/SmoothingLength"][:] * 0.5
+h_R = glass_R["/PartType0/SmoothingLength"][:] * 0.5
+
+# Merge things
+aa = pos_L - array([0.5, 0., 0.])
+pos_LL = append(pos_L, pos_L + array([0.5, 0., 0.]), axis=0)
+pos_RR = append(pos_R, pos_R + array([0.5, 0., 0.]), axis=0)
+pos = append(pos_LL - array([1.0, 0., 0.]), pos_RR, axis=0)
+h_LL = append(h_L, h_L)
+h_RR = append(h_R, h_R)
+h = append(h_LL, h_RR)
+
+numPart_L = size(h_LL)
+numPart_R = size(h_RR)
+numPart = size(h)
+
+vol_L = 0.25
+vol_R = 0.25
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+
+for i in range(numPart):
+ x = pos[i,0]
+
+ if x < 0: #left
+ u[i] = P_L / (rho_L * (gamma - 1.))
+ m[i] = rho_L * vol_L / numPart_L
+ v[i,0] = v_L
+ else: #right
+ u[i] = P_R / (rho_R * (gamma - 1.))
+ m[i] = rho_R * vol_R / numPart_R
+ v[i,0] = v_R
+
+# Shift particles
+pos[:,0] -= x_min
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
+grp.attrs["BoxSize"] = [boxSize, 0.5, 0.5]
grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
@@ -110,7 +99,7 @@ grp.attrs["Flag_Entropy_ICs"] = 0
#Runtime parameters
grp = file.create_group("/RuntimePars")
-grp.attrs["PeriodicBoundariesOn"] = periodic
+grp.attrs["PeriodicBoundariesOn"] = 1
#Units
grp = file.create_group("/Units")
@@ -122,7 +111,7 @@ grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group
grp = file.create_group("/PartType0")
-grp.create_dataset('Coordinates', data=coords, dtype='d')
+grp.create_dataset('Coordinates', data=pos, dtype='d')
grp.create_dataset('Velocities', data=v, dtype='f')
grp.create_dataset('Masses', data=m, dtype='f')
grp.create_dataset('SmoothingLength', data=h, dtype='f')
@@ -131,5 +120,3 @@ grp.create_dataset('ParticleIDs', data=ids, dtype='L')
file.close()
-
-
diff --git a/examples/SodShock_3D/plotSolution.py b/examples/SodShock_3D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..23a16e6aed73a7281cf78a215940ccdcff722a79
--- /dev/null
+++ b/examples/SodShock_3D/plotSolution.py
@@ -0,0 +1,288 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the Sod shock and plots the SPH answer
+
+
+# Generates the analytical solution for the Sod shock test case
+# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
+# This follows the solution given in (Toro, 2009)
+
+
+# Parameters
+gas_gamma = 5./3. # Polytropic index
+rho_L = 1. # Density left state
+rho_R = 0.125 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("sodShock_%03d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+x = sim["/PartType0/Coordinates"][:,0]
+v = sim["/PartType0/Velocities"][:,0]
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+N = 1000 # Number of points
+x_min = -1.
+x_max = 1.
+
+x += x_min
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave
+c_R = sqrt(gas_gamma * P_R / rho_R) # Speed of the shock front
+
+# Helpful variable
+Gama = (gas_gamma - 1.) / (gas_gamma + 1.)
+beta = (gas_gamma - 1.) / (2. * gas_gamma)
+
+# Characteristic function and its derivative, following Toro (2009)
+def compute_f(P_3, P, c):
+ u = P_3 / P
+ if u > 1:
+ term1 = gas_gamma*((gas_gamma+1.)*u + gas_gamma-1.)
+ term2 = sqrt(2./term1)
+ fp = (u - 1.)*c*term2
+ dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gas_gamma*(gas_gamma+1.)/P
+ else:
+ fp = (u**beta - 1.)*(2.*c/(gas_gamma-1.))
+ dfdp = 2.*c/(gas_gamma-1.)*beta*u**(beta-1.)/P
+ return (fp, dfdp)
+
+# Solution of the Riemann problem following Toro (2009)
+def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
+ P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gas_gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta)
+ P_3 = 0.5*(P_R + P_L)
+ f_L = 1.
+ while fabs(P_3 - P_new) > 1e-6:
+ P_3 = P_new
+ (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
+ (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
+ f = f_L + f_R + (v_R - v_L)
+ df = dfdp_L + dfdp_R
+ dp = -f/df
+ prnew = P_3 + dp
+ v_3 = v_L - f_L
+ return (P_new, v_3)
+
+
+# Solve Riemann problem for post-shock region
+(P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)
+
+# Check direction of shocks and wave
+shock_R = (P_3 > P_R)
+shock_L = (P_3 > P_L)
+
+# Velocity of shock front and and rarefaction wave
+if shock_R:
+ v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gas_gamma*(v_3-v_R))
+else:
+ v_right = c_R + 0.5*(gas_gamma+1.)*v_3 - 0.5*(gas_gamma-1.)*v_R
+
+if shock_L:
+ v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gas_gamma*(v_3-v_L))
+else:
+ v_left = c_L - 0.5*(gas_gamma+1.)*v_3 + 0.5*(gas_gamma-1.)*v_L
+
+# Compute position of the transitions
+x_23 = -fabs(v_left) * time
+if shock_L :
+ x_12 = -fabs(v_left) * time
+else:
+ x_12 = -(c_L - v_L) * time
+
+x_34 = v_3 * time
+
+x_45 = fabs(v_right) * time
+if shock_R:
+ x_56 = fabs(v_right) * time
+else:
+ x_56 = (c_R + v_R) * time
+
+
+# Prepare arrays
+delta_x = (x_max - x_min) / N
+x_s = arange(x_min, x_max, delta_x)
+rho_s = zeros(N)
+P_s = zeros(N)
+v_s = zeros(N)
+
+# Compute solution in the different regions
+for i in range(N):
+ if x_s[i] <= x_12:
+ rho_s[i] = rho_L
+ P_s[i] = P_L
+ v_s[i] = v_L
+ if x_s[i] >= x_12 and x_s[i] < x_23:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ else:
+ rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * time) + Gama * v_L/c_L + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = P_L*(rho_s[i] / rho_L)**gas_gamma
+ v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / time) + Gama*v_L
+ if x_s[i] >= x_23 and x_s[i] < x_34:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L)
+ else:
+ rho_s[i] = rho_L*(P_3 / P_L)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_34 and x_s[i] < x_45:
+ if shock_R:
+ rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R)
+ else:
+ rho_s[i] = rho_R*(P_3 / P_R)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_45 and x_s[i] < x_56:
+ if shock_R:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+ else:
+ rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*time) - Gama*v_R/c_R + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = p_R*(rho_s[i]/rho_R)**gas_gamma
+ v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/time) + Gama*v_R
+ if x_s[i] >= x_56:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(x, v, '.', color='r', ms=0.5)
+plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(-0.1, 0.95)
+
+# Density profile --------------------------------
+subplot(232)
+plot(x, rho, '.', color='r', ms=0.5)
+plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.05, 1.1)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(x, P, '.', color='r', ms=0.5)
+plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.01, 1.1)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(x, u, '.', color='r', ms=0.5)
+plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 2.2)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(x, S, '.', color='r', ms=0.5)
+plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 3.8)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Sod shock with $\\gamma=%.3f$ in 3D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$"%(P_L, rho_L, v_L), fontsize=10)
+text(-0.49, 0.7, "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$"%(P_R, rho_R, v_R), fontsize=10)
+plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.4, scheme, fontsize=10)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+
+savefig("SodShock.png", dpi=200)
diff --git a/examples/SodShock_3D/rhox.py b/examples/SodShock_3D/rhox.py
deleted file mode 100644
index 70493be3728cdeb27409a79f616fa3ec5bb9cdfd..0000000000000000000000000000000000000000
--- a/examples/SodShock_3D/rhox.py
+++ /dev/null
@@ -1,115 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- #
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- #
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- # GNU General Public License for more details.
- #
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program. If not, see <http://www.gnu.org/licenses/>.
- #
- ##############################################################################
-
-import h5py
-import random
-import sys
-import math
-from numpy import *
-
-# Reads the HDF5 output of SWIFT and generates a radial density profile
-# of the different physical values.
-
-# Input values?
-if len(sys.argv) < 3 :
- print "Usage: " , sys.argv[0] , " <filename> <nr. bins>"
- exit()
-
-# Get the input arguments
-fileName = sys.argv[1];
-nr_bins = int( sys.argv[2] );
-
-
-# Open the file
-fileName = sys.argv[1];
-file = h5py.File( fileName , 'r' )
-
-# Get the space dimensions.
-grp = file[ "/Header" ]
-boxsize = grp.attrs[ 'BoxSize' ]
-boxsize = boxsize[0]
-
-# Get the particle data
-grp = file.get( '/PartType0' )
-ds = grp.get( 'Coordinates' )
-coords = ds[...]
-ds = grp.get( 'Velocities' )
-v = ds[...]
-# ds = grp.get( 'Mass' )
-# m = ds[...]
-ds = grp.get( 'SmoothingLength' )
-h = ds[...]
-ds = grp.get( 'InternalEnergy' )
-u = ds[...]
-ds = grp.get( 'ParticleIDs' )
-ids = ds[...]
-ds = grp.get( 'Density' )
-rho = ds[...]
-
-# Get the maximum radius
-r_max = boxsize
-
-# Init the bins
-nr_parts = coords.shape[0]
-bins_v = zeros( nr_bins )
-bins_m = zeros( nr_bins )
-bins_h = zeros( nr_bins )
-bins_u = zeros( nr_bins )
-bins_rho = zeros( nr_bins )
-bins_count = zeros( nr_bins )
-bins_P = zeros( nr_bins )
-
-# Loop over the particles and fill the bins.
-for i in range( nr_parts ):
-
- # Get the box index.
- r = coords[i,0]
- ind = floor( r / r_max * nr_bins )
-
- # Update the bins
- bins_count[ind] += 1
- bins_v[ind] += v[i,0] # sqrt( v[i,0]*v[i,0] + v[i,1]*v[i,1] + v[i,2]*v[i,2] )
- # bins_m[ind] += m[i]
- bins_h[ind] += h[i]
- bins_u[ind] += u[i]
- bins_rho[ind] += rho[i]
- bins_P[ind] += (2.0/3)*u[i]*rho[i]
-
-# Loop over the bins and dump them
-print "# bucket left right count v m h u rho"
-for i in range( nr_bins ):
-
- # Normalize by the bin volume.
- r_left = r_max * i / nr_bins
- r_right = r_max * (i+1) / nr_bins
- vol = 4/3*math.pi*(r_right*r_right*r_right - r_left*r_left*r_left)
- ivol = 1.0 / vol
-
- print "%i %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e" % \
- ( i , r_left , r_right , \
- bins_count[i] * ivol , \
- bins_v[i] / bins_count[i] , \
- bins_m[i] * ivol , \
- bins_h[i] / bins_count[i] , \
- bins_u[i] / bins_count[i] , \
- bins_rho[i] / bins_count[i] ,
- bins_P[i] / bins_count[i] )
-
-
diff --git a/examples/SodShock_3D/run.sh b/examples/SodShock_3D/run.sh
index b8141e51543f348d6ec6be505d136aed7d803b2e..8ed85baf73425b75f402c491a3c66785f6c6fce0 100755
--- a/examples/SodShock_3D/run.sh
+++ b/examples/SodShock_3D/run.sh
@@ -1,10 +1,18 @@
#!/bin/bash
# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_64.hdf5 ]
+then
+ echo "Fetching initial glass file for the Sod shock example..."
+ ./getGlass.sh
+fi
if [ ! -e sodShock.hdf5 ]
then
- echo "Generating initial conditions for the SodShock example..."
+ echo "Generating initial conditions for the Sod shock example..."
python makeIC.py
fi
-../swift -s -t 16 sodShock.yml
+# Run SWIFT
+../swift -s -t 4 sodShock.yml 2>&1 | tee output.log
+
+python plotSolution.py 1
diff --git a/examples/SodShock_3D/sodShock.yml b/examples/SodShock_3D/sodShock.yml
index 19bacb57f3bfb173063dbac5d752929763dede29..1ab6eb626db09678f66322e8f0e8674c0931ddb6 100644
--- a/examples/SodShock_3D/sodShock.yml
+++ b/examples/SodShock_3D/sodShock.yml
@@ -9,15 +9,15 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1. # The end time of the simulation (in internal units).
+ time_end: 0.2 # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
- basename: sod # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
+ basename: sodShock # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.2 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
@@ -27,7 +27,7 @@ Statistics:
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
- max_smoothing_length: 0.01 # Maximal smoothing length allowed (in internal units).
+ max_smoothing_length: 0.05 # Maximal smoothing length allowed (in internal units).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
diff --git a/examples/SodShock_3D/solution.py b/examples/SodShock_3D/solution.py
deleted file mode 100644
index 39f25c625232eee9bae0300339955f775f3b46ed..0000000000000000000000000000000000000000
--- a/examples/SodShock_3D/solution.py
+++ /dev/null
@@ -1,186 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- #
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- #
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- # GNU General Public License for more details.
- #
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program. If not, see <http://www.gnu.org/licenses/>.
- #
- ##############################################################################
-
-import random
-from numpy import *
-
-# Generates the analytical solution for the Sod shock test case
-# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
-# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
-# entropic function on N points between x_min and x_max.
-# This follows the solution given in (Toro, 2009)
-
-
-# Parameters
-rho_L = 1
-P_L = 1
-v_L = 0.
-
-rho_R = 0.25
-P_R = 0.1795
-v_R = 0.
-
-gamma = 5./3. # Polytropic index
-
-t = 0.12 # Time of the evolution
-
-N = 1000 # Number of points
-x_min = -0.25
-x_max = 0.25
-
-
-# ---------------------------------------------------------------
-# Don't touch anything after this.
-# ---------------------------------------------------------------
-
-c_L = sqrt(gamma * P_L / rho_L) # Speed of the rarefaction wave
-c_R = sqrt(gamma * P_R / rho_R) # Speed of the shock front
-
-# Helpful variable
-Gama = (gamma - 1.) / (gamma + 1.)
-beta = (gamma - 1.) / (2. * gamma)
-
-# Characteristic function and its derivative, following Toro (2009)
-def compute_f(P_3, P, c):
- u = P_3 / P
- if u > 1:
- term1 = gamma*((gamma+1.)*u + gamma-1.)
- term2 = sqrt(2./term1)
- fp = (u - 1.)*c*term2
- dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gamma*(gamma+1.)/P
- else:
- fp = (u**beta - 1.)*(2.*c/(gamma-1.))
- dfdp = 2.*c/(gamma-1.)*beta*u**(beta-1.)/P
- return (fp, dfdp)
-
-# Solution of the Riemann problem following Toro (2009)
-def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
- P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta)
- P_3 = 0.5*(P_R + P_L)
- f_L = 1.
- while fabs(P_3 - P_new) > 1e-6:
- P_3 = P_new
- (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
- (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
- f = f_L + f_R + (v_R - v_L)
- df = dfdp_L + dfdp_R
- dp = -f/df
- prnew = P_3 + dp
- v_3 = v_L - f_L
- return (P_new, v_3)
-
-
-# Solve Riemann problem for post-shock region
-(P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)
-
-# Check direction of shocks and wave
-shock_R = (P_3 > P_R)
-shock_L = (P_3 > P_L)
-
-# Velocity of shock front and and rarefaction wave
-if shock_R:
- v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gamma*(v_3-v_R))
-else:
- v_right = c_R + 0.5*(gamma+1.)*v_3 - 0.5*(gamma-1.)*v_R
-
-if shock_L:
- v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gamma*(v_3-v_L))
-else:
- v_left = c_L - 0.5*(gamma+1.)*v_3 + 0.5*(gamma-1.)*v_L
-
-# Compute position of the transitions
-x_23 = -fabs(v_left) * t
-if shock_L :
- x_12 = -fabs(v_left) * t
-else:
- x_12 = -(c_L - v_L) * t
-
-x_34 = v_3 * t
-
-x_45 = fabs(v_right) * t
-if shock_R:
- x_56 = fabs(v_right) * t
-else:
- x_56 = (c_R + v_R) * t
-
-
-# Prepare arrays
-delta_x = (x_max - x_min) / N
-x_s = arange(x_min, x_max, delta_x)
-rho_s = zeros(N)
-P_s = zeros(N)
-v_s = zeros(N)
-
-# Compute solution in the different regions
-for i in range(N):
- if x_s[i] <= x_12:
- rho_s[i] = rho_L
- P_s[i] = P_L
- v_s[i] = v_L
- if x_s[i] >= x_12 and x_s[i] < x_23:
- if shock_L:
- rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L)
- P_s[i] = P_3
- v_s[i] = v_3
- else:
- rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * t) + Gama * v_L/c_L + (1.-Gama))**(2./(gamma-1.))
- P_s[i] = P_L*(rho_s[i] / rho_L)**gamma
- v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / t) + Gama*v_L
- if x_s[i] >= x_23 and x_s[i] < x_34:
- if shock_L:
- rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L)
- else:
- rho_s[i] = rho_L*(P_3 / P_L)**(1./gamma)
- P_s[i] = P_3
- v_s[i] = v_3
- if x_s[i] >= x_34 and x_s[i] < x_45:
- if shock_R:
- rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R)
- else:
- rho_s[i] = rho_R*(P_3 / P_R)**(1./gamma)
- P_s[i] = P_3
- v_s[i] = v_3
- if x_s[i] >= x_45 and x_s[i] < x_56:
- if shock_R:
- rho_s[i] = rho_R
- P_s[i] = P_R
- v_s[i] = v_R
- else:
- rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*t) - Gama*v_R/c_R + (1.-Gama))**(2./(gamma-1.))
- P_s[i] = p_R*(rho_s[i]/rho_R)**gamma
- v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/t) + Gama*v_R
- if x_s[i] >= x_56:
- rho_s[i] = rho_R
- P_s[i] = P_R
- v_s[i] = v_R
-
-
-# Additional arrays
-u_s = P_s / (rho_s * (gamma - 1.)) #internal energy
-s_s = P_s / rho_s**gamma # entropic function
-
-#---------------------------------------------------------------
-# Print arrays
-
-savetxt("rho.dat", column_stack((x_s, rho_s)))
-savetxt("P.dat", column_stack((x_s, P_s)))
-savetxt("v.dat", column_stack((x_s, v_s)))
-savetxt("u.dat", column_stack((x_s, u_s)))
-savetxt("s.dat", column_stack((x_s, s_s)))
diff --git a/examples/SquareTest_2D/makeIC.py b/examples/SquareTest_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..2cb624f9944fcd421d95934ffeded089613e8bc9
--- /dev/null
+++ b/examples/SquareTest_2D/makeIC.py
@@ -0,0 +1,127 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+from numpy import *
+
+# Generates a swift IC file for the Square test in a periodic box
+
+# Parameters
+L = 64 # Number of particles on the side
+gamma = 5./3. # Gas adiabatic index
+rho0 = 4 # Gas central density
+rho1 = 1 # Gas outskirt density
+P0 = 2.5 # Gas central pressure
+P1 = 2.5 # Gas central pressure
+vx = 142.3 # Random velocity for all particles
+vy = -31.4
+fileOutputName = "square.hdf5"
+#---------------------------------------------------
+
+vol = 1.
+
+numPart_out = L * L
+numPart_in = L * L * rho0 / rho1 / 4
+
+L_out = int(sqrt(numPart_out))
+L_in = int(sqrt(numPart_in))
+
+pos_out = zeros((numPart_out, 3))
+for i in range(L_out):
+ for j in range(L_out):
+ index = i * L_out + j
+ pos_out[index, 0] = i / (float(L_out)) + 1./(2. * L_out)
+ pos_out[index, 1] = j / (float(L_out)) + 1./(2. * L_out)
+h_out = ones(numPart_out) * (1. / L_out) * 1.2348
+m_out = ones(numPart_out) * vol * rho1 / numPart_out
+u_out = ones(numPart_out) * P1 / (rho1 * (gamma - 1.))
+
+pos_in = zeros((numPart_in, 3))
+for i in range(L_in):
+ for j in range(L_in):
+ index = i * L_in + j
+ pos_in[index, 0] = 0.25 + i / float(2. * L_in) + 1./(2. * 2. * L_in)
+ pos_in[index, 1] = 0.25 + j / float(2. * L_in) + 1./(2. * 2. * L_in)
+h_in = ones(numPart_in) * (1. / L_in) * 1.2348
+m_in = ones(numPart_in) * 0.25 * vol * rho0 / numPart_in
+u_in = ones(numPart_in) * P0 / (rho0 * (gamma - 1.))
+
+# Remove the central particles
+select_out = logical_or(logical_or(pos_out[:,0] < 0.25 , pos_out[:,0] > 0.75), logical_or(pos_out[:,1] < 0.25, pos_out[:,1] > 0.75))
+pos_out = pos_out[select_out, :]
+h_out = h_out[select_out]
+u_out = u_out[select_out]
+m_out = m_out[select_out]
+
+# Add the central region
+pos = append(pos_out, pos_in, axis=0)
+h = append(h_out, h_in, axis=0)
+u = append(u_out, u_in)
+m = append(m_out, m_in)
+numPart = size(h)
+ids = linspace(1, numPart, numPart)
+vel = zeros((numPart, 3))
+vel[:,0] = vx
+vel[:,1] = vy
+
+
+#File
+fileOutput = h5py.File(fileOutputName, 'w')
+
+# Header
+grp = fileOutput.create_group("/Header")
+grp.attrs["BoxSize"] = [vol, vol, 0.2]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+
+#Runtime parameters
+grp = fileOutput.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#Units
+grp = fileOutput.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+grp = fileOutput.create_group("/PartType0")
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = pos
+ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = vel
+ds = grp.create_dataset('Masses', (numPart, 1), 'f')
+ds[()] = m.reshape((numPart,1))
+ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+ds[()] = h.reshape((numPart,1))
+ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
+ds[()] = u.reshape((numPart,1))
+ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
+ds[()] = ids.reshape((numPart,1))
+
+fileOutput.close()
+
+
diff --git a/examples/SquareTest_2D/plotSolution.py b/examples/SquareTest_2D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..b9efe76de1e6c5993fa5333be76a13ba95bdab0f
--- /dev/null
+++ b/examples/SquareTest_2D/plotSolution.py
@@ -0,0 +1,186 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the square test
+
+# Parameters
+gas_gamma = 5./3. # Gas adiabatic index
+gamma = 5./3. # Gas adiabatic index
+rho0 = 4 # Gas central density
+rho1 = 1 # Gas outskirt density
+P0 = 2.5 # Gas central pressure
+P1 = 2.5 # Gas central pressure
+vx = 142.3 # Random velocity for all particles
+vy = -31.4
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+snap = int(sys.argv[1])
+
+# Read the simulation data
+sim = h5py.File("square_%03d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+# Analytical soltion
+centre_x = 0.5 + time * vx
+centre_y = 0.5 + time * vy
+
+while centre_x > 1.:
+ centre_x -= 1.
+while centre_x < 0.:
+ centre_x += 1.
+while centre_y > 1.:
+ centre_y -= 1.
+while centre_y < 0.:
+ centre_y += 1.
+
+pos = sim["/PartType0/Coordinates"][:,:]
+vel = sim["/PartType0/Velocities"][:,:]
+v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
+rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+x = pos[:,0] - centre_x
+y = pos[:,1] - centre_y
+
+# Box wrapping
+x[x>0.5] -= 1.
+x[x<-0.5] += 1.
+y[y>0.5] -= 1.
+y[y<-0.5] += 1.
+
+# Azimuthal velocity profile -----------------------------
+subplot(231)
+scatter(x, y, c=v_norm, cmap="PuBu", edgecolors='face', s=4, vmin=-1., vmax=1.)
+text(0.47, 0.47, "${\\rm{Velocity~norm}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+# Radial density profile --------------------------------
+subplot(232)
+scatter(x, y, c=rho, cmap="PuBu", edgecolors='face', s=4, vmin=0., vmax=4.)
+text(0.47, 0.47, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+# Radial pressure profile --------------------------------
+subplot(233)
+scatter(x, y, c=P, cmap="PuBu", edgecolors='face', s=4, vmin=2, vmax=4)
+text(0.47, 0.47, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+# Internal energy profile --------------------------------
+subplot(234)
+scatter(x, y, c=u, cmap="PuBu", edgecolors='face', s=4, vmin=0.5, vmax=4.)
+text(0.47, 0.47, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+# Radial entropy profile --------------------------------
+subplot(235)
+scatter(x, y, c=S, cmap="PuBu", edgecolors='face', s=4, vmin=0., vmax=3.)
+text(0.47, 0.47, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Square test with $\\gamma=%.3f$ at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Centre:~~~ $(P, \\rho) = (%.3f, %.3f)$"%(P0, rho0), fontsize=10)
+text(-0.49, 0.7, "Outskirts: $(P, \\rho) = (%.3f, %.3f)$"%(P1, rho1), fontsize=10)
+plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.4, scheme, fontsize=10)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+savefig("SquareTest.png", dpi=200)
diff --git a/examples/SquareTest_2D/run.sh b/examples/SquareTest_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..242f1b8b729979b026cfe31002e84cd9ef741129
--- /dev/null
+++ b/examples/SquareTest_2D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e square.hdf5 ]
+then
+ echo "Generating initial conditions for the square test ..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 square.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 40
diff --git a/examples/SquareTest_2D/square.yml b/examples/SquareTest_2D/square.yml
new file mode 100644
index 0000000000000000000000000000000000000000..4f39c6490899cfaafeda17fb0c28281cbadcbbea
--- /dev/null
+++ b/examples/SquareTest_2D/square.yml
@@ -0,0 +1,35 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 4. # The end time of the simulation (in internal units).
+ dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: square # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-1 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ max_smoothing_length: 0.02 # Maximal smoothing length allowed (in internal units).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./square.hdf5 # The file to read
diff --git a/examples/UniformBox_2D/run.sh b/examples/UniformBox_2D/run.sh
index fbabedb1cd564e0831228591ce32d9a73fe08d9a..ee3ef109968a65e2437ea17b42013266195d3314 100755
--- a/examples/UniformBox_2D/run.sh
+++ b/examples/UniformBox_2D/run.sh
@@ -7,4 +7,4 @@ then
python makeIC.py 100
fi
-../swift -s -t 16 uniformPlane.yml
+../swift -s -t 16 uniformPlane.yml 2>&1 | tee output.log
diff --git a/examples/UniformBox_3D/run.sh b/examples/UniformBox_3D/run.sh
index 0cb0a505915be47bafb99ed7531685bfeb3dc829..08891cdd08fccf8f43089951e94dddb33e162030 100755
--- a/examples/UniformBox_3D/run.sh
+++ b/examples/UniformBox_3D/run.sh
@@ -7,4 +7,4 @@ then
python makeIC.py 100
fi
-../swift -s -t 16 uniformBox.yml
+../swift -s -t 16 uniformBox.yml 2>&1 | tee output.log
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 755d4be527e4a81b2b4d2b3b829ea16a74ccc7c5..f9b67c96331e7572b0200093f0c32eee5d2391cd 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -17,6 +17,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_DEFAULT_GRAVITY_H
+#define SWIFT_DEFAULT_GRAVITY_H
#include <float.h>
#include "potentials.h"
@@ -151,3 +153,5 @@ __attribute__((always_inline)) INLINE static void external_gravity(
*/
__attribute__((always_inline)) INLINE static void gravity_kick_extra(
struct gpart* gp, float dt, float half_dt) {}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_H */
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index 7cf375a1fdf7bccc4131dc415ab2d4acbbf2d3bc..c284f543b3be06297600c010e302423eb683adc9 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_DEFAULT_GRAVITY_DEBUG_H
+#define SWIFT_DEFAULT_GRAVITY_DEBUG_H
__attribute__((always_inline)) INLINE static void gravity_debug_particle(
const struct gpart* p) {
@@ -27,3 +29,5 @@ __attribute__((always_inline)) INLINE static void gravity_debug_particle(
p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin,
p->ti_end);
}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_DEBUG_H */
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index aa285ce3245aeac608aa022924a613a7a7b00375..7c7c5b9d244534ef3f6b5f509062019bfcd5f9fb 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -17,8 +17,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_RUNNER_IACT_GRAV_H
-#define SWIFT_RUNNER_IACT_GRAV_H
+#ifndef SWIFT_DEFAULT_GRAVITY_IACT_H
+#define SWIFT_DEFAULT_GRAVITY_IACT_H
/* Includes. */
#include "const.h"
@@ -188,4 +188,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
#endif
}
-#endif /* SWIFT_RUNNER_IACT_GRAV_H */
+#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index e0fd15a91b6e20048f6844b03ea1dde40114d7ec..26eed8a0a680e13130dd257d8b31e6cd00392f6d 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_DEFAULT_GRAVITY_IO_H
+#define SWIFT_DEFAULT_GRAVITY_IO_H
#include "io_properties.h"
@@ -68,3 +70,5 @@ void darkmatter_write_particles(struct gpart* gparts, struct io_props* list,
list[4] = io_make_output_field("Acceleration", FLOAT, 3,
UNIT_CONV_ACCELERATION, gparts, a_grav);
}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_IO_H */
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 77e6543429993bea355c04dd28fcb5b04eee5519..1850ff0a1644d3593f78f150646eae8b2f074e1e 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -16,6 +16,9 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_DEFAULT_GRAVITY_PART_H
+#define SWIFT_DEFAULT_GRAVITY_PART_H
+
/* Some standard headers. */
#include <stdlib.h>
@@ -51,3 +54,5 @@ struct gpart {
long long id_or_neg_offset;
} __attribute__((aligned(gpart_align)));
+
+#endif /* SWIFT_DEFAULT_GRAVITY_PART_H */
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 51e09a5d943a7f3782e484289faddece2567f15d..021599cd2daf61ff35e5f29e3f13b2ad61c8947a 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_DEFAULT_HYDRO_H
+#define SWIFT_DEFAULT_HYDRO_H
#include "adiabatic_index.h"
#include "approx_math.h"
@@ -341,3 +343,5 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {}
+
+#endif /* SWIFT_DEFAULT_HYDRO_H */
diff --git a/src/hydro/Default/hydro_debug.h b/src/hydro/Default/hydro_debug.h
index 9b8136e8349398e4e9e804459a5de23d21043564..d02d3ef82c1b3d751731f49850c06df4b146b164 100644
--- a/src/hydro/Default/hydro_debug.h
+++ b/src/hydro/Default/hydro_debug.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_DEFAULT_HYDRO_DEBUG_H
+#define SWIFT_DEFAULT_HYDRO_DEBUG_H
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
const struct part* p, const struct xpart* xp) {
@@ -29,3 +31,5 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
p->ti_end);
}
+
+#endif /* SWIFT_DEFAULT_HYDRO_DEBUG_H */
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 4975cac35aaa3d4b4b9f99319dfd061998014453..51fa7d07229f86918ef2d7019a9708110cef02e3 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -17,8 +17,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_RUNNER_IACT_H
-#define SWIFT_RUNNER_IACT_H
+#ifndef SWIFT_DEFAULT_HYDRO_IACT_H
+#define SWIFT_DEFAULT_HYDRO_IACT_H
#include "adiabatic_index.h"
@@ -932,4 +932,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
#endif
}
-#endif /* SWIFT_RUNNER_IACT_H */
+#endif /* SWIFT_DEFAULT_HYDRO_IACT_H */
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index 1c3e70753d417eb59d0fda5a7acfd1c24ab93045..bb35c914bcab8787f609b4dd49acd0cc883b4263 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_DEFAULT_HYDRO_IO_H
+#define SWIFT_DEFAULT_HYDRO_IO_H
#include "io_properties.h"
#include "kernel_hydro.h"
@@ -112,3 +114,5 @@ void writeSPHflavour(hid_t h_grpsph) {
* @return 1 if entropy is in 'internal energy', 0 otherwise.
*/
int writeEntropyFlag() { return 0; }
+
+#endif /* SWIFT_DEFAULT_HYDRO_IO_H */
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index a3972dc7ebd1b04a18ccdb6f815d24e36aff270c..a2f4453dc69ed06ca4f315b6be29844c177d0435 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_DEFAULT_HYDRO_PART_H
+#define SWIFT_DEFAULT_HYDRO_PART_H
/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
@@ -117,3 +119,5 @@ struct part {
struct gpart* gpart;
} __attribute__((aligned(part_align)));
+
+#endif /* SWIFT_DEFAULT_HYDRO_PART_H */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 60f07c2fc31fd5f38d2929679d0e13beb1cc9131..e9d626cb8c147c0cf4fa8d27f8bab31d2471beae 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_GADGET2_HYDRO_H
+#define SWIFT_GADGET2_HYDRO_H
#include "adiabatic_index.h"
#include "dimension.h"
@@ -320,7 +322,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
p->force.h_dt *= p->h * hydro_dimension_inv;
- p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
+ p->entropy_dt *=
+ 0.5f * hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
}
/**
@@ -360,3 +363,5 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
/* We read u in the entropy field. We now get S from u */
p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
}
+
+#endif /* SWIFT_GADGET2_HYDRO_H */
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 0280ae7d7a0cdac4567b8e9d4c3a50faf20e93f6..7c8a6eaba96929b01f1901393d7d0498d58badf4 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_GADGET2_HYDRO_DEBUG_H
+#define SWIFT_GADGET2_HYDRO_DEBUG_H
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
const struct part* p, const struct xpart* xp) {
@@ -34,3 +36,5 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);
}
+
+#endif /* SWIFT_GADGET2_HYDRO_DEBUG_H */
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index d37ac491fe5d8ecdf127c217ca025080daf4bbfd..b9807b6332e08012d47ad63652377f4fe5337bf9 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -17,8 +17,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_RUNNER_IACT_LEGACY_H
-#define SWIFT_RUNNER_IACT_LEGACY_H
+#ifndef SWIFT_GADGET2_HYDRO_IACT_H
+#define SWIFT_GADGET2_HYDRO_IACT_H
/**
* @brief SPH interaction functions following the Gadget-2 version of SPH.
@@ -469,8 +469,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig);
/* Change in entropy */
- pi->entropy_dt += 0.5f * mj * visc_term * dvdr;
- pj->entropy_dt += 0.5f * mi * visc_term * dvdr;
+ pi->entropy_dt += mj * visc_term * dvdr;
+ pj->entropy_dt += mi * visc_term * dvdr;
}
/**
@@ -631,7 +631,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pjh_dt.v = mi.v * dvdr.v * ri.v / pirho.v * wj_dr.v;
/* Change in entropy */
- entropy_dt.v = vec_set1(0.5f) * visc_term.v * dvdr.v;
+ entropy_dt.v = visc_term.v * dvdr.v;
/* Store the forces back on the particles. */
for (k = 0; k < VEC_SIZE; k++) {
@@ -644,7 +644,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pi[k]->force.v_sig = fmaxf(pi[k]->force.v_sig, v_sig.f[k]);
pj[k]->force.v_sig = fmaxf(pj[k]->force.v_sig, v_sig.f[k]);
pi[k]->entropy_dt += entropy_dt.f[k] * mj.f[k];
- pj[k]->entropy_dt -= entropy_dt.f[k] * mi.f[k];
+ pj[k]->entropy_dt += entropy_dt.f[k] * mi.f[k];
}
#else
@@ -738,7 +738,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
/* Change in entropy */
- pi->entropy_dt += 0.5f * mj * visc_term * dvdr;
+ pi->entropy_dt += mj * visc_term * dvdr;
}
/**
@@ -894,7 +894,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
pih_dt.v = mj.v * dvdr.v * ri.v / pjrho.v * wi_dr.v;
/* Change in entropy */
- entropy_dt.v = vec_set1(0.5f) * mj.v * visc_term.v * dvdr.v;
+ entropy_dt.v = mj.v * visc_term.v * dvdr.v;
/* Store the forces back on the particles. */
for (k = 0; k < VEC_SIZE; k++) {
@@ -913,4 +913,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
#endif
}
-#endif /* SWIFT_RUNNER_IACT_LEGACY_H */
+#endif /* SWIFT_GADGET2_HYDRO_IACT_H */
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index b89b4f70aea1c29e2df4df3615e34552e4ea3f8d..433aef64c388c8bc4989e883f10a8f0d3eeb30e9 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_GADGET2_HYDRO_IO_H
+#define SWIFT_GADGET2_HYDRO_IO_H
#include "adiabatic_index.h"
#include "hydro.h"
@@ -107,9 +109,10 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
- "(No treatment) Legacy Gadget-2 as in Springel (2005)");
- writeAttribute_s(h_grpsph, "Viscosity Model",
- "Legacy Gadget-2 as in Springel (2005)");
+ "(No treatment) as in Springel (2005)");
+ writeAttribute_s(
+ h_grpsph, "Viscosity Model",
+ "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
}
@@ -120,3 +123,5 @@ void writeSPHflavour(hid_t h_grpsph) {
* @return 1 if entropy is in 'internal energy', 0 otherwise.
*/
int writeEntropyFlag() { return 0; }
+
+#endif /* SWIFT_GADGET2_HYDRO_IO_H */
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 0cd830c51a1836c0a56e7b4097990bf4d07a101f..484792438d2717413c1ca8d4f429eac2e6d21b20 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -16,6 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_GADGET2_HYDRO_PART_H
+#define SWIFT_GADGET2_HYDRO_PART_H
/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
@@ -109,3 +111,5 @@ struct part {
struct gpart* gpart;
} __attribute__((aligned(part_align)));
+
+#endif /* SWIFT_GADGET2_HYDRO_PART_H */
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 842a028b9031b3c1cf36e4735089d05e038abdca..a5d73aad02372aa22b840c2cb0d3100cd439e75d 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -16,10 +16,29 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_MINIMAL_HYDRO_H
+#define SWIFT_MINIMAL_HYDRO_H
+
+/**
+ * @file Minimal/hydro.h
+ * @brief Minimal conservative implementation of SPH (Non-neighbour loop
+ * equations)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
#include "adiabatic_index.h"
-#include "approx_math.h"
+#include "dimension.h"
#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "kernel_hydro.h"
/**
* @brief Returns the internal energy of a particle
@@ -34,7 +53,7 @@
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part *restrict p, float dt) {
- return p->u;
+ return p->u + p->u_dt * dt;
}
/**
@@ -46,7 +65,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
- return gas_pressure_from_internal_energy(p->rho, p->u);
+ const float u = p->u + p->u_dt * dt;
+
+ return gas_pressure_from_internal_energy(p->rho, u);
}
/**
@@ -62,7 +83,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part *restrict p, float dt) {
- return gas_entropy_from_internal_energy(p->rho, p->u);
+ const float u = p->u + p->u_dt * dt;
+
+ return gas_entropy_from_internal_energy(p->rho, u);
}
/**
@@ -74,7 +97,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) {
- return gas_soundspeed_from_internal_energy(p->rho, p->u);
+ const float u = p->u + p->u_dt * dt;
+
+ return gas_soundspeed_from_internal_energy(p->rho, u);
}
/**
@@ -150,7 +175,6 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
- xp->u_full = p->u;
}
/**
@@ -164,6 +188,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
*/
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p) {
+
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
@@ -227,7 +252,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) {
- p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
+ const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
+ const float pressure = hydro_get_pressure(p, half_dt);
+
+ p->force.pressure = pressure;
}
/**
@@ -268,10 +296,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part *restrict p, const struct xpart *restrict xp, int t0, int t1,
double timeBase) {
- p->u = xp->u_full;
-
- /* Need to recompute the pressure as well */
- p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
+ /* Drift the pressure */
+ const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
+ p->force.pressure = hydro_get_pressure(p, dt_entr);
}
/**
@@ -304,11 +331,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt,
float half_dt) {
- /* Kick in momentum space */
- xp->u_full += p->u_dt * dt;
+ /* Do not decrease the energy by more than a factor of 2*/
+ const float u_change = p->u_dt * dt;
+ if (u_change > -0.5f * p->u)
+ p->u += u_change;
+ else
+ p->u *= 0.5f;
- /* Get the predicted internal energy */
- p->u = xp->u_full - half_dt * p->u_dt;
+ /* Do not 'overcool' when timestep increases */
+ if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt;
}
/**
@@ -323,3 +354,5 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {}
+
+#endif /* SWIFT_MINIMAL_HYDRO_H */
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 7ab348d255173a33fc29e482eb6032c03c82d9eb..0a6e50da6051fd961f4bf35f32caee311888a13e 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -16,18 +16,36 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_MINIMAL_HYDRO_DEBUG_H
+#define SWIFT_MINIMAL_HYDRO_DEBUG_H
+
+/**
+ * @file Minimal/hydro_debug.h
+ * @brief Minimal conservative implementation of SPH (Debugging routines)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
const struct part* p, const struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
- "u_full=%.3e, u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
- "h=%.3e, dh/dt=%.3e "
- "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%d, t_end=%d\n",
+ "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+ "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, "
+ "t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
- xp->u_full, p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h,
- p->force.h_dt, (int)p->density.wcount, p->mass, p->rho_dh, p->rho,
- p->ti_begin, p->ti_end);
+ p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt,
+ (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
+ p->ti_end);
}
+
+#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 9eda45760248ebc0dde18e503abe0b1288659518..a8a855d9db81f6927c1d8b45410a57d50a8366de 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -16,18 +16,25 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_RUNNER_IACT_MINIMAL_H
-#define SWIFT_RUNNER_IACT_MINIMAL_H
-
-#include "adiabatic_index.h"
+#ifndef SWIFT_MINIMAL_HYDRO_IACT_H
+#define SWIFT_MINIMAL_HYDRO_IACT_H
/**
- * @brief Minimal conservative implementation of SPH
+ * @file Minimal/hydro_iact.h
+ * @brief Minimal conservative implementation of SPH (Neighbour loop equations)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
*
- * The thermal variable is the internal energy (u). No viscosity nor
- * thermal conduction terms are implemented.
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
*/
+#include "adiabatic_index.h"
+
/**
* @brief Density loop
*/
@@ -115,6 +122,8 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
__attribute__((always_inline)) INLINE static void runner_iact_force(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+ const float fac_mu = 1.f; /* Will change with cosmological integration */
+
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
@@ -143,34 +152,60 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
- const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
- const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
+ const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2];
- /* Compute sound speeds */
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = fminf(dvdr, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Compute sound speeds and signal velocity */
const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
- const float v_sig = ci + cj;
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+
+ /* Convolve with the kernel */
+ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* SPH acceleration term */
- const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+ const float sph_acc_term =
+ (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+ /* Assemble the acceleration */
+ const float acc = sph_acc_term + visc_acc_term;
/* Use the force Luke ! */
- pi->a_hydro[0] -= mj * sph_term * dx[0];
- pi->a_hydro[1] -= mj * sph_term * dx[1];
- pi->a_hydro[2] -= mj * sph_term * dx[2];
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
- pj->a_hydro[0] += mi * sph_term * dx[0];
- pj->a_hydro[1] += mi * sph_term * dx[1];
- pj->a_hydro[2] += mi * sph_term * dx[2];
+ pj->a_hydro[0] += mi * acc * dx[0];
+ pj->a_hydro[1] += mi * acc * dx[1];
+ pj->a_hydro[2] += mi * acc * dx[2];
/* Get the time derivative for u. */
- pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
- pj->u_dt += P_over_rho_j * mi * dvdr * r_inv * wj_dr;
+ const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+ const float sph_du_term_j = P_over_rho2_j * dvdr * r_inv * wj_dr;
+
+ /* Viscosity term */
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+
+ /* Assemble the energy equation term */
+ const float du_dt_i = sph_du_term_i + visc_du_term;
+ const float du_dt_j = sph_du_term_j + visc_du_term;
+
+ /* Internal energy time derivatibe */
+ pi->u_dt += du_dt_i * mj;
+ pj->u_dt += du_dt_j * mi;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
@@ -198,6 +233,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+ const float fac_mu = 1.f; /* Will change with cosmological integration */
+
const float r = sqrtf(r2);
const float r_inv = 1.0f / r;
@@ -226,29 +263,53 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
- const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
- const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
+ const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
(pi->v[1] - pj->v[1]) * dx[1] +
(pi->v[2] - pj->v[2]) * dx[2];
- /* Compute sound speeds */
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = fminf(dvdr, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Compute sound speeds and signal velocity */
const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
- const float v_sig = ci + cj;
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+
+ /* Convolve with the kernel */
+ const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* SPH acceleration term */
- const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+ const float sph_acc_term =
+ (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+ /* Assemble the acceleration */
+ const float acc = sph_acc_term + visc_acc_term;
/* Use the force Luke ! */
- pi->a_hydro[0] -= mj * sph_term * dx[0];
- pi->a_hydro[1] -= mj * sph_term * dx[1];
- pi->a_hydro[2] -= mj * sph_term * dx[2];
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
/* Get the time derivative for u. */
- pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
+ const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+
+ /* Viscosity term */
+ const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+
+ /* Assemble the energy equation term */
+ const float du_dt_i = sph_du_term_i + visc_du_term;
+
+ /* Internal energy time derivatibe */
+ pi->u_dt += du_dt_i * mj;
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
@@ -268,4 +329,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
"function does not exist yet!");
}
-#endif /* SWIFT_RUNNER_IACT_MINIMAL_H */
+#endif /* SWIFT_MINIMAL_HYDRO_IACT_H */
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 1a5b2938354fb8ee3ef936f86ff442c731e0d4c7..01a75b17fd5577cfcfb48d3afac22579f30fcf7a 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -16,7 +16,25 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_MINIMAL_HYDRO_IO_H
+#define SWIFT_MINIMAL_HYDRO_IO_H
+/**
+ * @file Minimal/hydro_io.h
+ * @brief Minimal conservative implementation of SPH (i/o routines)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
#include "io_properties.h"
#include "kernel_hydro.h"
@@ -51,6 +69,16 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
+float convert_S(struct engine* e, struct part* p) {
+
+ return hydro_get_entropy(p, 0);
+}
+
+float convert_P(struct engine* e, struct part* p) {
+
+ return hydro_get_pressure(p, 0);
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -61,7 +89,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 8;
+ *num_fields = 10;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -80,6 +108,11 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_ACCELERATION, parts, a_hydro);
list[7] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+ list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+ UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+ parts, rho, convert_S);
+ list[9] = io_make_output_field_convert_part(
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
}
/**
@@ -90,8 +123,9 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */
/* Nothing in this minimal model... */
- writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No model");
- writeAttribute_s(h_grpsph, "Viscosity Model", "No model");
+ writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
+ writeAttribute_s(h_grpsph, "Viscosity Model",
+ "Minimal treatment as in Monaghan (1992)");
/* Time integration properties */
writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
@@ -104,3 +138,5 @@ void writeSPHflavour(hid_t h_grpsph) {
* @return 1 if entropy is in 'internal energy', 0 otherwise.
*/
int writeEntropyFlag() { return 0; }
+
+#endif /* SWIFT_MINIMAL_HYDRO_IO_H */
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index f4feb2311ad2bfd8c7227f3b072c5971724a5815..ad65f8b44fc67f4aae6470246cbab91bc3710007 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -16,6 +16,22 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_MINIMAL_HYDRO_PART_H
+#define SWIFT_MINIMAL_HYDRO_PART_H
+
+/**
+ * @file Minimal/hydro_part.h
+ * @brief Minimal conservative implementation of SPH (Particle definition)
+ *
+ * The thermal variable is the internal energy (u). Simple constant
+ * viscosity term without switches is implemented. No thermal conduction
+ * term is implemented.
+ *
+ * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
+ * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
+ * pp. 759-794.
+ */
/**
* @brief Particle fields not needed during the SPH loops over neighbours.
@@ -31,8 +47,6 @@ struct xpart {
float v_full[3]; /*!< Velocity at the last full step. */
- float u_full; /*!< Thermal energy at the last full step. */
-
} __attribute__((aligned(xpart_align)));
/**
@@ -107,3 +121,5 @@ struct part {
struct gpart* gpart; /*!< Pointer to corresponding gravity part. */
} __attribute__((aligned(part_align)));
+
+#endif /* SWIFT_MINIMAL_HYDRO_PART_H */
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index c425ecce6de38b32645b367886e039bb950df68d..0e4eaf150c764e27156747bb01db20a71f03c7b5 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -70,7 +70,7 @@ void hydro_props_print(const struct hydro_props *p) {
message(
"Hydrodynamic integration: Max change of volume: %.2f "
"(max|dlog(h)/dt|=%f).",
- powf(expf(p->log_max_h_change), hydro_dimension), p->log_max_h_change);
+ pow_dimension(expf(p->log_max_h_change)), p->log_max_h_change);
if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
message("Maximal iterations in ghost task set to %d (default is %d)",
@@ -90,8 +90,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
writeAttribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
writeAttribute_f(h_grpsph, "Volume log(max(delta h))", p->log_max_h_change);
writeAttribute_f(h_grpsph, "Volume max change time-step",
- powf(expf(p->log_max_h_change), 3.f));
- writeAttribute_f(h_grpsph, "Max ghost iterations",
+ pow_dimension(expf(p->log_max_h_change)));
+ writeAttribute_i(h_grpsph, "Max ghost iterations",
p->max_smoothing_iterations);
}
#endif
diff --git a/tests/difffloat.py b/tests/difffloat.py
index 57707c5920997e3ef688606a0839f59a69d2e4f2..e0f0864372264899c6de1bf2f83ab678b7dd9ead 100644
--- a/tests/difffloat.py
+++ b/tests/difffloat.py
@@ -89,13 +89,12 @@ for i in range(n_lines_to_check):
abs_diff = abs(data1[i,j] - data2[i,j])
sum = abs(data1[i,j] + data2[i,j])
- if abs(data1[i,j]) + abs(data2[i,j]) < 2.5e-7: continue
if sum > 0:
rel_diff = abs(data1[i,j] - data2[i,j]) / sum
else:
rel_diff = 0.
- if( abs_diff > absTol[j]):
+ if( abs_diff > 1.1*absTol[j]):
print "Absolute difference larger than tolerance (%e) for particle %d, column %d:"%(absTol[j], i,j)
print "%10s: a = %e"%("File 1", data1[i,j])
print "%10s: b = %e"%("File 2", data2[i,j])
@@ -103,7 +102,9 @@ for i in range(n_lines_to_check):
print ""
error = True
- if( rel_diff > relTol[j]):
+ if abs(data1[i,j]) < 1e-6 and + abs(data2[i,j]) < 1e-6 : continue
+
+ if( rel_diff > 1.1*relTol[j]):
print "Relative difference larger than tolerance (%e) for particle %d, column %d:"%(relTol[j], i,j)
print "%10s: a = %e"%("File 1", data1[i,j])
print "%10s: b = %e"%("File 2", data2[i,j])
diff --git a/tests/test125cells.c b/tests/test125cells.c
index cbb3ec87c735164bb64a67d51e1f891a9030d50b..37d530fa7b3750a3304cb11868ac43e6df62781d 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -472,8 +472,8 @@ int main(int argc, char *argv[]) {
prog_const.const_newton_G = 1.f;
struct hydro_props hp;
- hp.target_neighbours = h * h * h * kernel_norm;
- hp.delta_neighbours = 1.;
+ hp.target_neighbours = pow_dimension(h) * kernel_norm;
+ hp.delta_neighbours = 2.;
hp.max_smoothing_iterations = 1;
hp.CFL_condition = 0.1;