Add raw results and all simulation scripts

A/plotPressure.sh

+#! /bin/bash
+gnuplot -p << EOF
+set grid
+set xlabel 'Time [s]' font ",12"
+set ylabel 'Pressure [Pa]' font ",12" offset -3
+set key center bottom
+set encoding utf8
+set xtics font ",12" 
+set ytics font ",12" 
+set key font ",14"
+set format y "%.1t*10^%T";
+set lmargin 15
+
+#set xrange [6.8:8.2]
+set yrange [-350000:350000]
+#set yrange [220000:320000]
+
+
+plot "results/watchpointI.txt" using 1:2 with lines dt 1 lc rgb "grey" lw 2 title "Crank-Nicolson θ = 0.5",\
+     "results/watchpointII.txt" using 1:2 with lines dt 3 lc rgb "black" lw 2 title "Crank-Nicolson θ = 0.55",\
+     "results/watchpointIII.txt" using 1:2 with lines dt 2 lc rgb "dark-grey" lw 2 title "Backward Euler"
+     
+#plot "results/watchpointI.txt" using 1:2 with linespoints pi 50 dt 1 pt 6 lt -1 title "Crank-Nicolson θ = 0.5",\
+#     "results/watchpointII.txt" using 1:2 with linespoints pi 50 dt 2 pt 6 lt -1title "Crank-Nicolson θ = 0.5",\
+#     "results/watchpointIII.txt" using 1:2 with linespoints pi 50 dt 3 pt 6 lt -1 title "Backward Euler"     
+     
+EOF
+

A/plotVelocity.sh

+#! /bin/bash
+gnuplot -p << EOF
+set grid
+set xlabel 'Time [s]'
+set ylabel 'Velocity'
+plot "watchpoint.txt" using 1:3 with lines title ""
+EOF
+

A/results/README.md

+I: normal CN, theta=0.5
+II: theta=0.55
+III: theta=1.0 -> BE

A/sonicLiquid.py

+#! /usr/bin/python3
+
+from nutils import *
+import numpy, itertools, treelog
+
+@function.replace
+def subs0(f):
+  if isinstance(f, function.Argument) and f._name == 'lhs':
+    return function.Argument(name='lhs0', shape=f.shape, nderiv=f._nderiv)
+
+def main(nelems=200, timestep=.005, refdensity=1e3, refpressure=101325., psi=1e-6, viscosity=1e-3, theta=1.0):
+
+  domain, geom = mesh.rectilinear([numpy.linspace(0, 1000, nelems+1)])
+  bezier = domain.sample('bezier', 2)
+
+  ns = function.Namespace()
+  ns._functions['t'] = lambda f: theta * f + (1-theta) * subs0(f)
+  ns._functions_nargs['t'] = 1
+  ns._functions['δt'] = lambda f: (f - subs0(f)) / timestep
+  ns._functions_nargs['δt'] = 1
+  ns.x = geom
+  ns.ρref = refdensity
+  ns.pref = refpressure
+  ns.pin = 98100
+  ns.μ = viscosity
+  ns.ψ = psi
+  ns.ρbasis, ns.ubasis = function.chain([domain.basis('std', degree=1), domain.basis('std', degree=2).vector(domain.ndims)])
+  ns.ρ = 'ρref + ρbasis_n ?lhs_n' # density is shifted by ref density
+  ns.u_i = 'ubasis_ni ?lhs_n'
+  ns.p = 'pref + (ρ - ρref) / ψ' # pressure-density connection, see equ (8)
+  ns.σ_ij = 'μ (u_i,j + u_j,i) - p δ_ij' # diffusive term and pressure gradient
+  ns.h = 1 / nelems
+  ns.k = 'ρ h / μ' # needs work, stabilization coeff
+
+  res = domain.integral('ρbasis_n (δt(ρ) + t((ρ u_k)_,k)) d:x' @ ns, degree=4) # mass balance
+  res += domain.integral('(ubasis_ni (δt(ρ u_i) + t((ρ u_i u_j)_,j)) + ubasis_ni,j t(σ_ij)) d:x' @ ns, degree=4) # momentum balance
+  res += domain.boundary['left'].integral('pin ubasis_ni n_i d:x' @ ns, degree=4) # pressure set at inlet
+  #res += domain.integral('k ubasis_ni,j (u_j / sqrt(u_k u_k)) (δρu_i / δt + (ρ u_i u_j - σ_ij)_,j) d:x' @ ns, degree=4) # SUPG stabilization
+
+  lhs = numpy.zeros(res.shape)
+  t = 0
+  
+  f = open("watchpoint.txt", "w")
+  
+  sqr = domain.boundary['right'].integral('(u_0 - 1)^2' @ ns, degree=4)
+  cons0 = solver.optimize('lhs', sqr, droptol=1e-14)
+
+  for istep in itertools.count():
+    with log.context('timestep', istep):
+   
+      x, p, ρ, u = bezier.eval(['x_i', 'p', 'ρ', 'u_i'] @ ns, lhs=lhs)
+      f.write("%e; %e; %e\n" % (t, p[399], u[199]))
+      f.flush()
+      
+      cons = min(.2 * t, 1) * cons0
+        
+      lhs = solver.newton('lhs', res, constrain=cons, arguments=dict(lhs0=lhs), lhs0=lhs).solve(1e-8)
+
+      t += timestep
+      
+      if t>=20:
+        break
+  f.close()
+
+
+if __name__ == '__main__':
+  cli.run(main)

B/Allclean

+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+echo "Cleaning..."
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+# Participant 1: Fluid1
+Participant1="Fluid1"
+cd ${Participant1}
+    # Clean the case
+    cleanCase
+    rm -rfv 0
+    # Create an empty .foam file for ParaView
+    # Note: ".foam" triggers the native OpenFOAM reader of ParaView.
+    # Change to ".OpenFOAM" to use the OpenFOAM reader provided with OpenFOAM.
+    touch ${Participant1}.foam
+cd ..
+# Remove the log files
+rm -fv ${Participant1}_blockMesh.log
+rm -fv ${Participant1}_checkMesh.log
+rm -fv ${Participant1}_potentialFoam.log
+rm -fv ${Participant1}_decomposePar.log
+rm -fv ${Participant1}.log
+rm -fv ${Participant1}_reconstructPar.log
+
+echo "Cleaning complete!"
+#------------------------------------------------------------------------------

B/Fluid1/0.orig/U

+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        outOfBounds     clamp;
+        // Time-varying outlet velocity
+        type            uniformFixedValue;
+        uniformValue    table
+                        (
+                            ( 0   ( 0 0 0 ) )
+                            ( 5   ( 0 1 0 ) )
+                            ( 100 ( 0 1 0 ) )
+                        );
+    }
+    xEmpty1
+    {
+       // type            fixedValue;
+       // value           uniform (0 0 0);
+
+     // type            empty;
+            type           symmetryPlane;
+    }
+
+    xEmpty2
+    {
+       // type            fixedValue;
+       // value           uniform (0 0 0);
+
+      // type            empty;
+          type           symmetryPlane;
+    }
+
+    yEmpty1
+    {
+       //type            fixedValue;
+      // value           uniform (0 0 0);
+
+       //type            empty;
+          type           symmetryPlane;
+    }
+
+    yEmpty2
+    {
+      //type            fixedValue;
+      // value           uniform (0 0 0);
+
+      // type            empty;
+          type           symmetryPlane;
+    }
+
+}
+
+
+// ************************************************************************* //

B/Fluid1/0.orig/p

+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 98100;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 98100;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+        value           uniform 98100;
+    }
+
+    xEmpty1
+    {
+       // type            zeroGradient;
+       // type            empty;
+        type           symmetryPlane;
+    }
+
+    xEmpty2
+    {
+      // type            zeroGradient;
+      // type            empty;
+          type           symmetryPlane;
+    }
+
+    yEmpty1
+    {
+         //type            zeroGradient;
+      // type            empty;
+           type           symmetryPlane;
+    }
+
+
+    yEmpty2
+    {
+         // type            zeroGradient;
+       //  type            empty;
+         type           symmetryPlane;
+    }
+
+}
+
+// ************************************************************************* //

B/Fluid1/constant/thermodynamicProperties

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermodynamicProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+rho0            rho0 [1 -3 0 0 0 0 0] 1000;
+
+p0              p0 [1 -1 -2 0 0 0 0] 101325;
+
+psi             psi [0 -2 2 0 0 0 0] 1e-6;
+
+
+// ************************************************************************* //

B/Fluid1/constant/transportProperties

+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+mu              mu  [1 -1 -1 0 0 0 0] 1e-3;
+
+// ************************************************************************* //

B/Fluid1/constant/turbulenceProperties

+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+// ************************************************************************* //

B/Fluid1/system/blockMeshDict

+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant/polyMesh";
+    object      blockMeshDict;
+}
+
+convertToMeters 1;
+
+vertices
+ (
+     ( -0.5  0  -0.5 )  //0
+     ( 0.5  0  -0.5 )   //1
+     ( 0.5 1000 -0.5 )   //2
+     ( -0.5 1000 -0.5 )  //3
+     ( -0.5  0  0.5 )  //4
+     ( 0.5  0  0.5 )    //5
+     ( 0.5 1000 0.5 )   //6
+     ( -0.5 1000 0.5 )  //7
+ );
+
+ blocks
+(
+     hex (0 1 2 3 4 5 6 7) (10 200 10) simpleGrading (1 1 1)  //1
+);
+
+
+  edges
+  (
+  );
+
+boundary
+(
+    outlet
+     {
+         type patch;
+             faces
+             (
+             (3 7 6 2)
+             );
+     }
+
+    inlet
+     {
+        type patch;
+             faces
+             (
+             (0 1 5 4)
+             );
+    }
+
+     xEmpty1
+     {
+     //   type empty;
+        //type wall;
+        type symmetryPlane;
+
+             faces
+             (
+             (4 7 3 0)
+             );
+     }
+
+     xEmpty2
+     {
+      //        type empty;
+        //type wall;
+        type symmetryPlane;
+
+             faces
+             (
+             (1 2 6 5)
+             );
+     }
+     yEmpty1
+     {
+        //  type empty;
+        //type wall;
+        type symmetryPlane;
+             faces
+             (
+             (0 3 2 1)
+             );
+     }
+
+     yEmpty2
+     {
+     // type empty;
+        //type wall;
+        type symmetryPlane;
+             faces
+             (
+             (5 6 7 4)
+             );
+     }
+
+);
+
+ mergePatchPairs
+ (
+ );
+
+ // ************************************************************************* //
+
+
+

B/Fluid1/system/controlDict

+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     sonicLiquidFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         20.0;
+
+deltaT          0.005;
+
+writeControl    adjustableRunTime;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+functions{#includeFunc  probes}
+
+
+/*functions
+{
+    preCICE_Adapter
+    {
+        type preciceAdapterFunctionObject;
+        libs ("libpreciceAdapterFunctionObject.so");
+    }
+}
+*/

B/Fluid1/system/fvSchemes

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    //default         Euler;
+    default         CrankNicolson 0.45;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss limitedLinearV 1;
+    div(phid,p)     Gauss limitedLinear 1;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //

B/Fluid1/system/fvSolution

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "p.*"
+    {
+        solver          PBiCG;
+        preconditioner  DILU
+        //smoother        symGaussSeidel;
+        tolerance       1e-06;
+        relTol          0.01;
+    }
+
+    "U.*"
+    {
+        $p;
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    "rho.*"
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+}
+
+PIMPLE
+{
+    nOuterCorrectors 4;
+    nCorrectors      2;
+    nNonOrthogonalCorrectors 0;
+}
+
+
+// ************************************************************************* //

B/Fluid1/system/probes

+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Web:      www.OpenFOAM.org
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+Description
+    Writes out values of fields from cells nearest to specified locations.
+
+\*---------------------------------------------------------------------------*/
+
+#includeEtc "caseDicts/postProcessing/probes/probes.cfg"
+
+fields (p U);
+probeLocations
+(
+    (0 999 0)
+    (0 500 0)
+);
+
+// ************************************************************************* //

B/Fluid1/system/probesDict

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  5                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;