file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/../resources DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
CreateTestExecAndInstall(CoupledTransportDiffusionEquations_1DHeatedChannel.cxx "${libs_for_tests}" )
+
CreateTestExecAndInstall(DiffusionEquation_1DHeatedRod.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(DiffusionEquation_1DHeatedRod_FE.cxx "${libs_for_tests}" )
+
CreateTestExecAndInstall(DriftModel_1DBoilingAssembly.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(DriftModel_1DBoilingChannel.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(DriftModel_1DChannelGravity.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(DriftModel_2DInclinedChannelGravity.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(DriftModel_2DInclinedChannelGravityBarriers.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(DriftModel_3DCanalCloison.cxx "${libs_for_tests}" )
+
CreateTestExecAndInstall(FiveEqsTwoFluid_1DBoilingChannel.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(FiveEqsTwoFluid_1DDepressurisation.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(FiveEqsTwoFluid_1DRiemannProblem.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(FiveEqsTwoFluid_2DInclinedBoilingChannel.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(FiveEqsTwoFluid_2DInclinedSedimentation.cxx "${libs_for_tests}" )
+
CreateTestExecAndInstall(IsothermalTwoFluid_1DDepressurisation.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(IsothermalTwoFluid_1DRiemannProblem.cxx "${libs_for_tests}" )
#CreateTestExecAndInstall(IsothermalTwoFluid_1DSedimentation.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(IsothermalTwoFluid_2DInclinedSedimentation.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(IsothermalTwoFluid_2DVidangeReservoir.cxx "${libs_for_tests}" )
+
CreateTestExecAndInstall(SinglePhase_1DDepressurisation.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(SinglePhase_1DHeatedChannel.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(SinglePhase_1DPorosityJump.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(SinglePhase_1DRiemannProblem.cxx "${libs_for_tests}" )
+CreateTestExecAndInstall(SinglePhase_1DRiemannProblem_Implicit.cxx "${libs_for_tests}" )
+CreateTestExecAndInstall(SinglePhase_1DRiemannProblem_Implicit_LineSearch.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(SinglePhase_2DHeatDrivenCavity.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(SinglePhase_2DHeatDrivenCavity_unstructured.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(SinglePhase_2DHeatedChannelInclined.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(SinglePhase_2DWallHeatedChannel.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(SinglePhase_3DHeatDrivenCavity.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(SinglePhase_HeatedWire_2Branches.cxx "${libs_for_tests}" )
+
CreateTestExecAndInstall(TransportEquation_1DHeatedChannel.cxx "${libs_for_tests}" )
+
CreateTestExecAndInstall(StationaryDiffusionEquation_2DEF_StructuredTriangles.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(StationaryDiffusionEquation_2DEF_StructuredTriangles_Neumann.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(StationaryDiffusionEquation_2DEF_UnstructuredTriangles.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(StationaryDiffusionEquation_2DFV_StructuredSquares.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(StationaryDiffusionEquation_3DEF_StructuredTetrahedra.cxx "${libs_for_tests}" )
CreateTestExecAndInstall(StationaryDiffusionEquation_3DFV_StructuredTetrahedra.cxx "${libs_for_tests}" )
+
CreateTestExecAndInstall(testEOS.cxx "${libs_for_tests}" )
--- /dev/null
+#include "SinglePhase.hxx"
+
+using namespace std;
+
+int main(int argc, char** argv)
+{
+ //Preprocessing: mesh and group creation
+ cout << "Building a cartesian mesh " << endl;
+ double xinf=0.0;
+ double xsup=1.0;
+ int nx=50;
+ Mesh M(xinf,xsup,nx);
+ double eps=1.E-8;
+ M.setGroupAtPlan(xsup,0,eps,"LeftBoundary");
+ M.setGroupAtPlan(xinf,0,eps,"RightBoundary");
+ int spaceDim = M.getSpaceDimension();
+
+ //initial data
+ double initialVelocity_Left=1;
+ double initialTemperature_Left=565;
+ double initialPressure_Left=155e5;
+ double initialVelocity_Right=1;
+ double initialTemperature_Right=565;
+ double initialPressure_Right=50e5;
+
+ SinglePhase myProblem(Liquid,around155bars600K,spaceDim);
+ // Prepare for the initial condition
+ int nVar = myProblem.getNumberOfVariables();
+ Vector VV_Left(nVar),VV_Right(nVar);
+ // left and right constant vectors
+ VV_Left[0] = initialPressure_Left;
+ VV_Left[1] = initialVelocity_Left;
+ VV_Left[2] = initialTemperature_Left ;
+ VV_Right[0] = initialPressure_Right;
+ VV_Right[1] = initialVelocity_Right;
+ VV_Right[2] = initialTemperature_Right ;
+
+ //Initial field creation
+ double discontinuity = (xinf+xsup)/2.;
+
+ cout << "Building initial data " << endl;
+ Field VV("Primitive", CELLS, M, nVar);
+
+ myProblem.setInitialFieldStepFunction(M,VV_Left,VV_Right,discontinuity);
+
+ //set the boundary conditions
+ myProblem.setNeumannBoundaryCondition("LeftBoundary");
+ myProblem.setNeumannBoundaryCondition("RightBoundary");
+
+ // set the numerical method
+ myProblem.setNumericalScheme(upwind, Implicit);
+
+ // name file save
+ string fileName = "1DRiemannProblem_Implicit";
+
+ // parameters calculation
+ unsigned MaxNbOfTimeStep = 3;
+ int freqSave = 1;
+ double cfl = 0.95;
+ double maxTime = 5;
+ double precision = 1e-6;
+
+ myProblem.setCFL(cfl);
+ myProblem.setPrecision(precision);
+ myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
+ myProblem.setTimeMax(maxTime);
+ myProblem.setFreqSave(freqSave);
+ myProblem.setFileName(fileName);
+ myProblem.saveConservativeField(true);
+ myProblem.setSaveFileFormat(CSV);
+
+ myProblem.setLinearSolver(GMRES, ILU);
+ myProblem.setNewtonSolver(precision,20, Newton_SOLVERLAB);
+ myProblem.usePrimitiveVarsInNewton(true);
+
+ // evolution
+ myProblem.initialize();
+
+ bool ok = myProblem.run();
+ if (ok)
+ cout << "Simulation "<<fileName<<" is successful !" << endl;
+ else
+ cout << "Simulation "<<fileName<<" failed ! " << endl;
+
+ cout << "------------ End of calculation !!! -----------" << endl;
+ myProblem.terminate();
+
+ return EXIT_SUCCESS;
+}
--- /dev/null
+#include "SinglePhase.hxx"
+
+using namespace std;
+
+int main(int argc, char** argv)
+{
+ //Preprocessing: mesh and group creation
+ cout << "Building a cartesian mesh " << endl;
+ double xinf=0.0;
+ double xsup=1.0;
+ int nx=2;
+ Mesh M(xinf,xsup,nx);
+ double eps=1.E-8;
+ M.setGroupAtPlan(xsup,0,eps,"LeftBoundary");
+ M.setGroupAtPlan(xinf,0,eps,"RightBoundary");
+ int spaceDim = M.getSpaceDimension();
+
+ //initial data
+ double initialVelocity_Left=1;
+ double initialTemperature_Left=565;
+ double initialPressure_Left=155e5;
+ double initialVelocity_Right=1;
+ double initialTemperature_Right=565;
+ double initialPressure_Right=50e5;
+
+ SinglePhase myProblem(Liquid,around155bars600K,spaceDim);
+ // Prepare for the initial condition
+ int nVar = myProblem.getNumberOfVariables();
+ Vector VV_Left(nVar),VV_Right(nVar);
+ // left and right constant vectors
+ VV_Left[0] = initialPressure_Left;
+ VV_Left[1] = initialVelocity_Left;
+ VV_Left[2] = initialTemperature_Left ;
+ VV_Right[0] = initialPressure_Right;
+ VV_Right[1] = initialVelocity_Right;
+ VV_Right[2] = initialTemperature_Right ;
+
+ //Initial field creation
+ double discontinuity = (xinf+xsup)/2.;
+
+ cout << "Building initial data " << endl;
+ Field VV("Primitive", CELLS, M, nVar);
+
+ myProblem.setInitialFieldStepFunction(M,VV_Left,VV_Right,discontinuity);
+
+ //set the boundary conditions
+ myProblem.setNeumannBoundaryCondition("LeftBoundary");
+ myProblem.setNeumannBoundaryCondition("RightBoundary");
+
+ // set the numerical method
+ myProblem.setNumericalScheme(upwind, Implicit);
+
+ // name file save
+ string fileName = "1DRiemannProblem_Implicit_LineSearch";
+
+ // parameters calculation
+ unsigned MaxNbOfTimeStep = 3;
+ int freqSave = 1;
+ double cfl = 1;
+ double maxTime = 5;
+ double precision = 1e-6;
+
+ myProblem.setCFL(cfl);
+ myProblem.setPrecision(precision);
+ myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
+ myProblem.setTimeMax(maxTime);
+ myProblem.setFreqSave(freqSave);
+ myProblem.setFileName(fileName);
+ myProblem.saveConservativeField(true);
+ //myProblem.setVerbose(true,false);
+ myProblem.setSaveFileFormat(CSV);
+
+ myProblem.setLinearSolver(GMRES, LU);
+ myProblem.setNewtonSolver(precision,1, Newton_PETSC_LINESEARCH);
+ myProblem.usePrimitiveVarsInNewton(false);
+
+ // evolution
+ myProblem.initialize();
+
+ bool ok = myProblem.run();
+ if (ok)
+ cout << "Simulation "<<fileName<<" is successful !" << endl;
+ else
+ cout << "Simulation "<<fileName<<" failed ! " << endl;
+
+ cout << "------------ End of calculation !!! -----------" << endl;
+ myProblem.terminate();
+
+ return EXIT_SUCCESS;
+}
CreatePythonTest(SinglePhase/SinglePhase_1DDepressurisation.py)
CreatePythonTest(SinglePhase/SinglePhase_1DHeatedAssembly.py)
CreatePythonTest(SinglePhase/SinglePhase_1DHeatedChannel.py)
+CreatePythonTest(SinglePhase/SinglePhase_1DHeatedChannel_Implicit.py)
CreatePythonTest(SinglePhase/SinglePhase_1DRiemannProblem.py)
+CreatePythonTest(SinglePhase/SinglePhase_1DRiemannProblem_Implicit.py)
+CreatePythonTest(SinglePhase/SinglePhase_1DRiemannProblem_Implicit_LineSearch.py)
CreatePythonTest(SinglePhase/SinglePhase_1DWaterHammer.py)
CreatePythonTest(SinglePhase/SinglePhase_2BranchesHeatedChannels.py)
CreatePythonTest(SinglePhase/SinglePhase_2DVidangeReservoir.py)
--- /dev/null
+#!/usr/bin/env python3
+# -*-coding:utf-8 -*
+
+import CoreFlows as cf
+import matplotlib.pyplot as plt
+import cdmath as cm
+import VTK_routines
+
+def SinglePhase_1DHeatedChannel_Implicit():
+
+ spaceDim = 1;
+ # Prepare for the mesh
+ print("Building mesh " );
+ xinf = 0 ;
+ xsup=4.2;
+ nx=50;
+
+ # set the limit field for each boundary
+
+ inletVelocityX=5;
+ inletTemperature=565;
+ outletPressure=155e5;
+
+ # physical parameters
+ heatPower=1e8;
+
+ myProblem = cf.SinglePhase(cf.Liquid,cf.around155bars600K,spaceDim);
+ nVar = myProblem.getNumberOfVariables();
+
+ # Prepare for the initial condition
+ VV_Constant =[0]*nVar;
+
+ # constant vector
+ VV_Constant[0] = outletPressure ;
+ VV_Constant[1] = inletVelocityX;
+ VV_Constant[2] = inletTemperature ;
+
+
+ #Initial field creation
+ print("Building initial data " );
+ myProblem.setInitialFieldConstant( spaceDim, VV_Constant, xinf, xsup, nx,"inlet","outlet");
+
+ # set the boundary conditions
+ myProblem.setInletBoundaryCondition("inlet",inletTemperature,inletVelocityX)
+ myProblem.setOutletBoundaryCondition("outlet", outletPressure,[xsup]);
+
+ # set physical parameters
+ myProblem.setHeatSource(heatPower);
+
+ # set the numerical method
+ myProblem.setNumericalScheme(cf.upwind, cf.Implicit);
+
+ # name of result file
+ fileName = "1DHeatedChannelUpwind_Implicit";
+
+ # simulation parameters
+ MaxNbOfTimeStep = 1000 ;
+ freqSave = 100;
+ cfl = 100;
+ maxTime = 500;
+ precision = 1e-7;
+
+ myProblem.setCFL(cfl);
+ myProblem.setPrecision(precision);
+ myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
+ myProblem.setTimeMax(maxTime);
+ myProblem.setFreqSave(freqSave);
+ myProblem.setFileName(fileName);
+ myProblem.setNewtonSolver(precision,20);
+ myProblem.saveConservativeField(True);
+ if(spaceDim>1):
+ myProblem.saveVelocity();
+ pass
+
+ myProblem.setLinearSolver(cf.GMRES, cf.ILU)
+ myProblem.setNewtonSolver(1e-3, 50, cf.Newton_PETSC_LINESEARCH)
+
+ # evolution
+ myProblem.initialize();
+
+ #Postprocessing
+ plt.xlabel('x')
+ plt.ylabel('Pressure')
+ plt.xlim(xinf,xsup)
+ plt.ylim( 0.999*outletPressure, 1.001*outletPressure )
+ plt.title('Solving Riemann problem for Euler equations\n with Finite volume schemes method')
+ dx=(xsup-xinf)/nx
+ x=[ i*dx for i in range(nx)] # array of cell center (1D mesh)
+
+ myPressureField = myProblem.getPressureField()
+ pressureArray=myPressureField.getFieldValues()
+ line_pressure, = plt.plot(x, pressureArray, label='Pressure time step 0')
+ plt.legend()
+ plt.savefig(fileName+".png")
+
+ ok = myProblem.run();
+
+ myPressureField = myProblem.getPressureField()
+ timeStep=myProblem.getNbTimeStep()#Final time step
+ pressureArray=myPressureField.getFieldValues()
+ line_pressure, = plt.plot(x, pressureArray, label='Pressure time step '+str(timeStep))
+ plt.legend()
+ plt.savefig(fileName+".png")
+
+ if (ok):
+ print( "Simulation python " + fileName + " is successful !" );
+ pass
+ else:
+ print( "Simulation python " + fileName + " failed ! " );
+ pass
+
+ print( "------------ End of calculation !!! -----------" );
+
+ myProblem.terminate();
+ return ok
+
+if __name__ == """__main__""":
+ SinglePhase_1DHeatedChannel_Implicit()
--- /dev/null
+#!/usr/bin/env python
+# -*-coding:utf-8 -*
+
+import CoreFlows as cf
+import cdmath as cm
+
+def SinglePhase_1DRiemannProblem_Implicit():
+
+ spaceDim = 1;
+ # Prepare for the mesh
+ print("Building mesh " );
+ xinf = 0 ;
+ xsup=4.2;
+ nx=100;
+ discontinuity=(xinf+xsup)/2
+ M=cm.Mesh(xinf,xsup,nx)
+ eps=1e-6
+ M.setGroupAtPlan(xsup,0,eps,"RightBoundary")
+ M.setGroupAtPlan(xinf,0,eps,"LeftBoundary")
+
+ # Prepare initial data
+ initialVelocity_Left=1;
+ initialTemperature_Left=565;
+ initialPressure_Left=155e5;
+ initialVelocity_Right=1;
+ initialTemperature_Right=565;
+ initialPressure_Right=1e5;
+
+ myProblem = cf.SinglePhase(cf.Liquid,cf.around155bars600K,spaceDim);
+ nVar = myProblem.getNumberOfVariables();
+
+ # Prepare for the initial condition
+ VV_Left = cm.Vector(nVar)
+ VV_Right = cm.Vector(nVar)
+
+ # left and right constant vectors
+ VV_Left[0] = initialPressure_Left;
+ VV_Left[1] = initialVelocity_Left;
+ VV_Left[2] = initialTemperature_Left ;
+ VV_Right[0] = initialPressure_Right;
+ VV_Right[1] = initialVelocity_Right;
+ VV_Right[2] = initialTemperature_Right ;
+
+
+ #Initial field creation
+ print("Building initial data " );
+ myProblem.setInitialFieldStepFunction(M,VV_Left,VV_Right,discontinuity);
+
+ # set the boundary conditions
+ myProblem.setNeumannBoundaryCondition("LeftBoundary");
+ myProblem.setNeumannBoundaryCondition("RightBoundary");
+
+ # set the numerical method
+ myProblem.setNumericalScheme(cf.upwind, cf.Implicit);
+
+ # name of result file
+ fileName = "1DRiemannProblem_Implicit";
+
+ # simulation parameters
+ MaxNbOfTimeStep = 3 ;
+ freqSave = 1;
+ cfl = 1;
+ maxTime = 500;
+ precision = 1e-6;
+
+ myProblem.setCFL(cfl);
+ myProblem.setPrecision(precision);
+ myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
+ myProblem.setTimeMax(maxTime);
+ myProblem.setFreqSave(freqSave);
+ myProblem.setFileName(fileName);
+ myProblem.setSaveFileFormat(cf.CSV)
+ myProblem.saveConservativeField(True);
+
+ myProblem.setLinearSolver(cf.GMRES, cf.LU)
+ myProblem.setNewtonSolver(precision,20, cf.Newton_SOLVERLAB)
+ myProblem.usePrimitiveVarsInNewton(False)
+
+ # evolution
+ myProblem.initialize();
+
+ ok = myProblem.run();
+ if (ok):
+ print( "Simulation python " + fileName + " is successful !" );
+ pass
+ else:
+ print( "Simulation python " + fileName + " failed ! " );
+ pass
+
+ print( "------------ End of calculation !!! -----------" );
+
+ myProblem.terminate();
+ return ok
+
+if __name__ == """__main__""":
+ SinglePhase_1DRiemannProblem_Implicit()
--- /dev/null
+#!/usr/bin/env python
+# -*-coding:utf-8 -*
+
+import CoreFlows as cf
+import cdmath as cm
+import matplotlib.pyplot as plt
+from numpy.linalg import norm
+import VTK_routines
+import exact_rs_stiffenedgas
+
+def SinglePhase_1DRiemannProblem_Implicit_LineSearch():
+
+ spaceDim = 1;
+ # Prepare for the mesh
+ print("Building mesh " );
+ xinf = 0 ;
+ xsup=4.2;
+ nx=50;
+ discontinuity=(xinf+xsup)/2
+ M=cm.Mesh(xinf,xsup,nx)
+ eps=1e-6
+ M.setGroupAtPlan(xsup,0,eps,"RightBoundary")
+ M.setGroupAtPlan(xinf,0,eps,"LeftBoundary")
+
+ # Prepare initial data
+ initialVelocity_Left=1;
+ initialTemperature_Left=565;
+ initialPressure_Left=155e5;
+ initialVelocity_Right=1;
+ initialTemperature_Right=565;
+ initialPressure_Right=1e5;
+
+ myProblem = cf.SinglePhase(cf.Liquid,cf.around155bars600K,spaceDim);
+ nVar = myProblem.getNumberOfVariables();
+
+ # Prepare for the initial condition
+ VV_Left = cm.Vector(nVar)
+ VV_Right = cm.Vector(nVar)
+
+ # left and right constant vectors
+ VV_Left[0] = initialPressure_Left;
+ VV_Left[1] = initialVelocity_Left;
+ VV_Left[2] = initialTemperature_Left ;
+ VV_Right[0] = initialPressure_Right;
+ VV_Right[1] = initialVelocity_Right;
+ VV_Right[2] = initialTemperature_Right ;
+
+
+ #Initial field creation
+ print("Building initial data " );
+ myProblem.setInitialFieldStepFunction(M,VV_Left,VV_Right,discontinuity);
+
+ # set the boundary conditions
+ myProblem.setNeumannBoundaryCondition("LeftBoundary");
+ myProblem.setNeumannBoundaryCondition("RightBoundary");
+
+ # set the numerical method
+ myProblem.setNumericalScheme(cf.upwind, cf.Implicit);
+
+ # name of result file
+ fileName = "1DRiemannProblem_Implicit_LineSearch";
+
+ # simulation parameters
+ MaxNbOfTimeStep = 3 ;
+ freqSave = 1;
+ cfl = 1;
+ maxTime = 500;
+ precision = 1e-6;
+
+ myProblem.setCFL(cfl);
+ myProblem.setPrecision(precision);
+ myProblem.setMaxNbOfTimeStep(MaxNbOfTimeStep);
+ myProblem.setTimeMax(maxTime);
+ myProblem.setFreqSave(freqSave);
+ myProblem.setFileName(fileName);
+ myProblem.setSaveFileFormat(cf.CSV)
+ myProblem.saveConservativeField(True);
+
+ myProblem.setLinearSolver(cf.GMRES, cf.ILU)
+ myProblem.setNewtonSolver(precision,20, cf.Newton_PETSC_LINESEARCH)
+ myProblem.usePrimitiveVarsInNewton(False)
+
+ # evolution
+ myProblem.initialize();
+
+ ### Postprocessing initial data
+ dx=(xsup-xinf)/nx
+ x=[ i*dx for i in range(nx+1)] # array of cell center (1D mesh)
+ fig, ([axDensity, axPressure], [axVelocity, axTemperature]) = plt.subplots(2, 2,sharex=True, figsize=(10,10))
+ plt.gcf().subplots_adjust(wspace = 0.5)
+
+ myEOS = myProblem.getFluidEOS()## Needed to retrieve gamma, pinfnity, convert (p,T) to density and (p, rho) to temperature
+
+ axPressure.set(xlabel='x (m)', ylabel='Pressure (bar)')
+ axPressure.set_xlim(xinf,xsup)
+ axPressure.set_ylim( initialPressure_Right - 0.1*(initialPressure_Left-initialPressure_Right), initialPressure_Left + 0.5*(initialPressure_Left-initialPressure_Right) )
+
+ myPressureField = myProblem.getPressureField()
+ myPressureField.writeVTK("PressureField")
+ pressureArray=VTK_routines.Extract_VTK_data_over_line_to_numpyArray("PressureField"+"_0.vtu", [xinf,0,0], [xsup,0,0],nx)
+ axPressure.plot(x, pressureArray, label='Pressure time step 0')
+
+ initialDensity_Left = myEOS.getDensity( initialPressure_Left, initialTemperature_Left)
+ initialDensity_Right = myEOS.getDensity( initialPressure_Right, initialTemperature_Right)
+
+ axDensity.set(xlabel='x (m)', ylabel='Density (Kg/m3)')
+ axDensity.set_xlim(xinf,xsup)
+ axDensity.set_ylim( initialDensity_Right - 0.1*(initialDensity_Left-initialDensity_Right), initialDensity_Left + 0.5*(initialDensity_Left-initialDensity_Right) )
+
+ myDensityField = myProblem.getDensityField()
+ myDensityField.writeVTK("DensityField")
+ densityArray=VTK_routines.Extract_VTK_data_over_line_to_numpyArray("DensityField"+"_0.vtu", [xinf,0,0], [xsup,0,0],nx)
+ axDensity.plot(x, densityArray, label='Density time step 0')
+
+ axVelocity.set(xlabel='x (m)', ylabel='Velocity (m/s)')
+ axVelocity.set_xlim(xinf,xsup)
+ axVelocity.set_ylim( 0.9*initialVelocity_Right - 0.1*(initialVelocity_Left-initialVelocity_Right), 22*initialVelocity_Left + 0.5*(initialVelocity_Left-initialVelocity_Right) )
+
+ myVelocityField = myProblem.getVelocityXField()
+ myVelocityField.writeVTK("VelocityField")
+ velocityArray=VTK_routines.Extract_VTK_data_over_line_to_numpyArray("VelocityField"+"_0.vtu", [xinf,0,0], [xsup,0,0],nx)
+ axVelocity.plot(x, velocityArray, label='Velocity time step 0')
+
+ axTemperature.set(xlabel='x (m)', ylabel='Temperature (K)')
+ axTemperature.set_xlim(xinf,xsup)
+ axTemperature.set_ylim( 0.999*initialTemperature_Right - 0.1*(initialTemperature_Left-initialTemperature_Right), 1.001*initialTemperature_Left + 0.5*(initialTemperature_Left-initialTemperature_Right) )
+
+ myTemperatureField = myProblem.getTemperatureField()
+ myTemperatureField.writeVTK("TemperatureField")
+ temperatureArray=VTK_routines.Extract_VTK_data_over_line_to_numpyArray("TemperatureField"+"_0.vtu", [xinf,0,0], [xsup,0,0],nx)
+ axTemperature.plot(x, temperatureArray, label='Temperature time step 0')
+
+ ### run simulation
+ ok = myProblem.run();
+
+ #Determine exact solution
+ exactDensity, exactVelocity, exactPressure = exact_rs_stiffenedgas.exact_sol_Riemann_problem(xinf, xsup, myProblem.presentTime(), myEOS.constante("gamma"), myEOS.constante("p0"), [ initialDensity_Left, initialVelocity_Left, initialPressure_Left ], [ initialDensity_Right, initialVelocity_Right, initialPressure_Right ], (xinf+xsup)/2, nx+1)
+
+ ### Plot curves
+ axPressure.plot(x, exactPressure, label='Exact Pressure ')
+ myPressureField = myProblem.getPressureField()
+ myPressureField.writeVTK("PressureField")
+ pressureArray=VTK_routines. Extract_VTK_data_over_line_to_numpyArray("PressureField_"+str(myProblem.getNbTimeStep())+".vtu", [xinf,0,0], [xsup,0,0],nx)
+ axPressure.plot(x, pressureArray, label='Pressure time step '+str(myProblem.getNbTimeStep()))
+ axPressure.legend()
+
+ axDensity.plot(x, exactDensity, label='Exact Density ')
+ myDensityField = myProblem.getDensityField()
+ myDensityField.writeVTK("DensityField")
+ densityArray=VTK_routines. Extract_VTK_data_over_line_to_numpyArray("DensityField_"+str(myProblem.getNbTimeStep())+".vtu", [xinf,0,0], [xsup,0,0],nx)
+ axDensity.plot(x, densityArray, label='Density time step '+str(myProblem.getNbTimeStep()))
+ axDensity.legend()
+
+ axVelocity.plot(x, exactVelocity, label='Exact Velocity ')
+ myVelocityField = myProblem.getVelocityXField()
+ myVelocityField.writeVTK("VelocityField")
+ velocityArray=VTK_routines. Extract_VTK_data_over_line_to_numpyArray("VelocityField_"+str(myProblem.getNbTimeStep())+".vtu", [xinf,0,0], [xsup,0,0],nx)
+ axVelocity.plot(x, velocityArray, label='Velocity time step '+str(myProblem.getNbTimeStep()))
+ axVelocity.legend()
+
+ exactTemperature = [0.]*(nx+1)
+ for i in range(nx+1):
+ exactTemperature[i] = myEOS.getTemperatureFromPressure(exactPressure[i], exactDensity[i])
+
+ axTemperature.plot(x, exactTemperature, label='Exact Temperature ')
+ myTemperatureField = myProblem.getTemperatureField()
+ myTemperatureField.writeVTK("TemperatureField")
+ temperatureArray=VTK_routines. Extract_VTK_data_over_line_to_numpyArray("TemperatureField_"+str(myProblem.getNbTimeStep())+".vtu", [xinf,0,0], [xsup,0,0],nx)
+ axTemperature.plot(x, temperatureArray, label='Temperature time step '+str(myProblem.getNbTimeStep()))
+ axTemperature.legend()
+
+ #plt.title('Solving Riemann problem for Euler equations\n with Finite volume schemes method')
+ plt.savefig(fileName+".png")
+
+ #Compute numerical error
+ error_pressure = norm( exactPressure - pressureArray )/norm( exactPressure )
+ print('relative error on pressure = ', error_pressure )
+
+ if (ok):
+ print( "Simulation python " + fileName + " is successful !" );
+ pass
+ else:
+ print( "Simulation python " + fileName + " failed ! " );
+ pass
+
+ print( "------------ End of calculation !!! -----------" );
+
+
+ myProblem.terminate();
+ return ok
+
+if __name__ == """__main__""":
+ SinglePhase_1DRiemannProblem_Implicit_LineSearch()
--- /dev/null
+#!/usr/bin/env python3
+# -*-coding:utf-8 -*
+
+######################################################################################################################
+# This file contains a class to solve for the exact solution of the Riemann Problem for the one dimensional Euler
+# equations with stiffened gas equation of state
+#
+# Author: Michael Ndjinga
+# Date: 18/02/2021
+# Description : Translated from C++ package developped by Murray Cutforth
+#######################################################################################################################
+
+from math import pow, fabs, sqrt
+
+def exact_sol_Riemann_problem(xmin, xmax, t, gamma, p0, WL, WR, offset, numsamples = 100):#offset= position of the initial discontinuity
+ print("")
+ print("Determination of the exact solution of the Riemann problem for the Euler equations, gamma=", gamma, ", p0= ", p0)
+
+ RS = exact_rs_stiffenedgas(gamma, gamma, p0, p0);
+ RS.solve_RP(WL,WR);
+
+ delx = (xmax - xmin)/numsamples;
+
+ density = [0.]*numsamples
+ velocity = [0.]*numsamples
+ pressure = [0.]*numsamples
+
+ for i in range(numsamples):
+ S = i*delx/t;
+ soln = RS.sample_solution(WL, WR, S - offset/t);
+ density[i] = soln[0]
+ velocity[i]= soln[1]
+ pressure[i]= soln[2]
+
+ return density, velocity, pressure
+
+class exact_rs_stiffenedgas :
+
+ def __init__(self, gamma_L, gamma_R, pinf_L, pinf_R, tol=1.e-6, max_iter=100):
+ self.TOL = tol
+ self.MAX_NB_ITER = max_iter
+
+ self.gamma_L = gamma_L
+ self.gamma_R = gamma_R
+ self.pinf_L = pinf_L
+ self.pinf_R = pinf_R
+
+ self.S_STAR = 0.
+ self.P_STAR = 0.
+ self.rho_star_L = 0.
+ self.rho_star_R = 0.
+
+ self.S_L = 0.
+ self.S_R = 0.
+ self.S_HL = 0.
+ self.S_TL = 0.
+ self.S_HR = 0.
+ self.S_TR = 0.
+
+
+
+ # Functions used to generate exact solutions to Riemann problems
+
+ def solve_RP (self, W_L, W_R):
+ assert len(W_L) == 3, "Left state should have three components (rho, u p)"
+ assert len(W_R) == 3, "Right state should have three components (rho, u p)"
+ assert W_L[0] >= 0.0, "Left density should be positive"
+ assert W_R[0] >= 0.0, "Right density should be positive"
+ # assert W_L[2] >= 0.0 # Since stiffened gases will often exhibit p<0..
+ # assert W_R[2] >= 0.0 #
+
+ print("")
+ print("Solving Riemann problem for left state W_L=", W_L, ", and right state W_R=",W_R)
+
+ # Calculate p_star
+
+ self.P_STAR = self.find_p_star_newtonraphson(W_L[0], W_L[1], W_L[2], W_R[0], W_R[1], W_R[2])
+
+
+ # Calculate u_star
+
+ self.S_STAR = 0.5*(W_L[1]+W_R[1]) + 0.5*(self.f(self.P_STAR,W_R[0],W_R[2],self.gamma_R,self.pinf_R) - self.f(self.P_STAR,W_L[0],W_L[2],self.gamma_L,self.pinf_L))
+
+
+ # Solution now depends on character of 1st and 3rd waves
+
+ if (self.P_STAR > W_L[2]):
+ # Left shock
+
+ self.rho_star_L = W_L[0]*((2.0*self.gamma_L*self.pinf_L + (self.gamma_L+1.0)*self.P_STAR + (self.gamma_L-1.0)*W_L[2])/(2.0*(W_L[2] + self.gamma_L*self.pinf_L) + (self.gamma_L-1.0)*self.P_STAR + (self.gamma_L-1.0)*W_L[2]))
+ self.S_L = W_L[1] - (self.Q_K(self.P_STAR,W_L[0],W_L[2],self.gamma_L,self.pinf_L)/W_L[0])
+ else:
+ # Left rarefaction
+
+ self.rho_star_L = W_L[0]*pow((self.P_STAR + self.pinf_L)/(W_L[2] + self.pinf_L), 1.0/self.gamma_L)
+
+ a_L = self.a(W_L[0], W_L[2], self.gamma_L, self.pinf_L)
+ a_star_L = a_L*pow((self.P_STAR + self.pinf_L)/(W_L[2] + self.pinf_L), (self.gamma_L-1.0)/(2.0*self.gamma_L))
+
+ self.S_HL = W_L[1] - a_L
+ self.S_TL = self.S_STAR - a_star_L
+
+ if (self.P_STAR > W_R[2]):
+ # Right shock
+
+ self.rho_star_R = W_R[0]*((2.0*self.gamma_R*self.pinf_R + (self.gamma_R+1.0)*self.P_STAR + (self.gamma_R-1.0)*W_R[2])/(2.0*(W_R[2] + self.gamma_R*self.pinf_R) + (self.gamma_R-1.0)*self.P_STAR + (self.gamma_R-1.0)*W_R[2]))
+
+ self.S_R = W_R[1] + (self.Q_K(self.P_STAR,W_R[0],W_R[2],self.gamma_R,self.pinf_R)/W_R[0])
+ else:
+ # Right rarefaction
+
+ self.rho_star_R = W_R[0]*pow((self.P_STAR + self.pinf_R)/(W_R[2] + self.pinf_R), 1.0/self.gamma_R)
+
+ a_R = self.a(W_R[0],W_R[2],self.gamma_R, self.pinf_R)
+ a_star_R = a_R*pow((self.P_STAR + self.pinf_R)/(W_R[2] + self.pinf_R), (self.gamma_R-1.0)/(2.0*self.gamma_R))
+
+ self.S_HR = W_R[1] + a_R
+ self.S_TR = self.S_STAR + a_star_R
+
+ def sample_solution (self, W_L, W_R, S):
+ W = [0.]*3
+
+ # Find appropriate part of solution and return primitives
+
+ if (S < self.S_STAR):
+ # To the left of the contact
+
+ if (self.P_STAR > W_L[2]):
+ # Left shock
+
+ if (S < self.S_L):
+ W = W_L
+ else:
+ W[0] = self.rho_star_L
+ W[1] = self.S_STAR
+ W[2] = self.P_STAR
+ else:
+ # Left rarefaction
+
+ if (S < self.S_HL):
+ W = W_L
+ else:
+ if (S > self.S_TL):
+ W[0] = self.rho_star_L
+ W[1] = self.S_STAR
+ W[2] = self.P_STAR
+ else:
+ self.set_left_rarefaction_fan_state(W_L, S, W)
+ else:
+ # To the right of the contact
+
+ if (self.P_STAR > W_R[2]):
+ # Right shock
+
+ if (S > self.S_R):
+ W = W_R
+ else:
+ W[0] = self.rho_star_R
+ W[1] = self.S_STAR
+ W[2] = self.P_STAR
+ else:
+ # Right rarefaction
+
+ if (S > self.S_HR):
+ W = W_R
+ else:
+ if (S < self.S_TR):
+ W[0] = self.rho_star_R
+ W[1] = self.S_STAR
+ W[2] = self.P_STAR
+ else:
+ self.set_right_rarefaction_fan_state(W_R, S, W)
+
+ return W
+
+ # Functions used to solve for p_star iteratively
+
+ def find_p_star_newtonraphson (self, rho_L, u_L, p_L, rho_R, u_R, p_R ):
+
+ # First we set the initial guess for p_star using a simple mean-value approximation
+
+ p_star_next = 0.5*(p_L+p_R)
+ n = 0
+
+
+ # Now use the Newton-Raphson algorithm
+
+ while True:#conversion of do ... while by while True... if (...) break
+ p_star = p_star_next
+
+ p_star_next = p_star - self.total_pressure_function(p_star,rho_L,u_L,p_L,rho_R,u_R,p_R)/self.total_pressure_function_deriv(p_star,rho_L,p_L,rho_R,p_R)
+
+ p_star_next = max(p_star_next, self.TOL)
+
+ n+=1
+
+ if not ((fabs(p_star_next - p_star)/(0.5*(p_star+p_star_next)) > self.TOL) and n < self.MAX_NB_ITER):
+ break
+
+ if (n == self.MAX_NB_ITER):
+ raise ValueError("!!!!!!!!!!Newton algorithm did not converge. Increase tolerance or maximum number of time steps. Current values : tol=" + str(self.TOL) + ", max_iter=" + str(self.MAX_NB_ITER) )
+ #p_star_next = 0.5*(p_L+p_R)
+
+ return p_star_next
+
+ def total_pressure_function (self, p_star, rho_L, u_L, p_L, rho_R, u_R, p_R ):
+
+ return self.f(p_star, rho_L, p_L, self.gamma_L, self.pinf_L) + self.f(p_star, rho_R, p_R, self.gamma_R, self.pinf_R) + u_R - u_L
+
+ def total_pressure_function_deriv (self, p_star, rho_L, p_L, rho_R, p_R ):
+
+ return self.f_deriv (p_star, rho_L, p_L, self.gamma_L, self.pinf_L) + self.f_deriv (p_star, rho_R, p_R, self.gamma_R, self.pinf_R)
+
+
+ def f (self, p_star, rho, p, gamma, pinf):
+ if (p_star > p):
+
+ return (p_star - p)/self.Q_K(p_star, rho, p, gamma, pinf)
+
+ else:
+
+ return (2.0*self.a(rho,p,gamma,pinf)/(gamma-1.0))*(pow((p_star + pinf)/(p + pinf), (gamma-1.0)/(2.0*gamma)) - 1.0)
+
+
+ def f_deriv (self, p_star, rho, p, gamma, pinf):
+ A = 2.0/((gamma+1.0)*rho)
+ B = (p+pinf)*(gamma-1.0)/(gamma+1.0)
+
+ if (p_star > p):
+
+ return sqrt(A/(B+p_star+pinf))*(1.0 - ((p_star-p)/(2.0*(B+p_star+pinf))))
+
+ else:
+
+ return (1.0/(rho*self.a(rho,p,gamma,pinf)))*pow((p_star+pinf)/(p+pinf), -(gamma+1.0)/(2.0*gamma))
+
+
+
+ # Functions to find the state inside a rarefaction fan
+
+ def set_left_rarefaction_fan_state (self, W_L, S, W):
+ a_L = self.a(W_L[0],W_L[2],self.gamma_L,self.pinf_L)
+ W[0] = W_L[0]*pow((2.0/(self.gamma_L+1.0)) + ((self.gamma_L-1.0)/(a_L*(self.gamma_L+1.0)))*(W_L[1] - S), 2.0/(self.gamma_L - 1.0))
+ W[1] = (2.0/(self.gamma_L+1.0))*(a_L + S + ((self.gamma_L-1.0)/2.0)*W_L[1])
+ W[2] = (W_L[2] + self.pinf_L)*pow((2.0/(self.gamma_L+1.0)) + ((self.gamma_L-1.0)/(a_L*(self.gamma_L+1.0)))*(W_L[1] - S), (2.0*self.gamma_L)/(self.gamma_L-1.0)) - self.pinf_L
+
+ def set_right_rarefaction_fan_state (self, W_R, S, W):
+ a_R = self.a(W_R[0],W_R[2],self.gamma_R,self.pinf_R)
+ W[0] = W_R[0]*pow((2.0/(self.gamma_R+1.0)) - ((self.gamma_R-1.0)/(a_R*(self.gamma_R+1.0)))*(W_R[1] - S), 2.0/(self.gamma_R - 1.0))
+ W[1] = (2.0/(self.gamma_R+1.0))*(- a_R + S + ((self.gamma_R-1.0)/2.0)*W_R[1])
+ W[2] = (W_R[2] + self.pinf_R)*pow((2.0/(self.gamma_R+1.0)) - ((self.gamma_R-1.0)/(a_R*(self.gamma_R+1.0)))*(W_R[1] - S), (2.0*self.gamma_R)/(self.gamma_R-1.0)) - self.pinf_R
+
+
+
+ # Misc functions
+
+ def Q_K (self, p_star, rho, p, gamma, pinf):
+ A = 2.0/((gamma+1.0)*rho)
+ B = (p+pinf)*(gamma-1.0)/(gamma+1.0)
+ return sqrt((p_star+pinf+B)/A)
+
+
+
+ # Equation of state functions
+
+ def a (self, rho, p, gamma, pinf):#sound speed
+ return sqrt(gamma*((p+pinf)/rho))
+
+
