From c5c161ee0e231de4d02237c96bfd71bc60a83de5 Mon Sep 17 00:00:00 2001
From: elemel <elise.lemeledo@math.uzh.ch>
Date: Wed, 12 Aug 2020 14:59:44 +0200
Subject: [PATCH] Added clean version with separated reconstructors
---
Pathes.py | 7 +-
...vection_Translation_2Fluids_FalseSystem.py | 41 --
.../PredefinedSettings/FirstOrderGalerkin.txt | 2 +-
.../SecondOrderGalerkin.txt | 6 +-
.../SecondOrderGalerkinStream.txt | 2 +-
.../Spatial1Time2_CGLFX.txt | 2 +-
RunPostProcess.py | 5 +-
RunProblem.py | 2 +-
SourceCode/Mappers.py | 146 ++---
SourceCode/PostProcess.py | 178 ++---
SourceCode/PostProcessing/Exports.py | 48 +-
SourceCode/PostProcessing/Plots_Matplotlib.py | 204 +++---
SourceCode/PostProcessing/Plots_Plotly.py | 231 ++++---
SourceCode/PostProcessing/Plots_Schlieren.py | 140 ++--
.../BoundaryConditions/Wall.py | 9 +-
SourceCode/ProblemDefinition/FluidModels.py | 37 +-
.../GoverningEquations/EulerEquations.py | 22 +-
.../LinearAdvectionRotation.py | 10 +-
.../LinearAdvectionTranslation.py | 12 +-
.../LinearAdvectionTranslationVectorTest.py | 492 --------------
SourceCode/Solve.py | 18 +-
.../_Solver/FluidSelectors/NaiveLS_CG.py | 383 ++++++-----
.../_Solver/FluidSelectors/NaiveLS_CG_LFX.py | 393 ++++++-----
.../FluidSelectors/NaiveLS_CG_LFX_Limited.py | 614 ++++++++++++++++++
SourceCode/_Solver/Solver_Initialisations.py | 144 ++--
SourceCode/_Solver/Solver_Spatial.py | 35 +-
SourceCode/_Solver/Solver_Temporal.py | 46 +-
.../SpatialModifiers/SpatialAmendements.py | 135 ++--
SourceCode/_Solver/SpatialSchemes/CG.py | 174 +++--
SourceCode/_Solver/SpatialSchemes/CG_LFX.py | 22 +-
.../_Solver/SpatialSchemes/CG_Primary.py | 26 +-
Templates/Template_Equations.py | 272 +++++---
Templates/Template_FluidSelector.py | 212 +++---
Templates/Template_Problem.py | 22 +-
Templates/Template_Settings.txt | 2 +-
Templates/Template_SpatialScheme.py | 84 +--
Templates/Template_TemporalScheme.py | 36 +-
_doc/source/ChildPages/DesignSettings.rst | 100 +--
_doc/source/ChildPages/SpatialModifiers.rst | 52 +-
_doc/source/index.rst | 58 +-
40 files changed, 2394 insertions(+), 2030 deletions(-)
delete mode 100644 PredefinedTests/PredefinedProblems/Test_Advection_Translation_2Fluids_FalseSystem.py
delete mode 100644 SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslationVectorTest.py
create mode 100644 SourceCode/_Solver/FluidSelectors/NaiveLS_CG_LFX_Limited.py
diff --git a/Pathes.py b/Pathes.py
index 7be8a76..272c4bd 100644
--- a/Pathes.py
+++ b/Pathes.py
@@ -1,7 +1,7 @@
"""
-This file can be located either in :code:`_SourceCode` or in LSRD folder,
-and specifies where to find the libraries containing the user-defined problems.
+This file can be located either in :code:`_SourceCode` or in LSRD folder,
+and specifies where to find the libraries containing the user-defined problems.
Edit it directly upon your local configuration. It has vocation to allow flexibility in where
your meshfiles, settings and problems are stored, and where you want to export the results to.
@@ -13,13 +13,14 @@ Just modify the following lines at your ease.
SettingPath = os.path.join(Dir, "Your wished location")
ProblemPath = os.path.join(Dir, "Your wished location")
SolutionPath = os.path.join(Dir, "Your wished location")
+
"""
# System import
import os
Dir = os.path.dirname(__file__)
-# Shortuct for the Path prefix that maps to the user's storage location, relatively to
+# Shortuct for the Path prefix that maps to the user's storage location, relatively to
# this file's storage location
MeshPath = os.path.join(Dir, "PredefinedTests", "PredefinedMeshes")
SettingPath = os.path.join(Dir, "PredefinedTests", "PredefinedSettings")
diff --git a/PredefinedTests/PredefinedProblems/Test_Advection_Translation_2Fluids_FalseSystem.py b/PredefinedTests/PredefinedProblems/Test_Advection_Translation_2Fluids_FalseSystem.py
deleted file mode 100644
index 291aee9..0000000
--- a/PredefinedTests/PredefinedProblems/Test_Advection_Translation_2Fluids_FalseSystem.py
+++ /dev/null
@@ -1,41 +0,0 @@
-################################################################################
-##
-## 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_2Fluids_FalseSystem" # Should match the module name
- self.MeshName = "BigSimpleSquare"
-
- # ******************************************************************************
- # Equation
- #*******************************************************************************
- self.GoverningEquationsIndex = 4
- 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 = [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,), (4,)]
-
- # 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/PredefinedSettings/FirstOrderGalerkin.txt b/PredefinedTests/PredefinedSettings/FirstOrderGalerkin.txt
index 5d997d6..69aeeef 100644
--- a/PredefinedTests/PredefinedSettings/FirstOrderGalerkin.txt
+++ b/PredefinedTests/PredefinedSettings/FirstOrderGalerkin.txt
@@ -17,4 +17,4 @@
# Export properties
m # Verbose mode (n=none, m=minimal, d=debug)
-1 # Time steps export intervals
+0.05 # Time steps export intervals
diff --git a/PredefinedTests/PredefinedSettings/SecondOrderGalerkin.txt b/PredefinedTests/PredefinedSettings/SecondOrderGalerkin.txt
index 134db27..89c8a6d 100644
--- a/PredefinedTests/PredefinedSettings/SecondOrderGalerkin.txt
+++ b/PredefinedTests/PredefinedSettings/SecondOrderGalerkin.txt
@@ -2,7 +2,7 @@
1 # Index of the scheme
2 B # Scheme properties (see doc for more details, if any, the order should be first)
Limiter_1 # Additives list Index of limiter (0=None)
-2 # Type of fluid spotter
+3 # Type of fluid spotter
2 B # Fluid spotter properties
# Quadrature properties
@@ -13,8 +13,8 @@ Limiter_1 # Additives list Index of limiter (0=None)
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
+3 # Tmax, in seconds
# Export properties
m # Verbose mode (n=none, m=minimal, d=debug)
-1 # Time steps export intervals
+0.05 # Time export intervals
diff --git a/PredefinedTests/PredefinedSettings/SecondOrderGalerkinStream.txt b/PredefinedTests/PredefinedSettings/SecondOrderGalerkinStream.txt
index 36ee022..4ab1ffd 100644
--- a/PredefinedTests/PredefinedSettings/SecondOrderGalerkinStream.txt
+++ b/PredefinedTests/PredefinedSettings/SecondOrderGalerkinStream.txt
@@ -17,4 +17,4 @@ Original+Filtering_1 # Additives list Index of limiter (0=None)
# Export properties
m # Verbose mode (n=none, m=minimal, d=debug)
-1 # Time steps export intervals
+0.05 # Time export intervals
diff --git a/PredefinedTests/PredefinedSettings/Spatial1Time2_CGLFX.txt b/PredefinedTests/PredefinedSettings/Spatial1Time2_CGLFX.txt
index 397f74e..b08ad4d 100644
--- a/PredefinedTests/PredefinedSettings/Spatial1Time2_CGLFX.txt
+++ b/PredefinedTests/PredefinedSettings/Spatial1Time2_CGLFX.txt
@@ -17,4 +17,4 @@ Limiter_1 # Additives list Index of limiter (0=None, see the docs for the strin
# Export properties
m # Verbose mode (n=none, m=minimal, d=debug) (yet inactive)
-1 # Time steps export intervals (yet inactive)
+0.05 # Time steps export intervals
diff --git a/RunPostProcess.py b/RunPostProcess.py
index 3942eaa..db3850b 100644
--- a/RunPostProcess.py
+++ b/RunPostProcess.py
@@ -73,6 +73,5 @@ if __name__ == '__main__':
#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)
+ #PostProcess.Plots_StaticSchlieren(ProblemDefinition, SettingsChoice, "all", MeshResolution)
+ #PostProcess.Plots_StaticPlotly(ProblemDefinition, SettingsChoice, [0.0], MeshResolution)
diff --git a/RunProblem.py b/RunProblem.py
index f2f8de6..3439f73 100644
--- a/RunProblem.py
+++ b/RunProblem.py
@@ -61,7 +61,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"
+ ProblemDefinition = "Euler_SOD_1Fluid"
SettingsChoice = "SecondOrderGalerkin"
MeshResolution = 20
diff --git a/SourceCode/Mappers.py b/SourceCode/Mappers.py
index 3288743..f176e37 100644
--- a/SourceCode/Mappers.py
+++ b/SourceCode/Mappers.py
@@ -10,7 +10,7 @@ Please then use this page as a reference when designing a new problem or setting
.. Note::
The mappers for the problem dependent equations of state, boundary and
- initial conditions should be provided in the ProblemFile module in the initialisation of their
+ initial conditions should be provided in the ProblemFile module in the initialisation of their
respective classes (see the problem definition's documentation).
|
@@ -51,32 +51,31 @@ import _Solver.FluidSelectors as FS
Mapper function for the governing equations
-----------------------------------------------------------------------------"""
def Governing_Equations(id):
- """
+ """
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
Mapper to the modules contained in the :code:`ProblemDefinition.GoverningEquations` library,
defining the governing equations.
-
+
.. rubric:: **Parameters**
-
- |bs| |bs| **id** *(integer)* -- the index corresponding to the investigated equation
+
+ |bs| |bs| **id** *(integer)* -- the index corresponding to the investigated equation
*Implemented equations given the id*
-
+
1. | Euler equations
2. | Linear advection (rotation)
3. | Linear advection (translation)
-
+
|
Returns:
- GoverningEquations.Module (module): the module itself,
+ GoverningEquations.Module (module): the module itself,
containing the classes associated to the given equation
(see the documentation of :code:`ProblemDefinition`)
----------------------------------------------------------------------------------------"""
- if id==0: return(GE.LinearAdvectionTranslationVectorTest)
- elif id==1: return(GE.EulerEquations)
+ if id==1: return(GE.EulerEquations)
elif id==2: return(GE.LinearAdvectionRotation)
elif id==3: return(GE.LinearAdvectionTranslation)
@@ -90,31 +89,31 @@ def Governing_Equations(id):
Mapper function for the pre-registered fluids
-----------------------------------------------------------------------------"""
def Fluid_Properties(id, Model):
- """
+ """
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
Mapper that defines the considered fluid properties to input in the model required by the
- GoverningEquations, according to the fluids templates (classes) defined in
+ GoverningEquations, according to the fluids templates (classes) defined in
:code:`ProblemDefinition.FluidModels` module.
-
+
.. rubric:: **Parameters**
-
+
|bs| |bs| **id** *(integer)* -- the index corresponding to the wished Fluid properties
-
+
|bs| |bs| **Model** *(string)* -- name of the fluid's model that one wants to use
-
+
*Implemented fluid types given the id*
-
+
0. | Air with properties :math:`\gamma=1.5` and :math:`p_\infty=0` and :math:`R=0.287`, with a blue representation color
1. | R22 fluid with properties :math:`\gamma=1.249` and :math:`p_\infty=0` and :math:`R=0.091`, with a red representation color
2. | Helium with properties :math:`\gamma=1.67` and :math:`p_\infty=0` and :math:`R=2.08`, with a green representation color
-
+
|
Returns:
Prop (Model instance): Returns an instance of the fluid model with filled properties
-
+
-------------------------------------------------------------------------------------------"""
if id==0: Properties = ['b', 1.4, 0, 0.287] # Air, ideal
@@ -138,23 +137,23 @@ def Fluid_Properties(id, Model):
Mapper function for the temporal schemes
-----------------------------------------------------------------------------"""
def Temporal_Scheme(id, params):
- """
+ """
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
Mapper that defines the temporal scheme to use, according to the predefined temporal schemes
given in :code:`SourceCode._Solver.TemporalSchemes` library.
-
+
.. rubric:: **Parameters**
-
+
|bs| |bs| **id** *(integer)* -- the index corresponding to the wished Fluid properties
-
+
|bs| |bs| **params** *(integers list)* -- list of the parameters to give to the time stepping functions as given by the user
-
+
*Implemented temporal schemes given the id*
-
+
1. | Deffered correction (DeC)
-
+
|
Returns:
@@ -162,7 +161,6 @@ def Temporal_Scheme(id, params):
---------------------------------------------------------------------------------------------"""
if id==1: return(fcts.partial(TI.DeC.DeC, *params))
- #elif id==2: return(fcts.partial(TI.RK.Iteration, *params))
else: raise ValueError("Error in the problem definition:\n\
Wrong index for the wished temporal scheme.\n\
@@ -174,34 +172,34 @@ def Temporal_Scheme(id, params):
Mapper function for the spatial schemes
-----------------------------------------------------------------------------"""
def Spatial_Scheme(id):
- """
-
+ """
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
Mapper that defines the temporal scheme to use, according to the predefined spatial schemes
given in :code:`SourceCode._Solver.SpatialSchemes` library.
-
+
.. rubric:: **Parameters**
-
+
|bs| |bs| **id** *(integer)* -- the index corresponding to the wished Fluid properties
-
+
*Implemented spatial schemes given the id*
-
+
0. | Continuous Galerkin
1. | Continuous Galerkin + Lax-Friedrichs flux
2. | Continuous Galerkin on Primary variables
-
+
|
Returns:
SpatialSchemes.Module (module): Returns the module associated to the wished spatial scheme
-
+
"""
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\
Check the Mapper and ProblemSettings files.\n\
@@ -215,18 +213,18 @@ def Limiter(id):
""" .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
Mapper that defines the limiting amendement to apply to the scheme, according to the given index.
- The user input for the admendements sequence is interpreted in the :code:`SourceCode._Solver.SpatialModifiers`
+ The user input for the admendements sequence is interpreted in the :code:`SourceCode._Solver.SpatialModifiers`
library. All the limiters should have a filename begining by "Limiter_".
-
+
.. rubric:: **Parameters**
-
+
|bs| |bs| **id** *(integer)* -- the index corresponding to the wished amendement to the scheme
-
+
*Implemented amendements given the id* and their type
-
- - |bs| Limiters
+
+ - |bs| Limiters
1. | Psi
-
+
|
Returns:
@@ -234,7 +232,7 @@ def Limiter(id):
----------------------------------------------------------------------------------"""
if id==1: return(SM.Limiter_Psi)
-
+
else: raise ValueError("Error in the problem definition:\n\
Wrong index for the wished limiter.\n\
Check the Mapper and ProblemSettings files.\n\
@@ -244,28 +242,28 @@ def Filtering(id):
""" .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
Mapper that defines the filtering amendement to apply to the scheme, according to the given index.
- The user input for the admendements sequence is interpreted in the :code:`SourceCode._Solver.SpatialModifiers`
+ The user input for the admendements sequence is interpreted in the :code:`SourceCode._Solver.SpatialModifiers`
library. All the filters should have a filename begining by "Filtering_".
-
+
.. rubric:: **Parameters**
-
+
|bs| |bs| **id** *(integer)* -- the index corresponding to the wished amendement to the scheme
-
+
*Implemented amendements given the id* and their type
- |bs| Filtering
- 1. | Streamline
-
+ 1. | Streamline (failing for systems)
+
|
Returns:
SpatialModifiers.Module (module): Returns the module associated to the wished limiter
----------------------------------------------------------------------------------"""
- if id==1:
+ 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\
Check the Mapper and ProblemSettings files.\n\
@@ -276,29 +274,31 @@ def Filtering(id):
Mapper function for the fluid selectors
-----------------------------------------------------------------------------"""
def Fluid_Selector(id):
- """
- Mapper that defines the the fluid selector to use, according to the ones defined in
+ """
+ Mapper that defines the the fluid selector to use, according to the ones defined in
:code:`SourceCode._Solver.FluidSelectors` module, containing its type, evolution equation,
fluid determination and possibly redistancing routines.
-
+
.. rubric:: **Parameters**
-
+
|bs| |bs| **id** *(integer)* -- the index corresponding to the wished Fluid selector
-
+
*Implemented fluid's selectors*
-
+
1. | NaiveLS_CG
2. | NaiveLS_CG_LFX
-
+ 3. | NaiveLS_CG_LFX_Limited
+
|
Returns:
FluidSelectors.Module (module): class containing all the fluid's selector tools
-
+
------------------------------------------------------------------------------------- """
if id == 1: return(FS.NaiveLS_CG)
elif id == 2: return(FS.NaiveLS_CG_LFX)
+ elif id == 3: return(FS.NaiveLS_CG_LFX_Limited)
else: raise ValueError("Error in the problem definition:\n\
Wrong index for the wished fluid selector.\n\
@@ -310,28 +310,28 @@ def Fluid_Selector(id):
Mapper function for the fluid selectors
-----------------------------------------------------------------------------"""
def Boundary_Conditions(id, *params):
- """
- Mapper that defines the the fluid selector to use, according to the ones defined in
+ """
+ Mapper that defines the the fluid selector to use, according to the ones defined in
:code:`SourceCode._Solver.FluidSelectors` module, containing its type, evolution equation,
fluid determination and possibly redistancing routines.
-
+
.. rubric:: **Parameters**
-
+
|bs| - |bs| **id** *(integer)* -- the index corresponding to the wished Fluid selector
- |bs| **id** *(params)* -- (optional) the parameters to pass to the selected boundary function routine
-
+
*Implemented boundary conditions*
-
+
0. | Outflow
1. | Inflow (see the inflow doc to know which inflow function and parameters are available)
2. | Wall
3. | Dirichlet
-
+
|
Returns:
FluidSelectors.Module (module): class containing all the fluid's selector tools
-
+
------------------------------------------------------------------------------------- """
if id == 0: return(BCs.Outflow.Outflow)
@@ -350,18 +350,18 @@ Mapper function for the quadratures
-----------------------------------------------------------------------------"""
def Quadrature_Rule(QuadratureType, QuadratureOrder):
"""
- Mapping to the wished quadrature rule
-
+ Mapping to the wished quadrature rule
+
.. note::
-
- So far this routine is inactive
+
+ So far this routine is inactive
"""
# Quadpy based quadrature
if QuadratureType==1:
if QuadratureOrder == 1: return(0)
elif QuadratureOrder == 2: return(0)
-
+
else: raise ValueError("Automatic quadrature: wrong index for the quadrature order.")
# By hand quadrature
diff --git a/SourceCode/PostProcess.py b/SourceCode/PostProcess.py
index f1dc1d4..34c0b62 100644
--- a/SourceCode/PostProcess.py
+++ b/SourceCode/PostProcess.py
@@ -1,6 +1,6 @@
"""
-The predefined set of exports plots are given in wrapper routines located in the file :code:`PostProcess` under the
+The predefined set of exports plots are given in wrapper routines located in the file :code:`PostProcess` under the
:code:`_Solver` folder. This is the main post-processing file that proposes predefined post-treatment tasks
|
@@ -44,24 +44,23 @@ import PostProcessing.Exports as EX
def Plots_StaticSchlieren(*argv):
""" Exports the usually welcome plots after computation, in Matplotlib. The inputs should be given in the
below given order.
-
- Args:
+
+ 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
time (float array-like or string): time points for which to process the plot, or "all"
Resolution (integer): mesh resolution
- Variable (string): the name of the variable one wants the image for
Returns:
- None: The solution's matplotlib plots for each time step given in :code:`time`.
-
+ None: The solution's matplotlib plots for each time step given in :code:`time`.
+
|
-
+
The automatically exported plots (called plotting routines) for this function are for the selected variable:
-
+
- | Compute2DSchlierenImages
-
-
+
+
---------------------------------------------------------------------------------------------------------"""
@@ -91,49 +90,52 @@ def Plots_StaticSchlieren(*argv):
if times.size==0: print("Error: No previous results exported. Pass.")
else:
times = argv[2]
-
+ lst = np.array([len(glob.glob(os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(argv[3]).replace(".","_"), "RawResults","{0:8.6f}.sol".format(tt).replace(".","_")))) for tt in times])
+ if np.any(lst==0): raise ValueError("Error: at least one of the wished time is not exported. Pass.")
+
# Loop on the wished times
for t in times:
# Load the solution
Mesh, ProblemData, Problem, Solution = EX.LoadSolution(ParametersId, argv[3], t)
# Export the schlieren images
- PSL.Compute2DSchlierenImages(Problem, Mesh, Solution, argv[4], ParametersId)
- PSL.Compute2DSchlierenImagesEos(Problem, Mesh, Solution, "p", ParametersId)
+ for VariableId in range(Problem.GoverningEquations.NbVariables):
+ PSL.Compute2DSchlierenImages(Problem, Mesh, Solution, VariableId, ParametersId)
+ PSL.Compute2DSchlierenImagesPrimary(Problem, Mesh, Solution, VariableId, ParametersId)
#### Plotting with Matplotlib static visualisation tools ####
def Plots_StaticMatplotlib(*argv):
""" Exports the usually welcome plots after computation, in Matplotlib. The inputs should be given in the
below given order.
-
- Args:
+
+ 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
time (float array-like or string): time points for which to process the plot, or "all"
Resolution (integer): mesh resolution
Returns:
- None: The solution's matplotlib plots for each time step given in :code:`time`.
-
+ None: The solution's matplotlib plots for each time step given in :code:`time`.
+
|
-
+
The automatically exported plots (called plotting routines) for this function are for each variable:
-
+
- | PlotVariableOnWholeDomainTri2D
- - | PlotVariableOnWholeDomainTri
- - | PlotVariableOnEachSubDomainTri2D
- - | PlotVariableOneEachSubDomainTri
- - | PlotVariableOnWholeDomain2D
- - | PlotVariableOnEachSubDomain2D
- - | PlotVariableOnWholeDomain
- - | PlotVariableOneEachSubDomain
-
+ - | PlotVariableOnWholeDomainTri
+ - | PlotVariableOnEachSubDomainTri2D
+ - | PlotVariableOneEachSubDomainTri
+ - | PlotVariableOnWholeDomain2D
+ - | PlotVariableOnEachSubDomain2D
+ - | PlotVariableOnWholeDomain
+ - | PlotVariableOneEachSubDomain
+
The automatically exported plots (called plotting routines) for the fluid's selector are:
-
+
- | PlotLSValuesOnWholeDomainTri
- | PlotLSValuesTri
- | PlotLSValuesOnWholeDomainTri
-
+
---------------------------------------------------------------------------------------------------------"""
@@ -163,7 +165,10 @@ def Plots_StaticMatplotlib(*argv):
if times.size==0: print("Error: No previous results exported. Pass.")
else:
times = argv[2]
-
+ lst = np.array([len(glob.glob(os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(argv[3]).replace(".","_"), "RawResults","{0:8.6f}.sol".format(tt).replace(".","_")))) for tt in times])
+ if np.any(lst==0): raise ValueError("Error: at least one of the wished time is not exported. Pass.")
+
+
# Loop on the wished times
for t in times:
# Load the solution
@@ -171,50 +176,49 @@ def Plots_StaticMatplotlib(*argv):
# Batch the plotting routines
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)
- 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)
- try: PMP.PlotLSValuesTri(Problem, Mesh, Solution, ParametersId)
- except: pass
- #PMP.PlotLSValuesOnWholeDomainTri(Problem, Mesh, Solution, ParametersId)
+ #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.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.PlotPrimaryVariableOnWholeDomain (Problem, Mesh, Solution, VariableId, ParametersId)
+ #PMP.PlotPrimaryVariableOnWholeDomainTri(Problem, Mesh, Solution, VariableId, ParametersId)
+ #PMP.PlotPrimaryVariableOneEachSubDomain(Problem, Mesh, Solution, VariableId, ParametersId)
+ PMP.PlotPrimaryVariableOneEachSubDomainTri(Problem, Mesh, Solution, VariableId, ParametersId)
+
+ PMP.PlotLSValuesOnWholeDomainTri(Problem, Mesh, Solution, ParametersId)
+ PMP.PlotLSValuesTri(Problem, Mesh, Solution, ParametersId)
#### Plotting with matpotlib some debug tools ####
def Plot_DebugMatplotlib(*argv):
"""
Exports the debug tools for multiphase computation, in Matplotlib (linked to the :code:`Plots_Matplotlib` module).
The inputs should be given in the below given order.
-
- Args:
+
+ 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
time (float array-like or string): time points for which to process the plot, or "all"
Resolution (integer): mesh resolution
Returns:
- None: The solution's matplotlib plots for each time step given in :code:`time`.
-
+ None: The solution's matplotlib plots for each time step given in :code:`time`.
+
|
-
+
The automatically debug (called routines) for this function are:
-
+
- | PlotFlags
- | PlotSubdomains
- | PlotLSubDomains
- | PlotLSValuesTri
-
+
---------------------------------------------------------------------------------------------------------"""
-
-
+
+
# ---------------------- Checking the user input -------------------------------------
# Checking the validity of the user input
if len(argv) < 4:
@@ -234,8 +238,18 @@ def Plot_DebugMatplotlib(*argv):
# Register the name of the test for export purposes
ParametersId = argv[0]+"_"+argv[1]
+ # Get the time stamps for which to generate the solutions' snapshots
+ if argv[2]=="all":
+ lst = glob.glob(os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(argv[3]).replace(".","_"), "RawResults","*.sol"))
+ times = np.sort([float(FileName[-12:-4].replace("_",".")) for FileName in lst])
+ if times.size==0: print("Error: No previous results exported. Pass.")
+ else:
+ times = argv[2]
+ lst = np.array([len(glob.glob(os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(argv[3]).replace(".","_"), "RawResults","{0:8.6f}.sol".format(tt).replace(".","_")))) for tt in times])
+ if np.any(lst==0): raise ValueError("Error: at least one of the wished time is not exported. Pass.")
+
# Loop on the wished times
- for t in argv[2]:
+ for t in times:
# Load the solution
Mesh, ProblemData, Problem, Solution = EX.LoadSolution(ParametersId, argv[3], t)
@@ -244,39 +258,39 @@ def Plot_DebugMatplotlib(*argv):
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 (linked to the :code:`Plots_Plotly` module).
The inputs should be given in the below given order.
-
- Args:
+
+ 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
time (float array-like or string): time points for which to process the plot, or "all"
Resolution (integer): mesh resolution
Returns:
- None: The solution's matplotlib plots for each time step given in :code:`time`.
-
+ None: The solution's matplotlib plots for each time step given in :code:`time`.
+
|
-
+
The automatically exported plots (called plotting routines) for this function are for each variable:
-
+
- | PlotSolution2DOnEachSubDomain
- | PlotSolution2D
- | PlotLevelSetOnEachSubDomain
- | PlotSolutionOnEachSubdomain
- | PlotSolution
-
-
+
+
The automatically exported plots (called plotting routines) for the fluid's selector are:
-
+
- | PlotEachLevelSet
- | PlotEachLevelSetTri
- | PlotAllLevelSets
-
+
---------------------------------------------------------------------------------------------------------"""
@@ -299,8 +313,18 @@ def Plots_StaticPlotly(*argv):
# Register the name of the test for export purposes
ParametersId = argv[0]+"_"+argv[1]
+ # Get the time stamps for which to generate the solutions' snapshots
+ if argv[2]=="all":
+ lst = glob.glob(os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(argv[3]).replace(".","_"), "RawResults","*.sol"))
+ times = np.sort([float(FileName[-12:-4].replace("_",".")) for FileName in lst])
+ if times.size==0: print("Error: No previous results exported. Pass.")
+ else:
+ times = argv[2]
+ lst = np.array([len(glob.glob(os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(argv[3]).replace(".","_"), "RawResults","{0:8.6f}.sol".format(tt).replace(".","_")))) for tt in times])
+ if np.any(lst==0): raise ValueError("Error: at least one of the wished time is not exported. Pass.")
+
# Loop on the wished times
- for t in argv[2]:
+ for t in times:
# Load the solution
Mesh, ProblemData, Problem, Solution = EX.LoadSolution(ParametersId, argv[3], t)
@@ -315,29 +339,29 @@ def Plots_StaticPlotly(*argv):
PLP.PlotLevelSetOnEachSubDomain(Problem, Mesh, Solution, ParametersId)
PLP.PlotSolutionOnEachSubdomain(Problem, Mesh, Solution, VariableId, ParametersId, [np.min(Solution.Sol[VariableId,:]), np.max(Solution.Sol[VariableId,:])])
PLP.PlotSolution(Problem, Mesh, Solution, VariableId, ParametersId)
-
+
#### Plotting with plotly static visualisation tools ####
def Plots_AnimationsPlotly(*argv):
- """
+ """
Exports the usually welcome animation after computation, in Plotly. The inputs should be given in the
below given order.
-
- Args:
+
+ 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): mesh resolution
Returns:
- None: The solution's plotly plot as html animation from all the raw files that have been exported in the results folder.
-
+ None: The solution's plotly plot as html animation from all the raw files that have been exported in the results folder.
+
|
-
+
The automatically exported animations (called plotting routines) for this function are for each variable:
-
+
- | Animate3DSolution
- | Animate3DLevelSet
- | Animate2DLevelSet (very slow, visualisation time slow as well, commented hardly by default)
-
+
---------------------------------------------------------------------------------------------------------"""
diff --git a/SourceCode/PostProcessing/Exports.py b/SourceCode/PostProcessing/Exports.py
index 8cfa711..285741a 100644
--- a/SourceCode/PostProcessing/Exports.py
+++ b/SourceCode/PostProcessing/Exports.py
@@ -2,7 +2,7 @@
The module :code:`Exports` located in the :code:`SourceCode/Postprocessing` furnishes the tools to export and load a solution and problem
that is stored in the solution and problem structures format. It contains binary exports/loads and text format export/loads, text files export compatible with the plt format
-.. note::
+.. note::
The binary exports use pickles, use it by knowing it and aknowledging the risks
@@ -37,14 +37,14 @@ import meshio
#### Clean the export environment ###
def CleanExportEnvironment(TestCase, Resolution):
""" Clean the export environment: removes all the files contained in the designated export folder.
-
- Args:
+
+ Args:
TestCase (string): name of the folder in which the solutions are stored
Resolution (integer): resolution of the associated mesh
-
+
Returns:
None: Files permanently removed from the designated folder
-
+
-----------------------------------------------------------------------------------"""
try:
Folder = os.path.join(Pathes.SolutionPath, TestCase, "Res"+str(Resolution).replace(".","_"))
@@ -55,7 +55,7 @@ def CleanExportEnvironment(TestCase, Resolution):
#### Load the file that has been exported in binary pickle format #####
def LoadSolution(TestCase, Resolution, t):
""" Load the file that has been exported in binary pickle format
-
+
Args:
TestCase (string): name of the folder in which the solutions are stored
Resolution (integer): resolution of the associated mesh
@@ -68,15 +68,15 @@ def LoadSolution(TestCase, Resolution, t):
Solution (Solution structure instance): the actual solution to the problem
.. note:
-
+
This routine loads a binary file that has been generated by pickles. To work, the used structures should already be known by the programm by the loading time.
-
+
----------------------------------------------------------------------------------------- """
# Selecting the file to load
- Folder = os.path.join(Pathes.SolutionPath, TestCase, "Res"+str(Resolution).replace(".","_"), "RawResults")
+ Folder = os.path.join(Pathes.SolutionPath, TestCase, "Res"+str(Resolution).replace(".","_"), "RawResults")
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)
@@ -85,8 +85,8 @@ def LoadSolution(TestCase, Resolution, t):
#### Export the solution in a binary pickle format #####
def ExportSolution(TestCase, Mesh, ProblemData, Problem, Solution):
- """ Export the solution in a binary pickle format
-
+ """ Export the solution in a binary pickle format
+
Args:
TestCase (string): name of the folder in which the solutions are stored
Mesh (MeshStructure instance): mesh on which the solution has been computed
@@ -95,14 +95,15 @@ def ExportSolution(TestCase, Mesh, ProblemData, Problem, Solution):
Solution (Solution structure instance): the actual solution to the problem
.. note::
-
- This routine loads a binary file that has been generated by pickles. To work, the used structures should already be known by the programm by the loading time.
-
+
+ - This routine loads a binary file that has been generated by pickles. To work, the used structures should already be known by the programm by the loading time.
+ - The time encoding format in the outfile is 8.6f, therefore any finer time export will overwrite the just exported results
+
----------------------------------------------------------------------------------------- """
# Creating the export environment
Folder = os.path.join(Pathes.SolutionPath, TestCase, "Res"+str(Mesh.MeshResolution).replace(".","_"), "RawResults")
- try: os.makedirs(Folder)
+ try: os.makedirs(Folder)
except: pass
# Opening the export file
@@ -111,14 +112,14 @@ def ExportSolution(TestCase, Mesh, ProblemData, Problem, Solution):
# Exporting the file in a binary structure
pickle.dump([Mesh, ProblemData, Problem, Solution], ToWrite)
-
+
# Closing the file
ToWrite.close()
#### Export the solution in VTK format through the meshio interface #####
def ExportSolutionVTK(TestCase, Mesh, ProblemData, Problem, Solution):
""" Export the solution in VTK format through the meshio interface.
-
+
Args:
TestCase (string): name of the folder in which the solutions are stored
Mesh (MeshStructure instance): mesh on which the solution has been computed
@@ -128,25 +129,22 @@ def ExportSolutionVTK(TestCase, Mesh, ProblemData, Problem, Solution):
.. note::
This routine is available only for triangulations
-
+
----------------------------------------------------------------------------------------- """
# Creating the export environment
Folder = os.path.join(Pathes.SolutionPath, TestCase, "Res"+str(Mesh.MeshResolution).replace(".","_"), "RawResults")
- try: os.makedirs(Folder)
+ try: os.makedirs(Folder)
except: pass
# Selecting the filename for the export
FileName = os.path.join(Folder, ("SolutionVTK_t{0:8.6f}".format(Solution.t)).replace(".","_")+".vtk")
-
+
# Extract the required information and convert them through meshio interface library
points = np.append(Mesh.CartesianPoints, np.zeros((Mesh.CartesianPoints.shape[0],1)), axis=1)
cells = [("triangle", np.array(Mesh.InnerElements))]
data = {Problem.GoverningEquations.VariablesNames[i]: Solution.RSol[i,:] for i in range(Problem.GoverningEquations.NbVariables)}
data.update({"Fluid#{0:d}".format(i): Solution.RLSValues[i,:] for i in range(Problem.NbFluids)})
-
+
# Exporting the solution in a VTK format
meshio.write_points_cells(FileName, points, cells, point_data = data)
-
-
-
diff --git a/SourceCode/PostProcessing/Plots_Matplotlib.py b/SourceCode/PostProcessing/Plots_Matplotlib.py
index 0347062..f8958b2 100644
--- a/SourceCode/PostProcessing/Plots_Matplotlib.py
+++ b/SourceCode/PostProcessing/Plots_Matplotlib.py
@@ -295,9 +295,10 @@ def PlotFlags(Problem, Mesh, Solution, ParametersId):
plt.plot(Mesh.CartesianPoints[Index,0], Mesh.CartesianPoints[Index,1],\
"+", color = Problem.FluidProp[i].color)
- # Plotting the fluids's contours
- 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])
+ # Plotting the fluids's contours in case of several fluids
+ if Problem.NbFluids>1:
+ 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])
# --------------- Export the figure at the right location -------------------------
@@ -333,12 +334,11 @@ def PlotSubdomains(Problem, Mesh, Solution, ParametersId):
# -------------- Initialising the figure ---------------------------------------
fig = plt.figure()
- ax = fig.add_subplot(1, 1, 1)
+ ax = fig.add_subplot(1, 1, 1)
# --------------- Creating the figure ------------------------------------------
Polys = Solution.Subdomains
- for i in range(len(Polys)):
- plotMultipoly(Polys[i], Problem.FluidProp[i].color, ax)
+ for i in range(len(Polys)): plotMultipoly(Polys[i], Problem.FluidProp[i].color, ax)
# --------------- Customisation properties -----------------------------------
# Customisation of labels and latex conversion
@@ -400,9 +400,13 @@ def PlotLSubDomains(Problem, Mesh, Solution, ParametersId):
Fluids = set(map(int, Solution.RFluidFlag))
Fluids.discard(-1)
for i in Fluids:
- ax.tricontour(Mesh.CartesianPoints[:,0], Mesh.CartesianPoints[:,1], Mesh.InnerElements, \
- Solution.RLSValues[i,:], [0],\
- colors=Problem.FluidProp[i].color, linewidths=len(Fluids)-i+2)
+ # Getting the unique values of the level set
+ UniqueValues = np.array(list(set(Solution.RLSValues[i,:])))
+
+ # Plot the global fluid-separation contour if relevant
+ if np.any(UniqueValues<=0):
+ ax.tricontour(Mesh.CartesianPoints[:,0], Mesh.CartesianPoints[:,1], Mesh.InnerElements, \
+ Solution.RLSValues[i,:], [0], colors=Problem.FluidProp[i].color, linewidths=len(Fluids)-i+2)
# --------------- Customisation properties -----------------------------------
@@ -462,21 +466,31 @@ def PlotLSValuesTri(Problem, Mesh, Solution, ParametersId):
# Span each fluid
for i in Fluids:
- # Plot the global fluid-separation contour
- ax.tricontour(Mesh.CartesianPoints[:,0], Mesh.CartesianPoints[:,1], Mesh.InnerElements, \
+ # Get the unique values of the fluid
+ UniqueValues = np.array(list(set(Solution.RLSValues[i,:])))
+
+ # Plot the global fluid-separation contour if relevant
+ if np.any(UniqueValues<=0):
+ ax.tricontour(Mesh.CartesianPoints[:,0], Mesh.CartesianPoints[:,1], Mesh.InnerElements, \
Solution.RLSValues[i,:], [0],
colors="k", linewidths=2)
- # Select where the values are positive and specify the contours to take into account
- ind = np.where(Solution.RLSValues[i,:]>0)[0]
+ # Select where the values are non-negative and specify the contours to take into account
+ ind = np.where(Solution.RLSValues[i,:]>=0)[0]
mm = np.max(Solution.RLSValues[i,ind])
lvls = np.linspace(0,mm,400)
- # Getting a surface from the vectors (hack from tricontourf to avoid patch tricontour and interpolation)
- Ct = ax.tricontourf(Mesh.CartesianPoints[:,0], Mesh.CartesianPoints[:,1], Mesh.InnerElements,\
- Solution.RLSValues[i,:],lvls,\
- cmap = ColorMapsMatplotlib[i])
-
+ # In the case of several fluids, getting a surface from the vectors (hack from tricontourf to avoid patch tricontour and interpolation)
+ if len(UniqueValues)>1:
+ Ct = ax.tricontourf(Mesh.CartesianPoints[:,0], Mesh.CartesianPoints[:,1], Mesh.InnerElements,\
+ Solution.RLSValues[i,:],lvls,\
+ cmap = ColorMapsMatplotlib[i])
+
+ # Else, just patch with a single color
+ else:
+ ax.tripcolor(Mesh.CartesianPoints[:,0], Mesh.CartesianPoints[:,1], Mesh.InnerElements,\
+ 0*Solution.RLSValues[i,:], cmap = ColorMapsMatplotlib[i])
+
# --------------- Customisation properties -----------------------------------
# Customisation of labels and latex conversion
rc('font',size = 12)
@@ -485,9 +499,9 @@ def PlotLSValuesTri(Problem, Mesh, Solution, ParametersId):
# 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_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_ha('left') for t in ax.get_xticklabels()]
# Set the grid and the pane colors
fig.suptitle(r'Solution for the variable at the time ${0:5.3f}$'.format(Solution.t))
@@ -497,7 +511,7 @@ def PlotLSValuesTri(Problem, Mesh, Solution, ParametersId):
# --------------- Export the figure at the right location -------------------------
FolderName = os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(Mesh.MeshResolution).replace(".","_"), "DomainPlots")
- try: os.makedirs(FolderName)
+ try: os.makedirs(FolderName)
except: pass
# Set the figure properties and shows it
@@ -539,7 +553,7 @@ def PlotLSValuesOnWholeDomainTri(Problem, Mesh, Solution, ParametersId):
ax = fig.gca(projection='3d')
ax.plot_trisurf(X, Y, Z, triangles=Mesh.InnerElements, cmap=ColorMapsMatplotlib[0],\
edgecolor=None ,linewidth=0, antialiased=False)
- #ax.tricontour(X,Y,Mesh.InnerElements , Z , [0], colors="r", offset=0.5 )
+ #ax.tricontour(X,Y,Mesh.InnerElements, Z , [0], colors="r", offset=0.5 )
# --------------- Customisation properties -----------------------------------
# Retrieve the variables names
@@ -612,10 +626,12 @@ def PlotVariableOnWholeDomainTri2D(Problem, Mesh, Solution, VariableId, Paramete
fig = plt.figure()
ax = fig.subplots()
- # Still plotting the contours of the LS
+ # Still plotting the contours of the LS in case of several fluids
for i in set(map(int, Solution.RFluidFlag)):
- ax.tricontour(X, Y, Mesh.InnerElements, Solution.RLSValues[i,:], \
- levels = [0], colors = "k", linewidths=2)
+ UniqueValues = np.array(list(set(Solution.RLSValues[i,:])))
+ if np.any(UniqueValues<=0):
+ ax.tricontour(X, Y, Mesh.InnerElements, Solution.RLSValues[i,:], \
+ levels = [0], colors = "k", linewidths=2)
# Plotting the solution homogeneously on the full subdomain
ax.tripcolor(X, Y, Z, triangles=Mesh.InnerElements, cmap=ColorMapsMatplotlib[0],\
@@ -623,9 +639,9 @@ def PlotVariableOnWholeDomainTri2D(Problem, Mesh, Solution, VariableId, Paramete
# --------------- Customisation properties -----------------------------------
# Retrieve the variables names
- variableName = Problem.GoverningEquations.VariablesNames[VariableId]
+ variableName = Problem.GoverningEquations.VariablesNames[VariableId]
variableLatexName = Problem.GoverningEquations.VariablesLatexNames[VariableId]
- variableLatexUnit = Problem.GoverningEquations.VariablesLatexUnits[VariableId]
+ variableLatexUnit = Problem.GoverningEquations.VariablesLatexUnits[VariableId]
# Customisation of labels and latex conversion
rc('font',size = 12)
@@ -645,7 +661,7 @@ def PlotVariableOnWholeDomainTri2D(Problem, Mesh, Solution, VariableId, Paramete
# --------------- Export the figure at the right location -------------------------
FolderName = os.path.join(Pathes.SolutionPath, ParametersId, "Res"+str(Mesh.MeshResolution).replace(".","_"), "2DPlots")
- try: os.makedirs(FolderName)
+ try: os.makedirs(FolderName)
except: pass
FileName = "Mplt_2DSol_OnWholeDomainTri_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
@@ -679,12 +695,15 @@ def PlotVariableOnEachSubDomainTri2D(Problem, Mesh, Solution, VariableId, Parame
fig = plt.figure()
ax = fig.subplots()
- # Still plotting the contours of the LS
+ # Spanning the fluids
for i in set(map(int, Solution.RFluidFlag)):
- ax.tricontour(X, Y, Mesh.InnerElements, LS[i,:], \
- colors = "k", levels = [0], linewidths=2)
+
+ # Still plotting the contours of the LS in case of several fluids
+ UniqueValues = np.array(list(set(Solution.RLSValues[i,:])))
+ if np.any(UniqueValues<=0): ax.tricontour(X, Y, Mesh.InnerElements, LS[i,:], colors = "k", levels = [0], linewidths=2)
- ind = np.where(LS[i,:]<=0)[0]
+ # Plotting the masked values for each fluid
+ ind = np.where(LS[i,:]<-1e-16)[0]
Ele = np.array(Mesh.InnerElements)
Mask = [bool(set(ele).intersection(ind)==set(ele)) for ele in Ele]
@@ -825,7 +844,7 @@ def PlotVariableOneEachSubDomainTri(Problem, Mesh, Solution, VariableId, Paramet
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]
+ ind = np.where(LS[i,:]<-1e-16)[0]
Ele = np.array(Mesh.InnerElements)
Mask = [set(ele).intersection(ind)==set(ele) for ele in Ele]
@@ -872,14 +891,14 @@ def PlotVariableOneEachSubDomainTri(Problem, Mesh, Solution, VariableId, Paramet
plt.close()
#### Plot the variable retrieved from the EoS independently of the subdomain ####
-def PlotVariableEoSOnWholeDomainTri(Problem, Mesh, Solution, Field, ParametersId):
+def PlotPrimaryVariableOnWholeDomainTri(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
+ Field (integer): the index of the primary variable to plot the solution for
ParametersId (string): name of the export folder (TestCase_TestSettings)
Returns:
@@ -893,23 +912,20 @@ def PlotVariableEoSOnWholeDomainTri(Problem, Mesh, Solution, Field, ParametersId
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)
+ # Computing the value from the conservative variables
+ Zbuf = Problem.GoverningEquations.ConservativeToPrimary(Z, Solution.RFluidFlag)[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)
+ 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
+ variableName = Problem.GoverningEquations.PrimaryVariablesNames[Field]
+ variableLatexName = Problem.GoverningEquations.PrimaryVariablesLatexNames[Field]
+ variableLatexUnit = Problem.GoverningEquations.PrimaryVariablesLatexUnits[Field]
# Customisation of labels and latex conversion
rc('font',size = 12)
@@ -941,19 +957,19 @@ def PlotVariableEoSOnWholeDomainTri(Problem, Mesh, Solution, Field, ParametersId
try: os.makedirs(FolderName)
except: pass
- FileName = "Mplt_3DSolEoS_OnWholeDomainTri_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
+ FileName = "Mplt_3DSol_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):
+def PlotPrimaryVariableOneEachSubDomainTri(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
+ Field (integer): the index of the primary variable one wants the plot for
ParametersId (string): name of the export folder (TestCase_TestSettings)
Returns:
@@ -968,12 +984,10 @@ def PlotVariableEoSOneEachSubDomainTri(Problem, Mesh, Solution, Field, Parameter
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)
+ # Computing the value from the conservative variables
+ Zbuf = Problem.GoverningEquations.ConservativeToPrimary(Z, Solution.RFluidFlag)[Field,:]
+
# --------------- Creating the figure ------------------------------------------
fig = plt.figure()
ax = fig.gca(projection='3d')
@@ -983,7 +997,7 @@ def PlotVariableEoSOneEachSubDomainTri(Problem, Mesh, Solution, Field, Parameter
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]
+ ind = np.where(LS[i,:]<-1e-16)[0]
Ele = np.array(Mesh.InnerElements)
Mask = [set(ele).intersection(ind)==set(ele) for ele in Ele]
@@ -992,10 +1006,10 @@ def PlotVariableEoSOneEachSubDomainTri(Problem, Mesh, Solution, Field, Parameter
# --------------- Customisation properties -----------------------------------
# Retrieve the variables names
- variableName = Field
- variableLatexName = Field
- variableLatexUnit = Field
-
+ variableName = Problem.GoverningEquations.PrimaryVariablesNames[Field]
+ variableLatexName = Problem.GoverningEquations.PrimaryVariablesLatexNames[Field]
+ variableLatexUnit = Problem.GoverningEquations.PrimaryVariablesLatexUnits[Field]
+
# Customisation of labels and latex conversion
rc('font',size = 12)
rc('font',family = 'serif')
@@ -1025,7 +1039,7 @@ def PlotVariableEoSOneEachSubDomainTri(Problem, Mesh, Solution, Field, Parameter
try: os.makedirs(FolderName)
except: pass
- FileName = "Mplt_3DSolEoS_OnSubDomainsTri_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
+ FileName = "Mplt_3DSol_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()
@@ -1051,9 +1065,9 @@ def PlotVariableOnWholeDomain2D(Problem, Mesh, Solution, VariableId, ParametersI
# --------------- Preparing the data to plot -----------------------------------
# Extracting the data of interest
- X = Mesh.CartesianPoints[:,0]
- Y = Mesh.CartesianPoints[:,1]
- Z = Solution.RSol[VariableId,:]
+ X = Mesh.CartesianPoints[:,0]
+ Y = Mesh.CartesianPoints[:,1]
+ Z = Solution.RSol[VariableId,:]
LS = Solution.RLSValues[:,:]
# --------------- Creating the figure ------------------------------------------
@@ -1073,9 +1087,14 @@ def PlotVariableOnWholeDomain2D(Problem, Mesh, Solution, VariableId, ParametersI
# Still plotting the contours of the LS
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,:])
- ax.contour(xx, yy, lzz, levels = [0], colors = "k", linewidths=2)
+ # Get the unique values of the fluid
+ UniqueValues = np.array(list(set(Solution.RLSValues[i,:])))
+
+ # Plot the global fluid-separation contour if relevant
+ if np.any(UniqueValues<=0):
+ # Interpolating the fluid selector on the visualisation grid
+ xx, yy, lzz = GetMeshgridFom3VectorsMaskX(Mesh, X, Y, LS[i,:])
+ ax.contour(xx, yy, lzz, levels = [0], colors = "k", linewidths=2)
# --------------- Customisation properties -----------------------------------
# Retrieve the variables names
@@ -1154,7 +1173,11 @@ def PlotVariableOnEachSubDomain2D(Problem, Mesh, Solution, VariableId, Parameter
# Plotting the surface
ax.pcolormesh(xx1, yy1, zz1, shading="gouraud", alpha=0.7, cmap=ColorMapsMatplotlib[i], \
vmin=np.min(Z)-0.2, vmax=np.max(Z)+0.2)
- ax.contour(xx, yy, lzz, levels = [0], colors = "k", linewidths=2)
+
+ # Get the unique values of the fluid
+ UniqueValues = np.array(list(set(Solution.RLSValues[i,:])))
+ # Plot the global fluid-separation contour if relevant
+ if np.any(UniqueValues<=0): ax.contour(xx, yy, lzz, levels = [0], colors = "k", linewidths=2)
# --------------- Customisation properties -----------------------------------
# Retrieve the variables names
@@ -1295,7 +1318,7 @@ def PlotVariableOneEachSubDomain(Problem, Mesh, Solution, VariableId, Parameters
xx, yy, lzz = GetMeshgridFom3VectorsMaskX(Mesh, X, Y, LS[i,:])
# Masking the points that lie outisde the fluid's area
- ind = np.where(lzz<=0)
+ ind = np.where(lzz<-1e-16)
xx[ind] = None
yy[ind] = None
@@ -1342,14 +1365,14 @@ def PlotVariableOneEachSubDomain(Problem, Mesh, Solution, VariableId, Parameters
plt.close()
#### Plot the variable coming from the EoS independently of the subdomain ####
-def PlotVariableEoSOnWholeDomain(Problem, Mesh, Solution, Field, ParametersId):
+def PlotPrimaryVariableOnWholeDomain(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
+ Field (integer): the index of the primary variable to plot
ParametersId (string): name of the export folder (TestCase_TestSettings)
Returns:
@@ -1361,15 +1384,10 @@ def PlotVariableEoSOnWholeDomain(Problem, Mesh, Solution, Field, ParametersId):
# 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]))
+ Z = Solution.RSol[:,:]
- 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)
-
+ # Computing the value from the conservative variables
+ Zbuf = Problem.GoverningEquations.ConservativeToPrimary(Z, Solution.RFluidFlag)[Field,:]
# --------------- Creating the figure ------------------------------------------
fig = plt.figure()
@@ -1385,9 +1403,9 @@ def PlotVariableEoSOnWholeDomain(Problem, Mesh, Solution, Field, ParametersId):
# --------------- Customisation properties -----------------------------------
# Retrieve the variables names
- variableName = Field
- variableLatexName = Field
- variableLatexUnit = Field
+ variableName = Problem.GoverningEquations.PrimaryVariablesNames[Field]
+ variableLatexName = Problem.GoverningEquations.PrimaryVariablesLatexNames[Field]
+ variableLatexUnit = Problem.GoverningEquations.PrimaryVariablesLatexUnits[Field]
# Customisation of labels and latex conversion
rc('font',size = 12)
@@ -1419,19 +1437,19 @@ def PlotVariableEoSOnWholeDomain(Problem, Mesh, Solution, Field, ParametersId):
try: os.makedirs(FolderName)
except: pass
- FileName = "Mplt_3DSolEoS_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()
#### Plot one variable coming from the EoS with a color scale depending on the subdomain ####
-def PlotVariableEoSOneEachSubDomain(Problem, Mesh, Solution, Field, ParametersId):
+def PlotPrimaryVariableOneEachSubDomain(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
+ Field (integer): the index of the primary variable to plot
ParametersId (string): name of the export folder (TestCase_TestSettings)
Returns:
@@ -1443,16 +1461,12 @@ def PlotVariableEoSOneEachSubDomain(Problem, Mesh, Solution, Field, ParametersId
# Extracting the data of interest
X = Mesh.CartesianPoints[:,0]
Y = Mesh.CartesianPoints[:,1]
- Z = Solution.RSol[VariableId,:]
+ Z = Solution.RSol[:,:]
LS = Solution.RLSValues[:,:]
- # Computing the value from the EoS
- Zbuf = np.zeros((Z.shape[1]))
+ # Computing the value from the conservative variables
+ Zbuf = Problem.GoverningEquations.ConservativeToPrimary(Z, Solution.RFluidFlag)[Field,:]
- 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()
@@ -1469,7 +1483,7 @@ def PlotVariableEoSOneEachSubDomain(Problem, Mesh, Solution, Field, ParametersId
xx, yy, lzz = GetMeshgridFom3VectorsMaskX(Mesh, X, Y, LS[i,:])
# Masking the points that lie outisde the fluid's area
- ind = np.where(lzz<=0)
+ ind = np.where(lzz<-1e-16)
xx[ind] = None
yy[ind] = None
@@ -1478,9 +1492,9 @@ def PlotVariableEoSOneEachSubDomain(Problem, Mesh, Solution, Field, ParametersId
# --------------- Customisation properties -----------------------------------
# Retrieve the variables names
- variableName = Field
- variableLatexName = Field
- variableLatexUnit = Field
+ variableName = Problem.GoverningEquations.PrimaryVariablesNames[Field]
+ variableLatexName = Problem.GoverningEquations.PrimaryVariablesLatexNames[Field]
+ variableLatexUnit = Problem.GoverningEquations.PrimaryVariablesLatexUnits[Field]
# Customisation of labels and latex conversion
rc('font',size = 12)
@@ -1511,7 +1525,7 @@ def PlotVariableEoSOneEachSubDomain(Problem, Mesh, Solution, Field, ParametersId
try: os.makedirs(FolderName)
except: pass
- FileName = "Mplt_3DSolEoS_OnSubDomains_{0!s}_t{1:5.3f}".format(variableName, Solution.t).replace(".","_")
+ FileName = "Mplt_3DSol_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()
diff --git a/SourceCode/PostProcessing/Plots_Plotly.py b/SourceCode/PostProcessing/Plots_Plotly.py
index 497a1e0..245b9c1 100644
--- a/SourceCode/PostProcessing/Plots_Plotly.py
+++ b/SourceCode/PostProcessing/Plots_Plotly.py
@@ -11,6 +11,12 @@ and animation routines.
The module is made to be used on-line, that is either while computing the solution (except videos) or after an exported solution has been loaded
+|
+
+.. warning::
+
+ When the level set is constant (i.e. 1 fluid), the 3D plot of the level set may fail
+
|
|
@@ -82,7 +88,7 @@ CSSColor = [ 'greenyellow', 'aliceblue', 'lightgreen', 'antiquewhite', 'aqua', '
#### Converts three vectors to a surface, masking the xy coodinates when holes ####
def GetMeshgridFom3VectorsMaskX(Mesh, X, Y, Z):
""" Interpolation routine in order to plot.
-
+
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
Args:
@@ -92,19 +98,19 @@ def GetMeshgridFom3VectorsMaskX(Mesh, X, Y, Z):
Z (float array like): the data's (scalar) value at the considered points
.. rubric:: Returns
-
+
- |bs| **xgrid** *(float 2D numpy array)* -- x meshgrid of the points where the solution Z has been interpolated
- |bs| **ygrid** *(float 2D numpy array)* -- y meshgrid of the points where the solution Z has been interpolated
- |bs| **zgrid** *(float 2D numpy array)* -- values of the interpolated Z at the generated grid points
-
+
|
-
+
.. note::
-
+
- When holes in the domains, the xy values are masked with None in both corresponding meshgrids
- The x resolution of the interpolation grid is set to 5*(xmax-xmin)/np.min(Mesh.InnerElementsDiameter), and respectively for the y axis
- The interpolation method is set to "linear", with a fill value of "-1e-15" in order to fit the 0 crossing at a physical boundary in case of interpolating a level set.
-
+
----------------------------------------------------------------------------------------"""
# ------------------------------- Initialisations -----------------------------------------
@@ -157,7 +163,7 @@ def GetMeshgridFom3VectorsMaskX(Mesh, X, Y, Z):
#### Converts three vectors to a surface, masking the z coodinates when holes ####
def GetMeshgridFom3VectorsMaskZ(Mesh, X, Y, Z):
""" Interpolation routine in order to plot.
-
+
Args:
Mesh (MeshStructure): the mesh in use
X (float array like): the x coordinates of the considered points
@@ -165,19 +171,19 @@ def GetMeshgridFom3VectorsMaskZ(Mesh, X, Y, Z):
Z (float array like): the data's (scalar) value at the considered points
.. rubric:: Returns
-
+
- |bs| **xgrid** *(float 2D numpy array)* -- x meshgrid of the points where the solution Z has been interpolated
- |bs| **ygrid** *(float 2D numpy array)* -- y meshgrid of the points where the solution Z has been interpolated
- |bs| **zgrid** *(float 2D numpy array)* -- values of the interpolated Z at the generated grid points
|
-
+
.. note::
-
+
- When holes in the domains, the z values are masked with None in both corresponding meshgrids
- The x resolution of the interpolation grid is set to 5*(xmax-xmin)/np.min(Mesh.InnerElementsDiameter), and respectively for the y axis
- The interpolation method is set to "linear", with a fill value of "-1e-15" in order to fit the 0 crossing at a physical boundary in case of interpolating a level set.
-
+
----------------------------------------------------------------------------------------"""
# ------------------------------- Initialisations -----------------------------------------
@@ -233,8 +239,8 @@ def GetMeshgridFom3VectorsMaskZ(Mesh, X, Y, Z):
#### Generic plotting routine from plotly that considers any 3D frame ####
def GenerateFigurePlotly3D(Frame, bounds, Title = "Plot Title", Zaxis = "Zaxis"):
""" A generic 3D plotly routine that encapsulates a plotted object in a figure.
-
- Args:
+
+ Args:
Frame (Plotly object): any plotly content that can be incorporated in a figure
bounds (float array-like): the xmin-xmax bounds of the plot (the same will be used for the y axis, guaranteeing the equal option: always use the larger data to not cut-off the picture unwinlingly)
title (string): (optional) the title of the figure
@@ -244,9 +250,9 @@ def GenerateFigurePlotly3D(Frame, bounds, Title = "Plot Title", Zaxis = "Zaxis")
fig (plotly figure instance): the figure in Plotly format
.. note::
-
+
The figure is only generated, not shown neither exported
-
+
----------------------------------------------------------------------------------------"""
# Creating a figure instance
@@ -273,19 +279,19 @@ def GenerateFigurePlotly3D(Frame, bounds, Title = "Plot Title", Zaxis = "Zaxis")
#### Generic plotting routine from plotly ####
def GenerateFigurePlotly2D(Frame, bounds, Title="Plot Title"):
""" A generic 2D plotly routine that encapsulates a plotted object in a figure.
-
- Args:
+
+ Args:
Frame (Plotly object): any plotly content that can be incorporated in a figure
bounds (float array-like): the xmin-xmax bounds of the plot (the same will be used for the y axis, guaranteeing the equal option: always use the larger data to not cut-off the picture unwinlingly)
title (string): (optional) the title of the figure
-
+
Returns:
fig (plotly figure instance): the figure in Plotly format
.. note::
-
+
The figure is only generated, not shown neither exported
-
+
----------------------------------------------------------------------------------------"""
# Creating a figure instance
@@ -312,7 +318,7 @@ def GenerateFigurePlotly2D(Frame, bounds, Title="Plot Title"):
def PlotAnimation3D(TestName, Resolution, FrameFunc, cmin, cmax, args, Title = "Plot Title", Zaxis = "Zaxis"):
""" Routine that creates a 3D plotly hmtl animation given the Frame function
framefunc and the solutions that are stored in the TestName folder.
-
+
Args:
TestName (string): the name of the folder where the results are stored
Resolution (integer): the resolution of the associated mesh
@@ -322,7 +328,7 @@ def PlotAnimation3D(TestName, Resolution, FrameFunc, cmin, cmax, args, Title = "
args (array-like): containing the further arguments to be passed to FrameFunc
Title (string): (optional) the title of the animation
Zaxis (string): (optional) the label of the z axis
-
+
Returns:
fig (plotly animation): a plotly animation
"""
@@ -347,7 +353,7 @@ def PlotAnimation3D(TestName, Resolution, FrameFunc, cmin, cmax, args, Title = "
# Setting the xy bounds upon it (assuming non-growing domain)
xmin, xmax = np.min(Mesh.CartesianPoints[:,0]), np.max(Mesh.CartesianPoints[:,0])
ymin, ymax = np.min(Mesh.CartesianPoints[:,1]), np.max(Mesh.CartesianPoints[:,1])
-
+
# Initialising figure and setting the first snapshot up
fig = go.Figure(frames=Frames)
for Frame in FrameContent: fig.add_trace(Frame)
@@ -403,7 +409,7 @@ def PlotAnimation3D(TestName, Resolution, FrameFunc, cmin, cmax, args, Title = "
def PlotAnimation2D(TestName, Resolution, FrameFunc, cmin, cmax, args, Title = "Plot Title"):
""" Routine that creates a 2D plotly hmtl animation given the Frame function
framefunc and the solutions that are stored in the TestName folder.
-
+
Args:
TestName (string): the name of the folder where the results are stored
Resolution (integer): the resolution of the associated mesh
@@ -412,10 +418,10 @@ def PlotAnimation2D(TestName, Resolution, FrameFunc, cmin, cmax, args, Title = "
cmax (float): the maximum value of the plotted quantities
args (array-like): containing the further arguments to be passed to FrameFunc
Title (string): (optional) the title of the animation
-
+
Returns:
fig (plotly animation): a plotly animation
-
+
"""
# ----------- Collecting the solution figures from the given data -------------------
@@ -496,8 +502,8 @@ def PlotAnimation2D(TestName, Resolution, FrameFunc, cmin, cmax, args, Title = "
#### Plot all the level sets on one single 3D plot #####
def PlotAllLevelSets(Mesh, Solution, ParametersId):
""" Plots all the level sets on the computational domain.
-
- Args:
+
+ Args:
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
ParametersId (string): the string forming the "ProblemName_ParameterName", mapping to the right solution export folder
@@ -506,9 +512,9 @@ def PlotAllLevelSets(Mesh, Solution, ParametersId):
None: The 3D plot of all the levels set on the considered domain are exported in the problem-dependent solution export folder.
.. note::
-
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -562,8 +568,8 @@ def PlotAllLevelSets(Mesh, Solution, ParametersId):
#### Generates the frame corresponding to a level set ####
def GenerateFrameEachLevelSet(Problem, Mesh, Solution, cbounds, VariableId, ParametersId):
""" Generates a figure's frame of the slevel set variableId on the full computational domain.
-
- Args:
+
+ Args:
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
@@ -575,9 +581,9 @@ def GenerateFrameEachLevelSet(Problem, Mesh, Solution, cbounds, VariableId, Para
None: The frame representing the solution on the considered domain with a subdomain-dependent color scale.
.. note::
-
+
The routine has been split from PlotSolutionOnEachSubdomain to ease the construction of a video frame later
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -601,7 +607,7 @@ def GenerateFrameEachLevelSet(Problem, Mesh, Solution, cbounds, VariableId, Para
# Getting a surface from the vectors
xx,yy,zz = GetMeshgridFom3VectorsMaskZ(Mesh,X,Y,Z[VariableId,:])
-
+
# Generating the surface and their 0 level
Frame = go.Surface(x = copy.deepcopy(xx),
y = copy.deepcopy(yy),
@@ -609,15 +615,15 @@ def GenerateFrameEachLevelSet(Problem, Mesh, Solution, cbounds, VariableId, Para
contours = {"z": {"show": True, "start": 0, "end": 0.1, "size": 0.1}},
colorscale = "haline",
name = Name)
-
+
# Returning the constructed frame
return([Frame])
#### Plot each level set on different 3D plots ####
def PlotEachLevelSet(Mesh, Solution, ParametersId):
""" Plots the level sets in independent figures, each on the full computational domain.
-
- Args:
+
+ Args:
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
ParametersId (string): the string forming the "ProblemName_ParameterName", mapping to the right solution export folder
@@ -625,10 +631,10 @@ def PlotEachLevelSet(Mesh, Solution, ParametersId):
Returns:
None: The 3D plot of all the levels set on the considered domain, exported in the problem-dependent solution export folder.
- .. note::
-
+ .. note::
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -676,20 +682,20 @@ def PlotEachLevelSet(Mesh, Solution, ParametersId):
#### Plot each level set on different 3D plots for triangular meshes ####
def PlotEachLevelSetTri(Mesh, Solution, ParametersId):
""" Plots the level sets in independent figures, each on the full computational domain.
-
- Args:
+
+ Args:
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
ParametersId (string): the string forming the "ProblemName_ParameterName", mapping to the right solution export folder
Returns:
-
+
None: The 3D plot of all the levels set on the considered domain, exported in the problem-dependent solution export folder.
- .. note::
-
+ .. note::
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -708,13 +714,13 @@ def PlotEachLevelSetTri(Mesh, Solution, ParametersId):
except: pass
# --------------- Creating the figure ------------------------------------------------
-
+
# Plotting the LS both one by one and on the same figure
for Id in range(np.shape(Z)[0]):
# Generating the name of the surface
Name = "LS n°{0:d}".format(Id+1)
-
+
# Generating the surface and their 0 level
fig = ff.create_trisurf(x=X, y=Y, z=Z,
simplices = Mesh.InnerElements,
@@ -734,7 +740,7 @@ def PlotEachLevelSetTri(Mesh, Solution, ParametersId):
font = dict(family = "Arial",
size = 14,
color = "#7f7f7f"))
-
+
# --------------- Export the figure at the right location -------------------------
# Creating the file name and exporting it in the right folder
try:
@@ -750,7 +756,7 @@ def PlotEachLevelSetTri(Mesh, Solution, ParametersId):
#### Plot a variable of the solution regardless the subdomain ####
def PlotSolution(Problem, Mesh, Solution, VariableId, ParametersId):
""" Plots the solution's variableId on the full computational domain.
-
+
Args:
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh
@@ -762,9 +768,9 @@ def PlotSolution(Problem, Mesh, Solution, VariableId, ParametersId):
The 3D plot of all the given solution set on the considered domain, exported in the problem-dependent solution export folder.
.. note::
-
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -812,22 +818,22 @@ def PlotSolution(Problem, Mesh, Solution, VariableId, ParametersId):
#### Generate frame of the solution with a colorscale per subdomain ####
def GenerateFrameSolutionOnEachSubdomain(Problem, Mesh, Solution, cbounds, VariableId, ParametersId):
""" Generates a figure's frame of the solution's variableId on the full computational domain with a color scale per subdomain.
-
- Args:
+
+ Args:
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
cbounds (float array-like): values determining the bounds of the color scale to use
VariableId, (integer): the index of the desired variable to be plot
ParametersId (string): the string forming the "ProblemName_ParameterName", mapping to the right solution export folder
-
+
Returns:
None: The frame representing the solution on the considered domain with a subdomain-dependent color scale.
- .. note::
-
+ .. note::
+
The routine has been split from PlotSolutionOnEachSubdomain to ease the construction of a video frame later
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -836,7 +842,7 @@ def GenerateFrameSolutionOnEachSubdomain(Problem, Mesh, Solution, cbounds, Varia
Y = Mesh.CartesianPoints[:,1]
Z = Solution.RSol[VariableId,:]
LS = Solution.RLSValues[:,:]
-
+
# -------------------- Customisation properties --------------------------------------
# Retrieve the variables names
variableName = Problem.GoverningEquations.VariablesNames[VariableId]
@@ -852,16 +858,16 @@ def GenerateFrameSolutionOnEachSubdomain(Problem, Mesh, Solution, cbounds, Varia
# Generating one surface per subfluid
FrameTot = []
-
+
for i in set(map(int,Solution.FluidFlag)):
# Selecting the part of the solution that lie in the fluid's subdomain
_,_,lzz = GetMeshgridFom3VectorsMaskX(Mesh,X,Y,LS[i,:])
ind = np.where(lzz < 0)
-
+
# Masking the disregarded fluid in the span
lz = copy.deepcopy(zz)
lz[ind] = None
-
+
# Generating the surface
Frame = go.Surface(x = copy.deepcopy(xx),
y = copy.deepcopy(yy),
@@ -878,8 +884,8 @@ def GenerateFrameSolutionOnEachSubdomain(Problem, Mesh, Solution, cbounds, Varia
#### Plot a variable of the solution with a color scale per subdomain ####
def PlotSolutionOnEachSubdomain(Problem, Mesh, Solution, VariableId, ParametersId, cbounds):
""" Generates a figure's frame of the solution's variableId on the full computational domain with a color scale per subdomain.
-
- Args:
+
+ Args:
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
@@ -891,9 +897,9 @@ def PlotSolutionOnEachSubdomain(Problem, Mesh, Solution, VariableId, ParametersI
None: The frame representing the solution on the considered domain with a subdomain-dependent color scale.
.. note::
-
+
The routine has been split from PlotSolutionOnEachSubdomain to ease the construction of a video frame later
-
+
--------------------------------------------------------------------------------------"""
# ----------------- Preparing the export at the right location -------------------------
@@ -912,9 +918,9 @@ def PlotSolutionOnEachSubdomain(Problem, Mesh, Solution, VariableId, ParametersI
# Visibility bounds for the figure
width = abs(np.min([*Mesh.CartesianPoints[:,0], *Mesh.CartesianPoints[:,1]])-np.max([*Mesh.CartesianPoints[:,0], *Mesh.CartesianPoints[:,1]]))
bounds = [np.min([*Mesh.CartesianPoints[:,0], *Mesh.CartesianPoints[:,1]]),np.max([*Mesh.CartesianPoints[:,0], *Mesh.CartesianPoints[:,1]])]
-
+
FrameTot = GenerateFrameSolutionOnEachSubdomain(Problem, Mesh, Solution, cbounds, VariableId, ParametersId)
- fig = GenerateFigurePlotly3D(FrameTot, bounds, variableName)
+ fig = GenerateFigurePlotly3D(FrameTot, bounds, variableName)
# --------------- Export the figure at the right location -------------------------
# Creating the results folder if not already existing
@@ -938,21 +944,21 @@ def PlotSolutionOnEachSubdomain(Problem, Mesh, Solution, VariableId, ParametersI
#### Plot the Level sets w.r.t. the subdomain ####
def PlotLevelSetOnEachSubDomain(Problem, Mesh, Solution, ParametersId):
""" Plots the positive part of each level set on their corresponding subdomain.
-
- Args:
+
+ Args:
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
ParametersId (string): the string forming the "ProblemName_ParameterName", mapping to the right solution export folder
Returns:
-
+
None: The 2D plot of all the levels set on their corresponding subdomains, exported in the problem-dependent solution export folder.
.. note::
-
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -1023,8 +1029,8 @@ def PlotLevelSetOnEachSubDomain(Problem, Mesh, Solution, ParametersId):
#### Generate a frame corresponding to a level set of the solution w.r.t. the subdomain ####
def GenerateFrameLevelSetOnEachSubDomain(Problem, Mesh, Solution, cbounds, VariableId, ParametersId):
""" Plots the positive part of each level set on their corresponding subdomain.
-
- Args:
+
+ Args:
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
@@ -1036,9 +1042,9 @@ def GenerateFrameLevelSetOnEachSubDomain(Problem, Mesh, Solution, cbounds, Varia
None: The 2D plot of all the levels set on their corresponding subdomains, exported in the problem-dependent solution export folder.
.. note::
-
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -1102,21 +1108,21 @@ def GenerateFrameLevelSetOnEachSubDomain(Problem, Mesh, Solution, cbounds, Varia
#### Plot a variable of the solution regardless the subdomain ####
def PlotSolution2D(Problem, Mesh, Solution, VariableId, ParametersId):
""" Plots the solution's variableId on the full computational domain in 2D.
-
- Args:
+
+ Args:
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
VariableId, (integer): the index of the desired variable to be plot
ParametersId (string): the string forming the "ProblemName_ParameterName", mapping to the right solution export folder
-
+
Returns:
None: The 2D plot of all the given solution set on the considered domain, exported in the problem-dependent solution export folder.
- .. note::
-
+ .. note::
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -1188,8 +1194,8 @@ def PlotSolution2D(Problem, Mesh, Solution, VariableId, ParametersId):
#### Plot a variable of the solution regardless the subdomain ####
def PlotSolution2DOnEachSubDomain(Problem, Mesh, Solution, VariableId, ParametersId):
""" Plots the solution's variableId on the full computational domain in 2D, with a subdomain-dependent colorscale.
-
- Args:
+
+ Args:
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh
Solution (Solution): the solution of the considered problem
@@ -1200,9 +1206,9 @@ def PlotSolution2DOnEachSubDomain(Problem, Mesh, Solution, VariableId, Parameter
None: The 2D plot of all the given solution set on the considered domain, exported in the problem-dependent solution export folder.
.. note::
-
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ----------------------- Retrieving the solution's values ---------------------------
@@ -1286,7 +1292,7 @@ def PlotSolution2DOnEachSubDomain(Problem, Mesh, Solution, VariableId, Parameter
### Animate each Level set of an exported solution accross the time ####
def Animate3DLevelSet(TestName, Resolution):
""" Animates the level set for each, on the full domain.
-
+
Args:
TestName (string): the string forming the "ProblemName_ParameterName", mapping to the right solution export folder
Resolution (integer): resolution of the considered mesh
@@ -1294,10 +1300,10 @@ def Animate3DLevelSet(TestName, Resolution):
Returns:
None: The 3D animation is exported an html format
- .. note::
-
+ .. note::
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ---------------- Creating the color bounds for the global plot -------------------
@@ -1325,8 +1331,11 @@ def Animate3DLevelSet(TestName, Resolution):
args = [i, TestName]
# Getting the bound range
- cmax = np.max([cmax, *Solution.RLSValues[i,:]])
- cmin = np.min([cmin, *Solution.RLSValues[i,:]])
+ rr = (np.max([*Solution.RLSValues[i,:]])-np.min([*Solution.RLSValues[i,:]]))
+ if np.isclose(rr,0): rr = 1
+
+ cmax = np.max([cmax, *Solution.RLSValues[i,:]])+0.1*rr
+ cmin = np.min([cmin, *Solution.RLSValues[i,:]])-0.1*rr
# Getting the animation
fig = PlotAnimation3D(TestName, Resolution, GenerateFrameEachLevelSet, cmin, cmax, args)
@@ -1347,19 +1356,19 @@ def Animate3DLevelSet(TestName, Resolution):
### Animate all the levels sets on the domain accross the time ####
def Animate2DLevelSet(TestName, Resolution):
""" Animates all the levels sets on the domain accross the time, on their respective positive parts.
-
+
Args:
TestName (string): the string forming the "ProblemName_ParameterName", mapping to the right solution export folder
Resolution (integer): resolution of the considered mesh
Returns:
-
+
None: Only exports the 2D animation in an html format
- .. note::
-
+ .. note::
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ---------------- Creating the color bounds for the global plot -------------------
@@ -1387,9 +1396,12 @@ def Animate2DLevelSet(TestName, Resolution):
args = [i, TestName]
# Getting the bound range
- cmax = np.max([cmax, *Solution.RLSValues[i,:]])
- cmin = np.min([cmin, *Solution.RLSValues[i,:]])
-
+ rr = (np.max([*Solution.RLSValues[i,:]])-np.min([*Solution.RLSValues[i,:]]))
+ if np.isclose(rr,0): rr = 0.1
+
+ cmax = np.max([cmax, *Solution.RLSValues[i,:]])+0.1*rr
+ cmin = np.min([cmin, *Solution.RLSValues[i,:]])-0.1*rr
+
# Getting the animation
fig = PlotAnimation2D(TestName, Resolution, GenerateFrameLevelSetOnEachSubDomain, cmin, cmax, args)
@@ -1412,7 +1424,7 @@ def Animate2DLevelSet(TestName, Resolution):
### Animate each variable of an exported solution ####
def Animate3DSolution(TestName, Resolution):
""" Animates the solution for each variable.
-
+
Args:
TestName (string): the string forming the "ProblemName_ParameterName", mapping to the right solution export folder
Resolution (integer): resolution of the considered mesh
@@ -1421,9 +1433,9 @@ def Animate3DSolution(TestName, Resolution):
None: Only exports the 3D animation in an html format
.. note::
-
+
The figure is only exported, not shown
-
+
--------------------------------------------------------------------------------------"""
# ---------------- Creating the color bounds for the global plot -------------------
@@ -1455,9 +1467,12 @@ def Animate3DSolution(TestName, Resolution):
args = [i, TestName]
# Getting the bound range
- cmax = np.max([cmax, *Solution.RSol[i,:]])+0.1
- cmin = np.min([cmin, *Solution.RSol[i,:]])-0.1
-
+ rr = (np.max([*Solution.RSol[i,:]])-np.min([*Solution.RSol[i,:]]))
+ if np.isclose(rr,0): rr = 0.1
+
+ cmax = np.max([cmax, *Solution.RSol[i,:]])+0.1*rr
+ cmin = np.min([cmin, *Solution.RSol[i,:]])-0.1*rr
+
# Getting the animation
fig = PlotAnimation3D(TestName, Resolution, GenerateFrameSolutionOnEachSubdomain, cmin, cmax, args)
diff --git a/SourceCode/PostProcessing/Plots_Schlieren.py b/SourceCode/PostProcessing/Plots_Schlieren.py
index b41bd2b..8d56e6e 100644
--- a/SourceCode/PostProcessing/Plots_Schlieren.py
+++ b/SourceCode/PostProcessing/Plots_Schlieren.py
@@ -5,7 +5,8 @@ The module :code:`Plots_Schileren` located in the :code:`SourceCode/Postprocessi
|
.. note::
- The module is made to be used on-line, that is either while computing the solution (except videos) or after the solution has been exported
+ - The module is made to be used on-line, that is either while computing the solution (except videos) or after the solution has been exported
+ - The plots will render properly only on decently fine meshes
|
|
@@ -201,18 +202,14 @@ def Compute2DSchlierenImages(Problem, Mesh, Solution, Field, ParametersId):
.. note::
- So far only avaialble for two fluids
-
- """
+ So far only available up to three fluids because of the color selection
- # -------------------- Check if the wished variable is indeed in the representable ones ------------------------
- try:
- FieldId = Problem.GoverningEquations.VariablesNames.index(Field)
- ok = 1
- except:
- ok = 0
- print("Warning: No Schlieren image exported, the variable name is unknown from the problem setting.")
+ .. warning::
+ On some variables for which the gradient's magnitude is null, it will fail
+
+
+ """
# ----------------------- Defining the static color plotting vectors -------------------------------------------
@@ -230,68 +227,70 @@ def Compute2DSchlierenImages(Problem, Mesh, Solution, Field, ParametersId):
cmap += [ListedColormap(vals)]
# ----------------------- Computing the Schlieren numerical images -----------------------------------------------
- # If we can export the plot for the wished variable
- if ok:
-
- # Getting the desired values and there spatial location
- X = Mesh.CartesianPoints[:,0]
- Y = Mesh.CartesianPoints[:,1]
- Z = Solution.RSol[FieldId,:]
- xx,yy,zz = GetMeshgridFom3VectorsMaskZ(Mesh, X, Y, Z)
-
- # 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 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)
+ # Getting the desired values and there spatial location
+ X = Mesh.CartesianPoints[:,0]
+ Y = Mesh.CartesianPoints[:,1]
+ Z = Solution.RSol[Field,:]
+ xx,yy,zz = GetMeshgridFom3VectorsMaskZ(Mesh, X, Y, Z)
+
+ # 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 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)
+
+ # 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 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)
-
- # 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])
+ # 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
- # --------------- 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
+ # Plot the quantity
+ plt.pcolor(xx,yy,Phi, cmap=cmap[i])
+
+
+ # --------------- Customize the figure --------------------------------------------
+ variableName = Problem.GoverningEquations.PrimaryVariablesNames[Field]
+ variableLatexName = Problem.GoverningEquations.PrimaryVariablesLatexNames[Field]
+ plt.title("Computational Schlieren image for the variable {0!s} at time {1:8.6f}".format(variableLatexName, Solution.t))
+
+ # --------------- 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_{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()
+ # Set the figure properties and shows it
+ FileName = "Mplt_Schlieren_{1!s}_t{0:5.3f}".format(Solution.t, variableName).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):
+def Compute2DSchlierenImagesPrimary(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
+ Field (integer): the index of the primary variable to derive the schlieren graph for
ParametersId (string): name of the export folder (TestCase_TestSettings)
Returns:
@@ -301,7 +300,7 @@ def Compute2DSchlierenImagesEos(Problem, Mesh, Solution, Field, ParametersId):
.. warning::
- Only valid for Euler equations, the other test case will be compatible on call but the results would be meaningless
+ On some variables for which the gradient's magnitude is null, it will fail
"""
@@ -326,12 +325,10 @@ def Compute2DSchlierenImagesEos(Problem, Mesh, Solution, Field, ParametersId):
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)
+ # Compute the primary variables
+ Zbuf = Problem.GoverningEquations.ConservativeToPrimary(Z, Solution.RFluidFlag)[Field,:]
+
# Mapping them to a cartesian grid
xx,yy,zz = GetMeshgridFom3VectorsMaskZ(Mesh, X, Y, Zbuf)
@@ -355,7 +352,7 @@ def Compute2DSchlierenImagesEos(Problem, Mesh, Solution, Field, ParametersId):
# 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)
@@ -367,6 +364,11 @@ def Compute2DSchlierenImagesEos(Problem, Mesh, Solution, Field, ParametersId):
# Plot the quantity
plt.pcolor(xx,yy,Phi, cmap=cmap[i])
+ # --------------- Customize the figure --------------------------------------------
+ variableName = Problem.GoverningEquations.PrimaryVariablesNames[Field]
+ variableLatexName = Problem.GoverningEquations.PrimaryVariablesLatexNames[Field]
+ plt.title("Computational Schlieren image for the variable {0!s} at time {1:8.6f}".format(variableLatexName, Solution.t))
+
# --------------- Export the figure at the right location -------------------------
# Creating the export environment
@@ -375,7 +377,7 @@ def Compute2DSchlierenImagesEos(Problem, Mesh, Solution, Field, ParametersId):
except: pass
# Set the figure properties and shows it
- FileName = "Mplt_Schlieren_Eos_{1!s}_t{0:5.3f}".format(Solution.t, Field).replace(".","_")
+ FileName = "Mplt_Schlieren_{1!s}_t{0:5.3f}".format(Solution.t, variableName).replace(".","_")
plt.savefig(os.path.join(FolderName, FileName)+".png", transparent=True, dpi=300)
plt.close()
diff --git a/SourceCode/ProblemDefinition/BoundaryConditions/Wall.py b/SourceCode/ProblemDefinition/BoundaryConditions/Wall.py
index e8ec676..65c55bf 100644
--- a/SourceCode/ProblemDefinition/BoundaryConditions/Wall.py
+++ b/SourceCode/ProblemDefinition/BoundaryConditions/Wall.py
@@ -13,7 +13,7 @@ def Wall(Problem, Mesh, Solution, Residuals, FluidId, Bd, BdIndex):
"""
Routine that does nothing, just lets outflow (there for compatibility)
- Args:
+ Args:
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh
Solution (Solution): the considered solution (as of the previous subtimestep (as used in the temporal scheme))
@@ -21,10 +21,9 @@ def Wall(Problem, Mesh, Solution, Residuals, FluidId, Bd, BdIndex):
FluidId (integer): the index of the considered fluid (warning: NOT its value in mappers)
Bd (integer): the tag of the boundary the considered points are lying on
BdIndex (integer list): the list of indices of the degrees of freedom to treat
- Constant (floar): the constant value to force on the degrees of freedom
-
+
Returns:
- res (2D float numpy array): the updated residuals
+ res (2D float numpy array): the updated residuals
|
@@ -36,7 +35,7 @@ def Wall(Problem, Mesh, Solution, Residuals, FluidId, Bd, BdIndex):
"""
-
+
# ------------------- Determine the wall conditions depending on the layout of the given Dofs -------------
# Retrieve the cartesian coordinates of the given Dofs
XY = Mesh.DofsXY[BdIndex,:]
diff --git a/SourceCode/ProblemDefinition/FluidModels.py b/SourceCode/ProblemDefinition/FluidModels.py
index 3cc3478..872f288 100644
--- a/SourceCode/ProblemDefinition/FluidModels.py
+++ b/SourceCode/ProblemDefinition/FluidModels.py
@@ -4,7 +4,7 @@
|
All the predefined Fluids' models (classes) that are casting
-the fluid's properties in order to be passed to a governing equations
+the fluid's properties in order to be passed to a governing equations
are given in the file :code:`SourceCode/ProblemDefinition/FluidModels.py`
|
@@ -19,7 +19,7 @@ are given in the file :code:`SourceCode/ProblemDefinition/FluidModels.py`
------------------------------------------
-Implemented models
+Implemented models
-------------------
|
@@ -32,18 +32,18 @@ Implemented models
# ==============================================================================
class Dummy():
- """ Dummy model that can be used to test the any simple problems or debug
- problem that does not rely on any Fluid's properties.
- This is a dummy model, it only contains a :code:`color` attribute for the visualisation process,
+ """ Dummy model that can be used to test the any simple problems or debug
+ problem that does not rely on any Fluid's properties.
+ This is a dummy model, it only contains a :code:`color` attribute for the visualisation process,
filled upon instance creation.
-
+
.. note::
-
+
- This is only a product-type class: no method is available
- At instance creation it takes as many input (whose values will be trashed) as wished for call time compatibility
-
+
|
-
+
"""
def __init__(self, color, *args):
@@ -51,21 +51,22 @@ class Dummy():
Args:
color (string or matplotlib color scheme): the color of the fluid for plotting
"""
-
+
self.color = color
+
class SimpleFluid():
""" SimpleFluid model that can be used to test the EulerEquations or Advection problems,
- for either the stiff or idel EOS. It fills the fluid's model fields :code:`color`,
+ for either the stiff or idel EOS. It fills the fluid's model fields :code:`color`,
:code:`gamma` and :code:`pinf` upon instance creation.
-
+
.. note::
-
+
- This is only a product-type class: no method is available
- It takes as many input (whose values will be trashed) as wished for call time compatibility
-
+
|
-
+
"""
def __init__(self, color, gamma, pinf, R, *args):
@@ -74,9 +75,9 @@ class SimpleFluid():
color (string or matplotlib color scheme): the color of the fluid for plotting
gamma (float): model attribute
pinf (float): model attribute
- """
-
+ """
+
self.gamma = gamma
self.pinf = pinf
- self.R = R
+ self.R = R
self.color = color
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py b/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py
index 729a4b0..f48c45d 100644
--- a/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py
+++ b/SourceCode/ProblemDefinition/GoverningEquations/EulerEquations.py
@@ -65,9 +65,13 @@ class Equations():
self.VariablesUnits = ["kg/m3", "m2/s", "m2/s", "J"]
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.VariablesLatexUnits = [r'$kg/m^3$', r'$kg/(m^2s)$', r'$kg/(m^2s)$', r'$J$']
self.XYLatexUnits = [r'$m$', r'$m$']
+ self.PrimaryVariablesLatexUnits = [r'$kg/m^3$', r'$m/s$', r'$m/s$', r'$Pa$']
+ self.PrimaryVariablesLatexNames = [r'$\rho$', r'$U_x$', r'$U_y$', r'$p$']
+ self.PrimaryVariablesNames = ["rho", "Ux", "Uy", "p"]
+
# Register the shortcut of the fluid's properties
self.FluidProp = FluidProp
self.EOs = EOs
@@ -97,7 +101,7 @@ class Equations():
# 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))
+ Fluids = map(int, list(set(FluidIndex)))
# Construct the fluidproperties depending o the right fluid selector's value and evaluating the
# EOS accordingly
@@ -116,7 +120,7 @@ class Equations():
return(PrimVars)
#### Retrieve the conservative variables from the primary ones ####
- def PrimaryToConservative(self, PrimVar, FluidIndex, *arrgs):
+ def PrimaryToConservative(self, PrimVar, FluidIndex, *args):
"""Converts the primary variables to the conservative ones
Args:
@@ -152,7 +156,7 @@ class Equations():
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 ----------------------------------------------------
+ # ------- Constructing the conservative variables ----------------------------------------------------
# Computing the primary variables
Vars[0,:] = PrimVar[0,:]
Vars[1,:] = PrimVar[1,:] # Directly from the previous semiconversion
@@ -287,7 +291,7 @@ class Equations():
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))
+ Fluids = map(int, list(set(FluidIndex)))
# Construct the fluidproperties depending o the right fluid selector's value and evaluating the
# EOS accordingly
@@ -345,7 +349,7 @@ class Equations():
c = np.zeros(len(Var[0,:]))
# List of fluids present in the list to be computed
- Fluids = list(set(FluidIndex))
+ Fluids = map(int,list(set(FluidIndex)))
# Getting the fluids properties depending on the fluid selector's value and evaluate the
# EOS accordingly
@@ -410,7 +414,7 @@ class Equations():
# --- EOS-dependent quantities
# List of fluids present in the list to be computed
- Fluids = list(set(FluidIndex))
+ Fluids = map(int, list(set(FluidIndex)))
# Getting the fluids properties depending on the fluid selector's value and evaluate the EOS accordingly
for Fluid in Fluids:
@@ -527,7 +531,7 @@ class Equations():
# 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
@@ -584,7 +588,7 @@ class Equations():
gmm = np.zeros(len(Var[0,:]))
# List of fluids present in the list to be computed
- Fluids = list(set(FluidIndex))
+ Fluids = map(int, list(set(FluidIndex)))
# Getting the fluids properties depending on the fluid selector's value and evaluate the
# EOS accordingly
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py
index c1d7f0a..546f541 100644
--- a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py
+++ b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionRotation.py
@@ -70,6 +70,10 @@ class Equations():
self.VariablesLatexUnits = [r'$m$']
self.XYLatexUnits = [r'$m$', r'$m$']
+ self.PrimaryVariablesLatexUnits = [r'$m$']
+ self.PrimaryVariablesLatexNames = [r'$h$']
+ self.PrimaryVariablesNames = ["h"]
+
# Register the shortcut of the fluid's properties
self.FluidProp = FluidProp
self.EOs = EOs
@@ -98,7 +102,7 @@ class Equations():
return(PrimVars)
#### Retrieve the conservative variables from the primary ones ####
- def PrimaryToConservative(self, PrimVar, FluidIndex, *arrgs):
+ def PrimaryToConservative(self, PrimVar, FluidIndex, *args):
"""Converts the primary variables to the conservative ones
Args:
@@ -113,7 +117,7 @@ class Equations():
- | *args is there only for compatibility reason at call time
"""
- # ------- Constructing the primary variables ----------------------------------------------------
+ # ------- Constructing the conservative variables ----------------------------------------------------
# Computing the primary variables
Vars = PrimVar[:]
@@ -273,7 +277,7 @@ class Equations():
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
+ lbd (numpy array): (NbGivenPoints x NbEigs) the eigenvalues at each given point
.. note::
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py
index 935c205..354d304 100644
--- a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py
+++ b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslation.py
@@ -71,6 +71,10 @@ class Equations():
self.VariablesLatexUnits = [r'$m$']
self.XYLatexUnits = [r'$m$', r'$m$']
+ self.PrimaryVariablesLatexUnits = [r'$m$']
+ self.PrimaryVariablesLatexNames = [r'$h$']
+ self.PrimaryVariablesNames = ["h"]
+
# Register the shortcut of the fluid's properties
self.FluidProp = FluidProp
self.EOs = EOs
@@ -99,7 +103,7 @@ class Equations():
return(PrimVars)
#### Retrieve the conservative variables from the primary ones ####
- def PrimaryToConservative(self, PrimVar, FluidIndex, *arrgs):
+ def PrimaryToConservative(self, PrimVar, FluidIndex, *args):
"""Converts the primary variables to the conservative ones
Args:
@@ -114,7 +118,7 @@ class Equations():
- | *args is there only for compatibility reason at call time
"""
- # ------- Constructing the primary variables ----------------------------------------------------
+ # ------- Constructing the conservative variables ----------------------------------------------------
# Computing the primary variables
Vars = PrimVar[:]
@@ -271,7 +275,7 @@ class Equations():
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
+ lbd (numpy array): (NbGivenPoints x NbEig) the four eigenvalues at each given point
.. note::
@@ -283,7 +287,7 @@ class Equations():
# --- 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)
diff --git a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslationVectorTest.py b/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslationVectorTest.py
deleted file mode 100644
index ba8ca1a..0000000
--- a/SourceCode/ProblemDefinition/GoverningEquations/LinearAdvectionTranslationVectorTest.py
+++ /dev/null
@@ -1,492 +0,0 @@
-"""
-
-Module that defines all the necessary routines that implements a
-simple 2D linear in a system shape for testing purposes
-
-.. rubric:: Contains (in scrolling order)
-
-- Directly accessible parameters: FluidModel, string, Name of the fluidModel to use
-- Equations class
-- EOS class
-- Initial conditions class
-
-"""
-
-# ==============================================================================
-# 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)
-
- .. rubric:: Parameters
-
- .. |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)
- - |bs| **EOS** *(list of EOS)* -- the list of Equation of state instances (NbFluids)
-
- It further contains the following attributes, filled upon initialisation
-
- - |bs| **NbVariables** *(integer)* -- number of variables of the model
- - |bs| **VariablesNames** *(list of strings)* -- ordered name of the variables
- - |bs| **VariablesUnits** *(list of strings)* -- symbols of the units of the variables
- - |bs| **VariablesLatexNames** *(list of strings)* -- ordered name of the variables, latex encoding
- - |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):
- """ Initialisation routine upon the considered fluid's properties (to have a direct access).
-
- Args:
- FluidProp (list of FluidModel): the list of fluid property instances (NbFluids)
- EOS (list of EOS): the list of Equation of state instances (NbFluids)
-
- ------------------------------------------------------------------------------------------"""
-
- # Defines the (conservative) variables number and names, in a latex typography
- self.NbVariables = 2
- self.VariablesNames = ["h1", "h2"]
- self.VariablesUnits = ["m", "m"]
-
- self.VariablesLatexNames = [r'$h1$', r'$h2$']
- self.VariablesLatexUnits = [r'$m$', 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)
-
- .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- 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:
- 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)
- - | *args is there only for compatibility reason at call time
-
- -----------------------------------------------------------------------------------------------------"""
-
- # --------- Computing the Euler flux directly in conservative variables ------------------------
- # Giving the flux furnishing the roation
- Fx = np.array([Var[1,:], Var[0,:]]) # Displacement in x direction
- Fy = 0*np.array([Var[0,:], Var[1,:]]) # 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, *args):
- """ 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:
- 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
-
- # Returning the value
- return(UV)
-
- #### Jacobian of the equations set ####
- def Jacobian(self, Var, x, *args):
- """ Computes the Jacobian of the flux for the LinearAdvection (rotation case)
-
- .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
- 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.
-
- |
-
- .. 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((2, 2, 2, NbPoints))
-
- # ------ Computing the actual Jacobian ---------------------------------
- # Creating the Jacobian for the first coordinate
- J[0, :, 0, :] = np.array([np.zeros(NbPoints),np.ones(NbPoints)])
- J[1, :, 0, :] = np.array([np.ones(NbPoints),np.zeros(NbPoints)])
-
- # Creating the Jacobian for the second coordinate
- J[0, :, 1, :] = np.array([np.zeros(NbPoints),np.zeros(NbPoints)])
- J[1, :, 1, :] = np.array([np.zeros(NbPoints),np.zeros(NbPoints)])
-
- # Returning the Jacobian
- return(J)
-
- #### 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)
-
- Output:
- 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)
-
- # 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 :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.
-
- 1. | Dummy
-
- .. Note::
-
- After having been initialised the instance, the EOS is accessible from the field :code:`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.
-
- 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.
-
- """
-
- # Mapping towards the right function
- if Id==1:
- self.EoS = self.Dummy
- self.EoSName = "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.
-
- 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:
- val: the variable retrieved from the EoS (here identity as dummy function)
-
- ----------------------------------------------------------------------------------------------------"""
-
- # 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.
-
- Access the routine computing the initial conditions through the field :code:`IC`
-
- --------------------------------------------------------------------------------------"""
-
- #### Automatic initialisation routine (mapper) ####
- def __init__(self, Id, *params):
- """ Automatic initialisation routine.
-
- 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
-
- Returns:
- None: fills the field :code:`IC` with the function callback corresponding to the wished initialisation function
-
- 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
-
- """
-
- # Mapping the initialisation routines
- if Id == 0: self.IC = self.ConstantState
- elif Id == 3: self.IC = self.ConstantState3
- elif Id == 1: self.IC = self.ConstantState2
- elif Id == 2: self.IC = self.Bump
- elif Id == 4: self.IC = self.BumpEdge
-
- 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.
-
- 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
- Init = np.ones((2, len(PointsID)))
-
- # Defining the raw initilisation variable
- Init[0,:] = 1#0.5
- Init[1,:] = 2
- # 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:
- 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
- Init = np.ones((2, len(PointsID)))
-
- # Defining the raw initilisation variable
- Init[0,:] = 0
- Init[1,:] = 0
-
- # 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:
- 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
- Init = np.ones((2, len(PointsID)))
-
- # Defining the raw initilisation variable
- Init[0,:] = 1
- Init[1,:] = 0
-
- # Returning the computed values
- return(Init)
-
- #### Initialise the solution to a Bump centred in the middle of the subdmain
- 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).
-
- 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((2, 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
-
- # 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
- Init[1,:] = -co
-
- # Returning the computed values
- return(Init)
-
- #### Initialise the solution to a Bump in the centre y of the subdomain, next to the max of the x span
- def BumpEdge(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:
- 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((2, len(PointsID)))
-
- # Retrieving the xy values of the points and the geometrical informations from the subdomain
- 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
- Init[1,:] = -co
-
- # Returning the computed values
- return(Init)
diff --git a/SourceCode/Solve.py b/SourceCode/Solve.py
index f710ada..6b12e2f 100644
--- a/SourceCode/Solve.py
+++ b/SourceCode/Solve.py
@@ -34,6 +34,7 @@ sys.tracebacklimit = None
# Import mathy library
import numpy as np
+import math
# Add to python path the location of the modules and meshes
import IncludePathes
@@ -136,21 +137,19 @@ def Solve(*argv):
\n Starting time dependent computations \
\n -------------------------------------")
- # Initialising the simple running time
- Solution.t=0
+ # Initialising the simple running time and the export periodicity index
+ Solution.t = 0
+ lastexport = 0
# Export the initial considition after having reconstructed the initial values at the physical vertices
- Solver.Reconstructor(Solution)
+ for Reconstructor in Solver.Reconstructor: Reconstructor(Solution)
EX.ExportSolution(TestName, Mesh, ProblemData, Problem, Solution)
# 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
@@ -164,12 +163,13 @@ def Solve(*argv):
# ----- 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)
+ for Reconstructor in Solver.Reconstructor: Reconstructor(Solution)
# ----- Inform the user and export the solution when relevant
- if np.mod(k,1)==0:
+ # May yield to duplicate if the time step is less than 1e-11
+ if (Solution.t-lastexport >= SolverParams.ExportInterval) or np.isclose(Solution.t, Solver.Tmax):
+ lastexport = Solution.t
PO.Printout_MinMaxSolVal(Solution.t, Problem, Solution)
EX.ExportSolutionVTK(TestName, Mesh, ProblemData, Problem, Solution)
EX.ExportSolution(TestName, Mesh, ProblemData, Problem, Solution)
diff --git a/SourceCode/_Solver/FluidSelectors/NaiveLS_CG.py b/SourceCode/_Solver/FluidSelectors/NaiveLS_CG.py
index 0e2b41c..5c78adc 100644
--- a/SourceCode/_Solver/FluidSelectors/NaiveLS_CG.py
+++ b/SourceCode/_Solver/FluidSelectors/NaiveLS_CG.py
@@ -23,20 +23,20 @@ import BarycentricTools as BT
# ==============================================================================
class FS():
""".. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Classical level set approach, initialised according to a signed distance and iterated
with a continuous galerkin - lfx scheme, without redistancing.
-
+
.. 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). It also contains the
+ the structure with the same names (see the parameters below). It also contains the
further field
-
+
- |bs| **FluidSelectorName** *(string)* -- the name of the level set
-
+
|
-
+
"""
#### Automatic initialisation routine ####
@@ -46,17 +46,17 @@ class FS():
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh instance
Params (string list): the level set's parameters as wished by the user
-
+
|
-
+
.. rubric:: Methods
"""
-
-
+
+
# Setting the name of the class to be accessible from outside
self.FluidSelectorName = "Naive Level Set -- Continuous Galerkin"
-
- # Setting the mesh and parameters for an internal access
+
+ # Setting the mesh and parameters for an internal access
self.Problem = Problem
self.Mesh = Mesh
self.Params = Params
@@ -68,13 +68,13 @@ class FS():
#### External initialisation routine ####
def Initialiser(self, Solution):
""" External initialisation routine.
-
+
Args:
Solution (Solution structure): the solution instance the scheme is working on
-
+
Returns:
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 ----------------------------------------------
@@ -101,28 +101,32 @@ class FS():
+[tuple(tuple(a) for a in np.array(interior.coords.xy).T) for interior in sub.interiors]
DomainBorder = ShapelyMultiLineString(Borders)
- # Looping on the intial number of fluids present in the computational domain
- for i in range(NbFluids):
- # Get the subdomain as a shapely multipolygon
- FluidArea = Subdomains[i]
-
- # Extract the boundary of the fluid
- Borders = []
- for sub in FluidArea:
- Borders += [ tuple(tuple(a) for a in np.array(sub.exterior.coords.xy).T)]\
- +[tuple(tuple(a) for a in np.array(interior.coords.xy).T) for interior in sub.interiors]
- FluidBorder = ShapelyMultiLineString(Borders)
-
- # Mask the parts that coincide with the domain's boundary itself
- FluidBorder = FluidBorder.difference(DomainBorder)
-
- # Computing the distance and location w.r.t. the subdomain through the shapely information
- Value = np.array([[ShapelyPoint(tuple(pt)).distance(FluidBorder),\
- 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)
+ # If we have only one fluid, set the LSValues everywhere equal to 1 for compatibility, otherwise fill the LSValues structure
+ if NbFluids==1: LSValues[:,:] = 1.
+ else:
+ # Looping on the intial number of fluids present in the computational domain
+ for i in range(NbFluids):
+
+ # Get the subdomain as a shapely multipolygon
+ FluidArea = Subdomains[i]
+
+ # Extract the boundary of the fluid
+ Borders = []
+ for sub in FluidArea:
+ Borders += [tuple(tuple(a) for a in np.array(sub.exterior.coords.xy).T)]\
+ +[tuple(tuple(a) for a in np.array(interior.coords.xy).T) for interior in sub.interiors]
+ FluidBorder = ShapelyMultiLineString(Borders)
+
+ # Mask the parts that coincide with the domain's boundary itself
+ FluidBorder = FluidBorder.difference(DomainBorder)
+
+ # Computing the distance and location w.r.t. the subdomain through the shapely information
+ Value = np.array([[ShapelyPoint(tuple(pt)).distance(FluidBorder),\
+ 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)
@@ -130,180 +134,180 @@ class FS():
#### Gives the dynamic of the Fluid's spotter through its flux ####
def Flux(self, Var, Motion, *args):
""" Flux function of the (conservative) EulerEquations
-
- Args:
+
+ Args:
Var (float numpy array): the value of the Level Set values (NbFluids x NbGivenPoints)
- Motion (float numpy array): the velocity array (2xNbPoints)
-
+ Motion (float numpy array): the velocity array (2xNbPoints)
+
Returns:
Flux (3D numpy array): the obtained flux (spatialdim x NbVars x NbPoints)
-
- .. note::
-
+
+ .. note::
+
*args is there for compatibility reason when e.g. a location xy is given when called
-
+
--------------------------------------------------------------------------------------------------------------------------------"""
# Determining the advection point-wise
Flux = np.array([Motion[0,:]*Var, Motion[1,:]*Var])
-
+
# Returning the values
return(Flux)
-
+
#### Jacobian of the equations set ####
def Jacobian(self, Var, Motion, *args):
""" .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Computes the Jacobian of the flux for the advection of the LevelSet
-
+
Args:
Var (float numpy array): the value of the Level Set values (NbFluids x NbGivenPoints)
- Motion (float numpy array): the velocity array (2xNbPoints)
-
- Returns:
- J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
-
+ Motion (float numpy array): the velocity array (2xNbPoints)
+
+ Returns:
+ J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
+
|
-
+
.. note::
-
+
- |bs| 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]
-
+
- |bs| *args is there for compatibility reason when e.g. a location xy is given when called the x-y locations at which the variables are given (NbGivenPoints x 2)
-
+
---------------------------------------------------------------------"""
-
+
# ------ Initialising the Jacobian structure ---------------------------
# Getting the number of furnished points
NbPoints = np.shape(Motion)[1]
J = np.zeros((self.Problem.NbFluids, self.Problem.NbFluids, 2, NbPoints))
-
+
for nf in range(self.Problem.NbFluids):
J[nf, nf, 0, :] = np.array(Motion[0,:])
J[nf, nf, 1, :] = np.array(Motion[1,:])
-
- # Returning the Jacobian
- return(J)
-
+
+ # Returning the Jacobian
+ return(J)
+
#### Spectral radius of the equations set ####
def SpectralRadius(self, Var, Motion, n, *args):
""" Computes the spectral radius associated to the motion flux.
-
- Args:
+
+ Args:
Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
- Motion (float numpy array): the velocity array (2xNbPoints)
+ Motion (float numpy array): the velocity array (2xNbPoints)
n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
-
- Returns:
+
+ Returns:
Lmb (numpy array): the array containing the spectral radius at each given point
-
+
.. note::
-
+
*args is there for compatibility reason when e.g. a the xy locations at which the variables are given (NbGivenPoints x 2)
-
+
---------------------------------------------------------------------"""
-
+
# Retrieving the Jacobian
J = self.Jacobian(Var, Motion)
-
+
# 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(self.Problem.NbFluids)]), axis=0)
-
+
# Returning the obtained values
return(Lmb)
-
+
# ************************************************************************
# * Definition of the Fluid spotter's spatial impact and rectification *
# ************************************************************************
-
+
#### 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
+ """ Flags the degrees of freedom according to the different values of the
associated level set.
-
- Args:
+
+ Args:
Solution (Solution structure): the solution instance the scheme is working on
-
- Returns:
+
+ Returns:
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
+ """ Method modifying the subdomains associated to the Solution upon the values of
LSValues.
-
- Args:
+
+ Args:
Solution (Solution): the solution instance the scheme is working on
-
- Returns:
+
+ Returns:
None: updating directly the subdomains in the solution instance
-
- .. warning:: Not safe yet if the first path given out is internal
-
+
+ .. warning:: Not safe yet if the first path given out is internal
+
---------------------------------------------------------------------------- """
-
+
# ------- Updating the subdomains from the countour values -----------------
# 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):
# ---------- Initialisations --------------------------------------------
# Buffering the LSvalues to modify them on the edges upon our wishes to infer closed contours
Buffer = Solution.LSValues[:,:]
-
+
# ----------- 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
+ # 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.
-
- Args:
+
+ Args:
Solution (Solution): the solution instance the scheme is working on
-
- Returns:
+
+ Returns:
None: updating directly the LSValues in the solution instance
-
+
----------------------------------------------------------------------------------"""
# Redistancing the points with the nearest
@@ -314,187 +318,187 @@ class FS():
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)
# ****************************************************************************
# * Definition of the Fluid spotter's evolution's process (e.g. advection) *
# ****************************************************************************
-
+
#### Emulates a flux computation as done in the iteration routine ####
def ComputeFlux(self, U, LSValues):
""" Emulates a flux computation as done in the iteration routine in order
to be used in the DeC subtimestep evaluation.
-
- Args:
+
+ Args:
U (numpy float array): buffer containing the current state values at each dof (NbVars x NbDofs)
LS (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 requires by the UpdateFSValues 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
fluxes = np.zeros((2, self.Mesh.NbInnerElements, 3, 2, self.Problem.NbFluids, 12))
-
- # 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]
-
+
# Retrieving the information on the local degree of freedom
dele = np.array(self.Mesh.ElementsDofs[i])
LS = LSValues[:, dele]
UU = U[:, dele]
-
+
# Get the coordinates of the vertices in the physical order
elevert = self.Mesh.CartesianPoints[ele,:]
-
- # ---- Computing the plain integral
- # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+
+ # ---- 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)
-
+
# Compute the flux according to the fluid's flag
UV = self.Problem.GoverningEquations.GetUV(Uloc, qp)
MidFluxes = self.Flux(LSloc, UV)
fluxes[0, i, 0, :, :, :] = MidFluxes
-
+
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
# (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(UU, EvaluatedBase.T)
LSloc = np.dot(LS, EvaluatedBase.T)
-
+
# Compute the flux
UV = self.Problem.GoverningEquations.GetUV(Uloc, bqp)
MidFluxes = self.Flux(LSloc, UV)
fluxes[1,i,face,:,:,0:3] = MidFluxes
-
+
# Returning the full array
return(fluxes)
#### (Crucial routine) Spatial scheme to update the level set ####
def Iteration(self, Solution, fluxes, i, du=0, dt=1):
""" Main iteration of the scheme, implementing the continuous galerkin scheme.
-
+
Args:
Solution (Solution): 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.LSValues)[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
coords = self.Mesh.DofsXY[dele,:]
Vars = Solution.Sol[:,dele]
LS = Solution.LSValues[:,dele]
dU = du[:, dele]
-
- # ---- Computing the plain integral
- # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+
+ # ---- 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)
-
+ EvaluatedGradients = BF.SimplicialGradients("B", Order, bqp, n, vol)
+
# Reconstructing the solution and its gradient 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(len(dele)):
# 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
-
+
# 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
@@ -503,16 +507,53 @@ class FS():
#### 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
-
+ """ Wrapper routine to update the quantitites that depend on the UpdateFSValues fresh output
+
Args:
Solution (Solution): the solution instance the scheme is working on
-
+
Returns:
None: updating directly the subvalues in the solution instance
-
+
--------------------------------------------------------------------------------- """
-
+
#self.UpdateSubdomains(Solution)
#self.Redistancing(Solution)
self.FlagDofs(Solution)
+
+ #### Routine that maps the values from the Dofs to the Physical vertices of the mesh ####
+ def ReconstructFSAtVertices(self, Solution):
+ """ Routine that maps the values from the Dofs to the Physical vertices of the mesh
+ of the FluidSpotter 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))
+
+ # ------------- 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)
+ 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/FluidSelectors/NaiveLS_CG_LFX.py b/SourceCode/_Solver/FluidSelectors/NaiveLS_CG_LFX.py
index d42cb7f..d55a5d2 100644
--- a/SourceCode/_Solver/FluidSelectors/NaiveLS_CG_LFX.py
+++ b/SourceCode/_Solver/FluidSelectors/NaiveLS_CG_LFX.py
@@ -23,20 +23,20 @@ import BarycentricTools as BT
# ==============================================================================
class FS():
""".. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Classical level set approach, initialised according to a signed distance and iterated
with a continuous galerkin - lfx scheme, without redistancing.
-
+
.. 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). It also contains the
+ the structure with the same names (see the parameters below). It also contains the
further field
-
+
- |bs| **FluidSelectorName** *(string)* -- the name of the level set
-
+
|
-
+
"""
#### Automatic initialisation routine ####
@@ -46,17 +46,17 @@ class FS():
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh instance
Params (string list): the level set's parameters as wished by the user
-
+
|
-
+
.. rubric:: Methods
"""
-
-
+
+
# Setting the name of the class to be accessible from outside
self.FluidSelectorName = "Naive Level Set -- Continuous Galerkin"
-
- # Setting the mesh and parameters for an internal access
+
+ # Setting the mesh and parameters for an internal access
self.Problem = Problem
self.Mesh = Mesh
self.Params = Params
@@ -68,13 +68,13 @@ class FS():
#### External initialisation routine ####
def Initialiser(self, Solution):
""" External initialisation routine.
-
+
Args:
Solution (Solution structure): the solution instance the scheme is working on
-
+
Returns:
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 ----------------------------------------------
@@ -86,7 +86,7 @@ class FS():
NbDofs = len(Solution.FluidFlag) # Only dofs are there
# Initialising the Level set values
- LSValues = np.zeros((len(Subdomains), len(Solution.FluidFlag)))
+ LSValues = np.zeros((NbFluids, NbDofs))
# ------------- Get the distance from the interfaces --------------------------------------------
# Evaluate the distance from Dofs to each boundary, and mark them by +-
@@ -101,28 +101,32 @@ class FS():
+[tuple(tuple(a) for a in np.array(interior.coords.xy).T) for interior in sub.interiors]
DomainBorder = ShapelyMultiLineString(Borders)
- # Looping on the intial number of fluids present in the computational domain
- for i in range(NbFluids):
- # Get the subdomain as a shapely multipolygon
- FluidArea = Subdomains[i]
-
- # Extract the boundary of the fluid
- Borders = []
- for sub in FluidArea:
- Borders += [ tuple(tuple(a) for a in np.array(sub.exterior.coords.xy).T)]\
- +[tuple(tuple(a) for a in np.array(interior.coords.xy).T) for interior in sub.interiors]
- FluidBorder = ShapelyMultiLineString(Borders)
-
- # Mask the parts that coincide with the domain's boundary itself
- FluidBorder = FluidBorder.difference(DomainBorder)
-
- # Computing the distance and location w.r.t. the subdomain through the shapely information
- Value = np.array([[ShapelyPoint(tuple(pt)).distance(FluidBorder),\
- 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)
+ # If we have only one fluid, set the LSValues everywhere equal to 1 for compatibility, otherwise fill the LSValues structure
+ if NbFluids==1: LSValues[:,:] = 1.
+ else:
+ # Looping on the intial number of fluids present in the computational domain
+ for i in range(NbFluids):
+
+ # Get the subdomain as a shapely multipolygon
+ FluidArea = Subdomains[i]
+
+ # Extract the boundary of the fluid
+ Borders = []
+ for sub in FluidArea:
+ Borders += [ tuple(tuple(a) for a in np.array(sub.exterior.coords.xy).T)]\
+ +[tuple(tuple(a) for a in np.array(interior.coords.xy).T) for interior in sub.interiors]
+ FluidBorder = ShapelyMultiLineString(Borders)
+
+ # Mask the parts that coincide with the domain's boundary itself
+ FluidBorder = FluidBorder.difference(DomainBorder)
+
+ # Computing the distance and location w.r.t. the subdomain through the shapely information
+ Value = np.array([[ShapelyPoint(tuple(pt)).distance(FluidBorder),\
+ 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)
@@ -130,180 +134,180 @@ class FS():
#### Gives the dynamic of the Fluid's spotter through its flux ####
def Flux(self, Var, Motion, *args):
""" Flux function of the (conservative) EulerEquations
-
- Args:
+
+ Args:
Var (float numpy array): the value of the Level Set values (NbFluids x NbGivenPoints)
- Motion (float numpy array): the velocity array (2xNbPoints)
-
+ Motion (float numpy array): the velocity array (2xNbPoints)
+
Returns:
Flux (3D numpy array): the obtained flux (spatialdim x NbVars x NbPoints)
-
- .. note::
-
+
+ .. note::
+
*args is there for compatibility reason when e.g. a location xy is given when called
-
+
--------------------------------------------------------------------------------------------------------------------------------"""
# Determining the advection point-wise
Flux = np.array([Motion[0,:]*Var, Motion[1,:]*Var])
-
+
# Returning the values
return(Flux)
-
+
#### Jacobian of the equations set ####
def Jacobian(self, Var, Motion, *args):
""" .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Computes the Jacobian of the flux for the advection of the LevelSet
-
+
Args:
Var (float numpy array): the value of the Level Set values (NbFluids x NbGivenPoints)
- Motion (float numpy array): the velocity array (2xNbPoints)
-
- Returns:
- J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
-
+ Motion (float numpy array): the velocity array (2xNbPoints)
+
+ Returns:
+ J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
+
|
-
+
.. note::
-
+
- |bs| 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]
-
+
- |bs| *args is there for compatibility reason when e.g. a location xy is given when called the x-y locations at which the variables are given (NbGivenPoints x 2)
-
+
---------------------------------------------------------------------"""
-
+
# ------ Initialising the Jacobian structure ---------------------------
# Getting the number of furnished points
NbPoints = np.shape(Motion)[1]
J = np.zeros((self.Problem.NbFluids, self.Problem.NbFluids, 2, NbPoints))
-
+
for nf in range(self.Problem.NbFluids):
J[nf, nf, 0, :] = np.array(Motion[0,:])
J[nf, nf, 1, :] = np.array(Motion[1,:])
-
- # Returning the Jacobian
- return(J)
-
+
+ # Returning the Jacobian
+ return(J)
+
#### Spectral radius of the equations set ####
def SpectralRadius(self, Var, Motion, n, *args):
""" Computes the spectral radius associated to the motion flux.
-
- Args:
+
+ Args:
Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
- Motion (float numpy array): the velocity array (2xNbPoints)
+ Motion (float numpy array): the velocity array (2xNbPoints)
n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
-
- Returns:
+
+ Returns:
Lmb (numpy array): the array containing the spectral radius at each given point
-
+
.. note::
-
+
*args is there for compatibility reason when e.g. a the xy locations at which the variables are given (NbGivenPoints x 2)
-
+
---------------------------------------------------------------------"""
-
+
# Retrieving the Jacobian
J = self.Jacobian(Var, Motion)
-
+
# 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(self.Problem.NbFluids)]), axis=0)
-
+
# Returning the obtained values
return(Lmb)
-
+
# ************************************************************************
# * Definition of the Fluid spotter's spatial impact and rectification *
# ************************************************************************
-
+
#### 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
+ """ Flags the degrees of freedom according to the different values of the
associated level set.
-
- Args:
+
+ Args:
Solution (Solution structure): the solution instance the scheme is working on
-
- Returns:
+
+ Returns:
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
+ """ Method modifying the subdomains associated to the Solution upon the values of
LSValues.
-
- Args:
+
+ Args:
Solution (Solution): the solution instance the scheme is working on
-
- Returns:
+
+ Returns:
None: updating directly the subdomains in the solution instance
-
- .. warning:: Not safe yet if the first path given out is internal
-
+
+ .. warning:: Not safe yet if the first path given out is internal
+
---------------------------------------------------------------------------- """
-
+
# ------- Updating the subdomains from the countour values -----------------
# 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):
# ---------- Initialisations --------------------------------------------
# Buffering the LSvalues to modify them on the edges upon our wishes to infer closed contours
Buffer = Solution.LSValues[:,:]
-
+
# ----------- 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
+ # 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.
-
- Args:
+
+ Args:
Solution (Solution): the solution instance the scheme is working on
-
- Returns:
+
+ Returns:
None: updating directly the LSValues in the solution instance
-
+
----------------------------------------------------------------------------------"""
# Redistancing the points with the nearest
@@ -314,205 +318,205 @@ class FS():
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)
# ****************************************************************************
# * Definition of the Fluid spotter's evolution's process (e.g. advection) *
# ****************************************************************************
-
+
#### Emulates a flux computation as done in the iteration routine ####
def ComputeFlux(self, U, LSValues):
""" Emulates a flux computation as done in the iteration routine in order
to be used in the DeC subtimestep evaluation.
-
- Args:
+
+ Args:
U (numpy float array): buffer containing the current state values at each dof (NbVars x NbDofs)
LS (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 requires by the UpdateFSValues 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
fluxes = np.zeros((2, self.Mesh.NbInnerElements, 3, 2, self.Problem.NbFluids, 12))
-
- # 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]
-
+
# Retrieving the information on the local degree of freedom
dele = np.array(self.Mesh.ElementsDofs[i])
LS = LSValues[:, dele]
UU = U[:, dele]
-
+
# Get the coordinates of the vertices in the physical order
elevert = self.Mesh.CartesianPoints[ele,:]
-
- # ---- Computing the plain integral
- # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+
+ # ---- 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)
-
+
# Compute the flux according to the fluid's flag
UV = self.Problem.GoverningEquations.GetUV(Uloc, qp)
MidFluxes = self.Flux(LSloc, UV)
fluxes[0, i, 0, :, :, :] = MidFluxes
-
+
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
# (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(UU, EvaluatedBase.T)
LSloc = np.dot(LS, EvaluatedBase.T)
-
+
# Compute the flux
UV = self.Problem.GoverningEquations.GetUV(Uloc, bqp)
MidFluxes = self.Flux(LSloc, UV)
fluxes[1,i,face,:,:,0:3] = MidFluxes
-
+
# Returning the full array
return(fluxes)
#### (Crucial routine) Spatial scheme to update the level set ####
def Iteration(self, Solution, fluxes, i, du=0, dt=1):
""" Main iteration of the scheme, implementing the continuous galerkin scheme.
-
+
Args:
Solution (Solution): 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.LSValues)[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
coords = self.Mesh.DofsXY[dele,:]
Vars = Solution.Sol[:,dele]
LS = Solution.LSValues[:,dele]
dU = du[:, dele]
-
- # ---- Computing the plain integral
- # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+
+ # ---- 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)
-
+ EvaluatedGradients = BF.SimplicialGradients("B", Order, bqp, n, vol)
+
# Reconstructing the solution and its gradient 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(len(dele)):
# 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
-
+
# 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 ------------------------
# Select the barycenter of the element as evaluation point
Bcoords = self.Mesh.DualCenters[i][-1]
-
+
# Getting the fluid's velocity quantities
UV = self.Problem.GoverningEquations.GetUV(Vars, coords)
-
+
# Computing the average state on the spanned element
LSb = np.array([np.mean(LS, axis=1)]).T
VAb = np.array([np.mean(UV, axis=1)]).T
-
+
# Computing the spectral radius at each vertex and extracting the maximum
# according to the average speed value
sprctrad = self.SpectralRadius(LSb, VAb, n)
-
+
# Modifying the residuals with the LFX flux
Resu[:,:] += dt*np.max(sprctrad)*(LS[:,:]-LSb[:])
@@ -521,16 +525,53 @@ class FS():
#### 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
-
+ """ Wrapper routine to update the quantitites that depend on the UpdateFSValues fresh output
+
Args:
Solution (Solution): the solution instance the scheme is working on
-
+
Returns:
None: updating directly the subvalues in the solution instance
-
+
--------------------------------------------------------------------------------- """
-
+
#self.UpdateSubdomains(Solution)
#self.Redistancing(Solution)
self.FlagDofs(Solution)
+
+ #### Routine that maps the values from the Dofs to the Physical vertices of the mesh ####
+ def ReconstructFSAtVertices(self, Solution):
+ """ Routine that maps the values from the Dofs to the Physical vertices of the mesh
+ of the FluidSpotter 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))
+
+ # ------------- 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)
+ 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/FluidSelectors/NaiveLS_CG_LFX_Limited.py b/SourceCode/_Solver/FluidSelectors/NaiveLS_CG_LFX_Limited.py
new file mode 100644
index 0000000..dccc184
--- /dev/null
+++ b/SourceCode/_Solver/FluidSelectors/NaiveLS_CG_LFX_Limited.py
@@ -0,0 +1,614 @@
+#"""
+#Module implementing a naive level set without redistancing
+#"""
+
+# ==============================================================================
+# Importing python libraries
+# ==============================================================================
+from shapely.geometry import MultiLineString as ShapelyMultiLineString
+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 numpy as np
+
+# Import custom modules
+import _Solver.QuadratureRules.HandBased as QR
+import BarycentricBasisFunctions as BF
+import BarycentricTools as BT
+
+# ==============================================================================
+# Defining the actuall class gathering the externally accessible functions
+# ==============================================================================
+class FS():
+ """.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
+
+ Classical level set approach, initialised according to a signed distance and iterated
+ with a continuous galerkin - lfx scheme, without redistancing.
+
+ .. 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). It also contains the
+ further field
+
+ - |bs| **FluidSelectorName** *(string)* -- the name of the level set
+
+ |
+
+ """
+
+ #### 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 level set's parameters as wished by the user
+
+ |
+
+ .. rubric:: Methods
+ """
+
+
+ # Setting the name of the class to be accessible from outside
+ self.FluidSelectorName = "Naive Level Set -- Continuous Galerkin"
+
+ # Setting the mesh and parameters for an internal access
+ self.Problem = Problem
+ self.Mesh = Mesh
+ self.Params = Params
+
+ # **********************************************************
+ # * Definition of the Fluid spotter's type and dynamic *
+ # **********************************************************
+
+ #### External initialisation routine ####
+ def Initialiser(self, Solution):
+ """ External initialisation routine.
+
+ Args:
+ Solution (Solution structure): the solution instance the scheme is working on
+
+ Returns:
+ 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((NbFluids, NbDofs))
+
+ # ------------- 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.)
+
+ # Extract the boundary of the domain
+ Borders = []
+ for sub in self.Problem.StudyDomain:
+ Borders += [tuple(tuple(a) for a in np.array(sub.exterior.coords.xy).T)]\
+ +[tuple(tuple(a) for a in np.array(interior.coords.xy).T) for interior in sub.interiors]
+ DomainBorder = ShapelyMultiLineString(Borders)
+
+ # If we have only one fluid, set the LSValues everywhere equal to 1 for compatibility, otherwise fill the LSValues structure
+ if NbFluids==1: LSValues[:,:] = 1.
+ else:
+ # Looping on the intial number of fluids present in the computational domain
+ for i in range(NbFluids):
+
+ # Get the subdomain as a shapely multipolygon
+ FluidArea = Subdomains[i]
+
+ # Extract the boundary of the fluid
+ Borders = []
+ for sub in FluidArea:
+ Borders += [ tuple(tuple(a) for a in np.array(sub.exterior.coords.xy).T)]\
+ +[tuple(tuple(a) for a in np.array(interior.coords.xy).T) for interior in sub.interiors]
+ FluidBorder = ShapelyMultiLineString(Borders)
+
+ # Mask the parts that coincide with the domain's boundary itself
+ FluidBorder = FluidBorder.difference(DomainBorder)
+
+ # Computing the distance and location w.r.t. the subdomain through the shapely information
+ Value = np.array([[ShapelyPoint(tuple(pt)).distance(FluidBorder),\
+ 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)
+
+ #### Gives the dynamic of the Fluid's spotter through its flux ####
+ def Flux(self, Var, Motion, *args):
+ """ Flux function of the (conservative) EulerEquations
+
+ Args:
+ Var (float numpy array): the value of the Level Set values (NbFluids x NbGivenPoints)
+ Motion (float numpy array): the velocity array (2xNbPoints)
+
+ Returns:
+ Flux (3D numpy array): the obtained flux (spatialdim x NbVars x NbPoints)
+
+ .. note::
+
+ *args is there for compatibility reason when e.g. a location xy is given when called
+
+ --------------------------------------------------------------------------------------------------------------------------------"""
+
+ # Determining the advection point-wise
+ Flux = np.array([Motion[0,:]*Var, Motion[1,:]*Var])
+
+ # Returning the values
+ return(Flux)
+
+ #### Jacobian of the equations set ####
+ def Jacobian(self, Var, Motion, *args):
+ """ .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
+
+ Computes the Jacobian of the flux for the advection of the LevelSet
+
+ Args:
+ Var (float numpy array): the value of the Level Set values (NbFluids x NbGivenPoints)
+ Motion (float numpy array): the velocity array (2xNbPoints)
+
+ Returns:
+ J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
+
+ |
+
+ .. note::
+
+ - |bs| 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]
+
+ - |bs| *args is there for compatibility reason when e.g. a location xy is given when called the x-y locations at which the variables are given (NbGivenPoints x 2)
+
+ ---------------------------------------------------------------------"""
+
+ # ------ Initialising the Jacobian structure ---------------------------
+ # Getting the number of furnished points
+ NbPoints = np.shape(Motion)[1]
+ J = np.zeros((self.Problem.NbFluids, self.Problem.NbFluids, 2, NbPoints))
+
+ for nf in range(self.Problem.NbFluids):
+ J[nf, nf, 0, :] = np.array(Motion[0,:])
+ J[nf, nf, 1, :] = np.array(Motion[1,:])
+
+ # Returning the Jacobian
+ return(J)
+
+ #### Spectral radius of the equations set ####
+ def SpectralRadius(self, Var, Motion, n, *args):
+ """ Computes the spectral radius associated to the motion flux.
+
+ Args:
+ Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
+ Motion (float numpy array): the velocity array (2xNbPoints)
+ n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
+
+ Returns:
+ Lmb (numpy array): the array containing the spectral radius at each given point
+
+ .. note::
+
+ *args is there for compatibility reason when e.g. a the xy locations at which the variables are given (NbGivenPoints x 2)
+
+ ---------------------------------------------------------------------"""
+
+ # Retrieving the Jacobian
+ J = self.Jacobian(Var, Motion)
+
+ # 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(self.Problem.NbFluids)]), axis=0)
+
+ # Returning the obtained values
+ return(Lmb)
+
+ # ************************************************************************
+ # * Definition of the Fluid spotter's spatial impact and rectification *
+ # ************************************************************************
+
+ #### 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.
+
+ Args:
+ Solution (Solution structure): the solution instance the scheme is working on
+
+ Returns:
+ 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.
+
+ Args:
+ Solution (Solution): the solution instance the scheme is working on
+
+ Returns:
+ 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 -----------------
+ # 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):
+ # ---------- Initialisations --------------------------------------------
+ # Buffering the LSvalues to modify them on the edges upon our wishes to infer closed contours
+ Buffer = Solution.LSValues[:,:]
+
+ # ----------- 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.
+
+ Args:
+ Solution (Solution): the solution instance the scheme is working on
+
+ Returns:
+ 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)
+
+ # ****************************************************************************
+ # * Definition of the Fluid spotter's evolution's process (e.g. advection) *
+ # ****************************************************************************
+
+ #### Emulates a flux computation as done in the iteration routine ####
+ def ComputeFlux(self, U, LSValues):
+ """ Emulates a flux computation as done in the iteration routine in order
+ to be used in the DeC subtimestep evaluation.
+
+ Args:
+ U (numpy float array): buffer containing the current state values at each dof (NbVars x NbDofs)
+ LS (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 requires by the UpdateFSValues 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
+ fluxes = np.zeros((2, self.Mesh.NbInnerElements, 3, 2, self.Problem.NbFluids, 12))
+
+ # 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]
+
+ # Retrieving the information on the local degree of freedom
+ dele = np.array(self.Mesh.ElementsDofs[i])
+ LS = LSValues[:, dele]
+ UU = U[:, dele]
+
+ # Get the coordinates of the vertices in the physical order
+ elevert = self.Mesh.CartesianPoints[ele,:]
+
+ # ---- 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)
+
+ # Compute the flux according to the fluid's flag
+ UV = self.Problem.GoverningEquations.GetUV(Uloc, qp)
+ MidFluxes = self.Flux(LSloc, UV)
+ fluxes[0, i, 0, :, :, :] = MidFluxes
+
+ 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
+ # (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
+
+ # Reconstruct the solution's value at this point
+ Uloc = np.dot(UU, EvaluatedBase.T)
+ LSloc = np.dot(LS, EvaluatedBase.T)
+
+ # Compute the flux
+ UV = self.Problem.GoverningEquations.GetUV(Uloc, bqp)
+ MidFluxes = self.Flux(LSloc, UV)
+ fluxes[1,i,face,:,:,0:3] = MidFluxes
+
+ # Returning the full array
+ return(fluxes)
+
+ #### (Crucial routine) Spatial scheme to update the level set ####
+ def Iteration(self, Solution, fluxes, i, du=0, dt=1):
+ """ Main iteration of the scheme, implementing the continuous galerkin scheme.
+
+ Args:
+ Solution (Solution): 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.LSValues)[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
+ coords = self.Mesh.DofsXY[dele,:]
+ Vars = Solution.Sol[:,dele]
+ LS = Solution.LSValues[:,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
+ 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(len(dele)):
+ # 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
+
+ # ------- Editing the CG residuals to become LFX ------------------------
+ # Select the barycenter of the element as evaluation point
+ Bcoords = self.Mesh.DualCenters[i][-1]
+
+ # Getting the fluid's velocity quantities
+ UV = self.Problem.GoverningEquations.GetUV(Vars, coords)
+
+ # Computing the average state on the spanned element
+ LSb = np.array([np.mean(LS, axis=1)]).T
+ VAb = np.array([np.mean(UV, axis=1)]).T
+
+ # Computing the spectral radius at each vertex and extracting the maximum
+ # according to the average speed value
+ sprctrad = self.SpectralRadius(LSb, VAb, n)
+
+ # Modifying the residuals with the LFX flux
+ Resu[:,:] += dt*np.max(sprctrad)*(LS[:,:]-LSb[:])
+
+ # ------ Editing the residuals to become limited -----------------------
+ # --- Performing the limiting initialisation
+ # Getting the number dofs contaned in the considered element
+ NbLS = np.shape(Solution.LSValues)[0]
+
+ # Initialising the mollification vector (residuals)
+ Lresu = np.zeros(np.shape(Resu))
+
+ # Setting the safety division trick
+ safe = 1e-42
+
+ # --- Limiting the given residual within the given element
+ # Checking if it is actually relevant and possible to limit, and if yes, select
+ # the variables on which to limit
+ Phi = np.array([np.sum(Resu, axis = 1)]).T
+ ModInd = np.where(Phi>safe)[0]
+ NodInd = np.setdiff1d(range(NbLS), 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)
+
+ #### 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
+
+ Args:
+ Solution (Solution): the solution instance the scheme is working on
+
+ Returns:
+ None: updating directly the subvalues in the solution instance
+
+ --------------------------------------------------------------------------------- """
+
+ #self.UpdateSubdomains(Solution)
+ #self.Redistancing(Solution)
+ self.FlagDofs(Solution)
+
+ #### Routine that maps the values from the Dofs to the Physical vertices of the mesh ####
+ def ReconstructFSAtVertices(self, Solution):
+ """ Routine that maps the values from the Dofs to the Physical vertices of the mesh
+ of the FluidSpotter 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))
+
+ # ------------- 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)
+ 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/Solver_Initialisations.py b/SourceCode/_Solver/Solver_Initialisations.py
index d5bcdaf..5610a07 100644
--- a/SourceCode/_Solver/Solver_Initialisations.py
+++ b/SourceCode/_Solver/Solver_Initialisations.py
@@ -4,8 +4,8 @@
The structures and classes that allow a shortcut of all the routines
used upon the user's choices (governing equations, spatial schemes, etc)
are defined in the module :code:`_Solver.Solver_Initialisations`.
-After shortuting the routine accordingly to the values registered in the file
-:code:`Mappers.py`, it automatically defines the spatial subdomains,
+After shortuting the routine accordingly to the values registered in the file
+:code:`Mappers.py`, it automatically defines the spatial subdomains,
initialises the solution and the fluid's selector.
The fields and methods of the three main instances representing the Problem, Solver and Solution
are given below.
@@ -19,7 +19,7 @@ are given below.
import ProblemDefinition.FluidSubdomains as SD
import _Solver.SpatialModifiers.SpatialAmendements as SA
import Mappers
-import Pathes
+import Pathes
# Import python and math libraries
import os, sys, re
@@ -30,20 +30,20 @@ import numpy as np
# Parameter (and reader) class of scheme properties and runtime informations
# ==============================================================================
class Parameters():
- """
+ """
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Parameter (and reader) class of scheme properties and runtime informations.
It furnishes all the information given from the user through the
setting file given in argument when running the script in a structure
containing the following fields:
-
+
|
-
+
.. rubric:: Parameters chosen by the user
-
+
The text content is exported into a ParamsReader structure variable containing:
-
+
- |bs| **ExportInterval** *(float)* -- time interval after which the computed results will be exported
- |bs| **Tmaw** *(float)* -- the (physical) end-time of the computations
- |bs| **CFL** *(float)* -- the CFl number
@@ -57,24 +57,24 @@ class Parameters():
- |bs| **QuadratureType** *(integer)* -- index of the wished quadrature (see the Mapper file)
- |bs| **QuadratureOrder** *(integer)* -- order of the wished quadrature
- |bs| **Verbose** *(string)* -- verbose mode: , "d": debug, "m": minimal, "n": none
-
+
|
-
+
.. rubric:: Automatically generated attributes:
-
+
*Retrieved directly at instance creation*
-
+
- |bs| **MeshType** *(string)* -- type of the mesh associated (CG, DG, RT, BDM, ...)
- |bs| **MeshOrder** *(integer)* -- order of the mesh
|
-
- .. note::
-
+
+ .. note::
+
This is only a product-type class: no method is available
-
+
|
-
+
--------------------------------------------------------------"""
#### Initialisation (reader) routine #####
@@ -82,7 +82,7 @@ class Parameters():
"""
Args:
SolverSpecs (string): the name of the setting file
-
+
"""
# Mapping to the file where the SolverData are stored and get all of its content
@@ -114,7 +114,7 @@ class Parameters():
self.QuadratureOrder = int(re.split(" |\t", Params[9])[0])
# Retrieve the runtime options
- self.ExportInterval = int(re.split(" |\t", Params[19])[0])
+ self.ExportInterval = float(re.split(" |\t", Params[19])[0])
self.Verbose = re.split(" |\t", Params[18])[0]
# Retrieve the type of required Mesh
@@ -128,7 +128,7 @@ class Parameters():
class Solution():
"""
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Class used as an product type object containing all the informations on
the computed solution, containing the following attributes:
@@ -146,43 +146,43 @@ class Solution():
- |bs| **t** *(float)* -- The time corresponding to the solution
|
-
+
Upon instance creation, itnitialises the solution structure upon the problem and the mesh (except LSValues).
-
+
|
.. note:
-
+
The class is not inherited from the mesh/problem in order to be able to keep the
solution as a subclass to be transported easily
-
+
|
-
+
"""
#### Actual initialisation routine upon the problem's parameters ####
def __init__(self, Problem, Mesh):
- """
-
-
+ """
+
+
Args:
Problem (Problem): structure containing all the problem info and shortcuted routines
Mesh (MeshStructure): the considered mesh
-
+
|
-
+
.. rubric:: Methods
-
+
"""
# Initialisation of the CFL related fix properties
self.t = 0
self.maxh = np.max(Mesh.InnerElementsDiameter)
self.minh = np.min(Mesh.InnerElementsDiameter)
-
+
# Initialisation of the Solution's fixed properties
self.Sol = np.zeros((Problem.GoverningEquations.NbVariables, Mesh.NbDofs))
self.FluidFlag = np.zeros(Mesh.NbDofs, dtype=int)
- self.LSValues = [] # Initialized to void.
+ self.LSValues = [] # Initialized to void.
# Will be initialized when the LS itlself will be set up.
# Initialisation of the Reconstructed solution's fixed properties
@@ -190,7 +190,7 @@ class Solution():
self.RXY = Mesh.CartesianPoints[:,:]
self.RFluidFlag = np.zeros(Mesh.NbMeshPoints, dtype=int)
self.RLSValues = [] # Initialized to void as for LSValues.
-
+
# Initialisation of the fluid-dependent parameters
self.Initialise_Flags(Problem, Mesh)
self.Initialise_SolutionValues(Problem, Mesh)
@@ -199,14 +199,14 @@ class Solution():
def Initialise_Flags(self, Problem, Mesh):
"""
Initialise the solution structure upon the problem and the mesh.
-
- Args:
+
+ Args:
Problem (Problem): structure containing all the problem info and shortcuted routines
Mesh (MeshStructure): the considered mesh
-
+
Returns:
None: filled FluidFlags vector within solution structure.
-
+
------------------------------------------------------------------------------------------------"""
# Exports the Flags
@@ -215,23 +215,23 @@ class Solution():
#### Initilising the actual values of the solutions ####
def Initialise_SolutionValues(self, Problem, Mesh):
- """
+ """
Initialise the solution structure upon the problem and the mesh.
-
- Args:
+
+ Args:
Problem (Problem): structure containing all the problem info and shortcuted routines
Mesh (MeshStructure): the considered mesh
-
- Returns:
+
+ Returns:
None: filled intial solution values within solution structure.
-
+
------------------------------------------------------------------------------------------------"""
# Initialise the solution by extracting the asked initial condition for each fluid aera
for i in range(Problem.NbFluids):
# Getting the index of the considered points
Index = np.where(self.FluidFlag[:]==i)[0]
-
+
# Getting the fluid's properties
EOS = Problem.FluidEOS[i]
FP = Problem.FluidProp[i]
@@ -248,14 +248,14 @@ class Solution():
# ==============================================================================
#### Defining the problem itself and the shortut to the relevant functions #####
class Problem():
- """
+ """
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Shortcutting class gathering all the methods and functions that the user targets
to solve the problem. It sets the problem up upon the solver parameters and the mesh,
and furnishes the following fields:
- - |bs| **NbFluids** *(integer)* -- number of fluids present initially in the problem's domain
+ - |bs| **NbFluids** *(integer)* -- number of fluids present initially in the problem's domain
- |bs| **FluidIndex** *(integer list)* -- the ordered list of fluid types corresponding to the fluids present initally in the problem's domain
- |bs| **InitialSubdomains** *(shapely multipolygon list)* -- list of the initial fluids subdomains, in the given order of the fluid lists
- |bs| **GoverningEquations** *(Equations class)* -- structure containing the methods of flux, jacobian, spectral radius, and a possible Roe matrix
@@ -263,21 +263,21 @@ class Problem():
- |bs| **FluidProp** *(Fluids structure list)* -- list of all the fluids properties correspoding to the list of fluids themselves
- |bs| **FluidEOS** *(function callback list)* -- list of all the equations of state in use for each fluid
- |bs| **EOSType** *(string)* -- the name of the EOS in use
-
+
|
------------------------------------------------------------------------------------------------------------ """
#### Automatic initialisation routine ####
def __init__(self, ProblemData, Mesh):
""" Automatic initialisation routine.
-
+
Args:
ProblemData (Parameters): the structure containing all the problem data
Mesh (MeshStructure): the mesh in use
-
- Returns:
+
+ Returns:
None: fills direclty the Problem structure
-
+
----------------------------------------------------------------------------------------"""
# Getting the fluids information from the ProblemData
@@ -291,7 +291,7 @@ class Problem():
# Setting the fluids properties and their equations of state
self.FluidProp = [Mappers.Fluid_Properties(i, FluidModel) for i in ProblemData.FluidIndex]
self.FluidEOS = [Setup.EOS(i) for i in ProblemData.EOSEquationsIndex]
-
+
# Getting the GoverningEquations themselves and retrieving the associated initial conditions and boundary conditions
self.GoverningEquations = Setup.Equations(self.FluidProp, self.FluidEOS)
self.InitialConditions = [Setup.InitialConditions(*id).IC for id in ProblemData.InitialConditions]
@@ -299,7 +299,7 @@ class Problem():
# Registering the Boundary conditions routine per boundary tag and per fluid
self.BoundaryConditions = np.zeros(ProblemData.BoundaryConditions.shape).tolist()
for i in range(ProblemData.BoundaryConditions.shape[0]):
- for j in range(ProblemData.BoundaryConditions.shape[1]):
+ for j in range(ProblemData.BoundaryConditions.shape[1]):
self.BoundaryConditions[i][j] = Mappers.Boundary_Conditions(*ProblemData.BoundaryConditions[i][j])
# Setting the initial subdomains where each of the initial fluid lie
@@ -307,24 +307,24 @@ class Problem():
#### Defining the solver class shortcuting all the routines that should be used #####
class Solver():
- """
+ """
.. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Shortcutting class gathering all the methods and functions that the solver
requires, upon the user desires and the inherent solver's properties,
furnishing the following fields:
- - |bs| **Quadrature** *(function callback)* -- quadrature rule
- - |bs| **FluidSelector** *(class instance)* -- FluidSelector instance containing all the required Fluid selector routines
- - |bs| **SpatialScheme** *(class instance)* -- Scheme instance, containing the mail "Iteration" method that corresponds to the spatial solver
- - |bs| **SPAdmendements** *(class instance)* -- Amendements instance, containing all the required routines to limit the scheme (if wished by the user)
- - |bs| **ControlVolumes** *(string)* -- The type of control volume required by the scheme: either "primal" or "dual"
- - |bs| **Reconstructor** *(function callback)* -- Function that reconstructs the solution at physical vertices from the degrees of freedom
- - |bs| **TimeStep** *(function callback)* -- Function implementing one time step iteration
- - |bs| **CFL** *(float)* -- A recall of the CFL number
- - |bs| **Tmax** *(float)* -- A recall of the final time
- - |bs| **t** *(float)* -- the current time
-
+ - |bs| **Quadrature** *(function callback)* -- quadrature rule
+ - |bs| **FluidSelector** *(class instance)* -- FluidSelector instance containing all the required Fluid selector routines
+ - |bs| **SpatialScheme** *(class instance)* -- Scheme instance, containing the mail "Iteration" method that corresponds to the spatial solver
+ - |bs| **SPAdmendements** *(class instance)* -- Amendements instance, containing all the required routines to limit the scheme (if wished by the user)
+ - |bs| **ControlVolumes** *(string)* -- The type of control volume required by the scheme: either "primal" or "dual"
+ - |bs| **Reconstructor** *(function callback list)* -- Functions that reconstructs the solution and Fluid selector values at physical vertices from the degrees of freedom
+ - |bs| **TimeStep** *(function callback)* -- Function implementing one time step iteration
+ - |bs| **CFL** *(float)* -- A recall of the CFL number
+ - |bs| **Tmax** *(float)* -- A recall of the final time
+ - |bs| **t** *(float)* -- the current time
+
|
"""
@@ -333,7 +333,7 @@ class Solver():
"""
Args:
SolverParams (Parameters): the structure containing all the problem data
-
+
"""
# Set up the quadrature type
@@ -347,8 +347,8 @@ class Solver():
self.SPAdmendements = SA.Amendements(SolverParams.ScModifiers, Problem, Mesh)
# Map the solution's reconstructor to the structure
- self.Reconstructor = self.SpatialScheme.ReconstructSolutionAtVertices
-
+ self.Reconstructor = [self.SpatialScheme.ReconstructSolutionAtVertices, self.FluidSelector.ReconstructFSAtVertices]
+
# Set up the time scheme and relocate the time parameters variables
self.TimeStep = Mappers.Temporal_Scheme(SolverParams.TimeSchemeId, SolverParams.TimerParams)
self.CFL = SolverParams.CFL
diff --git a/SourceCode/_Solver/Solver_Spatial.py b/SourceCode/_Solver/Solver_Spatial.py
index 9554615..35c78b2 100644
--- a/SourceCode/_Solver/Solver_Spatial.py
+++ b/SourceCode/_Solver/Solver_Spatial.py
@@ -7,7 +7,7 @@
# ==============================================================================
# Pythonic import
import numpy as np
-import copy
+import copy
# Custom import (has to be precisely imported, BCs are not included in the __all__)
import _Solver.Solver_ApplyBCs as BCs
@@ -20,13 +20,13 @@ import _Solver.Solver_ApplyBCs as BCs
#### Integrate the RD formulation in the Dec setting ####
def RD(Theta, m, M, dt, schemeSP, schemeFS, Problem, Mesh, Solution, up, lp):
""" Integrates the RD formulation in the Dec setting: interfaces from the scheme to the time stepper when used within a residual distribution framework
-
+
Args:
Theta (float numpy array): coefficients of the DeC routine
m (integer): current subtimestep of the DeC correction
M (integer): total number of DeC corrections
dt (float): the total current time step duration
-
+
schemeSP (scheme module list): handler of the spatial scheme module, followed (if any) by the list of its amendments (limiter, jumps etc)
schemeFS (scheme module list): handler of the FluidSpotter scheme module followed (if any) by the list of its amendments (limiter, jumps etc)
Problem (Problem): considered problem
@@ -37,37 +37,37 @@ def RD(Theta, m, M, dt, schemeSP, schemeFS, Problem, Mesh, Solution, up, lp):
Returns:
Residuals (3D numpy array): being of the form [upres (NbVarsxNbDofs), lsres(NbFluidsxNbDofs)], each term corresponds to the residuals of the state and the fluid's spotters, respectively.
-
+
"""
-
+
# --------------------- Initialisations -----------------------------------------------
# Intialising the value which will be considered at this time step
SolBuff = copy.deepcopy(Solution)
-
+
# Initialising the total residual vectors
upres = np.zeros(np.shape(Solution.Sol))
lsres = np.zeros(np.shape(Solution.LSValues))
-
+
# ------------------ Incrementing according to the DeC approximation ------------------
# Computing the temporal increment coming from dec
dup = up[:,:,m]-up[:,:,0]
dlp = lp[:,:,m]-lp[:,:,0]
-
-
+
+
# Computing the temporally local flux and the infinitesimal update to the solution
fincr = np.sum(np.array([Theta[r,m]*schemeSP[0].ComputeFlux(Solution.FluidFlag, up[:,:,r], lp[:,:,r]) for r in range(M+1)]), axis = 0)
sincr = np.sum(np.array([Theta[r,m]*up[:,:,r] for r in range(M+1)]), axis = 0)
-
+
# Computing the temporally local flux and the correction to the function
lincr = np.sum(np.array([Theta[r,m]*schemeFS[0].ComputeFlux(up[:,:,r], lp[:,:,r]) for r in range(M+1)]), axis = 0)
vincr = np.sum(np.array([Theta[r,m]*lp[:,:,r] for r in range(M+1)]), axis = 0)
-
+
# ------------------ Applying the dynamic to the iterated values ------------------
# Update fluid flag upon the new value of the level sets and filling a buffer structure
SolBuff.LSValues = vincr
SolBuff.Sol = sincr
schemeFS[0].FlagDofs(SolBuff)
-
+
# Applying the LS and Dynamic schemes upon the time-locally iterated solution, element-wise
# (the split on the elements already at this level is required by the possible mollification
# of the residual by limiters/streamline stabilisation)
@@ -75,26 +75,25 @@ def RD(Theta, m, M, dt, schemeSP, schemeFS, Problem, Mesh, Solution, up, lp):
# Computing the initial residual with the given scheme
resu = schemeSP[0].Iteration(SolBuff, fincr, i, dup, dt)[:]
resl = schemeFS[0].Iteration(SolBuff, lincr, i, dlp, dt)[:]
-
+
# Mollifying the obtained partial residuals by a the wished techniques if applicable
resu = schemeSP[1].Mollify(SolBuff, resu, fincr, i, dup, dt)
-
+
# Mapping the partial residuals to the total one
dele = np.array(Mesh.ElementsDofs[i])
upres[:, dele] += resu
lsres[:, dele] += resl
-
+
# ------------------ Applying the boundary conditions over the freshly computed residuals ------------------
upres = BCs.ApplyBCSState(Problem, Mesh, SolBuff, upres)
-
# ------------------ Applying the dynamic to the iterated values ------------------
# Updating each Dof (here representing the average) by the residual
# pondered (and lumped) by the control volume's area
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/Solver_Temporal.py b/SourceCode/_Solver/Solver_Temporal.py
index aadd84d..75b84e9 100644
--- a/SourceCode/_Solver/Solver_Temporal.py
+++ b/SourceCode/_Solver/Solver_Temporal.py
@@ -15,73 +15,73 @@ import numpy as np
#### Computing the mass lumping coefficients ####
def GetLumpingCoefficients(Mesh, Solver):
- """ Computes the mass lumping coefficients.
-
+ """ Computes the mass lumping coefficients.
+
Args:
Mesh (MeshStructure): the instance containing the mesh information
Solver (Solver): instance containing all the solver's parameters
-
+
Returns:
- None: fills directly the "Lumping" field of the mesh structure
-
+ None: fills directly the "Lumping" field of the mesh structure
+
---------------------------------------------"""
-
+
# ------------ Initialisations ------------------------------------
- # Getting the type of the control volume
+ # Getting the type of the control volume
typ = Solver.ControlVolumes
- # Initialising the
+ # Initialising the
Mesh.Lumping = np.zeros(Mesh.NbDofs)
-
+
# ------------ Actual mass lumping computations -------------------
# In the case of a dual mesh based control volume, compute the mass lumping coefficients
# by agglomerating the pondered volume of each cell (which is the information contained in the DualAreasTot)
if typ == "Dual":
for i in range(Mesh.NbInnerElements):
- Mesh.Lumping[Mesh.ElementsBoundaryDofs[i]] += Mesh.InnerElementsVolume[i]/len(Mesh.ElementsDofs[i])
+ Mesh.Lumping[Mesh.ElementsDofs[i]] += Mesh.InnerElementsVolume[i]/len(Mesh.ElementsDofs[i]) #Mesh.Lumping[Mesh.ElementsBoundaryDofs[i]]
Mesh.Lumping = 1/Mesh.Lumping
-
+
# In the case of a primal mesh based control volume, (warning: only cell-centred FV in the RD setting)
# register the elment's volume as the control one
elif typ == "Primal":
for i in range(Mesh.NbInnerElements):
Mesh.Lumping[Mesh.ElementsInnerDofs[i]] += Mesh.InnerElementsVolume[i]/len(Mesh.ElementsDofs[i])
Mesh.Lumping = 1/Mesh.Lumping
-
+
# Safe exception case
else:
raise ValueError("Error: Unknown type of control volume. Check you scheme definition. Abortion.")
#### Computing the time step upon the solution and the CFL ####
def GetTimeStep(CFL, Problem, Mesh, Solution):
- """ Computes the time step for each iteration upon the solution and the CFL.
-
+ """ Computes the time step for each iteration upon the solution and the CFL.
+
Args:
CFL (float): the CFL number
- Problem (Problem structure): the instance containing the problem information
+ Problem (Problem structure): the instance containing the problem information
Mesh (MeshStructure): the instance containing the mesh information
Solution (Solution): the solution of the problem that is being evolved in time
-
- Returns:
- timestep (float): the time step
-
+
+ Returns:
+ timestep (float): the time step
+
---------------------------------------------"""
-
+
# Preparing a huge increment, that is going to be reduced
delta = 1
-
+
# Looping on each element to find the smallest spectral radius
for i in range(Mesh.NbInnerElements):
# Extracting the physical element to get the number of edges
ele = Mesh.InnerElements[i]
dele = Mesh.ElementsInnerDofs[i]
-
+
# Looping on each edge
for f in range(len(ele)):
# Getting the variational information
dele += Mesh.ElementsBoundaryWiseDofs[i][f]
Flags = Solution.FluidFlag[dele]
Vars = Solution.Sol[:,dele]
-
+
# Getting the physical information
Coords = Mesh.DofsXY[dele, :]
nn = np.array([Mesh.ElementsBoundaryWiseNormals[i][f]])
diff --git a/SourceCode/_Solver/SpatialModifiers/SpatialAmendements.py b/SourceCode/_Solver/SpatialModifiers/SpatialAmendements.py
index f43fad8..11ea296 100644
--- a/SourceCode/_Solver/SpatialModifiers/SpatialAmendements.py
+++ b/SourceCode/_Solver/SpatialModifiers/SpatialAmendements.py
@@ -1,13 +1,13 @@
"""
-The module :code:`SpatialAmendements` located in the file :code:`_Solver/SpatialModifiers/SpatialAmendements.py`furnishes a single class
-:code:`Amendements` that interfaces the user instructions for limiting and amendments of the schemes to the actuals amemdement routine,
+The module :code:`SpatialAmendements` located in the file :code:`_Solver/SpatialModifiers/SpatialAmendements.py`furnishes a single class
+:code:`Amendements` that interfaces the user instructions for limiting and amendments of the schemes to the actuals amemdement routine,
respecting the sequence and order of mollifications.
-Upon instance creation, it converts the string-user input defining the sequence of modifications to a unique equivalent modification
-routine that is to apply on any freshly computed residual.
-To do so, it converts the string to a sequence of unitary application process and combines them according
+Upon instance creation, it converts the string-user input defining the sequence of modifications to a unique equivalent modification
+routine that is to apply on any freshly computed residual.
+To do so, it converts the string to a sequence of unitary application process and combines them according
to the given syntax. The method :code:`Mollify` is therefore the generic interface to apply any sequence of mollifyer to a spatial scheme.
|
@@ -35,7 +35,7 @@ class Amendements():
.. rubric:: Fields
- Upon instance creation, the given parameters are repeated within the structure
+ Upon instance creation, the given parameters are repeated within the structure
with identical names (see below). It also has for further fields:
- |bs| **MolliProcess** *(list of callbacks)* -- the list of amendements to apply to the designated sequential intermediate residuals
@@ -45,84 +45,85 @@ class Amendements():
#### Automatic initialisation routine ####
def __init__(self, AmendementsSequence, Problem, Mesh, Params = []):
- """
+ """
Args:
AmendementsSequence (string): the user-understandable list of amendements to apply to the designated sequential intermediate residuals
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.AmendementsSequence = AmendementsSequence
self.Problem = Problem
self.Params = Params
self.Mesh = Mesh
-
+
# Extracting the modifications sequence as tabular of unitary modifications
# (and mapping to the wished modifiers according to the mapper class)
# if any modifier is wished
if "0" != self.AmendementsSequence:
- self.StringToResu(self.ExtractSequence()) # Filling self.MolliProcess
+ # Filling self.MolliProcess
+ self.StringToResu(self.ExtractSequence())
else:
self.MolliProcess = 0
-
+
#### Converting the string user input into a sequential modification tabular ####
def ExtractSequence(self):
""" Converting the string user input into a tabular giving the sequence
of modification to apply
-
+
Args:
- None: self.AmendementsSequence of the structure should be filled
-
+ None: self.AmendementsSequence of the structure should be filled
+
Returns:
Tab (string list): the unitary modification list as strings, each in the format ModifyierType_Index(ApplicationTarget)
"""
-
+
# ------------------- Initialisations ---------------------------------------
string = "("+self.AmendementsSequence+")"
curr, Compact = [""], []
level = -1
-
+
# ------------------ Extracting the depth of each action --------------------
# (designed for e.g. a limiter applied on the top of a filter etc)
-
+
# Spanning each character of the string to find the max depth
for c in string:
# Detecting a new level
if c == "(":
- if curr != "":
- # Creating a new level in the level-list
- if len(Compact)<=level:
+ if curr != "":
+ # Creating a new level in the level-list
+ if len(Compact)<=level:
Compact += [[]]
curr += [""]
-
+
# Initialising the current level
level = level+1
curr[level] = ""
-
+
# Detecting an end of level
elif c==")":
# Safe expansion case if we went down the levl
- if len(Compact)<=level:
+ if len(Compact)<=level:
Compact += [[]]
curr += [""]
Compact[level] += [curr[level]]
-
+
# Expanding the upper levels by the one furnished there
for i in range(level-1, level):
curr[i] += "("+curr[level]+")"
-
+
# Reinitialising the level
curr[level] = ""
level -=1
-
+
# Otherwise, simply increment the considered string extraction
else: curr[level] += c
-
+
# ---- Converting the global modification as a sequence of unitary ones -------
# Initialise all the unitary modifiers
levels = len(Compact)
@@ -132,15 +133,15 @@ class Amendements():
# Span each level starting by the most nested one (to apply directly to the original residuals)
for lvl in range(levels)[::-1]:
-
+
# Get the iterative construction tabular for that level
for term in Compact[lvl]:
# Replace the terms that would already have been computed in the string litterals
for (i,te) in enumerate(TabTot[::-1]):
if te in term:
- # Spanning all occurences of the string and replacing the already computed
- # terms with its index while leaving the non-already computed terms with
- # similar syntax: e.g. if only Limiter_1 and Limiter_4 on the initial residual
+ # Spanning all occurences of the string and replacing the already computed
+ # terms with its index while leaving the non-already computed terms with
+ # similar syntax: e.g. if only Limiter_1 and Limiter_4 on the initial residual
# were already computed, in Limiter_4(Limiter_1)+Limiter_4, only the last term
# and Limiter_1 will be replaced by the corresponding intermediate residual indices.
search = te.replace("(","\(").replace(")","\)").replace("+","\+").replace("-","\-")
@@ -148,10 +149,10 @@ class Amendements():
k=0
for j in occurences:
idx = j-k*(len(te)-1)
- if (idx+len(te)>= len(term)) or not term[idx+len(te)]=="(":
+ if (idx+len(te)>= len(term)) or not term[idx+len(te)]=="(":
k=k+1
term = term[:idx]+"{0:d}".format(len(TabTot)-i-1)+term[idx+len(te):]
-
+
# Separate the different terms that are on the same level
uu = list(filter(None, re.split("[+\-*/]", term)))
@@ -169,10 +170,10 @@ class Amendements():
term2 = term2.replace("-{0:d}".format(i), "-{0!s}".format(TabTot[i])).replace("{0:d}-".format(i), "{0!s}-".format(TabTot[i]))
term2 = term2.replace("/{0:d}".format(i), "/{0!s}".format(TabTot[i])).replace("{0:d}/".format(i), "{0!s}/".format(TabTot[i]))
term2 = term2.replace("*{0:d}".format(i), "*{0!s}".format(TabTot[i])).replace("{0:d}*".format(i), "{0!s}*".format(TabTot[i]))
-
+
# Register each member of the same level as terms that have to be computed independently
if len(uu)>1:
- for u in uu:
+ for u in uu:
if not (u.isdigit() or u in Tab):
# Register the practical table
Tab.append(u)
@@ -181,7 +182,7 @@ class Amendements():
for i in range(len(Tab))[::-1]:
if "({0:d})".format(i) in littu: littu = littu.replace("({0:d})".format(i), "({0!s})".format(TabTot[i]))
TabTot.append(littu)
-
+
# Register in the practical table the aggregation term defining a same level
for (i,tt) in enumerate(Tab[::-1]):
if tt in term:
@@ -190,51 +191,51 @@ class Amendements():
k=0
for j in occurences:
idx = j-k*(len(tt)-1)
- if (idx+len(tt)>= len(term)) or not term[idx+len(tt)]=="(":
+ if (idx+len(tt)>= len(term)) or not term[idx+len(tt)]=="(":
k=k+1
term = term[:idx]+"{0:d}".format(len(TabTot)-i-1)+term[idx+len(tt):]
-
+
# If not already computed, register the new element
if not (term.isdigit() or term in Tab):
TabTot.append(term2)
Tab.append(term)
-
+
# Returning the unitary modifications sequence
- return(Tab)
-
+ return(Tab)
+
#### Interfacing the user wishes with the actual routines #####
def StringToResu(self, stringList):
""" Interfacing the user wishes with the actual routines, according to the unitary modification sequence list.
-
+
Args:
- StringList (list of strings): the unitary modification list, each in the format ModifyierType_Index(ApplicationTarget)
-
+ StringList (list of strings): the unitary modification list, each in the format ModifyierType_Index(ApplicationTarget)
+
Returns:
None: fills directly the self.MolliProcess of the structure, a tabular containing each step of modification [modifier, ApplicationTarget]
"""
-
+
# ------------------------- Initialisations -----------------------------------------
self.MolliProcess = [[]]*len(stringList)
# ------------- Converting the strings to callback depending on their type ----------
# Spanning each instruction line contained in the list
for (i, string) in enumerate(stringList):
-
+
# Modifiers of the type Modifyer_Index(Target)
if "(" in string and ")" in string:
# Extract the parameters if any
- if re.findall("\<(.*?)\>", string):
+ if re.findall("\<(.*?)\>", string):
Properties = list(filter(None, re.split(",| ", list(filter(None, re.findall("\<(.*?)\>", string)))[0])))
string = re.sub('\<(.*?)\>', '', string)
else: Properties = []
-
+
# Map to the right method
Type = (string.split("(")[0]).split("_")[0]
Index = (string.split("(")[0]).split("_")[1]
ResiInd = int(list(filter(None, re.findall("\((.*?)\)", string)))[0])
Modif = eval("Mappers."+Type+"("+Index+")."+Type+"(self.Problem, self.Mesh, *Properties)")
self.MolliProcess[i] = [copy.deepcopy(Modif.Mollify), ResiInd+1]
-
+
# Modifiers of the type Target1 operator Target2 operator Target3
elif "+" in string or "-" in string or "*" in string or "/" in string:
pieces = list(filter(None, re.split("[\*+\-/]", string)))
@@ -243,65 +244,61 @@ class Amendements():
if pieces[p].isdigit(): pieces[p] = "resu[int({0:d}),:,:]".format(p)
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))]
-
+
# 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):
+ if re.findall("\<(.*?)\>", string):
Properties = list(filter(None, re.split(",| ", list(filter(None, re.findall("\<(.*?)\>", string)))[0])))
string = re.sub('\<(.*?)\>', '', string)
else: Properties = []
-
+
# Map to the right method
Type = (string.split("(")[0]).split("_")[0]
Index = (string.split("(")[0]).split("_")[1]
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
-
+
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 modified residuals (NbVars x NbElementDofs)
-
+
.. warning::
the fluxes are not changing through the intermediary steps
-
+
---------------------------------------------------------------------"""
# Checking that there is indeed the wish of a mollifier
if self.MolliProcess == 0: return(resu)
-
+
# --------------- Initialisation ------------------------------------------
# Create a structure gathering all the intermediary steps' residuals
NbSteps = len(self.MolliProcess)+1
ModRes = np.zeros((NbSteps, np.shape(resu)[0], np.shape(resu)[1]))
ModRes[0,:,:] = copy.deepcopy(resu)
-
+
# --------------- Applying the sequencee of modifications -----------------
for (j,mol) in enumerate(self.MolliProcess):
ModRes[j+1,:,:] = mol[0](SolBuff, ModRes[mol[1],:,:], fluxes, i, du, dt)
-
+
# Return the last (and final) modified residual
return(ModRes[-1,:,:])
-
-
-
-
diff --git a/SourceCode/_Solver/SpatialSchemes/CG.py b/SourceCode/_Solver/SpatialSchemes/CG.py
index 20ce2a2..9709c16 100644
--- a/SourceCode/_Solver/SpatialSchemes/CG.py
+++ b/SourceCode/_Solver/SpatialSchemes/CG.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,71 +100,71 @@ class Scheme():
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
+
+ # 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]
-
+
# 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)
-
- # ---- Computing the plain integral
- # Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
+
+ # ---- 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)
-
+
# 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, :, :] = self.Problem.GoverningEquations.Flux(Uloc, qp, Flags)
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)
-
+
# Returning the full array
return(Flux)
@@ -172,99 +172,99 @@ 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
+
+ # ---- 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)
-
+ 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)
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
@@ -273,40 +273,34 @@ class Scheme():
#### 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.
-
+ """ Routine that maps the values from the Dofs to the Physical vertices of the mesh
+ of the solution.
+
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))
-
+ 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]))
-
+ RSO = np.mean(Solution.Sol[:,index], axis=1)
+
# Fill the reconstructed solution upon the found index
- Solution.RLSValues[:,i] = RS
- Solution.RFluidFlag[i] = FF
Solution.RSol[:,i] = RSO
diff --git a/SourceCode/_Solver/SpatialSchemes/CG_LFX.py b/SourceCode/_Solver/SpatialSchemes/CG_LFX.py
index 18d3ca4..e227efb 100644
--- a/SourceCode/_Solver/SpatialSchemes/CG_LFX.py
+++ b/SourceCode/_Solver/SpatialSchemes/CG_LFX.py
@@ -106,7 +106,7 @@ class Scheme():
# 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])
@@ -124,14 +124,14 @@ class Scheme():
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)
-
+
# ---- Computing the plain integral
# Retrieving the quadrature points as barycentric coordinates and convert them to cartesian
# in the Dof's vertex ordering referential
@@ -145,7 +145,7 @@ class Scheme():
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
Flux[0, i, 0, :, :, 0:12] = self.Problem.GoverningEquations.Flux(Uloc, qp, FlagsLoc)
@@ -291,7 +291,7 @@ class Scheme():
#### 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.
+ of the solution.
Args:
Solution (Solution): the currently being computed solution
@@ -306,9 +306,7 @@ class Scheme():
----------------------------------------------------------------------------------- """
# -------------- 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))
+ Solution.RSol = np.zeros((np.shape(Solution.Sol)[0], self.Mesh.NbMeshPoints))
# ------------- Reconstructing the values -----------------------------------------
# Extract the cartesian coordinates of the Dofs
@@ -319,11 +317,7 @@ class Scheme():
# 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]))
+ RSO = np.mean(Solution.Sol[:,index], axis=1)
# Fill the reconstructed solution upon the found index
- Solution.RLSValues[:,i] = RS
- Solution.RFluidFlag[i] = FF
Solution.RSol[:,i] = RSO
diff --git a/SourceCode/_Solver/SpatialSchemes/CG_Primary.py b/SourceCode/_Solver/SpatialSchemes/CG_Primary.py
index 348a474..b85cc4d 100644
--- a/SourceCode/_Solver/SpatialSchemes/CG_Primary.py
+++ b/SourceCode/_Solver/SpatialSchemes/CG_Primary.py
@@ -67,7 +67,7 @@ class Scheme():
.. note::
- This scheme is the CG one where the variables interpolated are rho, rhoU,
+ 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
@@ -112,7 +112,7 @@ class Scheme():
# 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])
@@ -130,19 +130,19 @@ class Scheme():
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
@@ -289,7 +289,7 @@ class Scheme():
#### 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.
+ of the solution.
Args:
Solution (Solution): the currently being computed solution
@@ -304,9 +304,7 @@ class Scheme():
----------------------------------------------------------------------------------- """
# -------------- 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))
+ Solution.RSol = np.zeros((np.shape(Solution.Sol)[0], self.Mesh.NbMeshPoints))
# ------------- Reconstructing the values -----------------------------------------
# Extract the cartesian coordinates of the Dofs
@@ -317,11 +315,7 @@ class Scheme():
# 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]))
+ RSO = np.mean(Solution.Sol[:,index], axis=1)
# Fill the reconstructed solution upon the found index
- Solution.RLSValues[:,i] = RS
- Solution.RFluidFlag[i] = FF
- Solution.RSol[:,i] = RSO
+ Solution.RSol[:,i] = RSO
diff --git a/Templates/Template_Equations.py b/Templates/Template_Equations.py
index 52cf70a..e203ee2 100644
--- a/Templates/Template_Equations.py
+++ b/Templates/Template_Equations.py
@@ -22,11 +22,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)
@@ -48,12 +48,13 @@ class Equations():
Args:
FluidProp (list of FluidModel): the list of fluid property instances (NbFluids)
EOS (list of callbacks): the list of Equation of state functions (NbFluids)
-
+
|
-
+
.. rubric:: Methods
- """
+ """
+
# Defines the (conservative) variables number and names, in a latex typography
self.NbVariables = 0
self.VariablesNames = ["The variables name in a string list"]
@@ -62,39 +63,83 @@ class Equations():
self.VariablesLatexNames = [r'$ The list of the variables name in a latex format $']
self.VariablesLatexUnits = [r'$ The list of the variables units in a latex format $']
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
+ """
+
+ # ------- Constructing the primary variables ----------------------------------------------------
+ # COMPUTE HERE THE PRIMARY VARIABLES
+
+ # Returning the Result
+ return(PrimVars)
+
+ #### Retrieve the conservative variables from the primary ones ####
+ def PrimaryToConservative(self, PrimVar, FluidIndex, *args):
+ """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 conservative variables ----------------------------------------------------
+ # COMPUTE HERE THE CONSERVATIVE VARIABLES
+
+ # 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:
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
"""
- # DEFINE THE FLUX FX (NbVars x NbPoints) AND FY (NbVars x NbPoints) OF
+ # DEFINE THE FLUX FX (NbVars x NbPoints) AND FY (NbVars x NbPoints) OF
# THE SPATIALLY 2D GOVERNING EQUATION THERE
# Returning the flux
@@ -105,121 +150,200 @@ class Equations():
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
"""
# COMPUTE THE FLUID VELOCITY VALUES VX AND VY FROM THE SET OF VARIABLES THERE
-
+
# 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]
-
- # EDIT THE INITIALISATION OF J BY THE NUMBER OF VARIABLES OF YOUR EQUATION SET
+
+ # EDIT THE INITIALISATION OF J BY THE NUMBER OF VARIABLES OF YOUR EQUATION SET
J = np.zeros((self.NbVariables, self.NbVariables, 2, NbPoints))
-
+
# FILL THE JACOBIAN IN X THERE (EXAMPLE OF FORMAT FOR FOUR EQUATIONS)
#J[0,:,0,:] = np.array([np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints)])
#J[1,:,0,:] = np.array([np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints)])
#J[2,:,0,:] = np.array([np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints)])
#J[3,:,0,:] = np.array([np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints)])
-
+
# FILL THE JACOBIAN IN Y THERE
#J[0,:,1,:] = np.array([np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints)])
#J[1,:,1,:] = np.array([np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints)])
#J[2,:,1,:] = np.array([np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints)])
#J[3,:,1,:] = np.array([np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints), np.zeros(NbPoints)])
-
- # Returning the Jacobian
+
+ # Returning the Jacobian
return(J)
#### 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)
"""
-
+
# DEFINE THE SPECTRAL RADIUS OF THE SYSTEM THERE
# 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 NbEigs) the 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
+ # COMPUTE HERE THE EIGENVALUES AT EACH GIVEN POINT
+
+ # 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]))
+ # COMPUTE HERE THE RIGHT EIGENVECTORS AT EACH GIVEN POINT
+
+ # 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]))
+ # COMPUTE HERE THE LEFT EIGENVECTORS AT EACH GIVEN POINT
+
+ # Returning the values
+ return(LEV)
+
+
# ==============================================================================
# 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 :code:`Chosen Fluid Model` definition.
+ used in combination with the :code:`Chosen Fluid Model` definition.
Upon creation, fills the field :code:`EOS` with the function callback
- corresponding to the wished initialisation function. The implemented EoS
+ corresponding to the wished initialisation function. The implemented EoS
are linked as follows. See their respective documentation for more information.
1. | EOS definition Template
@@ -231,27 +355,27 @@ 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
"""
# ADAPT HERE THE MAPPING TO THE EOS
if Id==1: self.EOS = self.EOSTemplate
-
+
else: raise NotImplementedError("Error in the selection of the Equation of State\n\
Unknown index. Check your Problem definition. Abortion.")
#### EOS for stiff fluids #####
def EOSTemplate(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:
ThewishedQuantity (float numpy array): the pressure values at the given points (NbGivenPoints)
"""
@@ -272,11 +396,11 @@ 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. | ICTemplate
-
+
"""
#### Automatic initialisation routine (mapper) ####
@@ -284,46 +408,46 @@ class InitialConditions():
"""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.ICTemplate
-
+
else: raise NotImplementedError("Error in the selection of the Initial conditions\
Unknown index. Check your Problem definition. Abortion.")
#### SOD test case ####
def ICTemplate(self, PointsID, Mesh, EOS, FluidProp, *argv):
- """Initialising the solution to a constant state that matches the SOD-inner value.
+ """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
"""
-
+
# CREATE THE INIT VECTOR ACCORDING TO THE NUMBER OF VARIABLES
NbVars = 0
Init = np.zeros((NbVars, len(PointsID)))
-
+
# DEFINE HERE THE INITIALISATION VALUES
-
+
# Returning the values
return(Init)
diff --git a/Templates/Template_FluidSelector.py b/Templates/Template_FluidSelector.py
index d70973e..f46ca7f 100644
--- a/Templates/Template_FluidSelector.py
+++ b/Templates/Template_FluidSelector.py
@@ -1,4 +1,4 @@
-# Template to copy in _Solver/FluidSelectors/
+# Template to copy in _Solver/FluidSelectors/
# Add the fluid slector from its module name under a new index in the file Mappers.py
@@ -23,20 +23,20 @@ import BarycentricTools as BT
# ==============================================================================
class FS():
""".. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Classical level set approach, initialised according to a signed distance and iterated
with a continuous galerkin - lfx scheme, without redistancing.
-
+
.. 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). It also contains the
+ the structure with the same names (see the parameters below). It also contains the
further field
-
+
- |bs| **FluidSelectorName** *(string)* -- the name of the level set
-
+
|
-
+
"""
#### Automatic initialisation routine ####
@@ -46,17 +46,17 @@ class FS():
Problem (Problem): the considered problem
Mesh (MeshStructure): the considered mesh instance
Params (string list): the level set's parameters as wished by the user
-
+
|
-
+
.. rubric:: Methods
"""
-
-
+
+
# Setting the name of the class to be accessible from outside
self.FluidSelectorName = "Naive Level Set -- Continuous Galerkin"
-
- # Setting the mesh and parameters for an internal access
+
+ # Setting the mesh and parameters for an internal access
self.Problem = Problem
self.Mesh = Mesh
self.Params = Params
@@ -68,10 +68,10 @@ class FS():
#### External initialisation routine ####
def Initialiser(self, Solution):
""" External initialisation routine.
-
+
Args:
Solution (Solution structure): the solution instance the scheme is working on
-
+
Returns:
LSValues (float numpy array): the level set values at the MeshVertex and Dofs, to be copied to the Solution structure (NbSubdomains x (NbMeshPoints+NbDofs))
"""
@@ -82,7 +82,7 @@ class FS():
# 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
+ NbDofs = len(Solution.FluidFlag) # Only dofs are there
# Initialising the fluid selector values
LSValues = np.zeros((len(Subdomains), len(Solution.FluidFlag)))
@@ -95,187 +95,187 @@ class FS():
#### Gives the dynamic of the Fluid's spotter through its flux ####
def Flux(self, Var, Motion, *args):
""" Flux function of the (conservative) EulerEquations
-
- Args:
+
+ Args:
Var (float numpy array): the value of the Level Set values (NbFluids x NbGivenPoints)
- Motion (float numpy array): the velocity array (2xNbPoints)
-
+ Motion (float numpy array): the velocity array (2xNbPoints)
+
Returns:
Flux (3D numpy array): the obtained flux (spatialdim x NbVars x NbPoints)
-
- .. note::
-
+
+ .. note::
+
*args is there for compatibility reason when e.g. a location xy is given when called
-
+
--------------------------------------------------------------------------------------------------------------------------------"""
# DEFINE THE EQUATION DRIVING YOUR FLUID SELECTOR HERE
-
+
# Returning the values
return(Flux)
-
+
#### Jacobian of the equations set ####
def Jacobian(self, Var, Motion, *args):
""" .. |bs| unicode:: U+00A0 .. NO-BREAK SPACE
-
+
Computes the Jacobian of the flux for the advection of the LevelSet
-
+
Args:
Var (float numpy array): the value of the Level Set values (NbFluids x NbGivenPoints)
- Motion (float numpy array): the velocity array (2xNbPoints)
-
- Returns:
- J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
-
+ Motion (float numpy array): the velocity array (2xNbPoints)
+
+ Returns:
+ J (3D numpy array): J[:,:,i] gives the jacobian of the flux taking care of the dynamic of the ith spatial coordinate.
+
|
-
+
.. note::
-
+
- |bs| 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]
-
+
- |bs| *args is there for compatibility reason when e.g. a location xy is given when called the x-y locations at which the variables are given (NbGivenPoints x 2)
-
+
---------------------------------------------------------------------"""
-
+
# DEFINE THE JACOBIAN ASSOCIATED TO THE ABOVE FLUX HERE
-
- # Returning the Jacobian
- return(J)
-
+
+ # Returning the Jacobian
+ return(J)
+
#### Spectral radius of the equations set ####
def SpectralRadius(self, Var, Motion, n, *args):
""" Computes the spectral radius associated to the motion flux.
-
- Args:
+
+ Args:
Var (2D numpy array): the variables of the problem (NbVars x NbGivenPoints)
- Motion (float numpy array): the velocity array (2xNbPoints)
+ Motion (float numpy array): the velocity array (2xNbPoints)
n (2D numpy array): the x-y values of the normal at each given point (NbGivenPoints x 2)
-
- Returns:
+
+ Returns:
Lmb (numpy array): the array containing the spectral radius at each given point
-
+
.. note::
-
+
*args is there for compatibility reason when e.g. a the xy locations at which the variables are given (NbGivenPoints x 2)
-
+
---------------------------------------------------------------------"""
-
+
# Retrieving the Jacobian
J = self.Jacobian(Var, Motion)
-
+
# Computing the spectral radius at each given point, shortcuted since here we are scalar
Lmb = np.abs([np.sum(J[i,i,:,:]*n.T, axis = 0) for i in range(self.Problem.NbFluids)])
-
+
# Returning the obtained values
return(Lmb)
-
+
# ************************************************************************
# * Definition of the Fluid spotter's spatial impact and rectification *
# ************************************************************************
-
+
#### 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
+ """ Flags the degrees of freedom according to the different values of the
associated level set.
-
- Args:
+
+ Args:
Solution (Solution structure): the solution instance the scheme is working on
-
- Returns:
+
+ Returns:
None: updating directly the value of FluidFlag in the solution instance
-
+
--------------------------------------------------------------------------------------- """
-
-
+
+
# FLAG THE DEGREES OF FREEDOM ACCORDING TO THE FLUID SELECTOR VALUES HERE
-
-
+
+
# ---------------- Updating the flags at Dofs ------------------------------
Solution.FluidFlag = 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
+ """ Method modifying the subdomains associated to the Solution upon the values of
LSValues.
-
- Args:
+
+ Args:
Solution (Solution): the solution instance the scheme is working on
-
- Returns:
+
+ Returns:
None: updating directly the subdomains in the solution instance
-
- .. warning:: Not safe yet if the first path given out is internal
-
+
+ .. warning:: Not safe yet if the first path given out is internal
+
---------------------------------------------------------------------------- """
-
+
# Update each of the fluid's locations
for i in range(self.Problem.NbFluids):
-
+
# UPDATE THE FLUID'S SUBDOMAINS HERE
-
+
Solution.Subdomains[i] = polygon
-
-
+
+
#### Redistancing the LSvalues according to the freshly computed level set values ####
def Redistancing(self, Solution):
""" Redistances the Level set values upon the computed subdomains.
-
- Args:
+
+ Args:
Solution (Solution): the solution instance the scheme is working on
-
- Returns:
+
+ Returns:
None: updating directly the LSValues in the solution instance
-
+
----------------------------------------------------------------------------------"""
for i in range(self.Problem.NbFluids):
-
+
# REDISTANCE THE FLUID SELECTOR VALUE IF NEEDED THERE
-
+
# Computing the signed distance
Solution.LSValues[i,:] = np.prod(Value, axis=0)
# ****************************************************************************
# * Definition of the Fluid spotter's evolution's process (e.g. advection) *
# ****************************************************************************
-
+
#### Emulates a flux computation as done in the iteration routine ####
def ComputeFlux(self, U, LSValues):
""" Emulates a flux computation as done in the iteration routine in order
to be used in the DeC subtimestep evaluation.
-
- Args:
+
+ Args:
U (numpy float array): buffer containing the current state values at each dof (NbVars x NbDofs)
LS (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 requires by the UpdateFSValues routine
-
+
-------------------------------------------------------------------------------------------------------- """
-
- # DEFINE HERE THE FLUXES INVOLVED IN THE FLUIDS'S SELECTOR SPATIAL ITERATION SCHEME
-
+
+ # DEFINE HERE THE FLUXES INVOLVED IN THE FLUIDS'S SELECTOR SPATIAL ITERATION SCHEME
+
# Returning the full array
return(fluxes)
#### (Crucial routine) Spatial scheme to update the level set ####
def Iteration(self, Solution, fluxes, i, du=0, dt=1):
""" Main iteration of the scheme, implementing the continuous galerkin scheme.
-
+
Args:
Solution (Solution): 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)
-
+
---------------------------------------------------------------------"""
# DEFINE HERE THE SPATIAL SCHEME TO EVOLVE THE FLUID SELECTOR VALUES
@@ -286,17 +286,17 @@ class FS():
#### 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
-
+ """ Wrapper routine to update the quantitites that depend on the UpdateFSValues fresh output
+
Args:
Solution (Solution): the solution instance the scheme is working on
-
+
Returns:
None: updating directly the subvalues in the solution instance
-
+
--------------------------------------------------------------------------------- """
-
+
# DEFINE HERE THE SEQUENCE OF ACTION POST-FLUID SPOTTER RESIDUAL COMPUTATION
# OF ONE (SUB)TIME-STEP
-
+
return(0)
diff --git a/Templates/Template_Problem.py b/Templates/Template_Problem.py
index 134b8e6..3b60ee1 100644
--- a/Templates/Template_Problem.py
+++ b/Templates/Template_Problem.py
@@ -10,8 +10,8 @@ class ProblemData():
#*******************************************************************************
# Define a name for the problem
self.ProblemID = "TheProblemName"
-
- # Select the geometry associated to the problem
+
+ # Select the geometry associated to the problem
# (the root name of the target mesh stored in the PredefinedTests/PredefinedMeshes, without the resolution)
self.MeshName = "AssociatedMesh"
@@ -20,7 +20,7 @@ class ProblemData():
#*******************************************************************************
# Index of the governing equation according to the Mappers.py file
self.GoverningEquationsIndex = 0
-
+
# Equation of state to consider for each fluid
# (the length of the list matches the number of fluids)
self.EOSEquationsIndex = [0,0,0]
@@ -36,9 +36,9 @@ class ProblemData():
# Define the subdomains on which the second and third fluids will live
# (see the documentation of the initial subdomains definition)
self.FluidInitialSubdomain = [0,0]
-
+
# Index of the initial conditions to apply on each subdomain
- # (the index is associated to the Governing equation. See the doc of the equation
+ # (the index is associated to the Governing equation. See the doc of the equation
# targetted by the GoverningEquationsIndex)
self.InitialConditions = [(0,),(0,),(0,)]
@@ -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,)],\
+ [(0,), (0,), (0,)],\
+ [(0,), (0,), (0,)]])
diff --git a/Templates/Template_Settings.txt b/Templates/Template_Settings.txt
index 8293b64..377e5d4 100644
--- a/Templates/Template_Settings.txt
+++ b/Templates/Template_Settings.txt
@@ -17,4 +17,4 @@ Limiter_1<0.2> # Additives list Index of limiter (0=None, see the docs for the s
# Export properties
m # Verbose mode (n=none, m=minimal, d=debug) (yet inactive)
-1 # Time steps export intervals (yet inactive)
+0.05 # Time export intervals
diff --git a/Templates/Template_SpatialScheme.py b/Templates/Template_SpatialScheme.py
index ffebca6..c11ed8c 100644
--- a/Templates/Template_SpatialScheme.py
+++ b/Templates/Template_SpatialScheme.py
@@ -21,7 +21,7 @@ SchemeName = "The scheme name"
# 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
@@ -30,25 +30,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 = "The mesh type in a Fenics format"
self.MeshOrder = VariationalOrder
@@ -60,16 +60,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 ####
@@ -79,13 +79,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
@@ -93,76 +93,82 @@ 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)
Returns:
fluxes (numpy float (multiD)array): fluxes in the format required by the Iteration routine
+ |
+
+ .. note::
+ If the scheme should be used with a filter, make sure that the last row of fluxes corresponds to the format
+ wished by the filters (fluxes[lastrowoffluxes, elementnumber, 0, :, :,1:NbDofs] = the fluxes computed at the dofs (dim x NbVars x NbDofs))
+
"""
-
-
- # DEFINE HERE THE FLUX TO BE A-PRIORI COMPUTED (ESP. FOR THE COMPATIBILITY WITH THE DEC ITERATION)
+
+
+ # DEFINE HERE THE FLUX TO BE A-PRIORI COMPUTED (ESP. FOR THE COMPATIBILITY WITH THE DEC ITERATION)
# AND PASSED TO THE MAIN ITERATION
-
-
+
+
# Returning the full array
return(Flux)
#### Main routine, definining the scheme ####
def Iteration(self, Solution, fluxes, i, du=0, dt=1):
"""Your scheme description
-
+
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)
"""
- # DEFINE HERE THE MAIN SPATIAL ITERATION SCHEME
+ # DEFINE HERE THE MAIN SPATIAL ITERATION SCHEME
# USE THE GIVEN PARAMETER FLUX FOR COMPATIBILITY WITH DEC REASONS
-
+
# 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
+ """ 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):
-
+
# RECONSTRUCT HERE YOUR SOLUTION AT THE SOLUTION VERTICES
-
+
# Fill the reconstructed solution upon the found index
Solution.RLSValues[:,i] = RS
Solution.RFluidFlag[i] = FF
diff --git a/Templates/Template_TemporalScheme.py b/Templates/Template_TemporalScheme.py
index fbf6e5d..041df6b 100644
--- a/Templates/Template_TemporalScheme.py
+++ b/Templates/Template_TemporalScheme.py
@@ -9,7 +9,7 @@
import copy
# Custom libraries import
-import _Solver.Solver_Spatial as SPL
+import _Solver.Solver_Spatial as SPL
# ==============================================================================
@@ -23,35 +23,35 @@ import _Solver.Solver_Spatial as SPL
def TheRoutineName(schemeSC, schemeFS, dt, Mesh, Sol, M, K, typpe = "Equispaced"):
""" Provides a full time iteration according to the DeC method
-
+
.. rubric:: Functional inputs
-
+
|bs| -- **schemeSP** *(scheme module list)* -- handler of the spatial scheme module, should be vectorised for using on several (spatial) points simultaneously followed (if any) by the list of its amendments (limiter, jumps etc)
-
+
|bs| -- **schemeFS** *(scheme module list)* -- handler of the FluidSpotter scheme module, should be vectorised for using on several (spatial) points simultaneously followed (if any) by the list of its amendments (limiter, jumps etc)
-
+
.. rubric:: Time-dependent inputs
-
+
|bs| -- **dt** *(float)* -- time span on which to evolve the solution
-
+
|bs| -- **Sol** *(Solution structure)* -- the initial condition, as a solution structure
-
+
.. rubric:: Input parameters:
-
+
|bs| -- **M** *(integer)* -- the number of subtimesteps
-
+
|bs| -- **K** *(integer)* -- the number of corrections
-
+
|bs| -- **typpe** *(string)* -- type of the interpolation (quadrature) nodes to be used: "Equispaced", "GaussLobatto". Default is set as Equispaced
.. rubric:: Returns
-
+
|bs| -- *Tspan (float)* -- the full time step interval
-
+
|bs| -- *U (float numpy array)* -- the updated solution (NbVars x (NbGivenPoints))
-
+
----------------------------------------------------------------------------------------"""
-
+
# --------- Initilising the variables --------------------------------------
# Initialising the variables of interest and buffering partially the quantities
# of interest in the solution (either the Solution itself or the level set only)
@@ -60,11 +60,11 @@ def TheRoutineName(schemeSC, schemeFS, dt, Mesh, Sol, M, K, typpe = "Equispaced"
LS = Sol.LSValues
# --------- Process the time iteration --------------------------------------
-
+
# DEFINE HERE YOUR TEMPORAL SCHEME
-
+
# Call the spatial scheme with a call analogous to
# SPL.RD(Theta, m, M, dt, schemeSC, schemeFS, Mesh, SolBuff, up, lp)
-
+
# Returning the full timestep
return([U_updated, LS_updated])
diff --git a/_doc/source/ChildPages/DesignSettings.rst b/_doc/source/ChildPages/DesignSettings.rst
index 78c1f01..c8bdfb6 100644
--- a/_doc/source/ChildPages/DesignSettings.rst
+++ b/_doc/source/ChildPages/DesignSettings.rst
@@ -19,13 +19,13 @@ in time and a second order Deferred Correction in time. We will name this settin
.. contents::
:local:
:depth: 1
-
+
|
----------------------------
-1. Create a new settings file
+1. Create a new settings file
---------------------------
Copy and paste the template inside the folder where the predefined settings are stored.
@@ -35,7 +35,7 @@ Copy and paste the template inside the folder where the predefined settings are
cp Templates/Template_Settings.txt PredefinedTests/PredefinedSettings/Spatial1Time2_CGLFX.txt
Open it with your favorite text editor to modify its values.
-
+
|
|
|
@@ -45,15 +45,15 @@ Open it with your favorite text editor to modify its values.
2. Set the spatial solver
---------------------------
-The spatial solver consists of a scheme, its properties and possibly a sequence of
+The spatial solver consists of a scheme, its properties and possibly a sequence of
modifications to apply to the achieved residuals.
|
.. rubric:: Select the spatial scheme
-To select the scheme one wants to use, get the list of implemented schemes and their
-associated index either in the file "Mappers.py" in the :code:`SourceCode`, or by using the
+To select the scheme one wants to use, get the list of implemented schemes and their
+associated index either in the file "Mappers.py" in the :code:`SourceCode`, or by using the
technical documentation of the class :code:`Spatial_Scheme` on the page :ref:`Mapping`.
|
@@ -63,7 +63,7 @@ In our case, we can use one of the following schemes.
.. code::
$ python -c "import SourceCode.Mappers as MP; print(MP.Spatial_Scheme.__doc__)"
-
+
# Find the lines similar to the following
*Implemented spatial schemes given the id*
1. | Continuous Galerkin
@@ -72,53 +72,53 @@ In our case, we can use one of the following schemes.
We select the scheme indexed by "2" and modify accordingly the first line of the given template.
.. code::
-
+
1 # Index of the scheme (see Mappers.py to know the indices meaning)
-
+
|
Specify then the properties of the scheme one line below, according to its given format. To find it, get the
scheme's argumenst list by querying the documentation of the scheme in the technical
-handbook (in the relevant rubric of :ref:`SpatialSchemes`), or in the module itself.
+handbook (in the relevant rubric of :ref:`SpatialSchemes`), or in the module itself.
In our case, the parameters that should be given are the wished order and the type of basis functions.
Indicate them on the second line, separated by spaces.
.. code::
- 1 B # Scheme properties
-
+ 1 B # Scheme properties
+
|
|
.. rubric:: Define a residual's amendement sequence
-It is possible to define a sequence of modifications to apply to the residuals computed by the selected
-spatial scheme. The specification has to be given in the third line of the given template,
-following the guidelines of the :code:`Usage` rubric of :ref:`SpatialModifiers`.
-The indices of the mollifiers can be looked for in the file "Mappers.py" located in the :code:`SourceCode`, or by using the
+It is possible to define a sequence of modifications to apply to the residuals computed by the selected
+spatial scheme. The specification has to be given in the third line of the given template,
+following the guidelines of the :code:`Usage` rubric of :ref:`SpatialModifiers`.
+The indices of the mollifiers can be looked for in the file "Mappers.py" located in the :code:`SourceCode`, or by using the
technical documentation of the classes :code:`Limiter` and :code:`Filtering` on the page :ref:`Mapping`.
In our case, we have for choices:
.. code::
$ python -c "import SourceCode.Mappers as MP; print(MP.Limiter.__doc__)"
-
+
# Find the lines similar to the following
*Implemented amendements given the id* and their type
- - Limiters
+ - Limiters
1. | Psi limiter
-
-
+
+
$ python -c "import SourceCode.Mappers as MP; print(MP.Filtering.__doc__)"
-
+
# Find the lines similar to the following
*Implemented amendements given the id* and their type
- - Filtering
+ - Filtering
1. | Streamline
Think to check each of the selected modifiers definition to check if they need arguments.
-There, we choose to apply a Psi limiter, not requiring any argument, on the original residual.
+There, we choose to apply a Psi limiter, not requiring any argument, on the original residual.
If you do not want to apply any modification, just write "0" instead.
.. code::
@@ -130,8 +130,8 @@ If you do not want to apply any modification, just write "0" instead.
.. rubric:: Choose a fluid selector
-To select the fluid selector one wants to use, get the list of implemented ones and their
-associated index either in the file "Mappers.py" located in :code:`SourceCode`, or by using the
+To select the fluid selector one wants to use, get the list of implemented ones and their
+associated index either in the file "Mappers.py" located in :code:`SourceCode`, or by using the
technical documentation of the class :code:`Fluid_Selector` on the page :ref:`Mapping`.
|
@@ -141,7 +141,7 @@ In our case, we can use one of the following schemes.
.. code::
$ python -c "import SourceCode.Mappers as MP; print(MP.Fluid_Selector.__doc__)"
-
+
# Find the lines similar to the following
*Implemented fluid's selectors*
1. | NaiveLS_CG
@@ -151,21 +151,21 @@ In our case, we can use one of the following schemes.
We select the scheme indexed by "2" and modify accordingly the fourth line of the given template.
.. code::
-
+
2 # Type of fluid spotter (1=LS)
|
Specify then the properties of the fluid selector one line below, according to its given format. To find it, get the
argumenst list by querying the documentation of the fluid selector in the technical
-handbook (in the relevant rubric of :ref:`FluidSelectors`), or in the module itself.
+handbook (in the relevant rubric of :ref:`FluidSelectors`), or in the module itself.
In our case, the parameters that should be given are the wished order and the type of basis functions.
Indicate them on the second line, separated by spaces.
.. code::
1 B # Fluid spotter properties
-
+
|
|
|
@@ -181,8 +181,8 @@ The temporal solver consists of a time scheme and time iteration constraints.
.. rubric:: Select the time scheme
-To select the scheme one wants to use in time, get the list of implemented schemes and their
-associated index either in the file "Mappers.py" in the :code:`SourceCode`, or by using the
+To select the scheme one wants to use in time, get the list of implemented schemes and their
+associated index either in the file "Mappers.py" in the :code:`SourceCode`, or by using the
technical documentation of the class :code:`Temporal_Scheme` on the page :ref:`Mapping`.
|
@@ -192,7 +192,7 @@ In our case, we only use (for now) the DeC time stepping.
.. code::
$ python -c "import SourceCode.Mappers as MP; print(MP.Temporal_Scheme.__doc__)"
-
+
# Find the lines similar to the following
*Implemented temporal schemes given the id*
1. | Deffered correction (DeC)
@@ -201,21 +201,21 @@ In our case, we only use (for now) the DeC time stepping.
We select the temporal scheme indexed by "1" and modify accordingly the line 13 of the given template.
.. code::
-
+
1 # Time scheme (see Mappers.py to know the indices meaning)
-
+
|
-Specify then the properties of the temporal scheme one line below, according to its given format.
+Specify then the properties of the temporal scheme one line below, according to its given format.
To find it, get the time scheme's argumenst list by querying the documentation in the technical
-handbook (in the relevant rubric of :ref:`TemporalSchemes`), or in the module itself.
+handbook (in the relevant rubric of :ref:`TemporalSchemes`), or in the module itself.
In our case, the parameters that should be given are the wished order and the number of corrections.
Indicate them on the second line, separated by spaces.
.. code::
- 2 2 # Properties
-
+ 2 2 # Properties
+
|
|
@@ -226,7 +226,7 @@ We consider that any problem always starts at "t=0". Modify the lines 14 and 15
to your wishes.
.. code::
-
+
0.4 # CFL
3 # Tmax, in seconds
@@ -240,8 +240,26 @@ to your wishes.
|
+.. rubric:: Select the result's export periodicity
+
+Select the interval time of result's export in the results folder. Whenever the elapsed solution time between the currently computed solution
+and the previously exported result is bigger than the specified interval, the current solution will be exported in the result folder.
+
+.. code::
+
+ 0.4 # CFL
+ 3 # Tmax, in seconds
+
+.. note::
+
+ Please do not use an interval smaller than 1e-7 seconds, the filenames' output format is not as precise and your exported solution may overwrite previous ones.
+
+|
+
+.. rubric:: Select the console verbose output mode
+
.. warning::
-
+
TODO
|
@@ -255,5 +273,5 @@ to your wishes.
|
.. warning::
-
+
TODO
diff --git a/_doc/source/ChildPages/SpatialModifiers.rst b/_doc/source/ChildPages/SpatialModifiers.rst
index 386d6be..1a1438d 100644
--- a/_doc/source/ChildPages/SpatialModifiers.rst
+++ b/_doc/source/ChildPages/SpatialModifiers.rst
@@ -7,16 +7,16 @@ Spatial Modifiers
-Modifications applied to the result of a classical scheme such as limiter or filters have to be
-implemented in :code:`_Solver/SpatialModifiers` as a module. The modifier has to be implemented within
-a class whose name should match its type (e.g "Limiter" or "Filtering"), and the main method that
+Modifications applied to the result of a classical scheme such as limiter or filters have to be
+implemented in :code:`_Solver/SpatialModifiers` as a module. The modifier has to be implemented within
+a class whose name should match its type (e.g "Limiter" or "Filtering"), and the main method that
is intended for modifying the residuals should be named :code:`Mollify` in any case.
|
From a user perspective, it is possible to apply several mollifiers to the original residual by specifying
-the wished sequence as a litteral string. The input will then be interpreted by the class Amendements upon
-the instance creation. In practice, the input is parsed by the method :code:`ExtractSequence` and
+the wished sequence as a litteral string. The input will then be interpreted by the class Amendements upon
+the instance creation. In practice, the input is parsed by the method :code:`ExtractSequence` and
converted to the actual modifyer by :code:`StringToResu`.
|
@@ -28,33 +28,36 @@ converted to the actual modifyer by :code:`StringToResu`.
:local:
:depth: 2
-
+
|
|
----------------------------------
-
+
Usage
----------------------------------
|
Assuming one has computed a residual :code:`resu` at some time setp :code:`n`, and likes to modify it by the means
-of some limiter or filter/streamlines dissipator. Then, in the settings file (see the :ref:`UserManual`) indicate
-for the amendements the composition of modifications in a natural way:
+of some limiter or filter/streamlines dissipator. Then, in the settings file (see the :ref:`UserManual`) indicate
+for the amendments the composition of modifications in a natural way:
.. code::
-
+
# No limiter
0
-
+
# A single modification
Limiter_$Limit
-
+
+ # A single modification that should be added to the original residual
+ Original+Filtering_$Filter
+
# A single modification with argument to give to the mollifier
Limiter_$Limit<listof,parameters>
-
+
# Or a combination
Limiter_$Limit+Filtering_$Filter<some,parameters,topass>
@@ -69,36 +72,36 @@ to the original residuals (those are the terms not followed by parenthesis).
.. code::
- Limiter_3<2>(Limiter_1<1.2>(Filtering_2)+Filtering_1<p>(Limiter_1<2.2>))*Limiter_2
+ Limiter_3<2>(Limiter_1<1.2>(Original+Filtering_2)+Filtering_1<p>(Limiter_1<2.2>))*Limiter_2
+
-
|
.. warning::
- There should be NO spurious spaces in the definition of the modification strings.
- The list of arguments to pass to the mollifier should be given ordered between the symbols :code:`<>` according to the targetted routine. They will be passed to the related function as strings.
-
+
|
|
|
----------------------------------
-
-Modifiers sequence generation
+
+Modifiers sequence generation
----------------------------------
|
.. automodule:: LSRD.SourceCode._Solver.SpatialModifiers.SpatialAmendements
:members:
- :member-order: bysource
+ :member-order: bysource
|
|
----------------------------------
-
+
Available limiters
----------------------------------
@@ -106,8 +109,8 @@ Available limiters
.. automodule:: LSRD.SourceCode._Solver.SpatialModifiers.Limiter_Psi
:members:
- :member-order: bysource
-
+ :member-order: bysource
+
|
|
@@ -122,9 +125,8 @@ Available filters
.. automodule:: LSRD.SourceCode._Solver.SpatialModifiers.Filtering_Streamline
:members:
- :member-order: bysource
-
+ :member-order: bysource
+
|
|
-
diff --git a/_doc/source/index.rst b/_doc/source/index.rst
index eaef8ce..1505121 100644
--- a/_doc/source/index.rst
+++ b/_doc/source/index.rst
@@ -6,10 +6,10 @@
<br />
-
-
-LSRD: Level Set Residual Distribution's documentation
+
+
+LSRD: Level Set Residual Distribution's documentation
=====================================================
|
@@ -23,9 +23,9 @@ All remarks, suggestions, comments, bugshooting and complains goes to elise(poin
.. rubric:: Documentation
-+------------------------------+--------------------------------------------+
-| :ref:`UserManual` | :ref:`TechnicalManual` |
-+------------------------------+--------------------------------------------+
++------------------------------+--------------------------------------------+
+| :ref:`UserManual` | :ref:`TechnicalManual` |
++------------------------------+--------------------------------------------+
|
@@ -54,7 +54,7 @@ A problem should always be defined by:
* the initial spatial distribution of the fluids coexisting in the computational domain
* the type of boundary conditions according to their label in the associated mesh
|
-
+
**2. Compute a numerical solution for the designed problem**
Once a problem is designed, the code allows to approximate a solution by choosing a solver and its properties (see the :ref:`UserManual` for more information).
@@ -67,26 +67,26 @@ Directly accessible from the solver settings's file, the user can:
- select of the type of the fluid's spotter and its associated solver (the fluid spotter is defined by its nature (e.g. level set), a spatial iterator (e.g. CG with psi-limiter), a redistancing process (e.g. HJB in pseudo-time), and a fluid's selection routine at degrees of freedom and quadrature points.)
|
-
+
**3. Post-process the solution**
-
+
Once a solution has been exported, the user can run through a semi-automated script post-processing tasks.
In particular, one can:
- plot the subdomains's evolution, either as static snapshots or as an html animation, in 2D or 3D
- plot the solution for each variable, either as static snapshots or as an html animation, in 2D or 3D
- plot the solution for each fluid's spotter, either as static snapshots or as an html animation, in 2D or 3D
- - export the solution as a binary file (with pickles), raw in a text-file or in a vtk format
-
+ - export the solution as a binary file (with pickles), raw in a text-file or in a vtk format
+
|
.. note::
- All the spatial schemes are recast into the Residual Distribution framework.
-
+
- The mass matrix is lumped
-
+
- Single-phase problems can also be studied in this framework.
-
+
|
|
|
@@ -121,7 +121,7 @@ Spatial schemes:
+----------------------------------+-----------------+--------------------+---------------------------------------------+
| GC + Lax-Friedrichs flux (LFX) | <=2 (bf) | Dual | Limiter_Psi |
| | | | |
-| | Reducing to 1 | | |
+| | Reducing to 1 | | |
+----------------------------------+-----------------+--------------------+---------------------------------------------+
|
@@ -162,7 +162,7 @@ Code Requirements
Basic requirement
^^^^^^^^^^^^^^^^^
-An installation of :code:`python3 > 3.5` and enough free space in the running folder to accept the export of the binary files of the solution.
+An installation of :code:`python3 > 3.5` and enough free space in the running folder to accept the export of the binary files of the solution.
(The export can be reduced/minimized/disabled manually in the file SourceCode/Solve.py).
A-priori it should be runable on any OS, though it has not been tested on W$ and Mac systems.
@@ -171,7 +171,7 @@ Just add/define your problem following the user manual and run the main file it
|
-Mandatory libraries furnished with the python package's distribution
+Mandatory libraries furnished with the python package's distribution
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+--------------------+---------------------+-----------------------+
@@ -196,7 +196,7 @@ Mandatory libraries furnished with the python package's distribution
|
-Mandatory external libraries
+Mandatory external libraries
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+--------------------+-------------------+-----------------------+
@@ -213,7 +213,7 @@ Mandatory external libraries
|
-Optional external libraries
+Optional external libraries
^^^^^^^^^^^^^^^^^^^^^^^^^^^
+----------------+-------------------+------------------------------------------------+
@@ -223,15 +223,23 @@ Optional external libraries
| | | possibly redistancing in LS, |
| | | plotmesh in primal mesh |
+----------------+-------------------+------------------------------------------------+
-| plotly | 0+untagged. | Plotly plots |
-| | 25.g3b07c4e | |
-+----------------+-------------------+------------------------------------------------+
-| scipy (▲▲) | 1.4.1 | Matplotlib plots, |
+| plotly | 0+untagged. | Plotly plots |
+| | 25.g3b07c4e | |
++----------------+-------------------+------------------------------------------------+
+| scipy (▲▲) | 1.4.1 | Matplotlib plots, |
| | | Plotly plots |
-+----------------+-------------------+------------------------------------------------+
++----------------+-------------------+------------------------------------------------+
| meshio | 4.0.15 | Export of the solution |
| | | in a VTK format |
-+----------------+-------------------+------------------------------------------------+
++----------------+-------------------+------------------------------------------------+
*(▲▲) If no wish to disable the associated feature, the dependency can be dropped upon significant modification of the code or code's specific routine removal*
+|
+
+Known bugs
+----------
+
+Few minor bugs have been spot and not (yet) corrected as non-prioritary. They may include:
+
+ - The plotting routine with Plotly applied to a constant Level Set (1Fluid case) returns an empty figure
--
GitLab