From 532a02f9c758b2c1278f0c8d6d5c9472dd7058cc Mon Sep 17 00:00:00 2001
From: michael <michael@localhost.localdomain>
Date: Sat, 12 Dec 2020 13:39:23 +0100
Subject: [PATCH] Now check maximum principle and conservativity after
 simulation

---
 ...inglePhase_2DSphericalExplosion_HEXAGON.py | 55 +++++++++++++------
 1 file changed, 39 insertions(+), 16 deletions(-)

diff --git a/CoreFlows/examples/Python/SinglePhase/SinglePhase_2DSphericalExplosion_HEXAGON.py b/CoreFlows/examples/Python/SinglePhase/SinglePhase_2DSphericalExplosion_HEXAGON.py
index cb7943b..89bf60c 100755
--- a/CoreFlows/examples/Python/SinglePhase/SinglePhase_2DSphericalExplosion_HEXAGON.py
+++ b/CoreFlows/examples/Python/SinglePhase/SinglePhase_2DSphericalExplosion_HEXAGON.py
@@ -7,35 +7,40 @@ import cdmath
 
 def SinglePhase_2DSphericalExplosion_HEXAGON():
 
-	inputfile="./resources/meshHexagonWithTriangles10.med";
+	inputfile="../resources/meshHexagonWithTriangles10.med";
 	my_mesh=cdmath.Mesh(inputfile);
 	spaceDim=2
 	
 	# Initial field data
 	nVar=2+spaceDim;
-	radius=0.5;
+	radius=0.35;
 	Center=cdmath.Vector(spaceDim);#default value is (0,0,0)
-	Vout=cdmath.Vector(nVar)
-	Vin =cdmath.Vector(nVar)
-	Vin[0]=155e5;
+	Vout = cdmath.Vector(nVar)
+	Vin  = cdmath.Vector(nVar)
+	
+	Pmin=100e5
+	Pmax=155e5
+	InitialTemperature = 563
+
+	Vin[0]=Pmax
 	Vin[1]=0;
 	Vin[2]=0;
-	Vin[3]=563;
-	Vout[0]=154e5;
+	Vin[3]=InitialTemperature
+	Vout[0]=Pmin;
 	Vout[1]=0;
 	Vout[2]=0;
-	Vout[3]=563;
+	Vout[3]=InitialTemperature
 
 	myProblem = cf.SinglePhase(cf.Liquid,cf.around155bars600K,spaceDim);
 
 	# Initial field creation
 	print ("Setting mesh and initial data" ) ;
-	myProblem.setInitialFieldSphericalStepFunction( my_mesh, Vout, Vin, radius, Center);
+	myProblem.setInitialFieldSphericalStepFunction( my_mesh, Vin, Vout, radius, Center);
 
 	# set the boundary conditions
 	wallVelocityX=0;
 	wallVelocityY=0;
-	wallTemperature=563;
+	wallTemperature=InitialTemperature;
 
 	myProblem.setWallBoundaryCondition("boundaries", wallTemperature, wallVelocityX, wallVelocityY);
 
@@ -47,8 +52,8 @@ def SinglePhase_2DSphericalExplosion_HEXAGON():
 
 	# parameters calculation
 	MaxNbOfTimeStep = 3 ;
-	freqSave = 5;
-	cfl = 0.49;
+	freqSave = 1;
+	cfl = 0.5;
 	maxTime = 5;
 	precision = 1e-6;
 
@@ -59,15 +64,15 @@ def SinglePhase_2DSphericalExplosion_HEXAGON():
 	myProblem.setFreqSave(freqSave);
 	myProblem.setFileName(fileName);
 	myProblem.setNewtonSolver(precision,20);
-	myProblem.saveConservativeField(False);
-	if(spaceDim>1):
-		myProblem.saveVelocity();
-		pass
+	myProblem.saveConservativeField(True);
 
 
 	# evolution
 	myProblem.initialize();
 
+	masseInitiale = abs( myProblem.getDensityField().integral()[0] )
+	initialTotalEnergy=abs(myProblem.getTotalEnergyField().integral()[0] )
+
 	ok = myProblem.run();
 	if (ok):
 		print( "Simulation python " + fileName + " is successful !" );
@@ -78,6 +83,24 @@ def SinglePhase_2DSphericalExplosion_HEXAGON():
 
 	print( "------------ End of calculation !!! -----------" );
 
+	# Control of the results
+	PressureField=myProblem.getPressureField()
+	print( "La pression est bornée par le maximum de pression initiale, à 10% près :\n (PressureField.max() - InitialPressure.max())/InitialPressure.max()= ", (PressureField.max() - Pmax)/Pmax )
+	assert PressureField.max() < Pmax
+
+	temperatureField=myProblem.getTemperatureField()
+	print("La température est constante à 0.1 % près : erreur relative = ", max(abs(temperatureField.max() - InitialTemperature),abs(temperatureField.min() - InitialTemperature))/InitialTemperature )
+	assert abs(temperatureField.max() - InitialTemperature)/InitialTemperature < 0.001
+	assert abs(temperatureField.min() - InitialTemperature)/InitialTemperature < 0.001
+
+	densityField=myProblem.getDensityField()
+	print("La masse totale est conservée au zero machine près : erreur relative = ", abs( (abs(densityField.integral()[0]) - masseInitiale)/masseInitiale ) )
+	assert abs( (abs(densityField.integral()[0]) - masseInitiale)/masseInitiale ) < 1e-14
+	
+	totalEnergyField=myProblem.getTotalEnergyField()
+	print("L'énergie totale est conservée au zero machine près : erreur relative = ", abs( (abs(totalEnergyField.integral()[0]) - initialTotalEnergy)/initialTotalEnergy ) )
+	assert abs( (abs(totalEnergyField.integral()[0]) - initialTotalEnergy)/initialTotalEnergy ) < 1e-13
+
 	myProblem.terminate();
 	return ok
 
-- 
2.39.2