From a2477507fe1404aacb8540caaa934f1043ddbc3a Mon Sep 17 00:00:00 2001
From: elemel <elise.lemeledo@math.uzh.ch>
Date: Mon, 27 Apr 2020 18:58:51 +0200
Subject: [PATCH] added predefined tests and tentative update subdomains
---
...1Fluid.py => Advection_Rotation_1Fluid.py} | 4 +-
.../Advection_Translation_1Fluid.py | 41 +++
...ds.py => Advection_Translation_2Fluids.py} | 8 +-
.../Advection_Translation_Obstacle_2Fluids.py | 41 +++
.../Euler_SODCste_2Fluid.py | 41 +++
.../PredefinedProblems/Euler_SODSOD_2Fluid.py | 41 +++
.../{SOD_1Fluid.py => Euler_SOD_1Fluid.py} | 2 +-
...id.py => Euler_StationaryVortex_1Fluid.py} | 2 +-
.../CustomSettingsTemplate.txt | 4 +-
RunPostProcess.py | 10 +-
RunProblem.py | 9 +-
SourceCode/Mappers.py | 30 +-
SourceCode/PostProcess.py | 58 +++-
SourceCode/PostProcessing/Exports.py | 4 +-
SourceCode/PostProcessing/Plots_Matplotlib.py | 22 +-
SourceCode/PostProcessing/Plots_Plotly.py | 24 +-
.../GoverningEquations/EulerEquations.py | 10 +-
...dvection.py => LinearAdvectionRotation.py} | 2 +-
.../LinearAdvectionTranslation.py | 327 ++++++++++++++++++
SourceCode/Solve.py | 21 +-
.../{NaiveLS.py => NaiveLS_CellCentredLFX.py} | 132 +++++--
.../NaiveLS_VertexCentredLFX.py | 235 +++++++++++++
22 files changed, 971 insertions(+), 97 deletions(-)
rename PredefinedTests/PredefinedProblems/{Advection_1Fluid.py => Advection_Rotation_1Fluid.py} (94%)
create mode 100644 PredefinedTests/PredefinedProblems/Advection_Translation_1Fluid.py
rename PredefinedTests/PredefinedProblems/{Advection_2Fluids.py => Advection_Translation_2Fluids.py} (90%)
create mode 100644 PredefinedTests/PredefinedProblems/Advection_Translation_Obstacle_2Fluids.py
create mode 100644 PredefinedTests/PredefinedProblems/Euler_SODCste_2Fluid.py
create mode 100644 PredefinedTests/PredefinedProblems/Euler_SODSOD_2Fluid.py
rename PredefinedTests/PredefinedProblems/{SOD_1Fluid.py => Euler_SOD_1Fluid.py} (96%)
rename PredefinedTests/PredefinedProblems/{StationaryVortex_1Fluid.py => Euler_StationaryVortex_1Fluid.py} (95%)
rename SourceCode/ProblemDefinition/GoverningEquations/{LinearAdvection.py => LinearAdvectionRotation.py} (99%)
create mode 100644 SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py
rename SourceCode/_Solver/FluidSelectors/{NaiveLS.py => NaiveLS_CellCentredLFX.py} (51%)
create mode 100644 SourceCode/_Solver/FluidSelectors/NaiveLS_VertexCentredLFX.py
diff --git a/PredefinedTests/PredefinedProblems/Advection_1Fluid.py b/PredefinedTests/PredefinedProblems/Advection_Rotation_1Fluid.py
similarity index 94%
rename from PredefinedTests/PredefinedProblems/Advection_1Fluid.py
rename to PredefinedTests/PredefinedProblems/Advection_Rotation_1Fluid.py
index bcf716d..eb62aa2 100644
--- a/PredefinedTests/PredefinedProblems/Advection_1Fluid.py
+++ b/PredefinedTests/PredefinedProblems/Advection_Rotation_1Fluid.py
@@ -11,7 +11,7 @@ class ProblemData():
# ******************************************************************************
# Meshing
#*******************************************************************************
- self.ProblemID = "StationaryVortex" # Should match the module name
+ self.ProblemID = "Advection_Rotation_1Fluid"
self.MeshName = "BigSimpleSquare"
# ******************************************************************************
@@ -31,7 +31,7 @@ class ProblemData():
self.FluidInitialSubdomain = [] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
# Index of the initial conditions to apply on each subdomain
- self.InitialConditions = [0]
+ self.InitialConditions = [2]
# 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/Advection_Translation_1Fluid.py b/PredefinedTests/PredefinedProblems/Advection_Translation_1Fluid.py
new file mode 100644
index 0000000..85b51fb
--- /dev/null
+++ b/PredefinedTests/PredefinedProblems/Advection_Translation_1Fluid.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 = "Advection_Translation_1Fluid"
+ self.MeshName = "BigSimpleSquare"
+
+ # ******************************************************************************
+ # Equation
+ #*******************************************************************************
+ self.GoverningEquationsIndex = 3
+ self.EOSEquationsIndex = [1]
+
+ # *******************************************************************************
+ # Fluids properties
+ #********************************************************************************
+ # Select the indices speciying the fluids properties living on each subdomain
+ self.FluidIndex = [0]
+
+ # The fluid at index 0 is the complementary to the specified ones.
+ # Only need to precise the location of other fluids
+ self.FluidInitialSubdomain = [] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
+
+ # Index of the initial conditions to apply on each subdomain
+ self.InitialConditions = [2]
+
+ # Row: BcTag, coklumn, index of the condition to apply depending on the fluid that comes next to it
+ self.BoundaryConditions = np.array(\
+ [[0],\
+ [0],\
+ [0],\
+ [0]])
diff --git a/PredefinedTests/PredefinedProblems/Advection_2Fluids.py b/PredefinedTests/PredefinedProblems/Advection_Translation_2Fluids.py
similarity index 90%
rename from PredefinedTests/PredefinedProblems/Advection_2Fluids.py
rename to PredefinedTests/PredefinedProblems/Advection_Translation_2Fluids.py
index 3c1a5b6..272f65c 100644
--- a/PredefinedTests/PredefinedProblems/Advection_2Fluids.py
+++ b/PredefinedTests/PredefinedProblems/Advection_Translation_2Fluids.py
@@ -11,27 +11,27 @@ class ProblemData():
# ******************************************************************************
# Meshing
#*******************************************************************************
- self.ProblemID = "Advection_2Fluids" # Should match the module name
+ self.ProblemID = "Advection_Translation_2Fluids" # Should match the module name
self.MeshName = "BigSimpleSquare"
# ******************************************************************************
# Equation
#*******************************************************************************
- self.GoverningEquationsIndex = 2
+ self.GoverningEquationsIndex = 3
self.EOSEquationsIndex = [1, 1]
# *******************************************************************************
# Fluids properties
#********************************************************************************
# Select the indices speciying the fluids properties living on each subdomain
- self.FluidIndex = [0, 0]
+ 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 = [5] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
# Index of the initial conditions to apply on each subdomain
- self.InitialConditions = [0, 0]
+ self.InitialConditions = [2, 0]
# 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/Advection_Translation_Obstacle_2Fluids.py b/PredefinedTests/PredefinedProblems/Advection_Translation_Obstacle_2Fluids.py
new file mode 100644
index 0000000..d515bc5
--- /dev/null
+++ b/PredefinedTests/PredefinedProblems/Advection_Translation_Obstacle_2Fluids.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 = "Advection_Translation_Obstacle_2Fluids" # Should match the module name
+ self.MeshName = "LShapeWithHoleAndRoundCorner"
+
+ # ******************************************************************************
+ # Equation
+ #*******************************************************************************
+ self.GoverningEquationsIndex = 3
+ self.EOSEquationsIndex = [1, 1]
+
+ # *******************************************************************************
+ # 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.18, 0.5, 0.14, 0.14, 0.14, 1, 1)] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
+
+ # Index of the initial conditions to apply on each subdomain
+ self.InitialConditions = [1, 0]
+
+ # Row: BcTag, coklumn, index of the condition to apply depending on the fluid that comes next to it
+ self.BoundaryConditions = np.array(\
+ [[0, 0],\
+ [0, 0],\
+ [0, 0],\
+ [0, 0]])
diff --git a/PredefinedTests/PredefinedProblems/Euler_SODCste_2Fluid.py b/PredefinedTests/PredefinedProblems/Euler_SODCste_2Fluid.py
new file mode 100644
index 0000000..5f1d7cb
--- /dev/null
+++ b/PredefinedTests/PredefinedProblems/Euler_SODCste_2Fluid.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_SODCste" # Should match the module name
+ self.MeshName = "BigSimpleRectangle"
+
+ # ******************************************************************************
+ # Equation
+ #*******************************************************************************
+ self.GoverningEquationsIndex = 1
+ self.EOSEquationsIndex = [1, 1]
+
+ # *******************************************************************************
+ # 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 = [5] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
+
+ # Index of the initial conditions to apply on each subdomain
+ self.InitialConditions = [3, 0]
+
+ # Row: BcTag, coklumn, index of the condition to apply depending on the fluid that comes next to it
+ self.BoundaryConditions = np.array(\
+ [[0],\
+ [0],\
+ [0],\
+ [0]])
diff --git a/PredefinedTests/PredefinedProblems/Euler_SODSOD_2Fluid.py b/PredefinedTests/PredefinedProblems/Euler_SODSOD_2Fluid.py
new file mode 100644
index 0000000..8b466cb
--- /dev/null
+++ b/PredefinedTests/PredefinedProblems/Euler_SODSOD_2Fluid.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_SODSOD" # Should match the module name
+ self.MeshName = "BigSimpleRectangle"
+
+ # ******************************************************************************
+ # Equation
+ #*******************************************************************************
+ self.GoverningEquationsIndex = 1
+ self.EOSEquationsIndex = [1, 1]
+
+ # *******************************************************************************
+ # 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 = [5] #(2,0.6,1.5,0,0.5),(4,0.5, 1.0, 0.8, 0.8, 0.05, 0.15)] #(4,0.5, 1.0, 0.8, 0.5, 0.05, 0.15) [[[0.5, 0],[0.5, 0.5],[0, 0.5]],[[0.75,0.75],[0.75,1],[1,1]], 2, (2,3), (1,3,5,6)] (0,9,0,0.5,-0.5,0.5,0.5)
+
+ # Index of the initial conditions to apply on each subdomain
+ self.InitialConditions = [3, 3]
+
+ # Row: BcTag, coklumn, index of the condition to apply depending on the fluid that comes next to it
+ self.BoundaryConditions = np.array(\
+ [[0],\
+ [0],\
+ [0],\
+ [0]])
diff --git a/PredefinedTests/PredefinedProblems/SOD_1Fluid.py b/PredefinedTests/PredefinedProblems/Euler_SOD_1Fluid.py
similarity index 96%
rename from PredefinedTests/PredefinedProblems/SOD_1Fluid.py
rename to PredefinedTests/PredefinedProblems/Euler_SOD_1Fluid.py
index e0854f5..1fc3a34 100644
--- a/PredefinedTests/PredefinedProblems/SOD_1Fluid.py
+++ b/PredefinedTests/PredefinedProblems/Euler_SOD_1Fluid.py
@@ -11,7 +11,7 @@ class ProblemData():
# ******************************************************************************
# Meshing
#*******************************************************************************
- self.ProblemID = "SOD" # Should match the module name
+ self.ProblemID = "Euler_SOD" # Should match the module name
self.MeshName = "SimpleSquare"
# ******************************************************************************
diff --git a/PredefinedTests/PredefinedProblems/StationaryVortex_1Fluid.py b/PredefinedTests/PredefinedProblems/Euler_StationaryVortex_1Fluid.py
similarity index 95%
rename from PredefinedTests/PredefinedProblems/StationaryVortex_1Fluid.py
rename to PredefinedTests/PredefinedProblems/Euler_StationaryVortex_1Fluid.py
index 87f3eeb..d3e6037 100644
--- a/PredefinedTests/PredefinedProblems/StationaryVortex_1Fluid.py
+++ b/PredefinedTests/PredefinedProblems/Euler_StationaryVortex_1Fluid.py
@@ -11,7 +11,7 @@ class ProblemData():
# ******************************************************************************
# Meshing
#*******************************************************************************
- self.ProblemID = "StationaryVortex" # Should match the module name
+ self.ProblemID = "Euler_StationaryVortex" # Should match the module name
self.MeshName = "BigSimpleRectangle"
# ******************************************************************************
diff --git a/PredefinedTests/PredefinedSettings/CustomSettingsTemplate.txt b/PredefinedTests/PredefinedSettings/CustomSettingsTemplate.txt
index 77c98b6..6aab559 100644
--- a/PredefinedTests/PredefinedSettings/CustomSettingsTemplate.txt
+++ b/PredefinedTests/PredefinedSettings/CustomSettingsTemplate.txt
@@ -1,5 +1,5 @@
# Spatial scheme properties
-2 # Index of the scheme
+1 # Index of the scheme
1 0 # Scheme properties (see doc for more details, if any, the order should be first)
0 # Index of limiter (0=None)
0 0 # Limiter properties (see doc for more details)
@@ -13,7 +13,7 @@
# 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
+0.1 # CFL
1 # Tmax, in seconds
# Export properties
diff --git a/RunPostProcess.py b/RunPostProcess.py
index e4efb8d..a893c6e 100644
--- a/RunPostProcess.py
+++ b/RunPostProcess.py
@@ -22,11 +22,19 @@ import SourceCode.PostProcess as PostProcess
if __name__ == '__main__':
# Point to the problem you would like to solve, and define how you would like to tackle it
- ProblemDefinition = "Advection_2Fluids"
+ # "Euler_StationaryVortex_1Fluid.py"
+ # "Euler_SODSOD_2Fluid.py"
+ # "Euler_SODCste_2Fluid.py" "Euler_SOD_1Fluid.py"
+ # "Advection_Translation_Obstacle_2Fluids.py"
+ # "Advection_Translation_2Fluids.py"
+ # "Advection_Translation_1Fluid.py"
+ # "Advection_Rotation_1Fluid.py"
+ ProblemDefinition = "Advection_Translation_1Fluid"
SettingsChoice = "CustomSettingsTemplate"
MeshResolution = 20
# Post - processing tasks
+ PostProcess.Plot_DebugMatplotlib(ProblemDefinition, SettingsChoice, [0.0], MeshResolution)
#PostProcess.Plots_StaticMatplotlib(ProblemDefinition, SettingsChoice, [0.0], MeshResolution)
PostProcess.Plots_AnimationsPlotly(ProblemDefinition, SettingsChoice, MeshResolution)
#PostProcess.Plots_StaticPlotly(ProblemDefinition, SettingsChoice, [0.0, 0.01, 0.02, 0.03], MeshResolution)
diff --git a/RunProblem.py b/RunProblem.py
index 4d50b1d..88d3619 100644
--- a/RunProblem.py
+++ b/RunProblem.py
@@ -19,7 +19,14 @@ 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 = "Advection_2Fluids"
+ # "Euler_StationaryVortex_1Fluid."
+ # "Euler_SODSOD_2Fluid.py"
+ # "Euler_SODCste_2Fluid.py" "Euler_SOD_1Fluid.py"
+ # "Advection_Translation_Obstacle_2Fluids.py"
+ # "Advection_Translation_2Fluids.py"
+ # "Advection_Translation_1Fluid.py"
+ # "Advection_Rotation_1Fluid.py"
+ ProblemDefinition = "Advection_Rotation_1Fluid"
SettingsChoice = "CustomSettingsTemplate"
MeshResolution = 20
diff --git a/SourceCode/Mappers.py b/SourceCode/Mappers.py
index 6150503..0faed09 100644
--- a/SourceCode/Mappers.py
+++ b/SourceCode/Mappers.py
@@ -43,7 +43,8 @@ def Governing_Equations(id):
It should return the module itself"""
if id==1: return(GE.EulerEquations)
- elif id==2: return(GE.LinearAdvection)
+ elif id==2: return(GE.LinearAdvectionRotation)
+ elif id==3: return(GE.LinearAdvectionTranslation)
else: raise ValueError("Error in the problem definition:\n\
Wrong index for the wished governing equation.\n\
@@ -124,7 +125,6 @@ def Limiter(id, params = []):
Abortion.")
-
"""-----------------------------------------------------------------------------
Mapper function for the fluid selectors
-----------------------------------------------------------------------------"""
@@ -132,8 +132,8 @@ def Fluid_Selector(id):
""" Mapper that defines the fluid selector to use.
Returns the class containing the LS tools"""
- if id == 1: return(FS.NaiveLS)
- elif id == 2: return(FS.VoF.VoF)
+ if id == 1: return(FS.NaiveLS_CellCentredLFX)
+ elif id == 2: return(FS.NaiveLS_VertexCentredLFX)
else: raise ValueError("Error in the problem definition:\n\
Wrong index for the wished fluid selector.\n\
@@ -145,16 +145,14 @@ def Fluid_Selector(id):
Mapper function for the quadratures
-----------------------------------------------------------------------------"""
def Quadrature_Rule(QuadratureType, QuadratureOrder):
+ # Quadpy based quadrature
if QuadratureType==1:
- # Quadpy based quadrature
- if QuadratureOrder == 1:
- return(0)
- elif QuadratureOrder==2:
- return(0)
- else:
- raise ValueError("Automatic quadrature: wrong index for the quadrature order.")
- elif QuadratureType == 2:
- return(0)
- # By hand quadrature
- else:
- raise ValueError("Wrong index for the quadrature type.")
+ if QuadratureOrder == 1: return(0)
+ elif QuadratureOrder == 2: return(0)
+
+ else: raise ValueError("Automatic quadrature: wrong index for the quadrature order.")
+
+ # By hand quadrature
+ elif QuadratureType == 2: return(0)
+
+ else: raise ValueError("Wrong index for the quadrature type.")
diff --git a/SourceCode/PostProcess.py b/SourceCode/PostProcess.py
index b2a2b43..62cc34e 100644
--- a/SourceCode/PostProcess.py
+++ b/SourceCode/PostProcess.py
@@ -50,9 +50,9 @@ import PostProcessing.Exports as EX
# ==============================================================================
# Post processing automatic plotting routines
# ==============================================================================
-#### Plotting with plotly static visualisation tools ####
+#### Plotting with Matplotlib static visualisation tools ####
def Plots_StaticMatplotlib(*argv):
- """ Exports the usually welcome plots after computation, in plotly.
+ """ Exports the usually welcome plots after computation, in Matplotlib.
Inputs (in that order): ProblemName: string, name of the problem module to investigate, stored
in the "PredefinedProblems" folder (see the problem definition's
instructions for more details)
@@ -93,11 +93,6 @@ def Plots_StaticMatplotlib(*argv):
Mesh, ProblemData, Problem, Solution = EX.LoadSolution(ParametersId, argv[3], t)
# Batch the plotting routines
- """PMP.PlotFlags(Problem, Mesh, Solution, ParametersId)
- PMP.PlotSubdomains(Problem, Mesh, Solution, ParametersId)
- PMP.PlotLSubDomains(Problem, Mesh, Solution, ParametersId)
- PMP.PlotLSValuesTri(Problem, Mesh, Solution, ParametersId)
- """
for VariableId in range(Problem.GoverningEquations.NbVariables):
PMP.PlotVariableOnWholeDomainTri2D(Problem, Mesh, Solution, VariableId, ParametersId)
PMP.PlotVariableOnEachSubDomainTri2D(Problem, Mesh, Solution, VariableId, ParametersId)
@@ -108,7 +103,54 @@ def Plots_StaticMatplotlib(*argv):
PMP.PlotVariableOnWholeDomain(Problem, Mesh, Solution, VariableId, ParametersId)
PMP.PlotVariableOneEachSubDomain(Problem, Mesh, Solution, VariableId, ParametersId)
-
+#### Plotting with matpotlib some debug tools ####
+def Plot_DebugMatplotlib(*argv):
+ """ Exports the debug tools for multiphase computation, in Matplotlib.
+ Inputs (in that order): 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
+ time, float array-like, time points for which to process the plot
+ Resolution: mesh resolution
+
+ Output: 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.
+ ---------------------------------------------------------------------------------------------------------"""
+
+
+ # ---------------------- Checking the user input -------------------------------------
+ # Checking the validity of the user input
+ if len(argv) < 4:
+ raise ValueError("""
+ Not enough arguments have been given.
+ The problem cannot be initialised.
+ Check the given problem inputs. Abortion.""")
+ elif not (type(argv[0])==str and type(argv[1])==str and type(argv[3])==int):
+ raise ValueError("""
+ Wrong type of input argument.
+ The problem cannot be initialised.
+ Check the given problem inputs. Abortion.
+ """)
+
+
+ # ----------- Load the problem's definition and solver's attributes ----------------
+ # Register the name of the test for export purposes
+ ParametersId = argv[0]+"_"+argv[1]
+
+ # Loop on the wished times
+ for t in argv[2]:
+ # Load the solution
+ Mesh, ProblemData, Problem, Solution = EX.LoadSolution(ParametersId, argv[3], t)
+
+ # Batch the plotting routines
+ PMP.PlotFlags(Problem, Mesh, Solution, ParametersId)
+ PMP.PlotSubdomains(Problem, Mesh, Solution, ParametersId)
+ PMP.PlotLSubDomains(Problem, Mesh, Solution, ParametersId)
+ PMP.PlotLSValuesTri(Problem, Mesh, Solution, ParametersId)
+
#### Plotting with plotly static visualisation tools ####
def Plots_StaticPlotly(*argv):
""" Exports the usually welcome plots after computation, in plotly.
diff --git a/SourceCode/PostProcessing/Exports.py b/SourceCode/PostProcessing/Exports.py
index d8bad20..6edd3b2 100644
--- a/SourceCode/PostProcessing/Exports.py
+++ b/SourceCode/PostProcessing/Exports.py
@@ -70,7 +70,7 @@ def LoadSolution(TestCase, Resolution, t):
# Selecting the file to load
Folder = os.path.join(Pathes.SolutionPath, TestCase, "Res"+str(Resolution).replace(".","_"), "RawResults")
- FileName = os.path.join(Folder, ("Solution_t{0:5.3f}".format(t)).replace(".","_")+".sol")
+ FileName = os.path.join(Folder, ("Solution_t{0:8.6f}".format(t)).replace(".","_")+".sol")
# Loading the file
Mesh, ProblemData, Problem, Solution = np.load(FileName, allow_pickle=True)
@@ -97,5 +97,5 @@ def ExportSolution(TestCase, Mesh, ProblemData, Problem, Solution):
except: pass
# Exporting the file in a binary structure
- BinaryFileName = os.path.join(Folder, ("Solution_t{0:5.3f}".format(Solution.t)).replace(".","_")+".sol")
+ BinaryFileName = os.path.join(Folder, ("Solution_t{0:8.6f}".format(Solution.t)).replace(".","_")+".sol")
pickle.dump([Mesh, ProblemData, Problem, Solution], open(BinaryFileName, "wb"))
diff --git a/SourceCode/PostProcessing/Plots_Matplotlib.py b/SourceCode/PostProcessing/Plots_Matplotlib.py
index 1922e54..f2f941a 100644
--- a/SourceCode/PostProcessing/Plots_Matplotlib.py
+++ b/SourceCode/PostProcessing/Plots_Matplotlib.py
@@ -259,14 +259,14 @@ def PlotFlags(Problem, Mesh, Solution, ParametersId):
for i in range(0, Problem.NbFluids):
# Selecting all the Dofs falling in the subdomains
Fluid = Problem.FluidIndex[i]
- Index = np.where(Solution.FluidFlag == i)[0]
+ Index = np.where(Solution.RFluidFlag == i)[0]
# Plotting the markers
plt.plot(Mesh.CartesianPoints[Index,0], Mesh.CartesianPoints[Index,1],\
- "+", color = Problem.FluidProp[Fluid].color)
+ "+", color = Problem.FluidProp[i].color)
# Plotting the fluids's contours
- xgrid, ygrid, zgrid = GetMeshgridFom3VectorsMaskX(Mesh, Mesh.CartesianPoints[:,0],Mesh.CartesianPoints[:,1],Solution.FluidFlag)
+ xgrid, ygrid, zgrid = GetMeshgridFom3VectorsMaskX(Mesh, Mesh.CartesianPoints[:,0],Mesh.CartesianPoints[:,1],Solution.RFluidFlag)
plt.contourf(xgrid, ygrid, zgrid, levels=range(0,Problem.NbFluids), alpha=0.4, colors = [Fluid.color for Fluid in Problem.FluidProp])
@@ -357,7 +357,7 @@ def PlotLSubDomains(Problem, Mesh, Solution, ParametersId):
for i in range(len(Solution.Subdomains)): plotMultipoly(Solution.Subdomains[i], Problem.FluidProp[i].color, ax)
# Plot the contour of each fluid that is still present in the domain (problems of undertermined fluid will then pop up)
- Fluids = set(map(int, Solution.FluidFlag))
+ Fluids = set(map(int, Solution.RFluidFlag))
Fluids.discard(-1)
for i in Fluids:
ax.tricontour(Mesh.CartesianPoints[:,0], Mesh.CartesianPoints[:,1], Mesh.InnerElements, \
@@ -412,7 +412,7 @@ def PlotLSValuesTri(Problem, Mesh, Solution, ParametersId):
for Poly in np.array(Mesh.Domain): plt.plot([*np.array(Poly)[:,0],np.array(Poly)[0,0]], [*np.array(Poly)[:,1],np.array(Poly)[0,1]], "--k")
# Plot the contour of each fluid that is still present in the domain (problems of undertermined fluid will then pop up)
- Fluids = set(map(int, Solution.FluidFlag))
+ Fluids = set(map(int, Solution.RFluidFlag))
Fluids.discard(-1)
# Span each fluid
@@ -487,7 +487,7 @@ def PlotVariableOnWholeDomainTri2D(Problem, Mesh, Solution, VariableId, Paramete
ax = fig.subplots()
# Still plotting the contours of the LS
- for i in set(map(int, Solution.FluidFlag)):
+ for i in set(map(int, Solution.RFluidFlag)):
ax.tricontour(X, Y, Mesh.InnerElements, Solution.RLSValues[i,:], \
levels = [0], colors = "k", linewidths=2)
@@ -549,7 +549,7 @@ def PlotVariableOnEachSubDomainTri2D(Problem, Mesh, Solution, VariableId, Parame
ax = fig.subplots()
# Still plotting the contours of the LS
- for i in set(map(int, Solution.FluidFlag)):
+ for i in set(map(int, Solution.RFluidFlag)):
ax.tricontour(X, Y, Mesh.InnerElements, LS[i,:], \
colors = "k", levels = [0], linewidths=2)
@@ -681,7 +681,7 @@ def PlotVariableOneEachSubDomainTri(Problem, Mesh, Solution, VariableId, Paramet
# 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.FluidFlag))
+ 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]
@@ -766,7 +766,7 @@ def PlotVariableOnWholeDomain2D(Problem, Mesh, Solution, VariableId, ParametersI
ax.pcolormesh(xx, yy, zz, shading="gouraud", alpha=0.7, cmap=ColorMapsMatplotlib[0], vmin=np.min(Z)-0.2, vmax=np.max(Z)+0.2)
# Still plotting the contours of the LS
- Fluids = set(map(int, Solution.FluidFlag))
+ 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,:])
@@ -831,7 +831,7 @@ def PlotVariableOnEachSubDomain2D(Problem, Mesh, Solution, VariableId, Parameter
# 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.FluidFlag))
+ 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,:])
@@ -969,7 +969,7 @@ def PlotVariableOneEachSubDomain(Problem, Mesh, Solution, VariableId, Parameters
# 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.FluidFlag))
+ 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,:])
diff --git a/SourceCode/PostProcessing/Plots_Plotly.py b/SourceCode/PostProcessing/Plots_Plotly.py
index d3d3a7d..3548cb0 100644
--- a/SourceCode/PostProcessing/Plots_Plotly.py
+++ b/SourceCode/PostProcessing/Plots_Plotly.py
@@ -302,7 +302,7 @@ def PlotAnimation3D(TestName, Resolution, FrameFunc, cmin, cmax, args, Title = "
# ----------- Collecting the solution figures from the given data -------------------
# Selecting all the available files w.r.t. to the exported time points
lst = glob.glob(os.path.join(Pathes.SolutionPath, TestName, "Res"+str(Resolution).replace(".","_"), "RawResults","*.sol"))
- Times = np.sort([float(FileName[-9:-4].replace("_",".")) for FileName in lst])
+ Times = np.sort([float(FileName[-12:-4].replace("_",".")) for FileName in lst])
# Creating one frame per time-point
Frames = []
@@ -389,7 +389,7 @@ def PlotAnimation2D(TestName, Resolution, FrameFunc, cmin, cmax, args, Title = "
# ----------- Collecting the solution figures from the given data -------------------
# Selecting all the available files w.r.t. to the exported time points
lst = glob.glob(os.path.join(Pathes.SolutionPath, TestName, "Res"+str(Resolution).replace(".","_"), "RawResults","*.sol"))
- Times = np.sort([float(FileName[-9:-4].replace("_",".")) for FileName in lst])
+ Times = np.sort([float(FileName[-12:-4].replace("_",".")) for FileName in lst])
# Creating one frame per time-point
Frames = []
@@ -517,7 +517,7 @@ def PlotAllLevelSets(Mesh, Solution, ParametersId):
# Creating the file name and exporting it in the right folder
try:
- FileName = ("Plty_LSValues3D_All_t{0:5.3f}".format(Solution.t)).replace(".","_")
+ FileName = ("Plty_LSValues3D_All_t{0:8.6f}".format(Solution.t)).replace(".","_")
fig.write_html(os.path.join(FolderName, FileName)+".html")
except:
print("WARNING: Error in the plot generation. Impossible export.")
@@ -621,7 +621,7 @@ def PlotEachLevelSet(Mesh, Solution, ParametersId):
# --------------- Export the figure at the right location -------------------------
# Creating the file name and exporting it in the right folder
try:
- FileName = ("Plty_LSValues3D_Fluid{0:02d}_t{1:5.3f}".format(Id, Solution.t)).replace(".","_")
+ FileName = ("Plty_LSValues3D_Fluid{0:02d}_t{1:8.6f}".format(Id, Solution.t)).replace(".","_")
fig.write_html(os.path.join(FolderName, FileName)+".html")
except:
print("WARNING: Error in the plot generation. Impossible export.")
@@ -680,7 +680,7 @@ def PlotSolution(Problem, Mesh, Solution, VariableId, ParametersId):
# --------------- Export the figure at the right location -------------------------
# Creating the file name and exporting it in the right folder
try:
- FileName = ("Plty_3DSol_OnWholeDomain_{0!s}_t{1:5.3f}".format(variableName, Solution.t)).replace(".","_")
+ FileName = ("Plty_3DSol_OnWholeDomain_{0!s}_t{1:8.6f}".format(variableName, Solution.t)).replace(".","_")
fig.write_html(os.path.join(FolderName, FileName)+".html")
except:
print("WARNING: Error in the plot generation. Impossible export.")
@@ -785,7 +785,7 @@ def PlotSolutionOnEachSubdomain(Problem, Mesh, Solution, VariableId, ParametersI
# Creating the file name and exporting it in the right folder
try:
- FileName = ("Plty_3DSol_OnSubDomains_{0!s}_t{1:5.3f}".format(variableName, Solution.t)).replace(".","_")
+ FileName = ("Plty_3DSol_OnSubDomains_{0!s}_t{1:8.6f}".format(variableName, Solution.t)).replace(".","_")
fig.write_html(os.path.join(FolderName, FileName)+".html")
except:
print("WARNING: Error in the plot generation. Impossible export.")
@@ -869,7 +869,7 @@ def PlotLevelSetOnEachSubDomain(Problem, Mesh, Solution, ParametersId):
# --------------- Export the figure at the right location -------------------------
# Creating the file name and exporting it in the right folder
try:
- FileName = ("Plty_LSValues2D_SubDomains_t{0:5.3f}".format(Solution.t)).replace(".","_")
+ FileName = ("Plty_LSValues2D_SubDomains_t{0:8.6f}".format(Solution.t)).replace(".","_")
fig.write_html(os.path.join(FolderName, FileName)+".html")
except:
print("WARNING: Error in the plot generation. Impossible export.")
@@ -1023,7 +1023,7 @@ def PlotSolution2D(Problem, Mesh, Solution, VariableId, ParametersId):
# --------------- Export the figure at the right location -------------------------
# Creating the file name and exporting it in the right folder
try:
- FileName = ("Plty_2DSol_OnWholeDomain_{0!s}_t{1:5.3f}".format(variableName, Solution.t)).replace(".","_")
+ FileName = ("Plty_2DSol_OnWholeDomain_{0!s}_t{1:8.6f}".format(variableName, Solution.t)).replace(".","_")
fig.write_html(os.path.join(FolderName, FileName)+".html")
except:
print("WARNING: Error in the plot generation. Impossible export.")
@@ -1111,7 +1111,7 @@ def PlotSolution2DOnEachSubDomain(Problem, Mesh, Solution, VariableId, Parameter
# --------------- Export the figure at the right location -------------------------
# Creating the file name and exporting it in the right folder
try:
- FileName = ("Plty_2DSol_OnSubDomains_{0!s}_t{1:5.3f}".format(variableName, Solution.t)).replace(".","_")
+ FileName = ("Plty_2DSol_OnSubDomains_{0!s}_t{1:8.6f}".format(variableName, Solution.t)).replace(".","_")
fig.write_html(os.path.join(FolderName, FileName)+".html")
except:
print("WARNING: Error in the plot generation. Impossible export.")
@@ -1137,7 +1137,7 @@ def Animate3DLevelSet(TestName, Resolution):
# ---------------- Creating the color bounds for the global plot -------------------
# Listing all the available files
lst = glob.glob(os.path.join(Pathes.SolutionPath, TestName, "Res"+str(Resolution).replace(".","_"), "RawResults","*.sol"))
- Times = np.sort([float(FileName[-9:-4].replace("_",".")) for FileName in lst])
+ Times = np.sort([float(FileName[-12:-4].replace("_",".")) for FileName in lst])
# Safety check
if Times.size==0: raise(ValueError("Error: No previous results exported. Abortion."))
@@ -1193,7 +1193,7 @@ def Animate2DLevelSet(TestName, Resolution):
# ---------------- Creating the color bounds for the global plot -------------------
# Listing all the available files
lst = glob.glob(os.path.join(Pathes.SolutionPath, TestName, "Res"+str(Resolution).replace(".","_"), "RawResults","*.sol"))
- Times = np.sort([float(FileName[-9:-4].replace("_",".")) for FileName in lst])
+ Times = np.sort([float(FileName[-12:-4].replace("_",".")) for FileName in lst])
# Safety check
if Times.size==0: raise(ValueError("Error: No previous results exported. Abortion."))
@@ -1251,7 +1251,7 @@ def Animate3DSolution(TestName, Resolution):
# ---------------- Creating the color bounds for the global plot -------------------
# Listing all the available files
lst = glob.glob(os.path.join(Pathes.SolutionPath, TestName, "Res"+str(Resolution).replace(".","_"), "RawResults","*.sol"))
- Times = np.sort([float(FileName[-9:-4].replace("_",".")) for FileName in lst])
+ Times = np.sort([float(FileName[-12:-4].replace("_",".")) for FileName in lst])
# Safety check
if Times.size==0: raise(ValueError("Error: No previous results exported. Abortion."))
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py b/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py
index a43c0fd..3efb3f2 100644
--- a/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py
+++ b/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py
@@ -77,9 +77,9 @@ class Equations():
Input: Var, float numpy array, the value of the variables (NbVariables x NbGivenPoints)
Optional arguments (in given order)
- FluidIndex, integer numpy array, the fluid present at each given point (NbGivenPoints)
x, float numpy array, (optional, not impacting here) x-y coordinates of the given points (NbGivenPoints x 2)
-
+ FluidIndex, integer numpy array, the fluid present at each given point (NbGivenPoints)
+
Output: None, fills direclty the Parameters data structure
Note: -> This function is Vectorised
@@ -89,7 +89,7 @@ class Equations():
# --------- Locating and extracting the right fluid's properties for each location ---------------
# Gettings the arguments
- FluidIndex = args[0]
+ FluidIndex = args[1]
# Initialising the pressure variable
p = np.zeros(len(Var[0,:]))
@@ -211,7 +211,7 @@ class Equations():
dfin/dx1, ....., dfin/dxn]
---------------------------------------------------------------------"""
# Retrieving the Jacobian
- J = self.Jacobian(Var, FluidIndex, x)
+ 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)
@@ -507,7 +507,7 @@ class InitialConditions():
Init[3,Ind1] = Init[0,Ind1]*e
# Computing the values that will be taken outside the SOD bump
- Ind2 = np.setdiff1d(PointsID, Ind1)
+ Ind2 = np.setdiff1d(range(len(PointsID)), Ind1)
Init[:,Ind2] = np.array([[0.125, 0., 0., 0.1]]*len(Ind2)).T
e = EOS(Init[:,Ind2], [FluidProp]*len(Ind2), "e")
Init[3,Ind2] = Init[0,Ind2]*e
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvection.py b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py
similarity index 99%
rename from SourceCode/ProblemDefinition/GoverningEquations/LinearAdvection.py
rename to SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py
index 98b0cca..c2d7151 100644
--- a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvection.py
+++ b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py
@@ -121,7 +121,7 @@ class Equations():
# Computing the velocity values from the conservative Euleur equations
UV = np.squeeze(self.Flux(Var, x))
-
+
# Returning the value
return(UV)
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py
new file mode 100644
index 0000000..27b4c89
--- /dev/null
+++ b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py
@@ -0,0 +1,327 @@
+"""
+################################################################################
+# Module that defines all the necessary routines that implements the #
+# 2D rotation, linear advection #
+################################################################################
+
+Contains (in scrolling order):
+ -> Directly accessible parameters:
+ FluidModel, string, Name of the fluidModel to use
+
+ -> Equations class
+
+ -> EOS class
+
+ -> Initial conditions class
+
+Usage:
+
+Note:
+-----------------------------------------------------------------------------"""
+
+# ==============================================================================
+# Preliminary imports
+# ==============================================================================
+import numpy as np
+
+
+# ==============================================================================
+# Parameter that should be accessible from outside the class
+# ==============================================================================
+FluidModel = "Dummy"
+IntegratedLS = False
+
+# ==============================================================================
+# Equation class giving the flux, jacobian and all tools required by schemes
+# ==============================================================================
+class Equations():
+ """ Class furnishing all the methods that are required to define numerical
+ schemes and evolve the solution according to the Linear Advection (rotation)
+
+ Initialisation Parameters (later available as well):
+ ----------------------------------------------------
+ FluidProp, list of FluidModel, the list of fluid property instances (NbFluids)
+ EOS, list of callbacks, the list of Equation of state functions (NbFluids)
+
+ Available Parameters:
+ ---------------------
+ NbVariables, integer, number of variables of the model
+ VariablesNames, list of strings, ordered name of the variables
+ VariablesUnits, list of strings, symbols of the units of the variables
+ VariablesLatexNames, list of strings, ordered name of the variables, latex encoding
+ VariablesLatexUnits, list of strings, symbols of the units of the variables, latex encoding
+ XYLatexUnits, list of strings, units of the x-y coordinates in latex encoding
+
+ Available Methods:
+ ------------------
+ Flux
+ Jacobian
+
+ --------------------------------------------------------------------------------------"""
+
+ #### Automatic initialisation routine ####
+ def __init__(self, FluidProp, EOs):
+ """ Initialisation routine upon the considered fluid's properties (to have a direct access).
+ Input: FluidProp, list of FluidModel, the list of fluid property instances (NbFluids)
+ EOS, list of callbacks, the list of Equation of state functions (NbFluids)
+ ------------------------------------------------------------------------------------------"""
+
+ # Defines the (conservative) variables number and names, in a latex typography
+ self.NbVariables = 1
+ self.VariablesNames = ["h"]
+ self.VariablesUnits = ["m"]
+
+ self.VariablesLatexNames = [r'$h$']
+ self.VariablesLatexUnits = [r'$m$']
+ self.XYLatexUnits = [r'$m$', r'$m$']
+
+ # Register the shortcut of the fluid's properties
+ self.FluidProp = FluidProp
+ self.EOs = EOs
+
+ #### Flux function of the equations set ####
+ def Flux(self, Var, x, *args):
+ """ Flux function of the LinearAdvection (rotation case)
+ Input: Var, float numpy array, the value of the variables (NbVariables x NbGivenPoints)
+
+ Optional arguments (in given order)
+ FluidIndex, integer numpy array, the fluid present at each given point (NbGivenPoints)
+ x, float numpy array, (generally optional, required for this problem)
+ x-y coordinates of the given points (NbGivenPoints x 2)
+
+ Output: 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)
+ -----------------------------------------------------------------------------------------------------"""
+
+ # --------- Computing the Euler flux directly in conservative variables ------------------------
+ # Giving the flux furnishing the roation
+ Fx = np.array([Var[0,:]]) # Displacement in x direction
+ Fy = 0*np.array([Var[0,:]]) # Displacement in y direction
+
+ # Returning the flux
+ return(np.array([Fx,Fy]))
+
+ #### Function that extracts the propagation speed from the variables to interface for the LevelSet ####
+ 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.
+ Input: Var, float numpy array, the value of the variables (NbVariables x NbGivenPoints)
+ x, float numpy array, (optional, not impacting here) x-y coordinates of the given points (NbGivenPoints x 2)
+
+ Optional arguments (in given order)
+ FluidIndex, integer numpy array, the fluid present at each given point (NbGivenPoints)
+
+ Output: UV, float numpy array, the velocities values at the points (2 x NbGivenPoints)
+
+ Note: -> This function is Vectorised
+ --------------------------------------------------------------------------------------------------------------------------------"""
+
+ # Computing the velocity values from the conservative Euleur equations
+ UV = np.squeeze(self.Flux(Var, x))
+
+ # Returning the value
+ return(UV)
+
+ #### Jacobian of the equations set ####
+ def Jacobian(self, Var, x, *argv):
+ """ Computes the Jacobian of the flux for the LinearAdvection (rotation case)
+ Input: Var, 2D numpy array, the variables of the problem (NbVars x NbGivenPoints)
+ x, 2D numpy array, (generally optional, required for this problem)
+ the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ Optional arguments (in given order)
+ FluidIndex, integer numpy array, the fluid present at each given point (NbGivenPoints)
+
+ Output: J, 3D numpy array, J[:,:,i] gives the jacobian of the flux taking care
+ of the dynamic of the ith spatial coordinate. For
+ each flux fi = (fi1,..., fin), the Jacobian reads:
+ J[:,:] = [dfi1/dx1, ....., dfi1/dxn
+ ....
+ dfin/dx1, ....., dfin/dxn]
+ ---------------------------------------------------------------------"""
+
+ # ------ Initialising the Jacobian structure ---------------------------
+ # Getting the number of furnished points
+ NbPoints = np.shape(Var)[1]
+ J = np.zeros((1, 1, 2, NbPoints))
+
+ # ------ Computing the actual Jacobian ---------------------------------
+ # Creating the Jacobian for the first coordinate
+ J[0, 0, 0, :] = 1+0*x[:,0]
+
+ # Creating the Jacobian for the second coordinate
+ J[0, 0, 1, :] = 0*x[:,0]
+
+ # Returning the Jacobian
+ return(J)
+
+ #### Spectral radius of the equations set ####
+ def SpectralRadius(self, Var, FluidIndex, n, x, *argv):
+ """ Computes the spectral radius associated to the flux.
+ Input: 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)
+
+ Output: J, 3D numpy array, J[:,:,i] gives the jacobian of the flux taking care
+ of the dynamic of the ith spatial coordinate. For
+ each flux fi = (fi1,..., fin), the Jacobian reads:
+ J[:,:] = [dfi1/dx1, ....., dfi1/dxn
+ ....
+ dfin/dx1, ....., dfin/dxn]
+ ---------------------------------------------------------------------"""
+ # Retrieving the Jacobian
+ J = self.Jacobian(Var, x)
+
+ # Computing the spectral radius at each given point, shortcuted since here we are scalar
+ Lmb = np.abs(np.sum(J[0,0,:,:]*n.T, axis = 0))
+
+ # Returning the spectral radius at each given point
+ return(Lmb)
+
+
+# ==============================================================================
+# Class containing all the methods for defining various Equation of states
+# ==============================================================================
+#### 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 SimpleFluid definition.
+
+ Note: After having initialised the instance, the EOS is accessible by the variable
+ EOS within this class.
+
+ WARNING: For advection, this class is void, just here for compatibility properties in the code
+ --------------------------------------------------------------------------------------"""
+
+ #### Automatic initialisation routine (mapper) ####
+ def __init__(self, Id):
+ """ Automatic initialisation routine."""
+
+ # Mapping towards the right function
+ if Id==1: self.EOS = self.Dummy
+
+ else: raise NotImplementedError("Error in the selection of the Equation of State\n\
+ Unknown index. Check your Problem definition. Abortion.")
+
+ #### Dummy EOS #####
+ def Dummy(self, Var, FluidProp, Type):
+ """ Dummy EOS returning the identity for the sake of the code compatibility.
+ Input: 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
+ Output: p, float numpy array, the pressure values at the given points (NbGivenPoints)
+ ----------------------------------------------------------------------------------------------------"""
+
+ # Void return
+ val = Var
+
+ # Returning the asked value
+ return(val)
+
+
+# ==============================================================================
+# Class containing all the predefined suitable initial conditions
+# ==============================================================================
+class InitialConditions():
+ """ Class furnishing the solution initialisation routines that are suitable
+ to study the the Euler Equations, defined for working on a subdomain
+
+ Note: After having initialised the instance, the Initialisation method is
+ accessible through the variable IC within this class
+ --------------------------------------------------------------------------------------"""
+
+ #### Automatic initialisation routine (mapper) ####
+ def __init__(self, Id):
+ """ Automatic initialisation routine."""
+
+ # Mapping the initialisation routines
+ if Id == 0: self.IC = self.ConstantState
+ elif Id == 1: self.IC = self.ConstantState2
+ elif Id == 2: self.IC = self.Bump
+
+ else: raise NotImplementedError("Error in the selection of the Initial conditions\
+ Unknown index. Check your Problem definition. Abortion.")
+
+ #### 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.
+ Input: 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 structure, (optional, not used here) the the properties of the fluid present where the given points are
+ Subdomain, shapely multipoly, (optional, not used here) the shapely polygon to which the given points belong
+ Output: 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.ones((1, len(PointsID)))
+
+ # Defining the raw initilisation variable
+ Init[0,:] = 0.5
+
+ # 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.
+ Input: 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 structure, (optional, not used here) the the properties of the fluid present where the given points are
+ Subdomain, shapely multipoly, (optional, not used here) the shapely polygon to which the given points belong
+ Output: 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.ones((1, len(PointsID)))
+
+ # Defining the raw initilisation variable
+ Init[0,:] = 1.5
+
+ # Returning the computed values
+ return(Init)
+
+ #### Initialise the solution to a StationaryVortex according to the given subdomain
+ def Bump(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).
+ Input: 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 structure, (optional, not used here) the the properties of the fluid present where the given points are
+ Subdomain, shapely multipoly, (optional, not used here) the shapely polygon to which the given points belong
+ Output: Init, float numpy array, The initialised values at the considered points (NbVariables x NbGivenPoints)
+
+ Notes: -> 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)))
+
+ # Retrieving the xy values of the points and the geometrical informations from the subdomain
+ xy = Mesh.DofsXY[PointsID,:]
+ xc = subdomain.centroid.x-3.5
+ 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
+ return(Init)
diff --git a/SourceCode/Solve.py b/SourceCode/Solve.py
index db39f01..f0d258a 100644
--- a/SourceCode/Solve.py
+++ b/SourceCode/Solve.py
@@ -139,7 +139,7 @@ def Solve(*argv):
# 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
print("Computing for the time {0!s}".format(t))
print(np.max(Solution.Sol))
@@ -147,18 +147,21 @@ def Solve(*argv):
# Main time loop
k=0
- while t<4:#Solver.Tmax:
+ while t<Solver.Tmax:
k=k+1
- # Get dt
+ # ----- Perform a time step
+ # Get the time step
dt = _Solver.Solver_Temporal.GetTimeStep(Solver.CFL, Problem, Mesh, Solution)
# One time iteration taking as parameter the spatial schemes
- #Solution.Sol = Solver.TimeStep(Solver.SpatialScheme, dt, Mesh, Solution, "Sol")
+ Solution.Sol = Solver.TimeStep(Solver.SpatialScheme, dt, Mesh, Solution, "Sol")
# If treated separately, update the level set upon the old solution
- Solution.LSValues = Solver.TimeStep(Solver.FluidSelector.Update, dt, Mesh, Solution, "LSValues")
+ Solution.LSValues = Solver.TimeStep(Solver.FluidSelector.UpdateFSValues, dt, Mesh, Solution, "LSValues")
+ Solver.FluidSelector.UpdateFSVAttributes(Solution)
+ # ----- Perform the subsidiary updates
# Increase the current time
t = t+dt
Solution.t = t
@@ -166,13 +169,13 @@ def Solve(*argv):
# Reconstruct the solution at the physical vertices of the mesh (just impacting the visualisation)
Solver.Reconstructor(Solution)
- # Inform the user
- if np.mod(k,10)==0:
- print("Computing for the time {0!s}".format(t))
+ # ----- Inform the user and export the solution when relevant
+ print("Computing for the time {0!s}".format(t))
+ if np.mod(k,1)==0:
print(np.max(Solution.Sol))
print(np.min(Solution.Sol))
EX.ExportSolution(TestName, Mesh, ProblemData, Problem, Solution)
-
+
# ------------------ Closing tasks -----------------------------------------
print("\n -------------------------------------\
diff --git a/SourceCode/_Solver/FluidSelectors/NaiveLS.py b/SourceCode/_Solver/FluidSelectors/NaiveLS_CellCentredLFX.py
similarity index 51%
rename from SourceCode/_Solver/FluidSelectors/NaiveLS.py
rename to SourceCode/_Solver/FluidSelectors/NaiveLS_CellCentredLFX.py
index d412eab..d69ce10 100644
--- a/SourceCode/_Solver/FluidSelectors/NaiveLS.py
+++ b/SourceCode/_Solver/FluidSelectors/NaiveLS_CellCentredLFX.py
@@ -13,13 +13,15 @@ Note:
# ==============================================================================
# Importing python libraries
# ==============================================================================
-from shapely.geometry import Polygon as ShapelyPolygon
-from shapely.geometry import Point as ShapelyPoint
+from shapely.geometry import MultiPolygon as ShapelyMultiPolygon
+from shapely.geometry import Polygon as ShapelyPolygon
+from shapely.geometry import Point as ShapelyPoint
+from shapely.ops import unary_union
import matplotlib.pyplot as plt
import matplotlib.path as PH
import numpy as np
-
+import Plots_Matplotlib as PP
# ==============================================================================
# Defining the actuall class gathering the externally accessible functions
# ==============================================================================
@@ -80,19 +82,101 @@ class FS():
# Returning the fresh LS values
return(LSValues)
- def Redistancing(self, Solution):
- # Redistances the Level set values and reassigns the new values to the fluidflags
- # This is useless if the surface is infered from the LSValues by the nearest interpolation,
- # use at least the linear interpolation
+ #### Determining the fluid flags depending on the various values of the level set #####
+ def FlagDofs(self, Solution):
+ """ Flags the degrees of freedom according to the different values of the
+ associated level set.
+
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: None, updating directly the value of FluidFlag in the
+ solution instance
+ --------------------------------------------------------------------------------------- """
+
+ # ---------------- Updating the flags at Dofs ------------------------------
+ Solution.FluidFlag = np.argmax(Solution.LSValues, axis=0)
+
+ #### Recomputing the subdomains according to the new values of the level set #####^
+ def UpdateSubdomains(self, Solution):
+ """ Method modifying the subdomains associated to the Solution upon the values of
+ LSValues.
+
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: None, updating directly the subdomains in the
+ solution instance
+
+ WARNING: Not safe yet if the first path given out is internal
+ ---------------------------------------------------------------------------- """
+
+ # ------- Updating the subdomains from the countour values -----------------
+ # Buffering the LSvalues to modify them on the edges upon our wishes to infer closed contours
+ Buffer = Solution.LSValues[:,:]
+
+ # Matplotlib initialisation to get back the countours
+ fig = plt.figure()
+ ax = fig.add_subplot(2, 1, 1)
+
+ # Update each of the fluid's locations
+ for i in range(self.Problem.NbFluids):
+ # ----------- Editing the LS value with a pythonic hack -----------------
+ # Force the boundary of the domain to be part of the contour if relevant
+ ind = np.where(Solution.LSValues[i,self.Mesh.BoundaryDofs]>=0)[0]
+ Buffer[i,np.array(self.Mesh.BoundaryDofs)[ind]]=0
+
+ cs = ax.tricontour(self.Mesh.DofsXY[:,0], self.Mesh.DofsXY[:,1],\
+ self.Mesh.ElementsBoundaryDofs, Buffer[i,:], [0])
+
+ # ----------- Extracting the contour of the level set -------------------
+ # And transforming into a shapely mutlipolygon
+
+ # Extracting the 0 level of the function and creating a polygon out of it,
+ # making sure that it is not self-intersecting
+ LV0 = cs.allsegs[0]
+ polygon = ShapelyPolygon([tuple(a) for a in np.array(LV0[0])]).buffer(0)
+
+ # Spanning over the possibly disconnected area encompassing the target fluid
+ for contour in LV0[1:]:
+ # Determining whether the patch is out or in the previously registered polygon
+ # and adding/removing it upon
+ patch = ShapelyPolygon([tuple(a) for a in np.array(contour)]).buffer(0)
+ if patch.difference(polygon).is_empty == False: polygon = polygon.union(patch)
+ else: polygon = polygon.difference(patch)
+
+ # Converting it to multipolygon for compatibility with other parts of the code (esp. plots)
+ if not (polygon.geom_type == 'MultiPolygon'): polygon = ShapelyMultiPolygon([polygon])
- # Finding the 0 level set
+ # Registering the new shape
+ Solution.Subdomains[i] = polygon
+
+ plt.close(fig)
- # Redistancing the points with the nearest
+ #### Redistancing the LSvalues according to the freshly computed level set values ####
+ def Redistancing(self, Solution):
+ """ Redistances the Level set values upon the computed subdomains.
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: None, updating directly the LSValues in the
+ solution instance
+ ----------------------------------------------------------------------------------"""
- pass
+ # Redistancing the points with the nearest
+ # Looping on the intial number of fluids present in the computational domain
+ for i in range(self.Problem.NbFluids):
+ # Computing the distance and location w.r.t. the subdomain through the shapely information
+ FluidArea = Solution.Subdomains[i]
+ Value = np.array([[FluidArea.boundary.distance(ShapelyPoint(tuple(pt))), \
+ int(2*FluidArea.intersects(ShapelyPoint(tuple(pt)).buffer(1e-14))-1)] \
+ for pt in self.Mesh.DofsXY]).T
+
+ # Computing the signed distance
+ Solution.LSValues[i,:] = np.prod(Value, axis=0)
- def Update(self, Solution):
- # Performs the level set update from stupid LFX centred finite volume scheme
+ #### (Crucial routine) Spatial scheme to update the level set ####
+ def UpdateFSValues(self, Solution):
+ """ Spatial scheme to update the level set:
+ Performs the level set update from stupid LFX centred finite volume scheme
+
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: Res, the residuals after one iteration
+ --------------------------------------------------------------------------------- """
# ---------------- Updating the level set values ---------------------------
# Initialising the mollification vector (residuals)
@@ -106,7 +190,8 @@ class FS():
Ele = np.array(ele + [ele[0]], dtype=int)
# Get the average state of the considered element
- LS1 = np.mean(Solution.LSValues[:,dele], axis = 1)
+ LS1 = np.mean(Solution.LSValues[:,dele], axis = 1)
+ Val1 = np.mean(Solution.Sol[:,dele], axis = 1)
res = np.zeros((1,self.Problem.NbFluids))
# Looping on each edge
@@ -124,11 +209,11 @@ class FS():
# Get the average state of the neighboring element
Dele2 = np.array(self.Mesh.ElementsBoundaryDofs[Nei])
LS2 = np.mean(Solution.LSValues[:,Dele2], axis=1)
+ Val2 = np.mean(Solution.Sol[:,Dele2], axis=1)
# Compute the input exchange from the neighboring element at the edge midpoint
EdgeMidPoint = np.array(2*[0.5*(np.sum(self.Mesh.CartesianPoints[Edge, :], axis=0))])
- UV = self.Problem.GoverningEquations.GetUV(np.array([LS1,LS2]).T, EdgeMidPoint)
-
+ UV = self.Problem.GoverningEquations.GetUV(np.array([Val1,Val2]).T, EdgeMidPoint)
MidFluxes = np.array([UV[0,:]*np.array([LS1,LS2]).T,\
UV[1,:]*np.array([LS1,LS2]).T])
@@ -154,12 +239,17 @@ class FS():
# Updating each nodal value (here representing the average)
Res[:, :] = Res[:, :]
+ # Returning the residuals
return(Res)
-
- # ---------------- Updating the flags at Dofs ------------------------------
-
-
- # ------- Updating the subdomains from the countour values -----------------
+ #### Wrapper routine to update the quantitites that depend on the UpdateFSValues output ####
+ def UpdateFSVAttributes(self, Solution):
+ """ Wrapper routine to update the quantitites that depend on the UpdateFSValues fresh output
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: None, updating directly the subvalues in the
+ solution instance
+ --------------------------------------------------------------------------------- """
-
+ #self.UpdateSubdomains(Solution)
+ #self.Redistancing(Solution)
+ self.FlagDofs(Solution)
diff --git a/SourceCode/_Solver/FluidSelectors/NaiveLS_VertexCentredLFX.py b/SourceCode/_Solver/FluidSelectors/NaiveLS_VertexCentredLFX.py
new file mode 100644
index 0000000..a5b2b90
--- /dev/null
+++ b/SourceCode/_Solver/FluidSelectors/NaiveLS_VertexCentredLFX.py
@@ -0,0 +1,235 @@
+"""
+################################################################################
+# Module implementing a naive level set with redistancing #
+################################################################################
+
+Contains (in scrolling order):
+
+Usage:
+
+Note:
+-----------------------------------------------------------------------------"""
+
+# ==============================================================================
+# Importing python libraries
+# ==============================================================================
+from shapely.geometry import Polygon as ShapelyPolygon
+from shapely.geometry import Point as ShapelyPoint
+import matplotlib.pyplot as plt
+import matplotlib.path as PH
+import numpy as np
+
+
+# ==============================================================================
+# Defining the actuall class gathering the externally accessible functions
+# ==============================================================================
+class FS():
+ """ Classical level set approach, initialised according to a signed distance
+ and benefiting from a redistancing/initialisation (one LS per fluid registered
+ in the domain)
+ ----------------------------------------------------------------------------"""
+
+ #### Automatic initialisation routine ####
+ def __init__(self, Problem, Mesh, Params = []):
+ """Input:
+ Mesh, MeshStructure, the mesh on which the problem is solved
+ Params, list of strings, list of parameters passed to the LS"""
+
+ # Setting the name of the class to be accessible from outside
+ self.FluidSelectorName = "Naive Level Set"
+ # Setting the mesh and parameters for an internal access
+ self.Problem = Problem
+ self.Mesh = Mesh
+ self.Params = Params
+
+ #### External initialisation routine ####
+ def Initialiser(self, Solution):
+ """ External initialisation routine.
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: LSValues, float numpy array, the level set values at the MeshVertex and Dofs,
+ to be copied to the Solution structure (NbSubdomains x (NbMeshPoints+NbDofs))
+ ------------------------------------------------------------------------------------------------------------------"""
+
+ # --------------- Shortcutting the required values ----------------------------------------------
+ # Getting the subdomains list
+ Subdomains = Solution.Subdomains
+
+ # Recompute the number of dofs to drop the dependency on the mesh
+ NbFluids = len(Subdomains) # Valid because straight after initialisation
+ NbDofs = len(Solution.FluidFlag) # Only dofs are there
+
+ # Initialising the Level set values
+ LSValues = np.zeros((len(Subdomains), len(Solution.FluidFlag)))
+
+ # ------------- Get the distance from the interfaces --------------------------------------------
+ # Evaluate the distance from Dofs to each boundary, and mark them by +-
+ # depending on their location with respect to the fluid's boundary
+ # Return the LS value for each fluid and for each Dof # may lead to non-precise determination
+ # in case of multiple junctions (side effect of the classicall method.)
+
+ # Looping on the intial number of fluids present in the computational domain
+ for i in range(NbFluids):
+ # Computing the distance and location w.r.t. the subdomain through the shapely information
+ FluidArea = Subdomains[i]
+ Value = np.array([[FluidArea.boundary.distance(ShapelyPoint(tuple(pt))), \
+ int(2*FluidArea.intersects(ShapelyPoint(tuple(pt)).buffer(1e-14))-1)] \
+ for pt in self.Mesh.DofsXY]).T
+ # Computing the signed distance
+ LSValues[i,:] = np.prod(Value, axis=0)
+
+ # Returning the fresh LS values
+ return(LSValues)
+
+ #### Determining the fluid flags depending on the various values of the level set #####
+ def FlagDofs(self, Solution):
+ """ Flags the degrees of freedom according to the different values of the
+ associated level set.
+
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: None, updating directly the value of FluidFlag in the
+ solution instance
+ --------------------------------------------------------------------------------------- """
+
+ # ---------------- Updating the flags at Dofs ------------------------------
+ Solution.FluidFlag = np.argmax(Solution.LSValues, axis=0)
+
+ #### Recomputing the subdomains according to the new values of the level set #####^
+ def UpdateSubdomains(self, Solution):
+ """ Method modifying the subdomains associated to the Solution upon the values of
+ LSValues.
+
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: None, updating directly the subdomains in the
+ solution instance
+
+ WARNING: Not safe yet if the first path given out is internal
+ ---------------------------------------------------------------------------- """
+
+ # ------- Updating the subdomains from the countour values -----------------
+ # Buffering the LSvalues to modify them on the edges upon our wishes to infer closed contours
+ Buffer = Solution.LSValues[:,:]
+
+ # Matplotlib initialisation to get back the countours
+ fig = plt.figure()
+ ax = fig.add_subplot(2, 1, 1)
+
+ # Update each of the fluid's locations
+ for i in range(self.Problem.NbFluids):
+ # ----------- Editing the LS value with a pythonic hack -----------------
+ # Force the boundary of the domain to be part of the contour if relevant
+ ind = np.where(Solution.LSValues[i,self.Mesh.BoundaryDofs]>=0)[0]
+ Buffer[i,np.array(self.Mesh.BoundaryDofs)[ind]]=0
+
+ cs = ax.tricontour(self.Mesh.DofsXY[:,0], self.Mesh.DofsXY[:,1],\
+ self.Mesh.ElementsBoundaryDofs, Buffer[i,:], [0])
+
+ # ----------- Extracting the contour of the level set -------------------
+ # And transforming into a shapely mutlipolygon
+
+ # Extracting the 0 level of the function and creating a polygon out of it,
+ # making sure that it is not self-intersecting
+ LV0 = cs.allsegs[0]
+ polygon = ShapelyPolygon([tuple(a) for a in np.array(LV0[0])]).buffer(0)
+
+ # Spanning over the possibly disconnected area encompassing the target fluid
+ for contour in LV0[1:]:
+ # Determining whether the patch is out or in the previously registered polygon
+ # and adding/removing it upon
+ patch = ShapelyPolygon([tuple(a) for a in np.array(contour)]).buffer(0)
+ if patch.difference(polygon).is_empty == False: polygon = polygon.union(patch)
+ else: polygon = polygon.difference(patch)
+
+ # Converting it to multipolygon for compatibility with other parts of the code (esp. plots)
+ if not (polygon.geom_type == 'MultiPolygon'): polygon = ShapelyMultiPolygon([polygon])
+
+ # Registering the new shape
+ Solution.Subdomains[i] = polygon
+
+ plt.close(fig)
+
+ #### Redistancing the LSvalues according to the freshly computed level set values ####
+ def Redistancing(self, Solution):
+ """ Redistances the Level set values upon the computed subdomains.
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: None, updating directly the LSValues in the
+ solution instance
+ ----------------------------------------------------------------------------------"""
+
+ # Redistancing the points with the nearest
+ # Looping on the intial number of fluids present in the computational domain
+ for i in range(self.Problem.NbFluids):
+ # Computing the distance and location w.r.t. the subdomain through the shapely information
+ FluidArea = Solution.Subdomains[i]
+ Value = np.array([[FluidArea.boundary.distance(ShapelyPoint(tuple(pt))), \
+ int(2*FluidArea.intersects(ShapelyPoint(tuple(pt)).buffer(1e-14))-1)] \
+ for pt in self.Mesh.DofsXY]).T
+
+ # Computing the signed distance
+ Solution.LSValues[i,:] = np.prod(Value, axis=0)
+
+ #### (Crucial routine) Spatial scheme to update the level set ####
+ def UpdateFSValues(self, Solution):
+ # Performs the level set update from stupid LFX centred finite volume scheme
+
+ # ---------------- Updating the level set values ---------------------------
+ # Initialising the mollification vector (residuals)
+ Res = np.zeros(np.shape(Solution.LSValues))
+
+ # Looping on each element
+ for i in range(self.Mesh.NbInnerElements):
+ ala = self.Mesh.InnerElements[i]
+
+ # Looping on each edge
+ for f in range(len(ala)):
+ # Getting the geometrical properties
+ edge = np.array(self.Mesh.ElementsBoundaryWiseDofs[i][f])
+ L = self.Mesh.DualEdgesWidths[i][f]
+ n = self.Mesh.DualNormals[i][f]
+
+ # Getting the variational quantitites
+ Vars = Solution.LSValues[:, edge]
+
+ # Creating the compatibility quantities to call the spectrum and flux functions
+ Coords = np.array(len(Vars[0,:])*[0.5*(self.Mesh.DualCenters[i][f] + self.Mesh.DualCenters[i][-1])])
+ nn = np.array(len(Vars[0,:])*[n])
+
+ # Getting the variational properties for creating the numerical flux
+ UV = self.Problem.GoverningEquations.GetUV(Vars, Coords)
+ MidFluxes = np.array([UV[0,:]*Vars,\
+ UV[1,:]*Vars])
+
+ # Computing the Jacobian for the second coordinate
+ J = np.zeros((self.Problem.NbFluids, self.Problem.NbFluids, 2, 2))
+ for nf in range(self.Problem.NbFluids):
+ J[nf, nf, 0, :] = np.array([UV[0,:]])
+ J[nf, nf, 1, :] = np.array([UV[1,:]])
+ sprctrad = np.max(np.abs([np.sum(J[i,i,:,:]*n.T, axis = 0) for i in range(self.Problem.NbFluids)]), axis=0)
+
+ # Compute the flux contribution
+ MidFluxes[:,:,1] = -MidFluxes[:,:,1]
+ avrg = 0.5*np.sum(MidFluxes, axis=2)
+ Fluxn = np.dot(avrg.T,n.T)
+
+ # Distibuting the residuals to the associated Dof
+ Quantity = L*(Fluxn + np.max(sprctrad)*(Vars[:,0]-Vars[:,1]))
+ Res[:,edge[0]] += Quantity
+ Res[:,edge[-1]] -= Quantity
+
+ # Updating each nodal value (here representing the average) by the residual
+ # pondered by the control volume's area
+ Res[:, :] = Res[:, :]/self.Mesh.DualAreasTot[:].T
+
+ # Return the computed residuals
+ return(Res)
+
+ #### Wrapper routine to update the quantitites that depend on the UpdateFSValues output ####
+ def UpdateFSVAttributes(self, Solution):
+ """ Wrapper routine to update the quantitites that depend on the UpdateFSValues fresh output
+ Input: Solution, Solution structure, the solution instance the scheme is working on
+ Output: None, updating directly the subvalues in the
+ solution instance
+ --------------------------------------------------------------------------------- """
+
+ self.UpdateSubdomains(Solution)
+ self.Redistancing(Solution)
+ self.FlagDofs(Solution)
--
GitLab