From 15b146554d7dae3d7bc2b69f0931aaebc68b087d Mon Sep 17 00:00:00 2001
From: elemel <elise.lemeledo@math.uzh.ch>
Date: Tue, 11 Aug 2020 08:01:05 +0200
Subject: [PATCH] Added an attemot of streamline
---
.../PredefinedProblems/Euler_KarniR22.py | 41 +
.../PredefinedProblems/Euler_KarniSquare.py | 41 +
.../Euler_SODCste_2Fluid.py | 2 +-
.../PredefinedProblems/Euler_ShockTube.py | 41 +
.../PredefinedProblems/Euler_StillBubble.py | 41 +
.../SecondOrderGalerkin.txt | 4 +-
.../SecondOrderGalerkinStream.txt | 20 +
RunPostProcess.py | 25 +-
RunProblem.py | 12 +-
SourceCode/Mappers.py | 14 +-
SourceCode/PostProcess.py | 13 +-
SourceCode/PostProcessing/Plots_Matplotlib.py | 341 +++++++-
SourceCode/PostProcessing/Plots_Schlieren.py | 99 ++-
.../GoverningEquations/EulerEquations.py | 755 ++++++++++++++----
.../LinearAdvectionRotation.py | 298 +++++--
.../LinearAdvectionTranslation.py | 333 +++++---
SourceCode/Solve.py | 46 +-
SourceCode/_Solver/Solver_Spatial.py | 1 +
.../SpatialModifiers/Filtering_Streamline.py | 206 ++---
.../SpatialModifiers/SpatialAmendements.py | 11 +-
SourceCode/_Solver/SpatialSchemes/CG.py | 8 +-
SourceCode/_Solver/SpatialSchemes/CG_LFX.py | 233 +++---
.../_Solver/SpatialSchemes/CG_Primary.py | 327 ++++++++
Templates/Template_Problem.py | 12 +-
24 files changed, 2281 insertions(+), 643 deletions(-)
create mode 100644 PredefinedTests/PredefinedProblems/Euler_KarniR22.py
create mode 100644 PredefinedTests/PredefinedProblems/Euler_KarniSquare.py
create mode 100644 PredefinedTests/PredefinedProblems/Euler_ShockTube.py
create mode 100644 PredefinedTests/PredefinedProblems/Euler_StillBubble.py
create mode 100644 PredefinedTests/PredefinedSettings/SecondOrderGalerkinStream.txt
create mode 100644 SourceCode/_Solver/SpatialSchemes/CG_Primary.py
diff --git a/PredefinedTests/PredefinedProblems/Euler_KarniR22.py b/PredefinedTests/PredefinedProblems/Euler_KarniR22.py
new file mode 100644
index 0000000..b7dd902
--- /dev/null
+++ b/PredefinedTests/PredefinedProblems/Euler_KarniR22.py
@@ -0,0 +1,41 @@
+################################################################################
+##
+## Template file to design the properties of the problem of interest ##
+##
+################################################################################
+import numpy as np
+
+class ProblemData():
+ def __init__(self):
+
+ # ******************************************************************************
+ # Meshing
+ #*******************************************************************************
+ self.ProblemID = "Euler_Karni"
+ self.MeshName = "Tube"
+
+ # ******************************************************************************
+ # Equation
+ #*******************************************************************************
+ self.GoverningEquationsIndex = 1
+ self.EOSEquationsIndex = [2, 2, 2]
+
+ # *******************************************************************************
+ # Fluids properties
+ #********************************************************************************
+ # Select the indices speciying the fluids properties living on each subdomain
+ self.FluidIndex = [0, 1]
+
+ # The fluid at index 0 is the complementary to the specified ones.
+ # Only need to precise the location of other fluids
+ self.FluidInitialSubdomain = [(4,0,0,0.5,0.5,0.5,4,4)]
+
+ # Index of the initial conditions to apply on each subdomain
+ self.InitialConditions = [(7, 1.22), (6,)]
+
+ # Row: BcTag, coklumn, index of the condition to apply depending on the fluid that comes next to it
+ self.BoundaryConditions = np.array(\
+ [[(0,), (0,)],\
+ [(2,), (2,)],\
+ [(0,), (0,)],\
+ [(2,), (2,)]])
diff --git a/PredefinedTests/PredefinedProblems/Euler_KarniSquare.py b/PredefinedTests/PredefinedProblems/Euler_KarniSquare.py
new file mode 100644
index 0000000..9d5cb62
--- /dev/null
+++ b/PredefinedTests/PredefinedProblems/Euler_KarniSquare.py
@@ -0,0 +1,41 @@
+################################################################################
+##
+## Template file to design the properties of the problem of interest ##
+##
+################################################################################
+import numpy as np
+
+class ProblemData():
+ def __init__(self):
+
+ # ******************************************************************************
+ # Meshing
+ #*******************************************************************************
+ self.ProblemID = "Euler_KarniSquare"
+ self.MeshName = "Tube"
+
+ # ******************************************************************************
+ # Equation
+ #*******************************************************************************
+ self.GoverningEquationsIndex = 1
+ self.EOSEquationsIndex = [2, 2, 2]
+
+ # *******************************************************************************
+ # Fluids properties
+ #********************************************************************************
+ # Select the indices speciying the fluids properties living on each subdomain
+ self.FluidIndex = [0, 2]
+
+ # The fluid at index 0 is the complementary to the specified ones.
+ # Only need to precise the location of other fluids
+ self.FluidInitialSubdomain = [[[-0.25, -0.250],[0.25, -0.25],[0.25, 0.25],[-0.25,0.25]]]
+
+ # Index of the initial conditions to apply on each subdomain
+ self.InitialConditions = [(7, 1.22), (6,)]
+
+ # Row: BcTag, coklumn, index of the condition to apply depending on the fluid that comes next to it
+ self.BoundaryConditions = np.array(\
+ [[(0,), (0,)],\
+ [(2,), (2,)],\
+ [(0,), (0,)],\
+ [(2,), (2,)]])
diff --git a/PredefinedTests/PredefinedProblems/Euler_SODCste_2Fluid.py b/PredefinedTests/PredefinedProblems/Euler_SODCste_2Fluid.py
index e61814b..5018023 100644
--- a/PredefinedTests/PredefinedProblems/Euler_SODCste_2Fluid.py
+++ b/PredefinedTests/PredefinedProblems/Euler_SODCste_2Fluid.py
@@ -31,7 +31,7 @@ class ProblemData():
self.FluidInitialSubdomain = [(4,0,0,0.5,0.5,0.5,4,4)]
# Index of the initial conditions to apply on each subdomain
- self.InitialConditions = [(1,), (3,)]
+ self.InitialConditions = [(5,), (4,)]
# Row: BcTag, coklumn, index of the condition to apply depending on the fluid that comes next to it
self.BoundaryConditions = np.array(\
diff --git a/PredefinedTests/PredefinedProblems/Euler_ShockTube.py b/PredefinedTests/PredefinedProblems/Euler_ShockTube.py
new file mode 100644
index 0000000..e7bf9b4
--- /dev/null
+++ b/PredefinedTests/PredefinedProblems/Euler_ShockTube.py
@@ -0,0 +1,41 @@
+################################################################################
+##
+## Template file to design the properties of the problem of interest ##
+##
+################################################################################
+import numpy as np
+
+class ProblemData():
+ def __init__(self):
+
+ # ******************************************************************************
+ # Meshing
+ #*******************************************************************************
+ self.ProblemID = "Euler_Karni"
+ self.MeshName = "Tube"
+
+ # ******************************************************************************
+ # Equation
+ #*******************************************************************************
+ self.GoverningEquationsIndex = 1
+ self.EOSEquationsIndex = [2]
+
+ # *******************************************************************************
+ # Fluids properties
+ #********************************************************************************
+ # Select the indices speciying the fluids properties living on each subdomain
+ self.FluidIndex = [2]
+
+ # The fluid at index 0 is the complementary to the specified ones.
+ # Only need to precise the location of other fluids
+ self.FluidInitialSubdomain = []
+
+ # Index of the initial conditions to apply on each subdomain
+ self.InitialConditions = [(7, 1.22)]
+
+ # Row: BcTag, coklumn, index of the condition to apply depending on the fluid that comes next to it
+ self.BoundaryConditions = np.array(\
+ [[(0,), (0,)],\
+ [(2,), (2,)],\
+ [(0,), (0,)],\
+ [(2,), (2,)]])
diff --git a/PredefinedTests/PredefinedProblems/Euler_StillBubble.py b/PredefinedTests/PredefinedProblems/Euler_StillBubble.py
new file mode 100644
index 0000000..c13500e
--- /dev/null
+++ b/PredefinedTests/PredefinedProblems/Euler_StillBubble.py
@@ -0,0 +1,41 @@
+################################################################################
+##
+## Template file to design the properties of the problem of interest ##
+##
+################################################################################
+import numpy as np
+
+class ProblemData():
+ def __init__(self):
+
+ # ******************************************************************************
+ # Meshing
+ #*******************************************************************************
+ self.ProblemID = "Euler_StillBubble"
+ self.MeshName = "Tube"
+
+ # ******************************************************************************
+ # Equation
+ #*******************************************************************************
+ self.GoverningEquationsIndex = 1
+ self.EOSEquationsIndex = [2,2]
+
+ # *******************************************************************************
+ # Fluids properties
+ #********************************************************************************
+ # Select the indices speciying the fluids properties living on each subdomain
+ self.FluidIndex = [0, 2]
+
+ # The fluid at index 0 is the complementary to the specified ones.
+ # Only need to precise the location of other fluids
+ self.FluidInitialSubdomain = [(4, 0,0,0.5,0.5,0.5,4,4)]
+
+ # Index of the initial conditions to apply on each subdomain
+ self.InitialConditions = [(4,), (6,)]
+
+ # Row: BcTag, coklumn, index of the condition to apply depending on the fluid that comes next to it
+ self.BoundaryConditions = np.array(\
+ [[(0,), (0,)],\
+ [(2,), (2,)],\
+ [(0,), (0,)],\
+ [(2,), (2,)]])
diff --git a/PredefinedTests/PredefinedSettings/SecondOrderGalerkin.txt b/PredefinedTests/PredefinedSettings/SecondOrderGalerkin.txt
index ef1f1ac..134db27 100644
--- a/PredefinedTests/PredefinedSettings/SecondOrderGalerkin.txt
+++ b/PredefinedTests/PredefinedSettings/SecondOrderGalerkin.txt
@@ -12,8 +12,8 @@ Limiter_1 # Additives list Index of limiter (0=None)
# Temporal scheme properties
1 # Time scheme
2 2 # Properties (RK order, Dec order, as many integers as needed, in the order required by the scheme (see docs))
-0.4 # CFL
-3 # Tmax, in seconds
+0.2 # CFL
+1.5 # Tmax, in seconds
# Export properties
m # Verbose mode (n=none, m=minimal, d=debug)
diff --git a/PredefinedTests/PredefinedSettings/SecondOrderGalerkinStream.txt b/PredefinedTests/PredefinedSettings/SecondOrderGalerkinStream.txt
new file mode 100644
index 0000000..36ee022
--- /dev/null
+++ b/PredefinedTests/PredefinedSettings/SecondOrderGalerkinStream.txt
@@ -0,0 +1,20 @@
+# Spatial scheme properties
+1 # Index of the scheme
+2 B # Scheme properties (see doc for more details, if any, the order should be first)
+Original+Filtering_1 # Additives list Index of limiter (0=None)
+2 # Type of fluid spotter
+2 B # Fluid spotter properties
+
+# Quadrature properties
+1 # 0: Quadpy 1: Custom
+1 # Order of exactness
+
+# Temporal scheme properties
+1 # Time scheme
+2 2 # Properties (RK order, Dec order, as many integers as needed, in the order required by the scheme (see docs))
+0.2 # CFL
+1.5 # Tmax, in seconds
+
+# Export properties
+m # Verbose mode (n=none, m=minimal, d=debug)
+1 # Time steps export intervals
diff --git a/RunPostProcess.py b/RunPostProcess.py
index 40174be..3942eaa 100644
--- a/RunPostProcess.py
+++ b/RunPostProcess.py
@@ -1,6 +1,6 @@
"""
-The file :code:`RunPostProcess.py` provides a simple interface to post-process a previously exported solution. To access the solution, one needs
+The file :code:`RunPostProcess.py` provides a simple interface to post-process a previously exported solution. To access the solution, one needs
- the problem's name
- the mesh's resolution that has been used to compute the solution
@@ -26,39 +26,39 @@ Select then the type of postprocessing routine you would like to execute by comm
following lines.
.. code::
-
+
# Plotting the solution's evolution in plotly considering all the exported timestamps
PostProcess.Plots_AnimationsPlotly(ProblemDefinition, SettingsChoice, MeshResolution)
-
- # Plotting graphics at given time steps. Furnish either the list of time step to export the plots
+
+ # Plotting graphics at given time steps. Furnish either the list of time step to export the plots
# at (should match the exported times, look within the folder before hand as only 7 digits are considered in the export file name)
# or the keyword "all" for which each exported result will lead to a plot.
PostProcess.Plots_StaticMatplotlib(ProblemDefinition, SettingsChoice, "all", MeshResolution)
PostProcess.Plots_StaticPlotly(ProblemDefinition, SettingsChoice, [0.0, 1.179924], MeshResolution)
-All those routine are gathering predefined sets of plotting functions whose list is available under the
+All those routine are gathering predefined sets of plotting functions whose list is available under the
routine documentation in the :ref:`TechnicalManual` under the rubric "Automatic post-processing routines".
Feel free to add other pre-defined plotting-routine wrap-up in the file :code:`SourceCode/PostProcess.py` at your wish.
|
-Running the postprocessing tasks is then done as simply by the two lines
+Running the postprocessing tasks is then done as simply by the two lines
.. code::
-
+
$ python RunPostProcess.py
|
|
.. note::
-
+
- | A debugging routine is also available, plotting raw information graphically
-
+
| :code:`PostProcess.Plot_DebugMatplotlib(ProblemDefinition, SettingsChoice, ExplicitTimeStampList, MeshResolution)`
-
+
- Make sure that you are using python3, otherwise some module-depency may fail (in particular)
-
+
"""
import SourceCode.PostProcess as PostProcess
@@ -67,11 +67,12 @@ if __name__ == '__main__':
# Point to the problem you would like to solve, and define how you would like to tackle it
ProblemDefinition = "Euler_Karni"
SettingsChoice = "SecondOrderGalerkin"
- MeshResolution = 50
+ MeshResolution = 20
# Post - processing tasks
#PostProcess.Plot_DebugMatplotlib(ProblemDefinition, SettingsChoice, [0.0], MeshResolution)
#PostProcess.Plots_AnimationsPlotly(ProblemDefinition, SettingsChoice, MeshResolution)
PostProcess.Plots_StaticMatplotlib(ProblemDefinition, SettingsChoice, "all", MeshResolution)
#PostProcess.Plots_StaticSchlieren(ProblemDefinition, SettingsChoice, "all", MeshResolution, "rho")
+
#PostProcess.Plots_StaticPlotly(ProblemDefinition, SettingsChoice, [0.0, 1.179924], MeshResolution)
diff --git a/RunProblem.py b/RunProblem.py
index 1f7ba66..f2f8de6 100644
--- a/RunProblem.py
+++ b/RunProblem.py
@@ -14,14 +14,14 @@ The numerical solution can be computed by editing the fields of the file accordi
on two fluids, using a coarse mesh and a first order solver one can set the fields of the file :code:`RunProblem.py` to
.. code::
-
+
ProblemDefinition = "Advection_Translation_2Fluids"
SettingsChoice = "FirstOrderGalerkin"
MeshResolution = 20
|
-Solving the problem is then done as simply by the two lines
+Solving the problem is then done as simply by the two lines
.. code::
@@ -30,10 +30,10 @@ Solving the problem is then done as simply by the two lines
|
-There is then just left to run the script by
+There is then just left to run the script by
.. code::
-
+
$ python RunProblem.py
|
@@ -50,7 +50,7 @@ the associated variables names are used.
|
.. note::
-
+
- Make sure that you are using python3, otherwise some module-depency may fail (in particular)
- The solution is only computed and periodically exported. No post-processing task is performed through this script.
@@ -62,7 +62,7 @@ import SourceCode.Solve as Solve
if __name__ == '__main__':
# Point to the problem you would like to solve, and define how you would like to tackle it
ProblemDefinition = "Euler_Karni"
- SettingsChoice = "FirstOrderGalerkin"
+ SettingsChoice = "SecondOrderGalerkin"
MeshResolution = 20
# Solve the problem
diff --git a/SourceCode/Mappers.py b/SourceCode/Mappers.py
index fbbc9ca..3288743 100644
--- a/SourceCode/Mappers.py
+++ b/SourceCode/Mappers.py
@@ -75,7 +75,7 @@ def Governing_Equations(id):
(see the documentation of :code:`ProblemDefinition`)
----------------------------------------------------------------------------------------"""
- if id==0: return(GE.LinearAdvectionTranslationVectorTest)
+ if id==0: return(GE.LinearAdvectionTranslationVectorTest)
elif id==1: return(GE.EulerEquations)
elif id==2: return(GE.LinearAdvectionRotation)
elif id==3: return(GE.LinearAdvectionTranslation)
@@ -119,7 +119,7 @@ def Fluid_Properties(id, Model):
if id==0: Properties = ['b', 1.4, 0, 0.287] # Air, ideal
elif id==1: Properties = ['r', 1.249, 0, 0.091] # R22, ideal
- elif id==2: Properties = ['g', 0.5, 0, 2.08] # Helium, ideal
+ elif id==2: Properties = ['g', 1.67, 0, 2.08] # Helium, ideal
try: Prop = (eval("FM."+ Model)(*Properties))
except: raise ValueError("Error in the problem definition:\n\
@@ -187,8 +187,9 @@ def Spatial_Scheme(id):
*Implemented spatial schemes given the id*
- 1. | Continuous Galerkin
- 2. | Continuous Galerkin + Lax-Friedrichs flux
+ 0. | Continuous Galerkin
+ 1. | Continuous Galerkin + Lax-Friedrichs flux
+ 2. | Continuous Galerkin on Primary variables
|
@@ -199,6 +200,7 @@ def Spatial_Scheme(id):
if id==0: return(SC.CG)
elif id==1: return(SC.CG_LFX)
+ elif id==2: return(SC.CG_Primary)
else: raise ValueError("Error in the problem definition:\n\
Wrong index for the wished spatial scheme.\n\
@@ -260,7 +262,9 @@ def Filtering(id):
SpatialModifiers.Module (module): Returns the module associated to the wished limiter
----------------------------------------------------------------------------------"""
- if id==1: return(SM.Filtering_Streamline)
+ if id==1:
+ print("WARNING: You selected an amendement to the scheme that is known for having bugs and is still under debug.")
+ return(SM.Filtering_Streamline)
else: raise ValueError("Error in the problem definition:\n\
Wrong index for the wished limiter.\n\
diff --git a/SourceCode/PostProcess.py b/SourceCode/PostProcess.py
index 1877d00..f1dc1d4 100644
--- a/SourceCode/PostProcess.py
+++ b/SourceCode/PostProcess.py
@@ -99,6 +99,7 @@ def Plots_StaticSchlieren(*argv):
# Export the schlieren images
PSL.Compute2DSchlierenImages(Problem, Mesh, Solution, argv[4], ParametersId)
+ PSL.Compute2DSchlierenImagesEos(Problem, Mesh, Solution, "p", ParametersId)
#### Plotting with Matplotlib static visualisation tools ####
def Plots_StaticMatplotlib(*argv):
@@ -172,15 +173,19 @@ def Plots_StaticMatplotlib(*argv):
for VariableId in range(Problem.GoverningEquations.NbVariables):
#PMP.PlotVariableOnWholeDomainTri2D (Problem, Mesh, Solution, VariableId, ParametersId)
#PMP.PlotVariableOnEachSubDomainTri2D (Problem, Mesh, Solution, VariableId, ParametersId)
- #PMP.PlotVariableOnWholeDomainTri (Problem, Mesh, Solution, VariableId, ParametersId)
- PMP.PlotVariableOneEachSubDomainTri (Problem, Mesh, Solution, VariableId, ParametersId)
+ PMP.PlotVariableOnWholeDomainTri (Problem, Mesh, Solution, VariableId, ParametersId)
+ try: PMP.PlotVariableOneEachSubDomainTri (Problem, Mesh, Solution, VariableId, ParametersId)
+ except: pass
#PMP.PlotVariableOnWholeDomain2D (Problem, Mesh, Solution, VariableId, ParametersId)
#PMP.PlotVariableOnEachSubDomain2D (Problem, Mesh, Solution, VariableId, ParametersId)
#PMP.PlotVariableOnWholeDomain (Problem, Mesh, Solution, VariableId, ParametersId)
#PMP.PlotVariableOneEachSubDomain (Problem, Mesh, Solution, VariableId, ParametersId)
-
+
+ PMP.PlotVariableEoSOnWholeDomainTri (Problem, Mesh, Solution, "p", ParametersId)
+
#PMP.PlotLSValuesOnWholeDomainTri(Problem, Mesh, Solution, ParametersId)
- PMP.PlotLSValuesTri(Problem, Mesh, Solution, ParametersId)
+ try: PMP.PlotLSValuesTri(Problem, Mesh, Solution, ParametersId)
+ except: pass
#PMP.PlotLSValuesOnWholeDomainTri(Problem, Mesh, Solution, ParametersId)
#### Plotting with matpotlib some debug tools ####
diff --git a/SourceCode/PostProcessing/Plots_Matplotlib.py b/SourceCode/PostProcessing/Plots_Matplotlib.py
index 50befa0..0347062 100644
--- a/SourceCode/PostProcessing/Plots_Matplotlib.py
+++ b/SourceCode/PostProcessing/Plots_Matplotlib.py
@@ -505,6 +505,7 @@ def PlotLSValuesTri(Problem, Mesh, Solution, ParametersId):
plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
plt.close()
+
# ==============================================================================
# Level Set tools only usable for triangles
# ==============================================================================
@@ -542,9 +543,9 @@ def PlotLSValuesOnWholeDomainTri(Problem, Mesh, Solution, ParametersId):
# --------------- Customisation properties -----------------------------------
# Retrieve the variables names
- variableName = "LevelSet_{0!s}".format(LevelSetId)
+ variableName = "LevelSet_{0!s}".format(LevelSetId)
variableLatexName = "Level Set #{0!s}".format(LevelSetId)
- variableLatexUnit = ""
+ variableLatexUnit = ""
# Customisation of labels and latex conversion
rc('font',size = 12)
@@ -580,6 +581,7 @@ def PlotLSValuesOnWholeDomainTri(Problem, Mesh, Solution, ParametersId):
plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
plt.close()
+
# ==============================================================================
# Soltuion tools only usable for triangles
# ==============================================================================
@@ -869,6 +871,165 @@ def PlotVariableOneEachSubDomainTri(Problem, Mesh, Solution, VariableId, Paramet
plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
plt.close()
+#### Plot the variable retrieved from the EoS independently of the subdomain ####
+def PlotVariableEoSOnWholeDomainTri(Problem, Mesh, Solution, Field, ParametersId):
+ """ Plot one variable independently of the subdomain, given a triangulation
+
+ Args:
+ Problem (Problem): the considered problem and its properties
+ Mesh (MeshStructure): the considered mesh
+ Solution (Solution): the solution of the considered problem
+ Field (string): the name to pass to the EoS in order to retrieve the desired variable
+ ParametersId (string): name of the export folder (TestCase_TestSettings)
+
+ Returns:
+ None: only exporting the plot
+
+ --------------------------------------------------------------------------------------------"""
+
+ # --------------- Preparing the data to plot -----------------------------------
+ # Extracting the data of interest
+ X = Mesh.CartesianPoints[:,0]
+ Y = Mesh.CartesianPoints[:,1]
+ Z = Solution.RSol[:,:]
+
+ # Computing the value from the EoS
+ Zbuf = np.zeros((Z.shape[1]))
+ for i in set(map(int,Solution.RFluidFlag)):
+ ind = np.where(Solution.RFluidFlag==i)[0]
+ Zbuf[ind] = Problem.FluidEOS[i].EoS(Z[:,ind], [Problem.FluidProp[i]]*len(ind), Field)
+
+ # --------------- Creating the figure ------------------------------------------
+ fig = plt.figure()
+ ax = fig.gca(projection='3d')
+ ax.plot_trisurf(X, Y, Zbuf, triangles=Mesh.InnerElements, cmap=ColorMapsMatplotlib[0],\
+ edgecolor=None ,linewidth=0, antialiased=False)
+
+ # --------------- Customisation properties -----------------------------------
+ # Retrieve the variables names
+ variableName = Field
+ variableLatexName = Field
+ variableLatexUnit = Field
+
+ # Customisation of labels and latex conversion
+ rc('font',size = 12)
+ rc('font',family = 'serif')
+ rc('axes',labelsize = 14)
+
+ # Customise the positioning of the labels
+ [t.set_va('center') for t in ax.get_yticklabels()]
+ [t.set_ha('right') for t in ax.get_yticklabels()]
+ [t.set_va('center') for t in ax.get_xticklabels()]
+ [t.set_ha('left') for t in ax.get_xticklabels()]
+ [t.set_va('center') for t in ax.get_zticklabels()]
+ [t.set_ha('right') for t in ax.get_zticklabels()]
+
+ # Set the grid and the pane colors
+ ax.grid(False)
+ ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+ ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+ ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+
+ fig.suptitle(r'Solution for the variable {0!s} at the time ${1:5.3f}$'.format(variableLatexName, Solution.t))
+ ax.set_xlabel(r'x [{0!s}]'.format(Problem.GoverningEquations.XYLatexUnits[0]))
+ ax.set_ylabel(r'y [{0!s}]'.format(Problem.GoverningEquations.XYLatexUnits[1]))
+ ax.set_zlabel(r'{0!s} [{1!s}]'.format(variableLatexName, variableLatexUnit))
+
+
+ # --------------- Export the figure at the right location -------------------------
+ FolderName = os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(Mesh.MeshResolution).replace(".","_"), "3DPlots")
+ try: os.makedirs(FolderName)
+ except: pass
+
+ FileName = "Mplt_3DSolEoS_OnWholeDomainTri_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
+ plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
+ plt.close()
+
+#### Plot the variable coming from the EoS with a color scale depending of the subdomain ####
+def PlotVariableEoSOneEachSubDomainTri(Problem, Mesh, Solution, Field, ParametersId):
+ """ Plot one variable with a color scale depending of the subdomain
+
+ Args:
+ Problem (Problem): the considered problem and its properties
+ Mesh (MeshStructure): the considered mesh
+ Solution (Solution): the solution of the considered problem
+ Field (string): the name to pass to the EoS to retrieve the desired variable to be plot
+ ParametersId (string): name of the export folder (TestCase_TestSettings)
+
+ Returns:
+ None: only exporting the plot
+
+ --------------------------------------------------------------------------------------------"""
+
+ # --------------- Preparing the data to plot -----------------------------------
+ # Extracting the data of interest
+ X = Mesh.CartesianPoints[:,0]
+ Y = Mesh.CartesianPoints[:,1]
+ Z = Solution.RSol[:,:]
+ LS = Solution.RLSValues[:,:]
+
+ # Computing the value from the EoS
+ Zbuf = np.zeros((Z.shape[1]))
+ for i in set(Solution.RFluidFlag):
+ ind = np.where(Solution.RFluidFlag==i)[0]
+ Zbuf[ind] = Problem.FluidEOS[i].EoS(Z[:,ind], [Problem.FluidProp[i]]*len(ind), Field)
+
+ # --------------- Creating the figure ------------------------------------------
+ fig = plt.figure()
+ ax = fig.gca(projection='3d')
+
+ # Span each fluid present in the domain and plot one variable surface per area and the countour of each area
+ # according to the interpolated value of the FluidSelector
+ Fluids = set(map(int, Solution.RFluidFlag))
+ for i in Fluids:
+ # Masking the points that lie outisde the fluid's area
+ ind = np.where(LS[i,:]<=0)[0]
+ Ele = np.array(Mesh.InnerElements)
+ Mask = [set(ele).intersection(ind)==set(ele) for ele in Ele]
+
+ ax.plot_trisurf(X, Y, Zbuf, triangles = Mesh.InnerElements, cmap=ColorMapsMatplotlib[i],\
+ mask = Mask, edgecolor=None ,linewidth=0, antialiased=False, alpha = 0.7)
+
+ # --------------- Customisation properties -----------------------------------
+ # Retrieve the variables names
+ variableName = Field
+ variableLatexName = Field
+ variableLatexUnit = Field
+
+ # Customisation of labels and latex conversion
+ rc('font',size = 12)
+ rc('font',family = 'serif')
+ rc('axes',labelsize =14)
+
+ # Customise the positioning of the labels
+ [t.set_va('center') for t in ax.get_yticklabels()]
+ [t.set_ha('right') for t in ax.get_yticklabels()]
+ [t.set_va('center') for t in ax.get_xticklabels()]
+ [t.set_ha('left') for t in ax.get_xticklabels()]
+ [t.set_va('center') for t in ax.get_zticklabels()]
+ [t.set_ha('right') for t in ax.get_zticklabels()]
+
+ # Set the grid and the pane colors
+ ax.grid(False)
+ ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+ ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+ ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+
+ fig.suptitle(r'Solution for the variable {0!s} at the time ${1:5.3f}$'.format(variableLatexName, Solution.t))
+ ax.set_xlabel(r'x [{0!s}]'.format(Problem.GoverningEquations.XYLatexUnits[0]))
+ ax.set_ylabel(r'y [{0!s}]'.format(Problem.GoverningEquations.XYLatexUnits[1]))
+ ax.set_zlabel(r'{0!s} [{1!s}]'.format(variableLatexName, variableLatexUnit))
+
+ # --------------- Export the figure at the right location -------------------------
+ FolderName = os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(Mesh.MeshResolution).replace(".","_"), "3DPlots")
+ try: os.makedirs(FolderName)
+ except: pass
+
+ FileName = "Mplt_3DSolEoS_OnSubDomainsTri_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
+ plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
+ plt.close()
+
+
# ==============================================================================
# Soltuion tools usable for any element type (less precise for contour plots)
# ==============================================================================
@@ -1092,7 +1253,7 @@ def PlotVariableOnWholeDomain(Problem, Mesh, Solution, VariableId, ParametersId)
try: os.makedirs(FolderName)
except: pass
- FileName = "Mplt_2DSol_OnWholeDomain_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
+ FileName = "Mplt_3DSol_OnWholeDomain_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
plt.close()
@@ -1180,6 +1341,180 @@ def PlotVariableOneEachSubDomain(Problem, Mesh, Solution, VariableId, Parameters
plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
plt.close()
+#### Plot the variable coming from the EoS independently of the subdomain ####
+def PlotVariableEoSOnWholeDomain(Problem, Mesh, Solution, Field, ParametersId):
+ """ Plot one variable on the whole subdomain
+
+ Args:
+ Problem (Problem): the considered problem and its properties
+ Mesh (MeshStructure): the considered mesh
+ Solution (Solution): the solution of the considered problem
+ Field (string): the name to give to the Eos to retrieeve the desired variable to be plot
+ ParametersId (string): name of the export folder (TestCase_TestSettings)
+
+ Returns:
+ None: only exporting the plot
+
+ --------------------------------------------------------------------------------------------"""
+
+ # --------------- Preparing the data to plot -----------------------------------
+ # Extracting the data of interest
+ X = Mesh.CartesianPoints[:,0]
+ Y = Mesh.CartesianPoints[:,1]
+ Z = Solution.RSol[VariableId,:]
+
+ # Computing the value from the EoS
+ Zbuf = np.zeros((Z.shape[1]))
+
+ for i in set(map(int,Solution.RFluidFlag)):
+ ind = np.where(Solution.RFluidFlag==i)[0]
+ Zbuf[ind] = Problem.FluidEOS[i].EoS(Z[:,ind], [Problem.FluidProp[i]]*len(ind), Field)
+
+
+ # --------------- Creating the figure ------------------------------------------
+ fig = plt.figure()
+ ax = fig.gca(projection='3d')
+
+ # Interpolating the data on the visualisation grid
+ xx, yy, zz = GetMeshgridFom3VectorsMaskX(Mesh, X, Y, Zbuf)
+
+ # --------------- Creating the figure ------------------------------------------
+ fig = plt.figure()
+ ax = fig.gca(projection='3d')
+ ax.plot_surface(xx, yy, zz, alpha=0.7, cmap=ColorMapsMatplotlib[0])
+
+ # --------------- Customisation properties -----------------------------------
+ # Retrieve the variables names
+ variableName = Field
+ variableLatexName = Field
+ variableLatexUnit = Field
+
+ # Customisation of labels and latex conversion
+ rc('font',size = 12)
+ rc('font',family = 'serif')
+ rc('axes',labelsize = 14)
+
+ # Customise the positioning of the labels
+ [t.set_va('center') for t in ax.get_yticklabels()]
+ [t.set_ha('right') for t in ax.get_yticklabels()]
+ [t.set_va('center') for t in ax.get_xticklabels()]
+ [t.set_ha('left') for t in ax.get_xticklabels()]
+ [t.set_va('center') for t in ax.get_zticklabels()]
+ [t.set_ha('right') for t in ax.get_zticklabels()]
+
+ # Set the grid and the pane colors
+ ax.grid(False)
+ ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+ ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+ ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+
+ fig.suptitle(r'Solution for the variable {0!s} at the time ${1:5.3f}$'.format(variableLatexName, Solution.t))
+ ax.set_xlabel(r'x [{0!s}]'.format(Problem.GoverningEquations.XYLatexUnits[0]))
+ ax.set_ylabel(r'y [{0!s}]'.format(Problem.GoverningEquations.XYLatexUnits[1]))
+ ax.set_zlabel(r'{0!s} [{1!s}]'.format(variableLatexName, variableLatexUnit))
+
+
+ # --------------- Export the figure at the right location -------------------------
+ FolderName = os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(Mesh.MeshResolution).replace(".","_"), "3DPlots")
+ try: os.makedirs(FolderName)
+ except: pass
+
+ FileName = "Mplt_3DSolEoS_OnWholeDomain_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
+ plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
+ plt.close()
+
+#### Plot one variable coming from the EoS with a color scale depending on the subdomain ####
+def PlotVariableEoSOneEachSubDomain(Problem, Mesh, Solution, Field, ParametersId):
+ """ Plot one variable with a color scale depending of the subdomain
+
+ Args:
+ Problem (Problem): the considered problem and its properties
+ Mesh (MeshStructure): the considered mesh
+ Solution (Solution): the solution of the considered problem
+ Field (string): the name to be passed to the EoS to retrieve the desired variable to be plot
+ ParametersId (string): name of the export folder (TestCase_TestSettings)
+
+ Returns:
+ None: only exporting the plot
+
+ --------------------------------------------------------------------------------------------"""
+
+ # --------------- Preparing the data to plot -----------------------------------
+ # Extracting the data of interest
+ X = Mesh.CartesianPoints[:,0]
+ Y = Mesh.CartesianPoints[:,1]
+ Z = Solution.RSol[VariableId,:]
+ LS = Solution.RLSValues[:,:]
+
+ # Computing the value from the EoS
+ Zbuf = np.zeros((Z.shape[1]))
+
+ for i in set(Solution.RFluidFlag):
+ ind = np.where(Solution.RFluidFlag==i)[0]
+ Zbuf[ind] = Problem.FluidEOS[i].EoS(Z[:,ind], [Problem.FluidProp[i]]*len(ind), Field)
+
+
+ # --------------- Creating the figure ------------------------------------------
+ fig = plt.figure()
+ ax = fig.gca(projection='3d')
+
+ # Interpolating the data on the visualisation grid
+ _, _, zz = GetMeshgridFom3VectorsMaskX(Mesh, X, Y, Zbuf)
+
+ # Span each fluid present in the domain and plot one variable surface per area and the countour of each area
+ # according to the interpolated value of the FluidSelector
+ Fluids = set(map(int, Solution.RFluidFlag))
+ for i in Fluids:
+ # Interpolating the fluid selector on the visualisation grid
+ xx, yy, lzz = GetMeshgridFom3VectorsMaskX(Mesh, X, Y, LS[i,:])
+
+ # Masking the points that lie outisde the fluid's area
+ ind = np.where(lzz<=0)
+ xx[ind] = None
+ yy[ind] = None
+
+ # Plotting the surface
+ ax.plot_surface(xx, yy, zz, alpha=0.7, cmap=ColorMapsMatplotlib[i])
+
+ # --------------- Customisation properties -----------------------------------
+ # Retrieve the variables names
+ variableName = Field
+ variableLatexName = Field
+ variableLatexUnit = Field
+
+ # Customisation of labels and latex conversion
+ rc('font',size = 12)
+ rc('font',family ='serif')
+ rc('axes',labelsize = 14)
+
+ # Customise the positioning of the labels
+ [t.set_va('center') for t in ax.get_yticklabels()]
+ [t.set_ha('right') for t in ax.get_yticklabels()]
+ [t.set_va('center') for t in ax.get_xticklabels()]
+ [t.set_ha('left') for t in ax.get_xticklabels()]
+ [t.set_va('center') for t in ax.get_zticklabels()]
+ [t.set_ha('right') for t in ax.get_zticklabels()]
+
+ # Set the grid and the pane colors
+ ax.grid(False)
+ ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+ ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+ ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.2))
+
+ fig.suptitle(r'Solution for the variable {0!s} at the time ${1:5.3f}$'.format(variableLatexName, Solution.t))
+ ax.set_xlabel(r'x [{0!s}]'.format(Problem.GoverningEquations.XYLatexUnits[0]))
+ ax.set_ylabel(r'y [{0!s}]'.format(Problem.GoverningEquations.XYLatexUnits[1]))
+ ax.set_zlabel(r'{0!s} [{1!s}]'.format(variableLatexName, variableLatexUnit))
+
+ # --------------- Export the figure at the right location -------------------------
+ FolderName = os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(Mesh.MeshResolution).replace(".","_"), "3DPlots")
+ try: os.makedirs(FolderName)
+ except: pass
+
+ FileName = "Mplt_3DSolEoS_OnSubDomains_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
+ plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
+ plt.close()
+
# ==============================================================================
# Checkup tools
diff --git a/SourceCode/PostProcessing/Plots_Schlieren.py b/SourceCode/PostProcessing/Plots_Schlieren.py
index 619c6f4..b41bd2b 100644
--- a/SourceCode/PostProcessing/Plots_Schlieren.py
+++ b/SourceCode/PostProcessing/Plots_Schlieren.py
@@ -183,7 +183,7 @@ def GetMeshgridFom3VectorsMaskZ(Mesh,X,Y,Z):
# Functions
# ==============================================================================
-#### Debug plotting routine that plots fluid's flag at each Dof ####
+#### Computing Schlieren images from the conservative variables ####
def Compute2DSchlierenImages(Problem, Mesh, Solution, Field, ParametersId):
""" Plots the numerical Schlieren images from the density field
@@ -282,4 +282,101 @@ def Compute2DSchlierenImages(Problem, Mesh, Solution, Field, ParametersId):
FileName = "Mplt_Schlieren_{1!s}_t{0:5.3f}".format(Solution.t, Field).replace(".","_")
plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
plt.close()
+
+#### Computing Schlieren images from the EOS linked pressure variable ####
+def Compute2DSchlierenImagesEos(Problem, Mesh, Solution, Field, ParametersId):
+ """ Plots the numerical Schlieren images from the value given out of the EoS (usually pressure)
+
+ Args:
+ Problem (Problem): the considered problem and its properties
+ Mesh (MeshStructure): the considered mesh
+ Solution (Solution): the solution of the considered problem
+ Field (string): the name of the variable to derive from the EOs
+ ParametersId (string): name of the export folder (TestCase_TestSettings)
+
+ Returns:
+ None: Only exporting the plot
+
+ |
+
+ .. warning::
+
+ Only valid for Euler equations, the other test case will be compatible on call but the results would be meaningless
+
+
+ """
+
+ # ----------------------- Defining the static color plotting vectors -------------------------------------------
+
+ # Defining the plotting representation quantities for several fluids
+ K = [60, 12, 40] # Smoothing factor in the shading function
+ start = [0, 230/255, 100/255] # Value indexing of the plotting (in the B scale)
+
+ # Creating the new colormaps for all fluids
+ cmap = []
+ for i in range(3):
+ vals = np.ones((256, 4))
+ vals[:, 0] = np.linspace(0, 200/255, 256)
+ vals[:, 1] = np.linspace(0, 230/255, 256)
+ vals[:, 2] = np.linspace(start[i], 1, 256)
+ cmap += [ListedColormap(vals)]
+
+ # ----------------------- Computing the Schlieren numerical images -----------------------------------------------
+ # Getting the desired values and there spatial location
+ X = Mesh.CartesianPoints[:,0]
+ Y = Mesh.CartesianPoints[:,1]
+ Z = Solution.RSol[:,:]
+ Zbuf = np.zeros((Z.shape[1]))
+
+ for i in set(map(int,Solution.RFluidFlag)):
+ ind = np.where(Solution.RFluidFlag==i)[0]
+ Zbuf[ind] = Problem.FluidEOS[i].EoS(Z[:,ind], [Problem.FluidProp[i]]*len(ind), Field)
+
+ # Mapping them to a cartesian grid
+ xx,yy,zz = GetMeshgridFom3VectorsMaskZ(Mesh, X, Y, Zbuf)
+
+ # Get the fluid's spotter markers on the interpolation nodes
+ lll = np.zeros(zz.shape+(Solution.RLSValues.shape[0],))
+ for i in range(Solution.RLSValues.shape[0]):
+ Z = Solution.RLSValues[i,:]
+ xx,yy,lll[:,:,i] = GetMeshgridFom3VectorsMaskZ(Mesh, X, Y, Z)
+ Fluid = np.argmax(lll,axis=2)
+
+ # Computing the shading functions for each fluid
+ plt.figure()
+ for i in range(Solution.RLSValues.shape[0]):
+ # Spotting the fluid's location
+ ind = np.where(Fluid==i)
+ ind2 = np.where(Fluid!=i)
+
+ # Computing the gradient field
+ Gradients = np.gradient(zz, yy[:,0], xx[0,:])
+
+ # Computing the gradient's magnitude
+ MDField = np.sqrt(Gradients[0]**2+Gradients[1]**2)
+ MmMDField = np.nanmax(MDField)
+
+ # Determining the shaded function
+ val = np.exp(-(K[i]*(Fluid[ind[0], ind[1]]+1)*MDField[ind[0], ind[1]])/MmMDField)
+
+ # Defining the plotting quantity
+ Phi = np.zeros(MDField.shape)
+ Phi[ind[0], ind[1]] = val
+ Phi[ind2[0], ind2[1]] = None
+
+ # Plot the quantity
+ plt.pcolor(xx,yy,Phi, cmap=cmap[i])
+
+
+ # --------------- Export the figure at the right location -------------------------
+ # Creating the export environment
+ FolderName = os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(Mesh.MeshResolution).replace(".","_"), "2DPlots")
+ try: os.makedirs(FolderName)
+ except: pass
+
+ # Set the figure properties and shows it
+ FileName = "Mplt_Schlieren_Eos_{1!s}_t{0:5.3f}".format(Solution.t, Field).replace(".","_")
+ plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
+ plt.close()
+
\ No newline at end of file
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py b/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py
index be866ea..729a4b0 100644
--- a/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py
+++ b/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py
@@ -1,6 +1,6 @@
"""
All the necessary routines that implement
-the 2D Euler Equations (given in conservative variables) are defined in
+the 2D Euler Equations (given in conservative variables) are defined in
:code:`GoverningEquations/EulerEquations.py`
"""
@@ -27,11 +27,11 @@ class Equations():
schemes and evolve the solution according to the Euler Equation
|
-
+
.. rubric:: Fields
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
It contains the following initialisation parameters (later available as well as attributes)
- |bs| **FluidProp** *(list of FluidModel)* -- the list of fluid property instances (NbFluids)
@@ -53,12 +53,12 @@ class Equations():
Args:
FluidProp (list of FluidModel): the list of fluid property instances (NbFluids)
EOs (list of EOS): the list of Equation of state instances (NbFluids)
-
+
|
-
+
.. rubric:: Methods
"""
-
+
# Defines the (conservative) variables number and names, in a latex typography
self.NbVariables = 4
self.VariablesNames = ["rho", "rhoUx", "rhoUy", "E"]
@@ -67,32 +67,122 @@ class Equations():
self.VariablesLatexNames = [r'$\rho$', r'$\rho U_x$', r'$\rho U_y$', r'$E$']
self.VariablesLatexUnits = [r'$kg/m^3$', r'$m^2/s$', r'$m^2/s$', r'$J$']
self.XYLatexUnits = [r'$m$', r'$m$']
-
+
# Register the shortcut of the fluid's properties
self.FluidProp = FluidProp
self.EOs = EOs
+ #### Retrieve the primary variables from the conservative ones ####
+ def ConservativeToPrimary(self, Var, FluidIndex, *args):
+ """Converts the conservative variables to the primary ones
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+
+ Returns:
+ Var (numpy array): (NbVars x NbGivenPoints) the corresponding primary variables
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # -------- Initialisations -----------------------------------------------------------------------
+ # Getting the number of furnished points
+ NbPoints = np.shape(Var)[1]
+ PrimVars = np.zeros(np.shape(Var))
+
+ # --------- Locating and extracting the right fluid's properties for each location ---------------
+ # Initialising the pressure and gamma variables
+ p = np.zeros(len(Var[0,:]))
+ # List of fluids present in the list to be computed
+ Fluids = list(set(FluidIndex))
+
+ # Construct the fluidproperties depending o the right fluid selector's value and evaluating the
+ # EOS accordingly
+ for Fluid in Fluids:
+ Index = np.where(Fluid==np.array(FluidIndex))[0]
+ p[Index] = self.EOs[Fluid].EoS(Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)), "p")
+
+ # ------- Constructing the primary variables ----------------------------------------------------
+ # Computing the primary variables
+ PrimVars[0,:] = Var[0,:]
+ PrimVars[1,:] = Var[1,:]/Var[0,:]
+ PrimVars[2,:] = Var[2,:]/Var[0,:]
+ PrimVars[3,:] = p[:]
+
+ # Returning the Result
+ return(PrimVars)
+
+ #### Retrieve the conservative variables from the primary ones ####
+ def PrimaryToConservative(self, PrimVar, FluidIndex, *arrgs):
+ """Converts the primary variables to the conservative ones
+
+ Args:
+ Var (2D numpy array): the primary variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+
+ Returns:
+ Var (numpy array): (NbVars x NbGivenPoints) the corresponding conservative variables
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # -------- Initialisations -----------------------------------------------------------------------
+ # Getting the number of furnished points
+ NbPoints = np.shape(PrimVar)[1]
+ Vars = np.zeros(np.shape(PrimVar))
+
+ # Semiconverting the variables to pass to the EOS
+ PrimVar[1,:] = PrimVar[0,:]*PrimVar[1,:]
+ PrimVar[2,:] = PrimVar[0,:]*PrimVar[2,:]
+
+ # --------- Locating and extracting the right fluid's properties for each location ---------------
+ # Initialising the pressure and gamma variables
+ e = np.zeros(len(PrimVar[0,:]))
+ # List of fluids present in the list to be computed
+ Fluids = list(set(FluidIndex))
+
+ # Construct the fluidproperties depending o the right fluid selector's value and evaluating the
+ # EOS accordingly
+ for Fluid in Fluids:
+ Index = np.where(Fluid==np.array(FluidIndex))[0]
+ e[Index] = self.EOs[Fluid].EoS(PrimVar[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)), "e")
+
+ # ------- Constructing the primary variables ----------------------------------------------------
+ # Computing the primary variables
+ Vars[0,:] = PrimVar[0,:]
+ Vars[1,:] = PrimVar[1,:] # Directly from the previous semiconversion
+ Vars[2,:] = PrimVar[2,:] # Directly from the previous semiconversion
+ Vars[3,:] = PrimVar[0,:]*e[:]+0.5*PrimVar[0,:]*(PrimVar[1,:]**2+PrimVar[2,:]**2)
+
+ # Returning the Result
+ return(Vars)
+
#### Flux function of the equations set ####
def Flux(self, Var, *args):
"""Flux function of the (conservative) EulerEquations
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Args:
+
+ Args:
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
-
+
.. rubric:: Optional argument
-
+
- |bs| **x** *(float numpy array)*, -- (optional, not impacting here) x-y coordinates of the given points (NbGivenPoints x 2)
- |bs| **FluidIndex** *(integer numpy array)* -- the fluid present at each given point (NbGivenPoints)
-
+
Returns:
None: fills direclty the Parameters data structure
|
-
+
.. note::
-
+
- | This function is Vectorised
- | Fluid index is the fluid index in the initially given list, not the fluid type
- | FluidProp is the list of all the fluid properties for each fluid (given the list, not the type)
@@ -100,9 +190,9 @@ class Equations():
"""
# --------- Locating and extracting the right fluid's properties for each location ---------------
- # Gettings the arguments
+ # Gettings the arguments
FluidIndex = args[1]
-
+
# Initialising the pressure variable
p = np.zeros(len(Var[0,:]))
# List of fluids present in the list to be computed
@@ -113,7 +203,7 @@ class Equations():
for Fluid in Fluids:
Index = np.where(Fluid==np.array(FluidIndex))[0]
p[Index] = self.EOs[Fluid].EoS(Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)), "p")
-
+
# --------- Computing the Euler flux directly in conservative variables ------------------------
Fx = np.array([Var[1,:], # RhoU
Var[1,:]**2/Var[0,:]+p, # RhoU**2+p
@@ -127,71 +217,71 @@ class Equations():
# Returning the flux
Flux = np.array([Fx, Fy])
-
+
return(Flux)
#### Function that extracts the propagation speed from the variables to interface for the LevelSet ####
def GetUV(self, Var, *args):
""" Function giving back the x-y velocity at the considered point, given the variational
values given for the the conservative Euler equation's variables.
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Args:
+
+ Args:
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
x (float numpy array): (generally optional, required for this problem), x-y coordinates of the given points (NbGivenPoints x 2)
.. rubric:: Optional argument
-
+
- |bs| **FluidIndex** *(integer numpy array)* -- the fluid present at each given point (NbGivenPoints)
-
+
Returns:
UV (float numpy array) -- the velocities values at the points (2 x NbGivenPoints)
.. note::
-
+
- | This function is Vectorised
- | *args is there only for compatibility reason at call time
"""
# Computing the velocity values from the conservative Euleur equations
UV = Var[1:3, :]/Var[0,:]
-
+
# Returning the value
return(UV)
#### Jacobian of the equations set ####
def Jacobian(self, Var, x, FluidIndex, *args):
"""Computes the Jacobian of the flux for the Euler equations
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Args:
+
+ Args:
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
x (float numpy array): (generally optional, required for this problem), x-y coordinates of the given points (NbGivenPoints x 2)
FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
-
- Returns:
- J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
-
+
+ Returns:
+ J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
+
|
-
+
.. note::
-
+
- | For each flux fi = (fi1,..., fin), the returned Jacobian reads:
-
+
| |bs| J[:,:] = [dfi1/dx1, ....., dfi1/dxn
| |bs| |bs| |bs| |bs| |bs| |bs| ....
| |bs| |bs| |bs| |bs| |bs| |bs| df in/dx1, ....., dfin/dxn]
- | *args is there only for compatibility reason at call time
-
+
"""
-
+
# ------ Initialising the Jacobian structure ---------------------------
# Getting the number of furnished points
NbPoints = np.shape(Var)[1]
J = np.zeros((4, 4, 2, NbPoints))
-
+
# --------- Locating and extracting the right fluid's properties for each location ---------------
# Initialising the pressure and gamma variables
p = np.zeros(len(Var[0,:]))
@@ -205,7 +295,7 @@ class Equations():
Index = np.where(Fluid==np.array(FluidIndex))[0]
p[Index] = self.EOs[Fluid].EoS(Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)), "p")
gmm[Index] = [self.FluidProp[Fluid].gamma]*len(Index)
-
+
# --------- Computing the Euler Jacobian directly from conservative variables ------------------------
# Getting back the values of speed from the conservative variables
u = Var[1,:]/Var[0,:] # get the x-speed
@@ -226,43 +316,273 @@ class Equations():
J[2,:,1,:] = np.array([-v**2 + 0.5*(gmm-1.)*uv2, u*(1.-gmm), v*(3.-gmm), gmm-1. ])
J[3,:,1,:] = np.array([v*(0.5*(gmm-1.)*uv2 - H), u*v*(1.-gmm), H-(gmm-1.)*v**2, gmm*v ])
- # Returning the Jacobian
+ # Returning the Jacobian
return(J)
+ #### Eigen values ####
+ def EigenValues(self, Var, FluidIndex, n, x, *args):
+ """Computes the eigenvalues associated to the flux.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+ x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Returns:
+ lbd (numpy array): (NbGivenPoints x 4) the four eigenvalues at each given point
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ---------- Ugly always working way ------------------------------
+
+ # --- Locating and extracting the right fluid's properties for each location
+ # Initialising the pressure and gamma variables
+ p = np.zeros(len(Var[0,:]))
+ c = np.zeros(len(Var[0,:]))
+
+ # List of fluids present in the list to be computed
+ Fluids = list(set(FluidIndex))
+
+ # Getting the fluids properties depending on the fluid selector's value and evaluate the
+ # EOS accordingly
+ for Fluid in Fluids:
+ # Extracting the points occupied by the selected fluid
+ Index = np.where(Fluid==np.array(FluidIndex))[0]
+
+ # Getting the EoS dependent quantitites
+ p[Index] = self.EOs[Fluid].EoS(Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)), "p")
+
+ # Getting the sound speed and velocity against normals
+ c[Index] = np.sqrt(self.FluidProp[Fluid].gamma*np.abs(p[Index])/Var[0,Index])
+
+ # --- Computing the spectral radius taking the sound speed into consideration
+ # Getting the velocity against normals
+ un = np.sum(Var[1:3,:].T*n[:], axis=1)/Var[0,:]
+ mod_n = np.sqrt(np.sum(n**2, axis=1))
+
+ # Getting the actual eigenvalues
+ lbd = np.zeros((len(Var[0,:]), 4))
+ lbd[:,0] = un - c*mod_n
+ lbd[:,1] = un
+ lbd[:,2] = un + c*mod_n
+ lbd[:,3] = un
+
+ return(lbd)
+
+ #### Right eigen vectors ####
+ def RightEigenVectors(self, Var, FluidIndex, n, x, *args):
+ """Computes the right eigenvectors associated to the eigenvalues.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+ x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Returns:
+ reg (numpy array): (NbVars x MbVars x NbGivenPoints) the matrix of eigenvectors for each given point
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ---------------- Initialisations -------------------------------------
+ p = np.zeros(len(Var[0,:]))
+ c = np.zeros(len(Var[0,:]))
+ dEoSRho = np.zeros(len(Var[0,:]))
+ dEoSE = np.zeros(len(Var[0,:]))
+
+ # ---------------- Precomputations -------------------------------------
+ # --- Variable-dependent quantities
+ # Normalising the normals
+ nn = n/(np.sqrt(np.sum(n**2, axis=1)).T)
+ # Retrieving the velocities
+ uv = Var[1:3]/Var[0]
+ # Retrieving the normal component of the velocity
+ un = np.sum(uv*n[:].T, axis=0)
+ # Retrieving the kinetic energy
+ ke = 0.5*np.sum(uv**2, axis=0)
+
+ # --- EOS-dependent quantities
+ # List of fluids present in the list to be computed
+ Fluids = list(set(FluidIndex))
+
+ # Getting the fluids properties depending on the fluid selector's value and evaluate the EOS accordingly
+ for Fluid in Fluids:
+ # Extracting the points occupied by the selected fluid
+ Index = np.where(Fluid==np.array(FluidIndex))[0]
+
+ # Getting the EoS dependent quantitites
+ p[Index] = self.EOs[Fluid].EoS (Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)), "p")
+ dEoSRho[Index] = self.EOs[Fluid].dEoSRho(Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)))
+ dEoSE[Index] = self.EOs[Fluid].dEoSE (Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)))
+
+ # Getting the sound speed and velocity against normals
+ c[Index] = np.sqrt(self.FluidProp[Fluid].gamma*np.abs(p[Index])/Var[0,Index])
+
+ # Computing the total Enthalpy
+ H = (Var[3,:]+p)/Var[0,:]
+
+ # -------------------- Computing the eigenvectors ----------------------
+ # Computing the tangents at each point
+ nt = np.zeros(nn.T.shape)
+ nt[0,:] = nn[:,1]
+ nt[1,:] = -nn[:,0]
+
+ # Computing the velocities against the tangent at each point and the quantities involved in the
+ uvnp = np.sum(uv*nt, axis=0)
+ xi = 2*ke - dEoSRho/dEoSE
+
+ # Creating right eigenvalues
+ REV = np.zeros((Var.shape[0], Var.shape[0], Var.shape[1]))
+ #print(c*nn.T)
+ REV[0, 0, :] = 1
+ REV[1:3, 0, :] = uv - c*nn.T
+ REV[3, 0, :] = H - c*un
+
+ REV[0, 1, :] = 1
+ REV[1:3, 1, :] = uv
+ REV[3, 1, :] = xi
+
+ REV[0, 2, :] = 1
+ REV[1:3, 2, :] = uv + c*nn.T
+ REV[3, 2, :] = H + c*un
+
+ REV[0, 3, :] = 0
+ REV[1:3, 3, :] = nt
+ REV[3, 3, :] = uvnp
+
+ # Returning the values
+ return(REV)
+
+ #### Right eigen vectors ####
+ def LeftEigenVectors(self, Var, FluidIndex, n, x, *args):
+ """Computes the left eigenvectors associated to the eigenvalues.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+ x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Returns:
+ reg (numpy array): (NbVars x MbVars x NbGivenPoints) the matrix of eigenvectors for each given point
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ---------------- Initialisations -------------------------------------
+ p = np.zeros(len(Var[0,:]))
+ c = np.zeros(len(Var[0,:]))
+ dEoSRho = np.zeros(len(Var[0,:]))
+ dEoSE = np.zeros(len(Var[0,:]))
+
+ # ---------------- Precomputations -------------------------------------
+ # --- Variable-dependent quantities
+ # Normalising the normals
+ nn = n/(np.sqrt(np.sum(n**2, axis=1)).T)
+ # Retrieving the velocities
+ uv = Var[1:3]/Var[0]
+ # Retrieving the normal component of the velocity
+ un = np.sum(uv*n.T, axis=0)
+ # Retrieving the kinetic energy
+ ke = 0.5*np.sum(uv**2, axis=0)
+
+ # --- EOS-dependent quantities
+ # List of fluids present in the list to be computed
+ Fluids = list(set(FluidIndex))
+
+ # Getting the fluids properties depending on the fluid selector's value and evaluate the EOS accordingly
+ for Fluid in Fluids:
+ # Extracting the points occupied by the selected fluid
+ Index = np.where(Fluid==np.array(FluidIndex))[0]
+
+ # Getting the EoS dependent quantitites
+ p[Index] = self.EOs[Fluid].EoS (Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)), "p")
+ dEoSRho[Index] = self.EOs[Fluid].dEoSRho(Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)))
+ dEoSE[Index] = self.EOs[Fluid].dEoSE (Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)))
+
+ # Getting the sound speed and velocity against normals
+ c[Index] = np.sqrt(self.FluidProp[Fluid].gamma*np.abs(p[Index])/Var[0,Index])
+
+ # Computing the total Enthalpy
+ H = (Var[3,:]+p)/Var[0,:]
+
+ # -------------------- Computing the eigenvectors ----------------------
+ # Computing the tangents at each point
+ nt = np.zeros(nn.T.shape)
+ nt[0,:] = nn[:,1]
+ nt[1,:] = -nn[:,0]
+
+ # Computing the velocities against the tangent at each point and the quantities involved in the
+ uvnp = np.sum(uv*nt, axis=0)
+ xi = 2*ke - dEoSRho/dEoSE
+
+ # Creating left eigenvalues
+ LEV = np.zeros((Var.shape[0], Var.shape[0], Var.shape[1]))
+
+ LEV[0, 0, :] = c*un + dEoSRho
+ LEV[1:3, 0, :] = -c*nn.T - dEoSE*uv
+ LEV[3, 0, :] = dEoSE
+
+ LEV[0, 1, :] = -2*(dEoSRho - c**2)
+ LEV[1:3, 1, :] = 2*dEoSE*uv
+ LEV[3, 1, :] = -2*dEoSE
+
+ LEV[0, 2, :] = -c*un + dEoSRho
+ LEV[1:3, 2, :] = c*nn.T - dEoSE*uv
+ LEV[3, 2, :] = dEoSE
+
+ LEV[0, 3, :] = -2*uvnp*c**2
+ LEV[1:3, 3, :] = 2*nt*c**2
+ LEV[3, 3, :] = 0
+
+ LEV = LEV/(2*(c**2))
+
+ # Returning the values
+ return(LEV)
+
#### Spectral radius of the equations set ####
def SpectralRadius(self, Var, FluidIndex, n, x, *args):
"""Computes the spectral radius associated to the flux.
-
+
Args:
Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
- Returns:
+ Returns:
Lmb (numpy array): the spectral radius computed at each given point
-
+
.. note::
-
+
- | *args is there only for compatibility reason at call time
"""
-
+
# ---------- Cheap way that may fail in case of zero speed ------------------------------
"""
# Retrieving the Jacobian
J = self.Jacobian(Var, x, FluidIndex)
-
+
# Computing the spectral radius at each given point, shortcuted since here we are scalar
Lmb = np.max(np.abs([np.sum(J[i,i,:,:]*n.T, axis = 0) for i in range(4)]), axis=0)
"""
-
+
# ---------- Ugly always working way ------------------------------
-
- # --- Locating and extracting the right fluid's properties for each location
+
+ # --- Locating and extracting the right fluid's properties for each location
# Initialising the pressure and gamma variables
p = np.zeros(len(Var[0,:]))
gmm = np.zeros(len(Var[0,:]))
-
+
# List of fluids present in the list to be computed
Fluids = list(set(FluidIndex))
@@ -270,23 +590,23 @@ class Equations():
# EOS accordingly
for Fluid in Fluids:
Index = np.where(Fluid==np.array(FluidIndex))[0]
- p[Index] = self.EOs[Fluid].EoS(Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)), "e")
+ p[Index] = self.EOs[Fluid].EoS(Var[:,Index], np.array([self.FluidProp[Fluid]]*len(Index)), "p")
gmm[Index] = [self.FluidProp[Fluid].gamma]*len(Index)
-
- # --- Computing the spectral radius taking the sound speed into consideration
- # Getting the sound speed
+
+ # --- Computing the spectral radius taking the sound speed into consideration
+ # Getting the sound speed
c = np.sqrt(gmm*np.abs(p)/Var[0,:])
vi = np.sqrt(np.sum(Var[1:3,:]**2, axis=0)/(Var[0,:]**2))
-
+
Lmb = np.sqrt(np.sum(n[:,:]**2, axis=1))*(vi+c)
# Returning the spectral radius at each given point
return(Lmb)
-
+
#### Roe matrix associated to the Euler problem ####
def RoeMatrix(self, Var, n, FluidIndex, *args):
"""A Roe matrix corresponding to the Euler Equations, accross one given edge
-
+
Args:
n (float numpy array): the outward normal of the considered element (nx, ny)
Var (float numpy array): the value of the variables (NbVariables x 2). Var[:,1] average value of the considered element, Var[:,2] average value of the neigborhing element.
@@ -297,12 +617,12 @@ class Equations():
RoeMatrix (float numpy array): a Roe matrix of the system
|
-
- .. note::
+
+ .. note::
- This function is NOT vectorised
- Fluid index is the fluid index in the initially given list, not the fluid type
- FluidProp is the list of all the fluid properties for each fluid (given the list, not the type)
-
+
--------------------------------------------------------------------------------------------------------------------------------"""
# --------- Locating and extracting the right fluid's properties for each location ---------------
@@ -354,10 +674,10 @@ class Equations():
#### Class containing the methods giving back the pressure ####
class EOS():
"""Class furnishing the methods giving the Equation of State that can be
- used in combination with the :code:`SimpleFluid` definition.
+ used in combination with the :code:`SimpleFluid` definition.
Upon creation, fills the field :code:`EOS` with the function callback
corresponding to the wished initialisation function and the field :code:`EOSName` with
- the associated label. The implemented EoS are linked as follows.
+ the associated label. The implemented EoS are linked as follows.
See their respective documentation for more information.
1. | Stiff fluid
@@ -370,19 +690,23 @@ class EOS():
"""
Args:
Id (integer): the index corresponding to the equation of fluid the user wants when considering the governing equations given in this module and the associated fluid.
-
+
|
-
+
.. rubric:: Methods
"""
# Mapping towards the right function
if Id==1:
self.EoS = self.Stiff
+ self.dEoSRho = self.dEoSRhoStiff
+ self.dEoSE = self.dEoSEStiff
self.EoSName = "Stiff"
-
+
elif Id==2:
self.EoS = self.Ideal
+ self.dEoSRho = self.dEoSRhoIdeal
+ self.dEoSE = self.dEoSEIdeal
self.EoSName = "Ideal"
else: raise NotImplementedError("Error in the selection of the Equation of State\n\
@@ -391,12 +715,12 @@ class EOS():
#### EOS for stiff fluids #####
def Stiff(self, Var, FluidProp, Type):
"""Equation of State for stiff fluids, giving back either the internal energy or pressure.
-
- Args:
+
+ Args:
Type (string): desired output: "e": internal energy, "p": pressure
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints). If type "e", the variables should be [rho, rhoU, rhoV, p].T. If type "p", the variables should be [rho, rhoU, rhoV, E].
FluidProp (FluidModel list): the list of the fluid properties associated to each given Point
-
+
Returns:
p (float numpy array): the pressure values at the given points (NbGivenPoints)
"""
@@ -415,15 +739,59 @@ class EOS():
# Returning the asked value
return(val)
+ #### EOS derivative with respect to rho variable ####
+ def dEoSRhoStiff(self, Var, FluidProp):
+ """Derivative of the Equation of State for stiff fluids along the density variable at the given points.
+
+ Args:
+ Type (string): desired output: "e": internal energy, "p": pressure
+ Var (float numpy array): the value of the conservative variables (NbVariables x NbGivenPoints).
+ FluidProp (FluidModel list): the list of the fluid properties associated to each given Point
+
+ Returns:
+ deos (float numpy array): derivative of the Equation of State for stiff fluids along the density variable (NbGivenPoints)
+ """
+
+ # Extracting the fluid properties
+ FluidSpecs = np.array([[Fluid.pinf, Fluid.gamma] for Fluid in FluidProp])
+
+ # Computing the derivative
+ deos = -0.5*(FluidSpecs[:,1].T-1)*np.sum(Var[1:3,:]**2, axis=0)/(Var[0,:]**2) # check minus
+
+ # Returning the value
+ return(deos)
+
+ #### EOS derivative with respect to rho variable ####
+ def dEoSEStiff(self, Var, FluidProp):
+ """Derivative of the Equation of State for stiff fluids along the Energy variable at the given points.
+
+ Args:
+ Type (string): desired output: "e": internal energy, "p": pressure
+ Var (float numpy array): the value of the conservative variables (NbVariables x NbGivenPoints).
+ FluidProp (FluidModel list): the list of the fluid properties associated to each given Point
+
+ Returns:
+ deos (float numpy array): derivative of the Equation of State for stiff fluids along the Energy variable (NbGivenPoints)
+ """
+
+ # Extracting the fluid properties
+ FluidSpecs = np.array([[Fluid.pinf, Fluid.gamma] for Fluid in FluidProp])
+
+ # Computing the derivative
+ deos = FluidSpecs[:,1].T-1
+
+ # Returning the value
+ return(deos)
+
#### EOS for ideal fluids ####
def Ideal(self, Var, FluidProp, Type):
"""Equation of State for ideal fluids, giving back either the internal energy or pressure.
-
- Args:
+
+ Args:
Type (string): desired output: "e": internal energy, "p": pressure
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints). If type "e", the variables should be [rho, rhoU, rhoV, p].T. If type "p", the variables should be [rho, rhoU, rhoV, E].
FluidProp (FluidModel list): the list of the fluid properties associated to each given Point
-
+
Returns:
p (float numpy array): the pressure values at the given points (NbGivenPoints)
"""
@@ -442,6 +810,51 @@ class EOS():
# Returning the asked value
return(val)
+ #### EOS derivative with respect to rho variable ####
+ def dEoSRhoIdeal(self, Var, FluidProp):
+ """Derivative of the Equation of State for ideal fluids along the density variable at the given points.
+
+ Args:
+ Type (string): desired output: "e": internal energy, "p": pressure
+ Var (float numpy array): the value of the conservative variables (NbVariables x NbGivenPoints).
+ FluidProp (FluidModel list): the list of the fluid properties associated to each given Point
+
+ Returns:
+ deos (float numpy array): derivative of the Equation of State for stiff fluids along the density variable (NbGivenPoints)
+ """
+
+ # Extracting the fluid properties
+ FluidSpecs = np.array([[Fluid.pinf, Fluid.gamma] for Fluid in FluidProp])
+
+ # Computing the derivative
+ deos = -0.5*(FluidSpecs[:,1].T-1)*np.sum(Var[1:3,:]**2, axis=0)/(Var[0,:]**2) # check minus
+
+ # Returning the value
+ return(deos)
+
+ #### EOS derivative with respect to rho variable ####
+ def dEoSEIdeal(self, Var, FluidProp):
+ """Derivative of the Equation of State for ideal fluids along the Energy variable at the given points.
+
+ Args:
+ Type (string): desired output: "e": internal energy, "p": pressure
+ Var (float numpy array): the value of the conservative variables (NbVariables x NbGivenPoints).
+ FluidProp (FluidModel list): the list of the fluid properties associated to each given Point
+
+ Returns:
+ deos (float numpy array): derivative of the Equation of State for stiff fluids along the Energy variable (NbGivenPoints)
+ """
+
+ # Extracting the fluid properties
+ FluidSpecs = np.array([[Fluid.pinf, Fluid.gamma] for Fluid in FluidProp])
+
+ # Computing the derivative
+ deos = FluidSpecs[:,1].T-1
+
+ # Returning the value
+ return(deos)
+
+
# ==============================================================================
# Class containing all the predefined suitable initial conditions
# ==============================================================================
@@ -450,9 +863,9 @@ class InitialConditions():
to study the the Euler Equations, defined for working on a subdomain.
Access the routine computing the initial conditions through the field :code:`IC`, filled upon creation
with the function callback corresponding to the index of the wished initialisation function.
- The implemented initialisation method are linked as follows.
+ The implemented initialisation method are linked as follows.
See their respective documentation for more information.
-
+
0. | ConstantState
1. | ConstantState2
2. | StationaryVortex
@@ -467,10 +880,10 @@ class InitialConditions():
def __init__(self, Id, *params):
"""Args:
Id (integer): the index corresponding to the initialisation method the user wants when considering the governing equations given in this module.
- params (list of arguments): the (fixed arguments to pass to the selected function)
-
+ params (list of arguments): the (fixed arguments to pass to the selected function)
+
|
-
+
.. rubric:: Methods"""
# Registering locally the parameters to give to the functions
@@ -490,91 +903,91 @@ class InitialConditions():
Unknown index. Check your Problem definition. Abortion.")
#### SOD test case ####
- def SOD_Fluid1(self, PointsID, Mesh, EoS, FluidProp):
- """Initialising the solution to a constant state that matches the SOD-inner value.
+ def SOD_Fluid1(self, PointsID, Mesh, EoS, FluidProp, *args):
+ """Initialising the solution to a constant state that matches the SOD-inner value.
(The given points should all belonging to a same subdomain).
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): (optional), the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): (optional, not used here) the the properties of the fluid present where the given points are
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
-
+
|
-
+
.. note::
-
+
- There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
- For running smothly the inner circle diameter of the given subdomain should be at least 2
-
+
"""
-
+
# Initialising the point-values
Init = np.zeros((4, len(PointsID)))
-
+
# Computing the values that will be taken with the SOD bump
Init[:,:] = np.array([[1., 0., 0., 1]]*len(PointsID)).T
e = EoS.EoS(Init[:,:], [FluidProp]*len(PointsID), "e")
Init[3,:] = Init[0,:]*e
-
+
# Returning the values
return(Init)
-
+
#### SOD test case ####
- def SOD_Fluid2(self, PointsID, Mesh, EoS, FluidProp):
- """Initialising the solution to a constant state that matches the SOD-outer value.
+ def SOD_Fluid2(self, PointsID, Mesh, EoS, FluidProp, *args):
+ """Initialising the solution to a constant state that matches the SOD-outer value.
(The given points should all belonging to a same subdomain).
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): the the properties of the fluid present where the given points are
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
-
+
|
-
+
.. note::
-
+
- There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
- For running smothly the inner circle diameter of the given subdomain should be at least 2
"""
-
+
# Initialising the point-values
Init = np.zeros((4, len(PointsID)))
-
+
# Computing the values that will be taken outside the SOD bump
Init[:,:] = np.array([[0.125, 0., 0., 0.1]]*len(PointsID)).T
e = EoS.EoS(Init[:,:], [FluidProp]*len(PointsID), "e")
Init[3,:] = Init[0,:]*e
-
+
# Returning the values
return(Init)
#### Initialise the solution to a constant state ####
- def ConstantState(self, PointsID, Mesh, EoS, FluidProp):
+ def ConstantState(self, PointsID, Mesh, EoS, FluidProp, *args):
"""Initialising the solution to a constant state on the given points, all belonging to a same subdomain.
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): the the properties of the fluid present where the given points are
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
-
+
|
-
+
.. note::
-
+
There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
-
+
"""
# Initialising the point-values
@@ -596,24 +1009,24 @@ class InitialConditions():
return(Init)
#### Initialise the solution to a constant state ####
- def ConstantState2(self, PointsID, Mesh, EoS, FluidProp):
+ def ConstantState2(self, PointsID, Mesh, EoS, FluidProp, *args):
"""Initialising the solution to a constant state on the given points, all belonging to a same subdomain.
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): (optional), the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): (optional, not used here) the the properties of the fluid present where the given points are
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
-
+
|
-
+
.. note::
-
+
- There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
-
+
"""
# Initialising the point-values
@@ -631,26 +1044,26 @@ class InitialConditions():
def StationaryVortex(self, PointsID, Mesh, EoS, FluidProp, subdomain):
""" Initialising the solution to a StationaryVortex centred at the center of mass
of the subdmain. (The given points should all belonging to a same subdomain).
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): the the properties of the fluid present where the given points are
subdomain (shapely multipolyon): (optional, not used here) the shapely polygon to which the given points belong
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
-
+
|
-
+
.. note::
-
+
- There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
- For running smothly the inner circle diameter of the given subdomain should be at least 2
-
+
"""
-
+
# Initialising the point-values
Init = np.ones((4, len(PointsID)))
@@ -692,85 +1105,85 @@ class InitialConditions():
# Returning the computed values
return(Init)
-
+
#### SOD test case ####
def SOD(self, PointsID, Mesh, EoS, FluidProp, subdomain):
""" Initialising the solution to a SOD centred at the center of mass, of width 0.5
of the subdmain. (The given points should all belonging to a same subdomain).
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): the the properties of the fluid present where the given points are
subdomain (shapely multipolyon): (optional, not used here) the shapely polygon to which the given points belong
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
-
+
|
-
+
.. note::
-
+
- There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
- For running smothly the inner circle diameter of the given subdomain should be at least 2
-
+
"""
-
+
# Initialising the point-values
Init = np.zeros((4, len(PointsID)))
# Retrieving the xy values of the points and the geometrical informations from the subdomain
xy = Mesh.DofsXY[PointsID,:]
xc = subdomain.centroid.x
yc = subdomain.centroid.y
-
+
# Computing the radius where the SOD-ready bump will be placed
r = np.sqrt(((xy[:,0]-xc)**2+(xy[:,1]-yc)**2))
-
+
# Computing the values that will be taken with the SOD bump
Ind1 = np.where(r<=0.5)[0]
Init[:,Ind1] = np.array([[1., 0., 0., 1]]*len(Ind1)).T
e = EoS.EoS(Init[:,Ind1], [FluidProp]*len(Ind1), "e")
Init[3,Ind1] = Init[0,Ind1]*e
-
+
# Computing the values that will be taken outside the SOD bump
Ind2 = np.setdiff1d(range(len(PointsID)), Ind1)
Init[:,Ind2] = np.array([[0.125, 0., 0., 0.1]]*len(Ind2)).T
e = EoS.EoS(Init[:,Ind2], [FluidProp]*len(Ind2), "e")
Init[3,Ind2] = Init[0,Ind2]*e
-
+
# Returning the values
return(Init)
#### Karni fluid 1 ####
def Bubble(self, PointsID, Mesh, EoS, FluidProp, subdomain):
- """Initialising the solution to a constant state corresponding to the bubble in the Karni test, where the
+ """Initialising the solution to a constant state corresponding to the bubble in the Karni test, where the
density is considered as Rconsidered/Rair fluid.
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): the properties of the fluid present where the given points are
subdomain (shapely multipolyon): (optional, not used here) the shapely polygon to which the given points belong
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
-
+
|
-
+
.. note::
-
+
There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
"""
-
+
# Initialising the point-values
Init = np.zeros((4, len(PointsID)))
-
+
# Computing the density value to use
Rair = 0.287
rho = Rair/FluidProp.R
-
+
# Computing the values that will be taken outside the SOD bump
Init[:,:] = np.array([[rho, 0., 0., 1]]*len(PointsID)).T
e = EoS.EoS(Init[:,:], [FluidProp]*len(PointsID), "e")
@@ -778,16 +1191,16 @@ class InitialConditions():
# Returning the values
return(Init)
-
-
+
+
#### SOD test case ####
-
+
#### Karni fluid 2 ####
def PlanarShock(self, PointsID, Mesh, EoS, FluidProp, subdomain, MachNum=[]):
"""Initialising the solution to two constante state distributed on either side of a shock for the Planar shock test.
The shock is set at 25% of the total x width of the domain, starting from the right. The (left) rested part is set to
be such that (p, u, v, p) = [1., 0., 0., 1]. The shock is set to move from right to left.
-
+
Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
@@ -795,28 +1208,28 @@ class InitialConditions():
FluidProp (FluidSpecs): the properties of the fluid present where the given points are
subdomain (shapely multipolyon): the shapely polygon to which the given points belong
MachNum (float): the desired Mach number (given on the still fluid)
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
-
+
|
-
+
.. warning::
-
+
The computation of the state assumes air as fluid component
-
+
|
-
+
.. note::
-
+
There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
"""
-
+
# ------------------ Plan the routine interface ------------------------------------------------
# Extract the parameters if called from the instance
if not MachNum: MachNum = self.Params[0]
-
+
# ------------------ Retrieving the spatial location of the shock ------------------------------
# Retrieving the xy values of the points and the geometrical informations from the subdomain
xy = Mesh.DofsXY[PointsID,:]
@@ -824,16 +1237,16 @@ class InitialConditions():
maxx = np.max(coords[0])
minx = np.min(coords[0])
shocklocation = maxx-0.25*abs(maxx-minx)
-
+
# ----------- Computing the quantities to assign to the point values ---------------------------
# Initialising the point-values
Init = np.zeros((4, len(PointsID)))
-
+
# Initialising the still part of the fluid
StillState = [1., 0., 0., 1]
- # Computing the associated reckless part of the fluid according to the Mach Number, the fluid
+ # Computing the associated reckless part of the fluid according to the Mach Number, the fluid
# properties and the direction of the shock's propagation.
gmm = FluidProp.gamma
Ms = MachNum
@@ -843,17 +1256,17 @@ class InitialConditions():
if EoS.EoSName == "Ideal":
h = (gmm-1)/(gmm+1)
c0 = np.sqrt(gmm*StillState[3]/StillState[0])
-
+
RecklessState = [0.]*4
RecklessState[0] = StillState[0]*(Ms**2)/(1-h+h*Ms**2)
RecklessState[1] = RecklessState[0]*(StillState[1]+direct*(1-h)*c0*(Ms-1/Ms))
RecklessState[2] = RecklessState[0]*(StillState[2])
RecklessState[3] = StillState[3]*((1+h)*Ms**2-h)
-
+
else:
raise(NotImplementedError("Error: the asked initial conditions (planar shock) are not yet implemented for \
the used Equation of state. Please check your problem definition."))
-
+
# ------------------- Assigning the right values to the right points ------------------------
# Computing the values that will be taken on the resting part of the fluid
@@ -861,12 +1274,12 @@ class InitialConditions():
Init[:,Ind1] = np.array([StillState]*len(Ind1)).T
e = EoS.EoS(Init[:,Ind1], [FluidProp]*len(Ind1), "e")
Init[3,Ind1] = Init[0,Ind1]*e
-
+
# Computing the values that will be taken after the moving shock of the fluid and filling the right values
Ind2 = np.setdiff1d(range(len(PointsID)), Ind1)
Init[:,Ind2] = np.array([RecklessState]*len(Ind2)).T
e = EoS.EoS(Init[:,Ind2], [FluidProp]*len(Ind2), "e")
- Init[3,Ind2] = Init[0,Ind2]*(e+Init[1,Ind2]**2)
-
+ Init[3,Ind2] = Init[0,Ind2]*e+0.5*Init[1,Ind2]**2/Init[0,Ind2]+0.5*Init[2,Ind2]**2/Init[0,Ind2]
+
# Returning the values
return(Init)
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py
index 13802a5..c1d7f0a 100644
--- a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py
+++ b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py
@@ -1,8 +1,8 @@
"""
-All the necessary routines that implement a
-simple 2D translation, linear advection are defined in
+All the necessary routines that implement a
+simple 2D translation, linear advection are defined in
:code:`GoverningEquations/LinearAdvectionRotation.py`
"""
@@ -27,11 +27,11 @@ class Equations():
schemes and evolve the solution according to the Linear Advection (translation)
|
-
+
.. rubric:: Fields
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
It contains the following initialisation parameters (later available as well as attributes)
- |bs| **FluidProp** *(list of FluidModel)* -- the list of fluid property instances (NbFluids)
@@ -51,13 +51,13 @@ class Equations():
def __init__(self, FluidProp, EOs):
"""
|
-
+
Args:
FluidProp (list of FluidModel): the list of fluid property instances (NbFluids)
EOs (list of EOS): the list of Equation of state instances (NbFluids)
-
+
|
-
+
.. rubric:: Methods
"""
@@ -74,35 +74,81 @@ class Equations():
self.FluidProp = FluidProp
self.EOs = EOs
+ #### Retrieve the primary variables from the conservative ones ####
+ def ConservativeToPrimary(self, Var, FluidIndex, *args):
+ """Converts the conservative variables to the primary ones
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+
+ Returns:
+ Var (numpy array): (NbVars x NbGivenPoints) the corresponding primary variables
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ------- Constructing the primary variables ----------------------------------------------------
+ # Computing the primary variables
+ PrimVars = Var[:]
+
+ # Returning the Result
+ return(PrimVars)
+
+ #### Retrieve the conservative variables from the primary ones ####
+ def PrimaryToConservative(self, PrimVar, FluidIndex, *arrgs):
+ """Converts the primary variables to the conservative ones
+
+ Args:
+ Var (2D numpy array): the primary variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+
+ Returns:
+ Var (numpy array): (NbVars x NbGivenPoints) the corresponding conservative variables
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ------- Constructing the primary variables ----------------------------------------------------
+ # Computing the primary variables
+ Vars = PrimVar[:]
+
+ # Returning the Result
+ return(Vars)
+
#### Flux function of the equations set ####
def Flux(self, Var, x, *args):
"""Flux function of the LinearAdvection (rotation case)
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Args:
+
+ Args:
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
x (float numpy array): (generally optional, required for this problem), x-y coordinates of the given points (NbGivenPoints x 2)
.. rubric:: Optional argument
-
+
-- |bs| **FluidIndex** *(integer numpy array)* -- the fluid present at each given point (NbGivenPoints)
-
+
|
-
+
Returns:
Flux (3D numpy array): the obtained flux (spatialdim x NbVars x NbPoints)
-
+
|
-
+
.. note::
-
+
- | This function is Vectorised
- | Fluid index is the fluid index in the initially given list, not the fluid type
- | FluidProp is the list of all the fluid properties for each fluid (given the list, not the type)
- | *args is there only for compatibility reason at call time
-
-
+
+
"""
# --------- Computing the Euler flux directly in conservative variables ------------------------
@@ -117,25 +163,25 @@ class Equations():
def GetUV(self, Var, x, *argv):
"""Function giving back the x-y velocity at the considered point, given the variational
values given for the advection.
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Args:
+
+ Args:
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
x (float numpy array): (generally optional, required for this problem), x-y coordinates of the given points (NbGivenPoints x 2)
.. rubric:: Optional argument
-
+
-- |bs| **FluidIndex** *(integer numpy array)* -- the fluid present at each given point (NbGivenPoints)
-
+
Returns:
UV (float numpy array) -- the velocities values at the points (2 x NbGivenPoints)
.. note::
-
+
- | This function is Vectorised
- | *args is there only for compatibility reason at call time
-
+
--------------------------------------------------------------------------------------------------------------------------------"""
# Computing the velocity values from the conservative Euleur equations
@@ -147,31 +193,31 @@ class Equations():
#### Jacobian of the equations set ####
def Jacobian(self, Var, x, *argv):
"""Computes the Jacobian of the flux for the LinearAdvection (rotation case)
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Args:
+
+ Args:
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
x (float numpy array): (generally optional, required for this problem), x-y coordinates of the given points (NbGivenPoints x 2)
.. rubric:: Optional argument
-
+
-- |bs| **FluidIndex** *(integer numpy array)* -- the fluid present at each given point (NbGivenPoints)
-
- Returns:
- J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
-
+
+ Returns:
+ J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
+
|
-
+
.. note::
-
+
- | For each flux fi = (fi1,..., fin), the returned Jacobian reads:
-
+
| |bs| J[:,:] = [dfi1/dx1, ....., dfi1/dxn
| |bs| |bs| |bs| |bs| |bs| |bs| ....
| |bs| |bs| |bs| |bs| |bs| |bs| df in/dx1, ....., dfin/dxn]
- | *args is there only for compatibility reason at call time
-
+
"""
# ------ Initialising the Jacobian structure ---------------------------
@@ -192,20 +238,20 @@ class Equations():
#### Spectral radius of the equations set ####
def SpectralRadius(self, Var, FluidIndex, n, x, *argv):
""" Computes the spectral radius associated to the flux.
-
+
Args:
Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
- Returns:
+ Returns:
Lmb (numpy array): the spectral radius computed at each given point
-
+
.. note::
-
+
- | *args is there only for compatibility reason at call time
-
+
---------------------------------------------------------------------"""
# Retrieving the Jacobian
J = self.Jacobian(Var, x)
@@ -216,6 +262,92 @@ class Equations():
# Returning the spectral radius at each given point
return(Lmb)
+ #### Eigen values ####
+ def EigenValues(self, Var, FluidIndex, n, x, *args):
+ """Computes the eigenvalues associated to the flux.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+ x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Returns:
+ lbd (numpy array): (NbGivenPoints x 4) the four eigenvalues at each given point
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ---------- Ugly always working way ------------------------------
+
+ # --- Locating and extracting the right fluid's properties for each location
+ # List of fluids present in the list to be computed
+ Fluids = list(set(FluidIndex))
+
+ # Retrieving the Jacobian
+ J = self.Jacobian(Var, x)
+
+ # Computing the spectral radius at each given point, shortcuted since here we are scalar
+ lbd = np.reshape(np.sum(J[0,0,:,:]*n.T, axis = 0), (1, 1, Var.shape[1]))
+
+ # Returning the values
+ return(lbd)
+
+ #### Right eigen vectors ####
+ def RightEigenVectors(self, Var, FluidIndex, n, x, *args):
+ """Computes the right eigenvectors associated to the eigenvalues.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+ x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Returns:
+ reg (numpy array): (NbVars x MbVars x NbGivenPoints) the matrix of eigenvectors for each given point
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ---------------- Right eigenvalues -------------------------------------
+ # Creating right eigenvalues
+ REV = np.zeros((Var.shape[0], Var.shape[0], Var.shape[1]))
+ REV[:, :, :] = 1
+
+ # Returning the values
+ return(REV)
+
+ #### Right eigen vectors ####
+ def LeftEigenVectors(self, Var, FluidIndex, n, x, *args):
+ """Computes the left eigenvectors associated to the eigenvalues.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+ x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Returns:
+ reg (numpy array): (NbVars x MbVars x NbGivenPoints) the matrix of eigenvectors for each given point
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ---------------- Creating left eigenvalues ---------------------------
+ # Creating left eigenvalues
+ LEV = np.zeros((Var.shape[0], Var.shape[0], Var.shape[1]))
+ LEV[:, :, :] = 1
+
+ # Returning the values
+ return(LEV)
+
+
# ==============================================================================
# Class containing all the methods for defining various Equation of states
@@ -223,43 +355,43 @@ class Equations():
#### Class containing the methods giving back the pressure ####
class EOS():
"""Class furnishing the methods giving the Equation of State that can be
- used in combination with the :code:`SimpleFluid` definition.
+ used in combination with the :code:`SimpleFluid` definition.
Upon creation, fills the field :code:`EOS` with the function callback
corresponding to the wished initialisation function and the field :code:`EOSName` with
- the associated label. The implemented EoS are linked as follows.
- See their respective documentation for more information.
-
+ the associated label. The implemented EoS are linked as follows.
+ See their respective documentation for more information.
+
1. | Dummy
-
+
|
-
+
.. Note::
-
+
After having been initialised the instance, the EOS is accessible from the field EOS within this class.
- .. warning::
-
+ .. warning::
+
For advection, this class is void, just here for compatibility properties in the code
-
+
|
-
+
"""
#### Automatic initialisation routine (mapper) ####
def __init__(self, Id):
- """
-
+ """
+
Args:
Id (integer): the index corresponding to the equation of fluid the user wants when considering the governing equations given in this module and the associated fluid.
-
+
|
-
+
.. rubric:: Methods
-
+
"""
-
+
# Mapping towards the right function
- if Id==1:
+ if Id==1:
self.EoS = self.Dummy
self.EoSName = "Dummy"
@@ -270,15 +402,15 @@ class EOS():
#### Dummy EOS #####
def Dummy(self, Var, FluidProp, Type):
"""Dummy EOS returning the identity for the sake of the code compatibility.
-
+
Args:
Type (string): desired output: "e": internal energy, "p": pressure
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
FluidProp (FluidModel list): the list of the fluid properties associated to each given Point
-
- Returns:
+
+ Returns:
val: the variable retrieved from the EoS (here identity as dummy function)
-
+
"""
# Void return
@@ -297,24 +429,24 @@ class InitialConditions():
Access the routine computing the initial conditions through the field :code:`IC`, filled upon creation
with the function callback corresponding to the index of the wished initialisation function.
The implemented initialisation method are linked as follows. See their respective documentation for more information.
-
+
0. | ConstantState
1. | Not implemented
2. | Bump
-
+
|
-
+
"""
#### Automatic initialisation routine (mapper) ####
def __init__(self, Id, *params):
- """
+ """
Args:
Id (integer): the index corresponding to the initialisation method the user wants when considering the governing equations given in this module.
- params (list of arguments): the (fixed) arguments to pass to the selected function
-
+ params (list of arguments): the (fixed) arguments to pass to the selected function
+
|
-
+
.. rubric:: Methods
"""
@@ -328,21 +460,21 @@ class InitialConditions():
#### Initialise the solution to a constant state ####
def ConstantState(self, PointsID, Mesh, EOS = [], FluidProp = [], subdomain = []):
""" Initialising the solution to a constant state on the given points, all belonging to a same subdomain.
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): (optional), the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): (optional, not used here) the the properties of the fluid present where the given points are
Subdomain (shapely multipolyon): (optional, not used here) the shapely polygon to which the given points belong
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
.. note::
-
+
There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
-
+
------------------------------------------------------------------------------------------------------------------"""
# Initialising the point-values
@@ -358,22 +490,22 @@ class InitialConditions():
def Bump(self, PointsID, Mesh, EOS = [], FluidProp = [], subdomain = []):
""" IInitialising the solution to a bump centred at the center of mass
of the subdmain. (The given points should all belonging to a same subdomain).
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): (optional), the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): (optional, not used here) the the properties of the fluid present where the given points are
Subdomain (shapely multipolyon): (optional, not used here) the shapely polygon to which the given points belong
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
.. note::
-
+
- | There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
- | For running smothly the inner circle diameter of the given subdomain should be at least 2
-
+
------------------------------------------------------------------------------------------------------------------"""
# Initialising the point-values
Init = np.ones((1, len(PointsID)))
@@ -386,7 +518,7 @@ class InitialConditions():
# Defining the space-dependent quantities
r2 = ((xy[:,0]-xc)**2+(xy[:,1]-yc)**2)
co = 3*np.exp(-1/(1-r2))*(r2<1)
-
+
Init[0,:] = co
# Returning the computed values
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py
index ed4b4a8..935c205 100644
--- a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py
+++ b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py
@@ -1,7 +1,7 @@
"""
-All the necessary routines that implements a
-simple 2D translation, linear advection are defined in
+All the necessary routines that implements a
+simple 2D translation, linear advection are defined in
:code:`GoverningEquations/LinearAdvectionTranslation.py`
"""
@@ -26,11 +26,11 @@ class Equations():
schemes and evolve the solution according to the Linear Advection (translation)
|
-
+
.. rubric:: Fields
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
It contains the following initialisation parameters (later available as well as attributes)
- |bs| **FluidProp** *(list of FluidModel)* -- the list of fluid property instances (NbFluids)
@@ -45,22 +45,21 @@ class Equations():
- |bs| **VariablesLatexUnits** *(list of strings)* -- symbols of the units of the variables, latex encoding
- |bs| **XYLatexUnits** *(list of strings)* -- units of the x-y coordinates in latex encoding
-
"""
#### Automatic initialisation routine ####
def __init__(self, FluidProp, EOs):
- """
+ """
|
-
+
Args:
FluidProp (list of FluidModel): the list of fluid property instances (NbFluids)
EOS (list of EOS): the list of Equation of state instances (NbFluids)
-
+
|
-
+
.. rubric:: Methods
-
+
"""
# Defines the (conservative) variables number and names, in a latex typography
@@ -76,32 +75,78 @@ class Equations():
self.FluidProp = FluidProp
self.EOs = EOs
+ #### Retrieve the primary variables from the conservative ones ####
+ def ConservativeToPrimary(self, Var, FluidIndex, *args):
+ """Converts the conservative variables to the primary ones
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+
+ Returns:
+ Var (numpy array): (NbVars x NbGivenPoints) the corresponding primary variables
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ------- Constructing the primary variables ----------------------------------------------------
+ # Computing the primary variables
+ PrimVars = Var[:]
+
+ # Returning the Result
+ return(PrimVars)
+
+ #### Retrieve the conservative variables from the primary ones ####
+ def PrimaryToConservative(self, PrimVar, FluidIndex, *arrgs):
+ """Converts the primary variables to the conservative ones
+
+ Args:
+ Var (2D numpy array): the primary variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+
+ Returns:
+ Var (numpy array): (NbVars x NbGivenPoints) the corresponding conservative variables
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ------- Constructing the primary variables ----------------------------------------------------
+ # Computing the primary variables
+ Vars = PrimVar[:]
+
+ # Returning the Result
+ return(Vars)
+
#### Flux function of the equations set ####
def Flux(self, Var, x, *args):
"""Flux function of the LinearAdvection (translation case)
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Args:
+
+ Args:
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
x (float numpy array): (generally optional, required for this problem), x-y coordinates of the given points (NbGivenPoints x 2)
.. rubric:: Optional argument
-
+
- |bs| **FluidIndex** *(integer numpy array)* -- the fluid present at each given point (NbGivenPoints)
-
+
Returns:
Flux (3D numpy array): the obtained flux (spatialdim x NbVars x NbPoints)
-
+
|
-
+
.. note::
-
+
- | This function is Vectorised
- | Fluid index is the fluid index in the initially given list, not the fluid type
- | FluidProp is the list of all the fluid properties for each fluid (given the list, not the type)
- | *args is there only for compatibility reason at call time
-
+
"""
# --------- Computing the Euler flux directly in conservative variables ------------------------
@@ -116,25 +161,25 @@ class Equations():
def GetUV(self, Var, x, *argv):
"""Function giving back the x-y velocity at the considered point, given the variational
values given for the advection.
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Args:
+
+ Args:
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
x (float numpy array): (generally optional, required for this problem), x-y coordinates of the given points (NbGivenPoints x 2)
.. rubric:: Optional argument
-
+
- |bs| **FluidIndex** *(integer numpy array)* -- the fluid present at each given point (NbGivenPoints)
-
+
Returns:
UV (float numpy array) -- the velocities values at the points (2 x NbGivenPoints)
.. note::
-
+
- | This function is Vectorised
- | *args is there only for compatibility reason at call time
-
+
"""
# Computing the velocity values from the conservative Euleur equations
@@ -146,31 +191,31 @@ class Equations():
#### Jacobian of the equations set ####
def Jacobian(self, Var, x, *argv):
"""Computes the Jacobian of the flux for the LinearAdvection (rotation case)
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Args:
+
+ Args:
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
x (float numpy array): (generally optional, required for this problem), x-y coordinates of the given points (NbGivenPoints x 2)
.. rubric:: Optional argument
-
+
- |bs| **FluidIndex** *(integer numpy array)* -- the fluid present at each given point (NbGivenPoints)
-
- Returns:
- J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
-
+
+ Returns:
+ J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
+
|
-
+
.. note::
-
+
- | For each flux fi = (fi1,..., fin), the returned Jacobian reads:
-
+
| |bs| J[:,:] = [dfi1/dx1, ....., dfi1/dxn
| |bs| |bs| |bs| |bs| |bs| |bs| ....
| |bs| |bs| |bs| |bs| |bs| |bs| df in/dx1, ....., dfin/dxn]
- | *args is there only for compatibility reason at call time
-
+
"""
# ------ Initialising the Jacobian structure ---------------------------
@@ -191,20 +236,20 @@ class Equations():
#### Spectral radius of the equations set ####
def SpectralRadius(self, Var, FluidIndex, n, x, *argv):
"""Computes the spectral radius associated to the flux.
-
+
Args:
Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
- Returns:
+ Returns:
Lmb (numpy array): the spectral radius computed at each given point
-
+
.. note::
-
+
- | *args is there only for compatibility reason at call time
-
+
---------------------------------------------------------------------"""
# Retrieving the Jacobian
J = self.Jacobian(Var, x)
@@ -215,6 +260,92 @@ class Equations():
# Returning the spectral radius at each given point
return(Lmb)
+ #### Eigen values ####
+ def EigenValues(self, Var, FluidIndex, n, x, *args):
+ """Computes the eigenvalues associated to the flux.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+ x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Returns:
+ lbd (numpy array): (NbGivenPoints x 4) the four eigenvalues at each given point
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ---------- Ugly always working way ------------------------------
+
+ # --- Locating and extracting the right fluid's properties for each location
+ # List of fluids present in the list to be computed
+ Fluids = list(set(FluidIndex))
+
+ # Retrieving the Jacobian
+ J = self.Jacobian(Var, x)
+
+ # Computing the spectral radius at each given point, shortcuted since here we are scalar
+ lbd = np.reshape(np.sum(J[0,0,:,:]*n.T, axis = 0), (1, 1, Var.shape[1]))
+
+ # Returning the values
+ return(lbd)
+
+ #### Right eigen vectors ####
+ def RightEigenVectors(self, Var, FluidIndex, n, x, *args):
+ """Computes the right eigenvectors associated to the eigenvalues.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+ x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Returns:
+ reg (numpy array): (NbVars x MbVars x NbGivenPoints) the matrix of eigenvectors for each given point
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ---------------- Right eigenvalues -------------------------------------
+ # Creating right eigenvalues
+ REV = np.zeros((Var.shape[0], Var.shape[0], Var.shape[1]))
+ REV[:, :, :] = 1
+
+ # Returning the values
+ return(REV)
+
+ #### Right eigen vectors ####
+ def LeftEigenVectors(self, Var, FluidIndex, n, x, *args):
+ """Computes the left eigenvectors associated to the eigenvalues.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ FluidIndex (integer numpy array): the fluid present at each given point (NbGivenPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+ x (2D numpy array): (generally optional, required for this problem) the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Returns:
+ reg (numpy array): (NbVars x MbVars x NbGivenPoints) the matrix of eigenvectors for each given point
+
+ .. note::
+
+ - | *args is there only for compatibility reason at call time
+ """
+
+ # ---------------- Creating left eigenvalues ---------------------------
+ # Creating left eigenvalues
+ LEV = np.zeros((Var.shape[0], Var.shape[0], Var.shape[1]))
+ LEV[:, :, :] = 1
+
+ # Returning the values
+ return(LEV)
+
+
# ==============================================================================
# Class containing all the methods for defining various Equation of states
@@ -222,36 +353,36 @@ class Equations():
#### Class containing the methods giving back the pressure ####
class EOS():
"""Class furnishing the methods giving the Equation of State that can be
- used in combination with the :code:`SimpleFluid` definition.
+ used in combination with the :code:`SimpleFluid` definition.
Upon creation, fills the field :code:`EOS` with the function callback
corresponding to the wished initialisation function and the field :code:`EOSName` with
- the associated label. The implemented EoS are linked as follows.
+ the associated label. The implemented EoS are linked as follows.
See their respective documentation for more information.
1. | Dummy
-
+
|
-
+
.. Note::
-
+
After having been initialised the instance, the EOS is accessible from the field EOS within this class.
- .. warning::
-
+ .. warning::
+
For advection, this class is void, just here for compatibility properties in the code
-
+
|
-
+
"""
#### Automatic initialisation routine (mapper) ####
def __init__(self, Id):
- """
+ """
Args:
Id (integer): the index corresponding to the equation of fluid the user wants when considering the governing equations given in this module and the associated fluid.
-
+
|
-
+
.. rubric:: Methods
"""
@@ -266,15 +397,15 @@ class EOS():
#### Dummy EOS #####
def Dummy(self, Var, FluidProp, Type):
"""Dummy EOS returning the identity for the sake of the code compatibility.
-
+
Args:
Type (string): (dummy here)
Var (float numpy array): the value of the variables (NbVariables x NbGivenPoints)
FluidProp (FluidModel list): (dummy here) the list of the fluid properties associated to each given Point
-
- Returns:
+
+ Returns:
val: the variable retrieved from the EoS (here identity as dummy function)
-
+
----------------------------------------------------------------------------------------------------"""
# Void return
@@ -292,32 +423,32 @@ class InitialConditions():
to study the the Euler Equations, defined for working on a subdomain.
Access the routine computing the initial conditions through the field :code:`IC`, filled upon creation
with the function callback corresponding to the index of the wished initialisation function.
- The implemented initialisation method are linked as follows.
+ The implemented initialisation method are linked as follows.
See their respective documentation for more information.
-
+
0. | ConstantState
1. | ConstantState2
2. | Bump
3. | ConstantState3
4. | BumpEdge
-
+
|
-
+
"""
#### Automatic initialisation routine (mapper) ####
def __init__(self, Id, *params):
- """
+ """
Args:
Id (integer): the index corresponding to the initialisation method the user wants when considering the governing equations given in this module.
params (list of arguments): the (fixed) arguments to pass to the selected function
-
+
|
-
+
.. rubric:: Methods
-
+
"""
-
+
# Mapping the initialisation routines
if Id == 0: self.IC = self.ConstantState
elif Id == 3: self.IC = self.ConstantState3
@@ -330,25 +461,25 @@ class InitialConditions():
#### Initialise the solution to a constant state ####
def ConstantState(self, PointsID, Mesh, EOS = [], FluidProp = [], subdomain = []):
- """
-
+ """
+
Initialising the solution to a constant state on the given points, all belonging to a same subdomain.
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): (optional), the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): (optional, not used here) the the properties of the fluid present where the given points are
Subdomain (shapely multipolyon): (optional, not used here) the shapely polygon to which the given points belong
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
.. note::
-
+
There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
-
-
+
+
------------------------------------------------------------------------------------------------------------------"""
# Initialising the point-values
@@ -359,25 +490,25 @@ class InitialConditions():
# Returning the computed values
return(Init)
-
+
#### Initialise the solution to a constant state ####
def ConstantState3(self, PointsID, Mesh, EOS = [], FluidProp = [], subdomain = []):
""" Initialising the solution to a constant state on the given points, all belonging to a same subdomain.
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): (optional), the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): (optional, not used here) the the properties of the fluid present where the given points are
Subdomain (shapely multipolyon): (optional, not used here) the shapely polygon to which the given points belong
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
.. note::
-
+
There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
-
+
------------------------------------------------------------------------------------------------------------------"""
# Initialising the point-values
@@ -388,25 +519,25 @@ class InitialConditions():
# Returning the computed values
return(Init)
-
+
#### Initialise the solution to a constant state ####
def ConstantState2(self, PointsID, Mesh, EOS = [], FluidProp = [], subdomain = []):
""" Initialising the solution to a constant state on the given points, all belonging to a same subdomain.
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): (optional), the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): (optional, not used here) the the properties of the fluid present where the given points are
Subdomain (shapely multipolyon): (optional, not used here) the shapely polygon to which the given points belong
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
.. note::
-
+
There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
-
+
------------------------------------------------------------------------------------------------------------------"""
# Initialising the point-values
@@ -422,7 +553,7 @@ class InitialConditions():
def BumpEdge(self, PointsID, Mesh, EOS = [], FluidProp = [], subdomain = []):
"""
Initialising the solution to a Bump centred next to the righter boundary of the subdmain. (The given points should all belonging to a same subdomain).
-
+
Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
@@ -439,7 +570,7 @@ class InitialConditions():
- | For running smothly the inner circle diameter of the given subdomain should be at least 2
"""
-
+
# Initialising the point-values
Init = np.ones((1, len(PointsID)))
@@ -447,11 +578,11 @@ class InitialConditions():
xy = Mesh.DofsXY[PointsID,:]
xc = np.max([np.max(np.array(subdo.exterior.coords.xy).T[:,0]) for subdo in subdomain])-1
yc = subdomain.centroid.y
-
+
# Defining the space-dependent quantities
r2 = ((xy[:,0]-xc)**2+(xy[:,1]-yc)**2)
co = 3*np.exp(-1/(1-r2))*(r2<1)
-
+
Init[0,:] = co
# Returning the computed values
@@ -461,24 +592,24 @@ class InitialConditions():
def Bump(self, PointsID, Mesh, EOS = [], FluidProp = [], subdomain = []):
"""Initialising the solution to a bump centred at the center of mass
of the subdmain. (The given points should all belonging to a same subdomain).
-
- Args:
+
+ Args:
PointsID (integer array-like): the index of the points to consider
Mesh (MeshStructure): the considered mesh
EOS (function callback): (optional), the equation of state given by the model of the fluid that is present at the given points
FluidProp (FluidSpecs): (optional, not used here) the the properties of the fluid present where the given points are
Subdomain (shapely multipolyon): (optional, not used here) the shapely polygon to which the given points belong
-
+
Returns:
Init (float numpy array): The initialised values at the considered points (NbVariables x NbGivenPoints)
.. note::
-
+
- | There is usually one initialisation routine per subdomain, and the PointsID are not necessarily contigous, be careful when assigning the values back in the regular solution.
- | For running smothly the inner circle diameter of the given subdomain should be at least 2
-
+
------------------------------------------------------------------------------------------------------------------"""
-
+
# Initialising the point-values
Init = np.ones((1, len(PointsID)))
@@ -486,12 +617,12 @@ class InitialConditions():
xy = Mesh.DofsXY[PointsID,:]
xc = subdomain.centroid.x
yc = subdomain.centroid.y
-
+
# Defining the space-dependent quantities
r2 = ((xy[:,0]-xc)**2+(xy[:,1]-yc)**2)
co = 3*np.exp(-1/(1-r2))*(r2<1)
#Init[0,:] = co-2.
-
+
#co = (r2<1)
Init[0,:] = co
diff --git a/SourceCode/Solve.py b/SourceCode/Solve.py
index 6ef6384..f710ada 100644
--- a/SourceCode/Solve.py
+++ b/SourceCode/Solve.py
@@ -8,9 +8,9 @@ Main file that solves a prefefined mutliphase problem
- Create first a problem in a given module (follow the problem generation docs)
- Create or use a prefefined mesh for which this problem should be solved (follow the meshing docs)
- Solve your problem in a python environment by loading this module and simply running
-
+
..code ::
-
+
Main.Solve(ProblemDefinition, SettingsChoice, MeshResolution)
.. note::
@@ -53,17 +53,17 @@ import _Solver
def Solve(*argv):
"""
Solving routine that considers a problem, a mesh resolution and a solver settings set.
- Inputs (in that order):
-
+ Inputs (in that order):
+
Args:
-
+
ProblemName (string): name of the problem module to investigate, stored in the "PredefinedProblems" folder (see the problem definition's instructions for more details)
Settings (string): name of the settings file specifying the scheme, time options etc to use for the study, stored in the "PredefinedSettings" folder
Resolution (integer): the resolution of the mesh to be used (defined as in the Fenics' library. A corresponding list of max element's size can be found in the root folder of this code)
Returns:
None: The solution's file in both Numpy binary format and text format for each time step at the intervals given in the Settings' file. The post-processing of the solution is not automatically performed.
-
+
---------------------------------------------------------------------------------------------------------"""
@@ -95,8 +95,8 @@ def Solve(*argv):
ProblemData = eval("Prob."+argv[0]+".ProblemData()")
SolverParams = _Solver.Solver_Initialisations.Parameters(argv[1])
print("--> Done.")
-
- # Clean all previously computed solution within the destination folder
+
+ # Clean all previously computed solution within the destination folder
print("Cleaning the export environment")
EX.CleanExportEnvironment(TestName, argv[2])
print("--> Done.")
@@ -106,7 +106,7 @@ def Solve(*argv):
Mesh = Meshing.Meshing_DualMesh.DualMesh(ProblemData.MeshName, SolverParams.MeshOrder,\
SolverParams.MeshType, argv[2]).Mesh
print("--> Done.")
-
+
# ----------- Perform the initialisation tasks ------------------------------------------------
# Initialise the variables to work with shortcut the mapped scheme routines
print("Initialising the problem")
@@ -117,7 +117,7 @@ def Solve(*argv):
print("--> Initialising the solution")
Solution = _Solver.Solver_Initialisations.Solution(Problem, Mesh)
print("--> Done.")
-
+
# Generating the mass lumping coefficients
print("Generating the mass lumping coefficients")
_Solver.Solver_Temporal.GetLumpingCoefficients(Mesh, Solver)
@@ -138,20 +138,20 @@ def Solve(*argv):
# Initialising the simple running time
Solution.t=0
-
+
# Export the initial considition after having reconstructed the initial values at the physical vertices
Solver.Reconstructor(Solution)
EX.ExportSolution(TestName, Mesh, ProblemData, Problem, Solution)
-
- # Inform the user
+
+ # Inform the user
PO.Printout_MinMaxSolVal(0.0, Problem, Solution)
-
+
# Main time loop
k=0
while Solution.t < Solver.Tmax:
# ----- Update step informations
k=k+1
-
+
# ----- Perform a time step
# Get the time step
dt = _Solver.Solver_Temporal.GetTimeStep(Solver.CFL, Problem, Mesh, Solution)
@@ -160,23 +160,23 @@ def Solve(*argv):
schemeSC = [Solver.SpatialScheme, Solver.SPAdmendements]
schemeFS = [Solver.FluidSelector]
Solution.Sol, Solution.LSValues = Solver.TimeStep(schemeSC, schemeFS, dt, Problem, Mesh, Solution)
-
+
# ----- Perform the subsidiary updates
# Increase the current time
Solution.t += dt
-
+
# Reconstruct the solution at the physical vertices of the mesh (just impacting the visualisation)
Solver.Reconstructor(Solution)
-
+
# ----- Inform the user and export the solution when relevant
- if np.mod(k,10)==0:
+ if np.mod(k,1)==0:
PO.Printout_MinMaxSolVal(Solution.t, Problem, Solution)
EX.ExportSolutionVTK(TestName, Mesh, ProblemData, Problem, Solution)
EX.ExportSolution(TestName, Mesh, ProblemData, Problem, Solution)
- else: print("Computing for the time {0!s} in progress".format(Solution.t))
-
-
-
+ else: print("Computing for the time {0!s}".format(Solution.t))
+
+
+
# ------------------ Closing tasks -----------------------------------------
print("\n -------------------------------------\
\n Performing post-computations tasks \
diff --git a/SourceCode/_Solver/Solver_Spatial.py b/SourceCode/_Solver/Solver_Spatial.py
index 88c4295..9554615 100644
--- a/SourceCode/_Solver/Solver_Spatial.py
+++ b/SourceCode/_Solver/Solver_Spatial.py
@@ -95,5 +95,6 @@ def RD(Theta, m, M, dt, schemeSP, schemeFS, Problem, Mesh, Solution, up, lp):
upres = upres*Mesh.Lumping
lsres = lsres*Mesh.Lumping
+
# Retuning the residuals'increments (not the new value of the solution)
return([upres, lsres])
diff --git a/SourceCode/_Solver/SpatialModifiers/Filtering_Streamline.py b/SourceCode/_Solver/SpatialModifiers/Filtering_Streamline.py
index 81b2583..3dcb27d 100644
--- a/SourceCode/_Solver/SpatialModifiers/Filtering_Streamline.py
+++ b/SourceCode/_Solver/SpatialModifiers/Filtering_Streamline.py
@@ -10,6 +10,8 @@ import copy
# Import custom modules
import _Solver.QuadratureRules.HandBased as QR
+import BarycentricBasisFunctions as BF
+import BarycentricTools as BT
# ==============================================================================
# Parameter (and reader) class of scheme properties and runtime informations
@@ -22,17 +24,23 @@ LimiterName = "Streamline"
# ==============================================================================
class Filtering():
""" .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Class furnishing all the schemes tools and registers the scheme's wished
parameters. The main iteration is accessible with the method "Mollify".
-
+
.. rubric:: Fields
-
+
All the fields given as arguments at the instance's creation are repeated within
the structure with the same names (see the parameters below).
-
+
|
+
+ .. warning::
+ The streamline stabilisation only seems to work for advection so far, it crashes on euler
+
+ |
+
"""
#### Automatic initialisation routine ####
@@ -42,142 +50,158 @@ class Filtering():
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh instance
Params (string list): the limiters's parameters as wished by the user
-
+
|
-
+
.. rubric:: Methods
-
+
"""
-
+
self.Problem = Problem
self.Params = Params
self.Mesh = Mesh
#### Main routine, definining the limiting mollifier ####
- def Mollify(self, SolBuff, resu, fluxes, i, du=0, dt=1):
+ def Mollify(self, Solution, resu, fluxes, i, du=0, dt=1):
""" .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- Streamline dissipation (not implemented yet)
-
+
+ Streamline dissipation, that has to be ADDED to the residual. In the input line, do not
+ simply write Filtering_1, but write 0+Filtering_1 so that it gets added to the initial residual
+
Args:
Solution (solution structure): structure containing the current solution's values to iterate
resu (float numpy array): previoulsy computed residuals that have to be modified (NbVars x NbElementDofs)
fluxes (numpy (multiD)array): pre-computed fluxes at the points of interest. For this scheme, access with fluxes[element, face, coordinate(fx or fy), variable, pointindex]
i (integer): the index of the considered element within which the partial residuals will be computed
- du (float numpy array): (optional) when using DeC, the difference in the time iteration
+ du (float numpy array): (optional) when using DeC, the difference in the time iteration
dt (float): (optional) when using DeC, the full time step
-
+
Returns:
Lresu (float numpy array): the limited residuals (NbVars x NbElementDofs)
-
+
---------------------------------------------------------------------"""
- # -------- Initialisation ------------------------------------------------
+ # -------- Initialisation ----------------------------------------------
# Getting the number dofs contaned in the considered element
NbVars, NbDofs = np.shape(resu)
-
+
# Initialising the mollification vector (residuals)
Lresu = np.zeros(np.shape(resu))
-
+
# Setting the safety division trick
safe = 1e-42
-
- # -------- Computations of the actual streamline ------------------------
- # Retrieving the scheme order
- Order = int(self.Params[0])
-
+
+ # -------- Precomputations of the actual streamline --------------------
+ # Retrieving the mesh order, furnishing the variatianal order to set as streamline
+ # according to the number of Dofs per element
+ Order = int(self.Mesh.MeshOrder)
+
# Getting the dofs contaned in the considered element
dele = np.array(self.Mesh.ElementsDofs[i])
NbDofs = len(dele)
# Initialising the mollification vector (residuals)
Resu = np.zeros([np.shape(Solution.Sol)[0], NbDofs])
-
- # Retrieving the inner and boundary fluxes
- Flx = fluxes[0, :, 0, :, :]
- Fln = fluxes[1,:,:,:,:,0:3]
-
+
+ # Retrieving the fluxes at the Dofs
+ Fld = fluxes[-1,:,:,:,:,0:len(dele)]
+
# ---- Getting the element and state informations
- # Retrieving the physical information on the spanned element's vertices
+ # Retrieving the physical information on the spanned element's vertices
ele = self.Mesh.InnerElements[i]
vol = self.Mesh.InnerElementsVolume[i]
n = self.Mesh.ElementsBoundaryWiseNormals[i]
-
+
# Retrieving the information on the local degree of freedom
LS = Solution.LSValues[:, dele]
U = Solution.Sol[:, dele]
dU = du[:, dele]
-
- # ---- Computing the plain integral
- # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+
+ # Get the coordinates of the vertices in the physical order
+ elevert = self.Mesh.CartesianPoints[ele,:]
+
+ # ---- Computing the evolution information
+ # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
# in the Dof's vertex ordering referential
bqp, weights = QR.InnerQuadrature()
-
+ qp = BT.GetCartesianCoordinates(bqp, elevert)
+
# Evaluate the basis function in the order of the Dofs
EvaluatedBase = BF.SimplicialBasis("B", Order, bqp)
- EvaluatedGradients = BF.SimplicialGradients("B", Order, bqp, n, vol)
-
+ EvaluatedGradients = BF.SimplicialGradients("B", Order, bqp, n, vol)
+
# Reconstructing the solution and its gradient at the quadrature point
Uloc = np.dot(U, EvaluatedBase.T)
dUloc = np.dot(dU, EvaluatedBase.T)
+ GUloc = np.matmul(EvaluatedGradients.T, U.T)
LSloc = np.dot(LS, EvaluatedBase.T)
-
- # Compute the flux according to the fluid's flag
- Flux = Flx[i]
-
- """
-
- u_loc=0._dp
-
- DO l=1, e%nvertex
- u_loc=u_loc+u(l)
- ENDDO
-
- u_loc=u_loc/REAL(e%nvertex,dp)
- Nmat=u_loc%Nmat(e%nvertex, e%n/e%volume ,xx, dt)
- DO i=1, e%nsommets
- DO iq=1, e%nquad
+ # Reconstructing back the conservative variables at the wished points
+ # (helps to preserve the pressure contacts, comment if E has to be preserved as well)
+ FlagsLoc = np.argmax(LSloc, axis=0)
- DO l=1, e%nsommets
- base(l )=e%base (l, e%quad(:,iq))
- grad(:,l)=e%gradient(l, e%quad(:,iq))
- ENDDO
-
-
-
- u_loc=SUM(base*u)
- du_loc=SUM(base*du)
- gru_loc(1)=SUM( grad(1,:)*u(:))
- gru_loc(2)=SUM( grad(2,:)*u(:))
- Jac=u_loc%Jacobian(xx)
-
- divflux_loc=0._dp
- DO l=1, n_dim
- divflux_loc = divflux_loc + SUM(grad(l,:)*flux(l,:))
- !divflux_loc%u=divflux_loc%u+ matmul(Jac(:,:,l), gru_loc(l)%u )
- ENDDO
-
-
- divflux_loc%u=MATMUL(nmat, du_loc%u(:)+ dt* divflux_loc%u(:) )
-
-
- ! res(i)%u = res(i)%u + MATMUL(TRANSPOSE( Jac(:,:,1)*grad(1,i)+Jac(:,:,2)*grad(2,i) ), divflux_loc%u )*e%weight(iq)
- res(i)%u = res(i)%u + MATMUL( Jac(:,:,1)*grad(1,i)+Jac(:,:,2)*grad(2,i) , divflux_loc%u )*e%weight(iq)
-
- ENDDO ! iq
-
- res(i)= res(i)* e%volume
- ENDDO ! i
+ # Compute the flux according to the fluid's flag
+ Flux = Fld[i,0]
+
+ # Compute the jacobian at the local points
+ Jac = self.Problem.GoverningEquations.Jacobian(Uloc, qp, FlagsLoc)
+
+
+ # --------------- Computing the stabilisation material --------------------
+ # --- Extract the required quantities
+ # Retrieving the Dofs that are associated to the vertices and averaging the solution value
+ verts = np.array(self.Mesh.Vertex2Dofs)[ele].flatten()
+ # Extracting the Dofs of the considered triangle (required in case of DG only)
+ verts = [v for v in verts if v in dele]
+ # Averaging the solution from the considered dofs
+ unmat = np.array([np.mean(Solution.Sol[:,verts], axis=1)]).T
+ # Considering the cell center as a point of evaluation
+ xymat = np.array([self.Mesh.DualCenters[i][-1]])
+ # Getting the fluid on which it is
+ lsmat = np.array([np.argmax(np.mean(Solution.LSValues[:,verts], axis=1))])
+
+ # Computing a first order gradient per given submit
+ if len(ele)==3: Grad = np.roll(self.Mesh.ElementsBoundaryWiseNormals[i], -1, axis=0)/self.Mesh.InnerElementsVolume[i]
+ else: raise NotImplementedError("No gradient has been implemented correctly for the filtering streamline yet.")
+
+ # Initialising Nmat
+ Nmat = np.zeros((self.Problem.GoverningEquations.NbVariables, self.Problem.GoverningEquations.NbVariables))
+
+ # Computing Nmat
+ for v in range(len(ele)):
+ # Getting the eigenvalues at the center of mass of the element
+ eig = self.Problem.GoverningEquations.EigenValues(unmat, lsmat, Grad[v:v+1,:], xymat)[0,:]
+ # Getting the safety net ratio
+ beta = np.max(np.abs(eig))*0.1
+ dd = np.zeros(eig.shape)
+ for j in range(len(eig)):
+ if np.abs(eig[j]) >= beta: dd[j] = np.abs(eig[j])
+ else: dd[j] = ((eig[j])**2+beta**2)/(2*beta)
+
+ # Constructing the matrix
+ dia = np.diag(dd)
+ R = self.Problem.GoverningEquations.RightEigenVectors(unmat, lsmat, Grad[v:v+1,:], xymat)[:,:,0]
+ L = self.Problem.GoverningEquations.LeftEigenVectors (unmat, lsmat, Grad[v:v+1,:], xymat)[:,:,0]
+ A = np.dot(R, np.dot(dia, L))
+ Nmat += A
+
+ # Inverting the obtained term
+ Nmat = np.linalg.inv(Nmat)
+
+ # ---- Computing the local divergence flux
+ # Getting the flux values at the
+ divflux = np.sum(np.matmul(np.swapaxes(EvaluatedGradients, 0,2),np.swapaxes(Flux, 2, 1)), axis=0)
+ divflux = np.matmul(Nmat, dUloc + dt*divflux.T)
+
+ # ---- Computing the stabilisation term
+ # Spanning the dofs to mollify
+ for dof in range(NbDofs):
+ # Performing the quadrature
+ stream = np.sum([weights[q]*np.matmul(EvaluatedGradients[dof,q,0]*Jac[:,:,0,q]+EvaluatedGradients[dof,q,1]*Jac[:,:,1,q], divflux[:,q]) for q in range(len(weights))], axis=0)
+
+ # Updating the residual at the spanned Dof by the local contribution
+ Lresu[:,dof] += (stream)*vol
-
- """
-
-
# -------- Returning the additive ----------------------------------------
-
- print("WARNING:THIS ROUTINE IS VOID AND HAS TO BE FURTHER IMPLEMENTED")
- Lresu = resu
-
- # Returning the partial residuals
+ # Returning the stream term that should be added to the original resuduals
return(Lresu)
diff --git a/SourceCode/_Solver/SpatialModifiers/SpatialAmendements.py b/SourceCode/_Solver/SpatialModifiers/SpatialAmendements.py
index b3bbde6..f43fad8 100644
--- a/SourceCode/_Solver/SpatialModifiers/SpatialAmendements.py
+++ b/SourceCode/_Solver/SpatialModifiers/SpatialAmendements.py
@@ -244,7 +244,13 @@ class Amendements():
string = "".join([piec for piec2 in [[pieces[k], operator[k]] for k in range(len(operator))] for piec in piec2]+[pieces[-1]])
self.MolliProcess[i] = [copy.deepcopy(lambda dumb, resu, *argv, string=string: eval(string)), np.arange(len(pieces))]
- # Modifiers of the type Modifyer_Index acting on the original residual
+ # Case of untouched original residual (used for simple additions of e.g. filters to the original residual)
+ elif "Original" in string:
+ # Map to the identity method
+ Modif = lambda dumb, resu, *argv: resu[:,:]
+ self.MolliProcess[i] = [copy.deepcopy(Modif), 0]
+
+ # Modifiers of the type Modifyer_Index, acting on the original residual
else:
# Extract the parameters if any
if re.findall("\<(.*?)\>", string):
@@ -258,7 +264,8 @@ class Amendements():
ResiInd = 0
Modif = eval("Mappers."+Type+"("+Index+")."+Type+"(self.Problem, self.Mesh, *Properties)")
self.MolliProcess[i] = [copy.deepcopy(Modif.Mollify), ResiInd]
-
+
+
#### Main routine, definining the limiting mollifier ####
def Mollify(self, SolBuff, resu, fluxes, i, du=0, dt=1):
""" Limting routine according to the given sequence stored in the structure
diff --git a/SourceCode/_Solver/SpatialSchemes/CG.py b/SourceCode/_Solver/SpatialSchemes/CG.py
index 1f8c250..20ce2a2 100644
--- a/SourceCode/_Solver/SpatialSchemes/CG.py
+++ b/SourceCode/_Solver/SpatialSchemes/CG.py
@@ -105,7 +105,7 @@ class Scheme():
# Defining the output fluxes (inner and boundary ones)
# The very first coordinate corresponds to the inner flux. Only Flux[0, :, 0, :, :, :] is filled
# The second coordinate corresponds to the inner flux. Only Flux[1, :, :, :, :, 0:3] is filled
- Flux = np.zeros((2, self.Mesh.NbInnerElements, 3, 2, self.Problem.GoverningEquations.NbVariables, 12))
+ Flux = np.zeros((3, self.Mesh.NbInnerElements, 3, 2, self.Problem.GoverningEquations.NbVariables, np.max([12, *[len(dele) for dele in self.Mesh.ElementsDofs]])))
# Retrieving the scheme order
Order = int(self.Params[0])
@@ -125,8 +125,12 @@ class Scheme():
UU = U[:, dele]
# Get the coordinates of the vertices in the physical order
- elevert = self.Mesh.CartesianPoints[ele,:]
+ elevert = self.Mesh.CartesianPoints[ele,:]
+ delevert = self.Mesh.DofsXY[dele,:]
+ # ---- Computing the flux at the Dofs directly
+ Flux[2, i, 0, :, :, 0:len(dele)] = self.Problem.GoverningEquations.Flux(UU, delevert, Flags)
+
# ---- Computing the plain integral
# Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
# in the Dof's vertex ordering referential
diff --git a/SourceCode/_Solver/SpatialSchemes/CG_LFX.py b/SourceCode/_Solver/SpatialSchemes/CG_LFX.py
index 3079798..18d3ca4 100644
--- a/SourceCode/_Solver/SpatialSchemes/CG_LFX.py
+++ b/SourceCode/_Solver/SpatialSchemes/CG_LFX.py
@@ -18,7 +18,7 @@ SchemeName = "Continuous Galerkin"
# Mapping to the mesh that should be loaded
class AssociatedMeshType():
""" .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Class that only gives out the mesh type that should be loaded when tackling the problem
with this scheme, given the wished variational order. Access the associated mesh properties
through the inner variables
@@ -27,25 +27,25 @@ class AssociatedMeshType():
- |bs| **MeshOrder** *(integer)* -- the mesh order
|
-
+
.. Note::
-
+
This is only a product-type class: no method is available
-
+
|
"""
#### Automatic initialisation routine ####
def __init__(self, SchemeParams = [1]):
""" .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
- Args:
+ Args:
SchemeParams (string list, optional): the parameters of the spatial scheme wished by the user
-
+
"""
-
+
# Extracting the schemes properties that are required to determine the Mesh
VariationalOrder = int(SchemeParams[0])
-
+
# Selecting the corresponding mesh type
self.MeshType = "CG"
self.MeshOrder = VariationalOrder
@@ -57,16 +57,16 @@ class AssociatedMeshType():
# ==============================================================================
class Scheme():
""".. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
-
+
+
Class furnishing all the schemes tools and registers the scheme's wished
parameters. The main iteration is accessible with the method "Iteration".
-
+
Repeats the instance's initialisation arguments as fields with identical names
within the structure.
|
-
+
"""
#### Automatic initialisation routine ####
@@ -76,13 +76,13 @@ class Scheme():
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh instance
Params (string list): the scheme's parameters as wished by the user
-
+
|
-
+
.. rubric:: Methods
-
+
"""
-
+
self.Problem = Problem
self.Params = Params
self.Mesh = Mesh
@@ -90,8 +90,8 @@ class Scheme():
#### Emulates a flux computation as done in the iteration routine ####
def ComputeFlux(self, FluidFlag, U, LSValues):
""" Emulates a flux computation as done in the iteration routine in order to be used in the DeC subtimestep evaluation.
-
- Args:
+
+ Args:
FluidFlag (numpy float array): buffer vector containing the current fluid flags at each dof (NbDofs)
U (numpy float array): buffer containing the current state values at each dof (NbVars x NbDofs)
LSValues (numpy float array): buffer containing the current FluidSpotters values at each dof(NbFluids x NbDofs)
@@ -100,67 +100,72 @@ class Scheme():
fluxes (numpy float (multiD)array): fluxes in the format required by the Iteration routine
---------------------------------------------------- """
-
- # -------------------------- Initialisations -------------------------------------------------
+
+ # -------------------------- Initialisations -------------------------------------------------
# Defining the output fluxes (inner and boundary ones)
# The very first coordinate corresponds to the inner flux. Only Flux[0, :, 0, :, :, :] is filled
# The second coordinate corresponds to the inner flux. Only Flux[1, :, :, :, :, 0:3] is filled
- Flux = np.zeros((2, self.Mesh.NbInnerElements, 3, 2, self.Problem.GoverningEquations.NbVariables, 12))
+ Flux = np.zeros((3, self.Mesh.NbInnerElements, 3, 2, self.Problem.GoverningEquations.NbVariables, np.max([12, *[len(dele) for dele in self.Mesh.ElementsDofs]])))
- # Retrieving the scheme order
+ # Retrieving the scheme order
Order = int(self.Params[0])
-
+
# ------------- Flux computation at the control's volume interface ----------------------------
# Spanning on each element
for i in range(self.Mesh.NbInnerElements):
# ---- Getting the element and state informations
- # Retrieving the physical information on the spanned element's vertices
+ # Retrieving the physical information on the spanned element's vertices
ele = self.Mesh.InnerElements[i]
vol = self.Mesh.InnerElementsVolume[i]
n = self.Mesh.ElementsBoundaryWiseNormals[i]
-
+
# Retrieving the information on the local degree of freedom
- dele = np.array(self.Mesh.ElementsDofs[i])
- LS = LSValues[:, dele]
- UU = U[:, dele]
+ dele = np.array(self.Mesh.ElementsDofs[i])
+ LS = LSValues[:, dele]
+ UU = U[:, dele]
+ Flags = np.argmax(LS, axis=0)
# Get the coordinates of the vertices in the physical order
- elevert = self.Mesh.CartesianPoints[ele,:]
+ elevert = self.Mesh.CartesianPoints[ele,:]
+ delevert = self.Mesh.DofsXY[dele,:]
- # ---- Computing the plain integral
- # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+ # ---- Computing the flux at the Dofs directly
+ Flux[2, i, 0, :, :, 0:len(dele)] = self.Problem.GoverningEquations.Flux(UU, delevert, Flags)
+
+ # ---- Computing the plain integral
+ # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
# in the Dof's vertex ordering referential
bqp, weights = QR.InnerQuadrature()
qp = BT.GetCartesianCoordinates(bqp, elevert)
-
+
# Evaluate the basis function in the order of the Dofs
- EvaluatedBase = BF.SimplicialBasis("B", Order, bqp)
-
+ EvaluatedBase = BF.SimplicialBasis("B", Order, bqp)
+
# Reconstructing the solution and its gradient at the quadrature point
- Uloc = np.dot(UU, EvaluatedBase.T)
- LSloc = np.dot(LS, EvaluatedBase.T)
+ Uloc = np.dot(UU, EvaluatedBase.T)
+ LSloc = np.dot(LS, EvaluatedBase.T)
+ FlagsLoc = np.argmax(LSloc, axis=0)
# Compute the flux according to the fluid's flag
- Flags = np.argmax(LSloc, axis=0)
- Flux[0, i, 0, :, :] = self.Problem.GoverningEquations.Flux(Uloc, qp, Flags)
+ Flux[0, i, 0, :, :, 0:12] = self.Problem.GoverningEquations.Flux(Uloc, qp, FlagsLoc)
for face in range(len(ele)):
- # Get the coordinates of the quadrature points on the edge, flipping the
+ # Get the coordinates of the quadrature points on the edge, flipping the
# points to match the Dof's order with the physical vertices
qbbp, weights = QR.BoundaryQuadrature(face)
bqp = BT.GetCartesianCoordinates(qbbp, elevert)
-
+
# Compute each basis function value at those points
- EvaluatedBase = BF.SimplicialBasis("B", Order, qbbp)
-
+ EvaluatedBase = BF.SimplicialBasis("B", Order, qbbp)
+
# Reconstruct the solution's value at this point
Uloc = np.dot(UU, EvaluatedBase.T)
LSloc = np.dot(LS, EvaluatedBase.T)
-
+
# Evaluate the flux there
- Flags = np.argmax(LSloc, axis=0)
- Flux[1,i,face,:,:,0:3] = self.Problem.GoverningEquations.Flux(Uloc, bqp, Flags)
-
+ FlagsLoc = np.argmax(LSloc, axis=0)
+ Flux[1,i,face,:,:,0:3] = self.Problem.GoverningEquations.Flux(Uloc, bqp, FlagsLoc)
+
# Returning the full array
return(Flux)
@@ -168,188 +173,156 @@ class Scheme():
def Iteration(self, Solution, fluxes, i, du=0, dt=1):
""" Main iteration of the scheme, implementing the most stupid scheme
you can think of.
-
+
Args:
Solution (solution structure): structure containing the current solution's values to iterate
fluxes (numpy (multiD)array): pre-computed fluxes at the points of interest. For this scheme, access with fluxes[element, face, coordinate(fx or fy), variable, pointindex]
i (integer): the index of the considered element within which the partial residuals will be computed
- du (float numpy array): (optional) when using DeC, the difference in the time iteration
+ du (float numpy array): (optional) when using DeC, the difference in the time iteration
dt (float): (optional) when using DeC, the full time step
-
+
Returns:
Resu (float numpy array): the computed residuals (NbVars x NbDofs)
-
+
---------------------------------------------------------------------"""
# -------- Initialisation ------------------------------------------------
- # Retrieving the scheme order
+ # Retrieving the scheme order
Order = int(self.Params[0])
-
+
# Getting the dofs contaned in the considered element
dele = np.array(self.Mesh.ElementsDofs[i])
NbDofs = len(dele)
# Initialising the mollification vector (residuals)
Resu = np.zeros([np.shape(Solution.Sol)[0], NbDofs])
-
+
# Retrieving the inner and boundary fluxes
Flx = fluxes[0, :, 0, :, :]
Fln = fluxes[1,:,:,:,:,0:3]
-
+
# ------- Getting the CG residuals for each element -------------------------
-
+
# ---- Getting the element and state informations
- # Retrieving the physical information on the spanned element's vertices
+ # Retrieving the physical information on the spanned element's vertices
ele = self.Mesh.InnerElements[i]
vol = self.Mesh.InnerElementsVolume[i]
n = self.Mesh.ElementsBoundaryWiseNormals[i]
-
+
# Retrieving the information on the local degree of freedom
- LS = Solution.LSValues[:, dele]
- U = Solution.Sol[:, dele]
- dU = du[:, dele]
-
- # ---- Computing the plain integral
- # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+ LS = Solution.LSValues[:, dele]
+ U = Solution.Sol[:, dele]
+ dU = du[:, dele]
+
+ # ---- Computing the plain integral
+ # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
# in the Dof's vertex ordering referential
bqp, weights = QR.InnerQuadrature()
-
+
# Evaluate the basis function in the order of the Dofs
EvaluatedBase = BF.SimplicialBasis("B", Order, bqp)
- EvaluatedGradients = BF.SimplicialGradients("B", Order, bqp, n, vol)
-
- # Reconstructing the solution and its gradient at the quadrature point
- Uloc = np.dot(U, EvaluatedBase.T)
+ EvaluatedGradients = BF.SimplicialGradients("B", Order, bqp, n, vol)
+
+ # Reconstructing the time iteration buffer at the quadrature point
dUloc = np.dot(dU, EvaluatedBase.T)
- LSloc = np.dot(LS, EvaluatedBase.T)
-
+
# Compute the flux according to the fluid's flag
Flux = Flx[i]
-
+
# Complete the resodual at each Dof by the relevant quadrature value
for dof in range(NbDofs):
# Evaluating the flux at the quadrature point for the reconstructed flux variable
qt = np.sum([weights[q]*(EvaluatedGradients[dof,q,0]*Flux[0,:,q]+EvaluatedGradients[dof,q,1]*Flux[1,:,q]) for q in range(len(weights))], axis=0)
DEC = np.sum( weights[:]*EvaluatedBase[:,dof]*dUloc, axis=1)
-
+
# Updating the residual at the spanned Dof by the local contribution
Resu[:,dof] += (DEC - dt*qt)*vol
- # ---- Computing the boundary integral
- # Spanning each edge in the physical element
+ # ---- Computing the boundary integral
+ # Spanning each edge in the physical element
for face in range(len(ele)):
- # Get the coordinates of the quadrature points on the edge, flipping the
+ # Get the coordinates of the quadrature points on the edge, flipping the
# points to match the Dof's order with the physical vertices
qbbp, weights = QR.BoundaryQuadrature(face)
-
+
# Compute each basis function value at those points
# (get the order of the basis in the order of the physical vertices:
# the dele have to be ordered in the same output as the evaluated base:
- # check also for higher order the order of the dele given in the mesh and
+ # check also for higher order the order of the dele given in the mesh and
# the order of the basis function given here)
EvaluatedBase = BF.SimplicialBasis("B", Order, qbbp) # To move as qbbp are not tabulated at the right dof location,, tabultaed on physical and not dof ordered
-
- # Reconstruct the solution's value at this point
- Uloc = np.dot(U, EvaluatedBase.T)
- LSloc = np.dot(LS, EvaluatedBase.T)
-
+
# Evaluate the flux there
Flux = Fln[i,face]
-
+
# Compute the flux's normal contribution to each dof
Qt = (Flux[0,:,:]*n[face,0]+Flux[1,:,:]*n[face,1])
WeigthedBaseWiseQt = np.array([weights[pt]*np.array([Qt[:,pt]]).T*np.array([EvaluatedBase[pt,:]]) for pt in range(len(weights))])
PonderedQt = np.sum(WeigthedBaseWiseQt, axis=0)
-
+
# Adjust the residuals
Resu[:,:] += dt*PonderedQt
# ------- Editing the CG residuals to become LFX ------------------------
# Selecting the barycenter of the element as evaluation point
Bcoords = self.Mesh.DualCenters[i][-1]
-
+
# Computing the average state on the spanned element
Ub = np.array([np.mean(U, axis=1)]).T
-
+
# Preparing the evaluation of the spectral radius at the centre of gravity of the element
# for each dof value and each edge's normal direction
Uu = np.tile(U, [1,np.shape(n)[0]])
nn = np.array([nc for nc in n for cp in range(NbDofs)])
Coor = np.zeros((np.shape(Uu)[1],2))
- Coor[:,:] = Bcoords
-
+ Coor[:,:] = Bcoords
+
# Computing the spectral radius at each vertex and extracting the maximum
FluidIndex = np.argmax(LS, axis=0)
sprctrad = self.Problem.GoverningEquations.SpectralRadius(Uu, FluidIndex, nn, Coor)
+
+ # Dirty hack to preserve the pressure contact
+
Resu[:,:] += dt*np.max(sprctrad)*(U[:,:]-Ub[:])
-
-
- # -------- Limitting the Fluid's spotter residuals ------------------------
- NbVars, NbDofs = np.shape(Resu)
- Lresu = np.zeros(np.shape(Resu))
- safe = 1e-42
-
- Phi = np.array([np.sum(Resu, axis = 1)]).T
- ModInd = np.where(Phi>safe)[0]
- NodInd = np.setdiff1d(range(NbVars), ModInd)
-
- # Modify the residuals for the variables corresponding to the indices of ModInd
- if ModInd.size>0:
- # Using the definition of the Psi limiter and retrieving the quantities from the
- # current residuals. Note that in mol, the abs is here to compensate
- # the negative zero python convention used in this case
- buff = Resu[ModInd, :]/Phi[ModInd]
- mol = np.abs(buff*(buff>0))
- den = np.array([np.sum(mol, axis=1)+safe]).T
-
- # Getting the new distribution coefficients and the blending ratio
- bet = mol/den
- tet = np.abs(Phi[ModInd])/(np.array([np.sum(np.abs(Resu[ModInd, :]), axis=1)]).T+safe)
-
- # Blending the scheme
- Lresu[ModInd,:] = (1-tet)*bet*Phi[ModInd] + tet*Resu[ModInd, :]
- # Not touching the residuals for the variables corresponding to the indices of ModInd
- if NodInd.size>0: Lresu[NodInd, :] = Resu[NodInd, :]
-
# Returning the partial residuals
- return(Lresu)
+ return(Resu)
#### Routine that maps the values from the Dofs to the Physical vertices of the mesh ####
def ReconstructSolutionAtVertices(self, Solution):
- """ Routine that maps the values from the Dofs to the Physical vertices of the mesh
+ """ Routine that maps the values from the Dofs to the Physical vertices of the mesh
of the FluidSpotter, solution and flags.
-
+
Args:
Solution (Solution): the currently being computed solution
-
+
Returns:
None: fills direclty the RSol and RFluidFlag values in the data structure
-
+
.. Note::
-
+
This is only for later plotting purposes.
-
+
----------------------------------------------------------------------------------- """
-
+
# -------------- Initialising the reconstructed values ----------------------------
Solution.RFluidFlag = np.zeros(self.Mesh.NbMeshPoints)
Solution.RLSValues = np.zeros((np.shape(Solution.LSValues)[0], self.Mesh.NbMeshPoints))
Solution.RSol = np.zeros((np.shape(Solution.Sol)[0], self.Mesh.NbMeshPoints))
-
+
# ------------- Reconstructing the values -----------------------------------------
# Extract the cartesian coordinates of the Dofs
xyDofs = self.Mesh.DofsXY[:,:]
-
+
# Spanning all the mesh points where to fill the values
for i in range(self.Mesh.NbMeshPoints):
-
- # Averaging all the Dofs values located at the corresponding physical vertex
+
+ # Averaging all the Dofs values located at the corresponding physical vertex
index = self.Mesh.Vertex2Dofs[i]
RS = np.mean(Solution.LSValues[:,index], axis=1)
RSO = np.mean(Solution.Sol[:,index], axis=1)
FF = int(np.mean(Solution.FluidFlag[index]))
-
+
# Fill the reconstructed solution upon the found index
Solution.RLSValues[:,i] = RS
Solution.RFluidFlag[i] = FF
diff --git a/SourceCode/_Solver/SpatialSchemes/CG_Primary.py b/SourceCode/_Solver/SpatialSchemes/CG_Primary.py
new file mode 100644
index 0000000..348a474
--- /dev/null
+++ b/SourceCode/_Solver/SpatialSchemes/CG_Primary.py
@@ -0,0 +1,327 @@
+# ==============================================================================
+# Preliminary imports
+# ==============================================================================
+import numpy as np
+import copy
+
+# Import custom modules
+import _Solver.QuadratureRules.HandBased as QR
+import BarycentricBasisFunctions as BF
+import BarycentricTools as BT
+
+# ==============================================================================
+# Parameter (and reader) class of scheme properties and runtime informations
+# ==============================================================================
+# Name of the module to be externally accessible
+SchemeName = "Continuous Galerkin"
+
+# Mapping to the mesh that should be loaded
+class AssociatedMeshType():
+ """ .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
+
+ Class that only gives out the mesh type that should be loaded when tackling the problem
+ with this scheme, given the wished variational order. Access the associated mesh properties
+ through the inner variables
+
+ - |bs| **MeshType** *(string)* -- the mesh type according to fenics's classification
+ - |bs| **MeshOrder** *(integer)* -- the mesh order
+
+ |
+
+ .. Note::
+
+ This is only a product-type class: no method is available
+
+ |
+ """
+
+ #### Automatic initialisation routine ####
+ def __init__(self, SchemeParams = [1]):
+ """ .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
+ Args:
+ SchemeParams (string list, optional): the parameters of the spatial scheme wished by the user
+
+ """
+
+ # Extracting the schemes properties that are required to determine the Mesh
+ VariationalOrder = int(SchemeParams[0])
+
+ # Selecting the corresponding mesh type
+ self.MeshType = "CG"
+ self.MeshOrder = VariationalOrder
+ self.ControlVolumes = "Dual"
+
+
+# ==============================================================================
+# Class furnishing all the scheme tools
+# ==============================================================================
+class Scheme():
+ """.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
+
+
+ Class furnishing all the schemes tools and registers the scheme's wished
+ parameters. The main iteration is accessible with the method "Iteration".
+
+ Repeats the instance's initialisation arguments as fields with identical names
+ within the structure.
+
+ .. note::
+
+ This scheme is the CG one where the variables interpolated are rho, rhoU,
+ rhoV, p instead of the conservative ones, in order to preserve the pressure
+ contacts if any
+
+ |
+
+ """
+
+ #### Automatic initialisation routine ####
+ def __init__(self, Problem, Mesh, Params = []):
+ """
+ Args:
+ Problem (Problem): the considered problem
+ Mesh (MeshStructure): the considered mesh instance
+ Params (string list): the scheme's parameters as wished by the user
+
+ |
+
+ .. rubric:: Methods
+
+ """
+
+ self.Problem = Problem
+ self.Params = Params
+ self.Mesh = Mesh
+
+ #### Emulates a flux computation as done in the iteration routine ####
+ def ComputeFlux(self, FluidFlag, U, LSValues):
+ """ Emulates a flux computation as done in the iteration routine in order to be used in the DeC subtimestep evaluation.
+
+ Args:
+ FluidFlag (numpy float array): buffer vector containing the current fluid flags at each dof (NbDofs)
+ U (numpy float array): buffer containing the current state values at each dof (NbVars x NbDofs)
+ LSValues (numpy float array): buffer containing the current FluidSpotters values at each dof(NbFluids x NbDofs)
+
+ Returns:
+ fluxes (numpy float (multiD)array): fluxes in the format required by the Iteration routine
+
+ ---------------------------------------------------- """
+
+ # -------------------------- Initialisations -------------------------------------------------
+ # Defining the output fluxes (inner and boundary ones)
+ # The very first coordinate corresponds to the inner flux. Only Flux[0, :, 0, :, :, :] is filled
+ # The second coordinate corresponds to the inner flux. Only Flux[1, :, :, :, :, 0:3] is filled
+ Flux = np.zeros((3, self.Mesh.NbInnerElements, 3, 2, self.Problem.GoverningEquations.NbVariables, np.max([12, *[len(dele) for dele in self.Mesh.ElementsDofs]])))
+
+ # Retrieving the scheme order
+ Order = int(self.Params[0])
+
+ # ------------- Flux computation at the control's volume interface ----------------------------
+ # Spanning on each element
+ for i in range(self.Mesh.NbInnerElements):
+ # ---- Getting the element and state informations
+ # Retrieving the physical information on the spanned element's vertices
+ ele = self.Mesh.InnerElements[i]
+ vol = self.Mesh.InnerElementsVolume[i]
+ n = self.Mesh.ElementsBoundaryWiseNormals[i]
+
+ # Retrieving the information on the local degree of freedom
+ dele = np.array(self.Mesh.ElementsDofs[i])
+ LS = LSValues[:, dele]
+ UU = U[:, dele]
+ Flags = np.argmax(LS, axis=0)
+
+ # Get the coordinates of the vertices in the physical order
+ elevert = self.Mesh.CartesianPoints[ele,:]
+ delevert = self.Mesh.DofsXY[dele,:]
+
+ # Computing the flux at the Dofs directly
+ Flux[2, i, 0, :, :, 0:len(dele)] = self.Problem.GoverningEquations.Flux(UU, delevert, Flags)
+
+ # Converting the conservative variables to the primary ones before the interpolation
+ # (helps to preserve the pressure contacts, comment if E has to be preserved as well)
+ UU = self.Problem.GoverningEquations.ConservativeToPrimary(UU, Flags)
+
+
+ # ---- Computing the plain integral
+ # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+ # in the Dof's vertex ordering referential
+ bqp, weights = QR.InnerQuadrature()
+ qp = BT.GetCartesianCoordinates(bqp, elevert)
+
+ # Evaluate the basis function in the order of the Dofs
+ EvaluatedBase = BF.SimplicialBasis("B", Order, bqp)
+
+ # Reconstructing the solution and its gradient at the quadrature point
+ Uloc = np.dot(UU, EvaluatedBase.T)
+ LSloc = np.dot(LS, EvaluatedBase.T)
+
+ # Reconstructing back the conservative variables at the wished points
+ # (helps to preserve the pressure contacts, comment if E has to be preserved as well)
+ FlagsLoc = np.argmax(LSloc, axis=0)
+ Uloc = self.Problem.GoverningEquations.PrimaryToConservative(Uloc, FlagsLoc)
+
+ # Compute the flux according to the fluid's flag
+ Flux[0, i, 0, :, :, 0:12] = self.Problem.GoverningEquations.Flux(Uloc, qp, FlagsLoc)
+
+ for face in range(len(ele)):
+ # Get the coordinates of the quadrature points on the edge, flipping the
+ # points to match the Dof's order with the physical vertices
+ qbbp, weights = QR.BoundaryQuadrature(face)
+ bqp = BT.GetCartesianCoordinates(qbbp, elevert)
+
+ # Compute each basis function value at those points
+ EvaluatedBase = BF.SimplicialBasis("B", Order, qbbp)
+
+ # Reconstruct the solution's value at this point
+ Uloc = np.dot(UU, EvaluatedBase.T)
+ LSloc = np.dot(LS, EvaluatedBase.T)
+
+ # Reconstructing back the conservative variables at the wished points
+ # (helps to preserve the pressure contacts, comment if E has to be preserved as well)
+ FlagsLoc = np.argmax(LSloc, axis=0)
+ Uloc = self.Problem.GoverningEquations.PrimaryToConservative(Uloc, FlagsLoc)
+
+ # Evaluate the flux there
+ Flux[1,i,face,:,:,0:3] = self.Problem.GoverningEquations.Flux(Uloc, bqp, FlagsLoc)
+
+ # Returning the full array
+ return(Flux)
+
+ #### Main routine, definining the scheme ####
+ def Iteration(self, Solution, fluxes, i, du=0, dt=1):
+ """ Main iteration of the scheme, implementing the most stupid scheme
+ you can think of.
+
+ Args:
+ Solution (solution structure): structure containing the current solution's values to iterate
+ fluxes (numpy (multiD)array): pre-computed fluxes at the points of interest. For this scheme, access with fluxes[element, face, coordinate(fx or fy), variable, pointindex]
+ i (integer): the index of the considered element within which the partial residuals will be computed
+ du (float numpy array): (optional) when using DeC, the difference in the time iteration
+ dt (float): (optional) when using DeC, the full time step
+
+ Returns:
+ Resu (float numpy array): the computed residuals (NbVars x NbDofs)
+
+ ---------------------------------------------------------------------"""
+
+ # -------- Initialisation ------------------------------------------------
+ # Retrieving the scheme order
+ Order = int(self.Params[0])
+
+ # Getting the dofs contaned in the considered element
+ dele = np.array(self.Mesh.ElementsDofs[i])
+ NbDofs = len(dele)
+
+ # Initialising the mollification vector (residuals)
+ Resu = np.zeros([np.shape(Solution.Sol)[0], NbDofs])
+
+ # Retrieving the inner and boundary fluxes
+ Flx = fluxes[0, :, 0, :, :]
+ Fln = fluxes[1,:,:,:,:,0:3]
+
+ # ------- Getting the CG residuals for each element -------------------------
+
+ # ---- Getting the element and state informations
+ # Retrieving the physical information on the spanned element's vertices
+ ele = self.Mesh.InnerElements[i]
+ vol = self.Mesh.InnerElementsVolume[i]
+ n = self.Mesh.ElementsBoundaryWiseNormals[i]
+
+ # Retrieving the information on the local degree of freedom
+ LS = Solution.LSValues[:, dele]
+ U = Solution.Sol[:, dele]
+ dU = du[:, dele]
+
+ # ---- Computing the plain integral
+ # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+ # in the Dof's vertex ordering referential
+ bqp, weights = QR.InnerQuadrature()
+
+ # Evaluate the basis function in the order of the Dofs
+ EvaluatedBase = BF.SimplicialBasis("B", Order, bqp)
+ EvaluatedGradients = BF.SimplicialGradients("B", Order, bqp, n, vol)
+
+ # Reconstructing the time iteration buffer at the quadrature point
+ dUloc = np.dot(dU, EvaluatedBase.T)
+
+ # Compute the flux according to the fluid's flag
+ Flux = Flx[i]
+
+ # Complete the resodual at each Dof by the relevant quadrature value
+ for dof in range(NbDofs):
+ # Evaluating the flux at the quadrature point for the reconstructed flux variable
+ qt = np.sum([weights[q]*(EvaluatedGradients[dof,q,0]*Flux[0,:,q]+EvaluatedGradients[dof,q,1]*Flux[1,:,q]) for q in range(len(weights))], axis=0)
+ DEC = np.sum( weights[:]*EvaluatedBase[:,dof]*dUloc, axis=1)
+
+ # Updating the residual at the spanned Dof by the local contribution
+ Resu[:,dof] += (DEC - dt*qt)*vol
+
+ # ---- Computing the boundary integral
+ # Spanning each edge in the physical element
+ for face in range(len(ele)):
+ # Get the coordinates of the quadrature points on the edge, flipping the
+ # points to match the Dof's order with the physical vertices
+ qbbp, weights = QR.BoundaryQuadrature(face)
+
+ # Compute each basis function value at those points
+ # (get the order of the basis in the order of the physical vertices:
+ # the dele have to be ordered in the same output as the evaluated base:
+ # check also for higher order the order of the dele given in the mesh and
+ # the order of the basis function given here)
+ EvaluatedBase = BF.SimplicialBasis("B", Order, qbbp) # To move as qbbp are not tabulated at the right dof location,, tabultaed on physical and not dof ordered
+
+ # Evaluate the flux there
+ Flux = Fln[i,face]
+
+ # Compute the flux's normal contribution to each dof
+ Qt = (Flux[0,:,:]*n[face,0]+Flux[1,:,:]*n[face,1])
+ WeigthedBaseWiseQt = np.array([weights[pt]*np.array([Qt[:,pt]]).T*np.array([EvaluatedBase[pt,:]]) for pt in range(len(weights))])
+ PonderedQt = np.sum(WeigthedBaseWiseQt, axis=0)
+
+ # Adjust the residuals
+ Resu[:,:] += dt*PonderedQt
+
+
+ # Returning the partial residuals
+ return(Resu)
+
+ #### Routine that maps the values from the Dofs to the Physical vertices of the mesh ####
+ def ReconstructSolutionAtVertices(self, Solution):
+ """ Routine that maps the values from the Dofs to the Physical vertices of the mesh
+ of the FluidSpotter, solution and flags.
+
+ Args:
+ Solution (Solution): the currently being computed solution
+
+ Returns:
+ None: fills direclty the RSol and RFluidFlag values in the data structure
+
+ .. Note::
+
+ This is only for later plotting purposes.
+
+ ----------------------------------------------------------------------------------- """
+
+ # -------------- Initialising the reconstructed values ----------------------------
+ Solution.RFluidFlag = np.zeros(self.Mesh.NbMeshPoints)
+ Solution.RLSValues = np.zeros((np.shape(Solution.LSValues)[0], self.Mesh.NbMeshPoints))
+ Solution.RSol = np.zeros((np.shape(Solution.Sol)[0], self.Mesh.NbMeshPoints))
+
+ # ------------- Reconstructing the values -----------------------------------------
+ # Extract the cartesian coordinates of the Dofs
+ xyDofs = self.Mesh.DofsXY[:,:]
+
+ # Spanning all the mesh points where to fill the values
+ for i in range(self.Mesh.NbMeshPoints):
+
+ # Averaging all the Dofs values located at the corresponding physical vertex
+ index = self.Mesh.Vertex2Dofs[i]
+ RS = np.mean(Solution.LSValues[:,index], axis=1)
+ RSO = np.mean(Solution.Sol[:,index], axis=1)
+ FF = int(np.mean(Solution.FluidFlag[index]))
+
+ # Fill the reconstructed solution upon the found index
+ Solution.RLSValues[:,i] = RS
+ Solution.RFluidFlag[i] = FF
+ Solution.RSol[:,i] = RSO
diff --git a/Templates/Template_Problem.py b/Templates/Template_Problem.py
index 9264323..134b8e6 100644
--- a/Templates/Template_Problem.py
+++ b/Templates/Template_Problem.py
@@ -47,9 +47,9 @@ class ProblemData():
# Column: Index of the defined fluid
# Content: Index of the boundary condition to apply if the considered fluid touches the given boundary segment
self.BoundaryConditions = np.array(\
- [[0, 0],\
- [0, 0],\
- [0, 0],\
- [0, 0],\
- [0, 0],\
- [0, 0]])
+ [[(0,), (0,)],\
+ [(0,), (0,)],\
+ [(0,), (0,)],\
+ [(0,), (0,)],\
+ [(0,), (0,)],\
+ [(0,), (0,)]])
--
GitLab