From: michael Date: Mon, 7 Jun 2021 18:27:56 +0000 (+0200) Subject: Added test for line search methods in Newton algorithm X-Git-Tag: V9_8_0~97 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=50aeb479f4dfe239966c029191ef0612efd74900;p=tools%2Fsolverlab.git Added test for line search methods in Newton algorithm --- diff --git a/CoreFlows/examples/C/CMakeLists.txt b/CoreFlows/examples/C/CMakeLists.txt index 339b28f..873ca92 100755 --- a/CoreFlows/examples/C/CMakeLists.txt +++ b/CoreFlows/examples/C/CMakeLists.txt @@ -47,8 +47,10 @@ set( libs_for_tests CoreFlowsLibs ) 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}" ) @@ -61,20 +63,25 @@ CreateTestExecAndInstall(DriftModel_2DInclinedBoilingChannel.cxx "${libs_for_te 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}" ) @@ -86,7 +93,9 @@ CreateTestExecAndInstall(SinglePhase_2DWallHeatedChannel_ChangeSect.cxx "${libs 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}" ) @@ -95,6 +104,7 @@ CreateTestExecAndInstall(StationaryDiffusionEquation_2DFV_StructuredTriangles_Ne 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}" ) diff --git a/CoreFlows/examples/C/SinglePhase_1DRiemannProblem_Implicit.cxx b/CoreFlows/examples/C/SinglePhase_1DRiemannProblem_Implicit.cxx new file mode 100755 index 0000000..e520c4c --- /dev/null +++ b/CoreFlows/examples/C/SinglePhase_1DRiemannProblem_Implicit.cxx @@ -0,0 +1,89 @@ +#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 "<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() diff --git a/CoreFlows/examples/Python/SinglePhase/SinglePhase_1DRiemannProblem_Implicit.py b/CoreFlows/examples/Python/SinglePhase/SinglePhase_1DRiemannProblem_Implicit.py new file mode 100755 index 0000000..7a49ad2 --- /dev/null +++ b/CoreFlows/examples/Python/SinglePhase/SinglePhase_1DRiemannProblem_Implicit.py @@ -0,0 +1,96 @@ +#!/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() diff --git a/CoreFlows/examples/Python/SinglePhase/SinglePhase_1DRiemannProblem_Implicit_LineSearch.py b/CoreFlows/examples/Python/SinglePhase/SinglePhase_1DRiemannProblem_Implicit_LineSearch.py new file mode 100644 index 0000000..ecd8c2a --- /dev/null +++ b/CoreFlows/examples/Python/SinglePhase/SinglePhase_1DRiemannProblem_Implicit_LineSearch.py @@ -0,0 +1,193 @@ +#!/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() diff --git a/CoreFlows/examples/Python/SinglePhase/exact_rs_stiffenedgas.py b/CoreFlows/examples/Python/SinglePhase/exact_rs_stiffenedgas.py new file mode 100644 index 0000000..7a6037e --- /dev/null +++ b/CoreFlows/examples/Python/SinglePhase/exact_rs_stiffenedgas.py @@ -0,0 +1,269 @@ +#!/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)) + +