Commit c16ab640 authored by Élise Le Mélédo's avatar Élise Le Mélédo
Browse files

Cleaned code and added extensive documentation

parent a62bfbe9
find . -name "*.pyc" -exec rm -f {} \;
find . -type d -name __pycache__ -exec rm -r {} \+
"""
################################################################################
# Module that defines a Vertex centred finite volume scheme #
################################################################################
Contains (in scrolling order):
-> User defined data reader tools
-> Structure definition
-> Initialisation routines
Usage:
Note:
-----------------------------------------------------------------------------"""
# ==============================================================================
# Preliminary imports
# ==============================================================================
import numpy as np
import copy
# ==============================================================================
# Parameter (and reader) class of scheme properties and runtime informations
# ==============================================================================
# Name of the module to be externally accessible
SchemeName = "Finite Volumes -- Vertex centred with natural flux"
# Mapping to the mesh that should be loaded
class AssociatedMeshType():
""" 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
MeshType, string, the mesh type according to fenics's classification
MeshOrder, integer, the mesh order
Note: This is only a product-type class: no method is available
--------------------------------------------------------------------------------------"""
#### Automatic initialisation routine ####
def __init__(self, VariationalOrder = 0):
""" Automatic initialisation routine defining the type of mesh
to be loaded for this solver."""
self.MeshType = "DG"
self.MeshOrder = 0 # Or dependent on VariationalOrder
self.ControlVolumes = "Primal"
# ==============================================================================
# Class furnishing all the scheme tools
# ==============================================================================
class Scheme():
""" Class furnishing all the schemes tools and registers the scheme's wished
parameters. The main iteration is accessible with the method "Iteration".
--------------------------------------------------------------------------------------"""
#### Automatic initialisation routine ####
def __init__(self, Problem, Mesh, Params = []):
""" Automatic initialisation routine. Just registers the user defined scheme parameters
within the instance"""
self.Problem = Problem
self.Params = Params
self.Mesh = Mesh
#### Emulates a flux computation as done in the iteration routine ####
def ComputeFlux(self, FluidFlag, U, LSValues):
""" Emulates a flux computation as done in the iteration routine in order
to be used in the DeC subtimestep evaluation.
Input:
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)
Output:
fluxes, numpy float (multiD)array, fluxes in the format requires by the Iteration routine
---------------------------------------------------- """
# -------------------------- Initialisations -------------------------------------------------
fluxes = np.zeros((self.Mesh.NbInnerElements, 3, 2, self.Problem.GoverningEquations.NbVariables,2))
# ------------- Flux computation at the control's volume interface ----------------------------
for i in range(self.Mesh.NbInnerElements):
# Retrieving the information on the local degree of freedom
dele = self.Mesh.ElementsInnerDofs[i][0]
ele = self.Mesh.InnerElements[i]
Eele = self.Mesh.InnerElements[i] + [self.Mesh.InnerElements[i][0]]
# Looping on each edge
for f in range(len(ele)):
# Getting the variational info of the immediate neighbor
ala = self.Mesh.Neighborhoods[i][f]
dala = self.Mesh.ElementsInnerDofs[ala][0]
# Reconstructing the solution's value at the considered point (usually quadrature point)
Flags = FluidFlag[np.array([dele,dala])]
Var = np.array([U[:,dele], U[:,dala]])
# Compute the input exchange from the neighboring element at the edge midpoint
EdgeMidPoint = np.array(2*[0.5*np.sum(self.Mesh.CartesianPoints[Eele[f:f+2],:], axis=0)])
fluxes[i,f,:,:,:] = self.Problem.GoverningEquations.Flux(Var, EdgeMidPoint, Flags)
# Returning the full array
return(fluxes)
#### Main routine, definining the scheme ####
def Iteration(self, Solution, fluxes, 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]
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)
---------------------------------------------------------------------"""
# -------- Initialising the mollification vector (residuals) -------------
Resu = np.zeros(np.shape(Solution.Sol))
# ------- Getting the residuals for each element -------------------------
for i in range(self.Mesh.NbInnerElements):
# Retrieving the information on the local degree of freedom
dele = self.Mesh.ElementsInnerDofs[i][0]
ele = self.Mesh.InnerElements[i]
# Looping on each edge
for f in range(len(ele)):
ala = self.Mesh.Neighborhoods[i][f]
dala = self.Mesh.ElementsInnerDofs[ala][0]
# Reconstructing the solution's value at the considered point (usually quadrature point)
Flags = FluidFlag[np.array([dele,dala])]
Var = np.array([U[:,dele], U[:,dala]])
# Get the physical properties of the considered edge
L = self.Mesh.ElementsBoundaryWiseWidth[i][f]
n = self.Mesh.ElementsBoundaryWiseNormals[i][f,:]
EdgeMidPoint = np.array(2*[0.5*(np.sum(self.Mesh.CartesianPoints[Edge, :], axis=0))])
# Retrieving the flux computed through DeC
MidFluxes = fluxes[i,f,:,:,:]
MidFluxes[:,:,1] = -MidFluxes[:,:,1]
sprctrad = self.Problem.GoverningEquations.SpectralRadius(Var.T, Flags, np.array(2*[n]), EdgeMidPoint)
# Compute the pondered contribution
FluxN = np.dot(0.5*np.sum(MidFluxes, axis=2).T,n)
quant = L*FluxN - L*np.max(sprctrad)*(U[:,dele]-U[:,dala])
# Distribute it equally to each submit
Resu[:,dele] = Resu[:,dele] + (quant)
#----- Modifying the residual on each Dof by addin the DeC increment --------
# Adapting the residual to Dec
Resu += du[:,:]/self.Mesh.Lumping
# Returning the result
return(Resu)
#### Routine that maps the values from the Dofs to the Physical vertices of the mesh ####
def ReconstructSolutionAtVertices(self, Solution):
""" Routine that maps the values from the Dofs to the Physical vertices of the mesh
of the FluidSpotter, solution and flags.
Input: SolverSpecs, string, the name of the setting file
Output: None, fills direclty the Parameters data structure
Note: This is only for later plotting purposes.
----------------------------------------------------------------------------------- """
# -------------- Initialising the reconstructed values ----------------------------
Solution.RFluidFlag = np.zeros(np.shape(Solution.FluidFlag))
Solution.RLSValues = np.zeros(np.shape(Solution.LSValues))
Solution.RSol = np.zeros(np.shape(Solution.Sol))
# ------------- 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):
# Retrieve the catesian coordinates of the considered point
xy = self.Mesh.CartesianPoints[i,:]
# Stupid reconstruction at the closest value for Dofs #TODO better
norm = np.linalg.norm(xy-xyDofs, 2, axis=1)
index = np.argsort(norm)[0]
# Fill the reconstructed solution upon the found index
Solution.RLSValues[:,i] = Solution.LSValues[:,index]
Solution.RFluidFlag[i] = Solution.FluidFlag[index]
Solution.RSol[:,i] = Solution.Sol[:,index]
"""
################################################################################
# Module that defines a Vertex centred finite volume scheme #
################################################################################
Contains (in scrolling order):
-> User defined data reader tools
-> Structure definition
-> Initialisation routines
Usage:
Note:
-----------------------------------------------------------------------------"""
# ==============================================================================
# Preliminary imports
# ==============================================================================
import numpy as np
# ==============================================================================
# Parameter (and reader) class of scheme properties and runtime informations
# ==============================================================================
# Name of the module to be externally accessible
SchemeName = "Finite Volumes -- Cell centred with LFX flux"
# Mapping to the mesh that should be loaded
class AssociatedMeshType():
""" 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
MeshType, string, the mesh type according to fenics's classification
MeshOrder, integer, the mesh order
Note: This is only a product-type class: no method is available
--------------------------------------------------------------------------------------"""
#### Automatic initialisation routine ####
def __init__(self, VariationalOrder = 0):
""" Automatic initialisation routine defining the type of mesh
to be loaded for this solver."""
self.MeshType = "CG"
self.MeshOrder = 1 # Or dependent on VariationalOrder
self.ControlVolumes = "Dual"
# ==============================================================================
# Class furnishing all the scheme tools
# ==============================================================================
class Scheme():
""" Class furnishing all the schemes tools and registers the scheme's wished
parameters. The main iteration is accessible with the method "Iteration".
--------------------------------------------------------------------------------------"""
#### Automatic initialisation routine ####
def __init__(self, Problem, Mesh, Params = []):
""" Automatic initialisation routine. Just registers the user defined scheme parameters
within the instance"""
self.Problem = Problem
self.Params = Params
self.Mesh = Mesh
def ComputeFlux(self, FluidFlag, U, LSValues):
fluxes = np.zeros((self.Mesh.NbInnerElements, 3, 2, self.Problem.GoverningEquations.NbVariables,2))
# Initialising the mollification vector (residuals)
Res = np.zeros(np.shape(U))
# Looping on each element
for i in range(self.Mesh.NbInnerElements):
dele = self.Mesh.ElementsBoundaryDofs[i]
ele = self.Mesh.InnerElements[i]
Dele = np.array(dele + [dele[0]], dtype=int)
Ele = np.array(ele + [ele[0]], dtype=int)
# Get the average state of the considered element
Var1 = np.mean(U[:,dele], axis = 1)
res = np.zeros((1,self.Problem.GoverningEquations.NbVariables))
# Looping on each edge
for f in range(len(ele)):
# Get the physical vertices of the considered edge
n = self.Mesh.ElementsBoundaryWiseNormals[i][f]
L = self.Mesh.ElementsBoundaryWiseWidth[i][f]
Nei = self.Mesh.Neighborhoods[i][f]
Edge = Ele[f:f+2]
# Get the variational information at each Dof lying on the edge
DEdge = Dele[f:f+2]
Flags = FluidFlag[DEdge]
# Get the average state of the neighboring element
Dele2 = np.array(self.Mesh.ElementsBoundaryDofs[Nei])
Var2 = np.mean(U[:,Dele2], axis=1)
# Compute the input exchange from the neighboring element at the edge midpoint
EdgeMidPoint = np.array(2*[0.5*(np.sum(self.Mesh.CartesianPoints[Edge, :], axis=0))])
fluxes[i,f,:,:] = self.Problem.GoverningEquations.Flux(np.array([Var1, Var2]).T, EdgeMidPoint, Flags)
return(fluxes)
#### Main routine, defininin the scheme ####
def Iteration(self, Solution, fflx, du=0, dt=1):
""" Main iteration of the scheme, implementing the most stupid scheme
you can think of.
Args:
SolverSpecs (string): the name of the setting file
Returns:
None : fills direclty the Parameters data structure
---------------------------------------------------------------------"""
# Initialising the mollification vector (residuals)
Res = np.zeros(np.shape(Solution.Sol))
# Looping on each element
for i in range(self.Mesh.NbInnerElements):
dele = self.Mesh.ElementsBoundaryDofs[i]
ele = self.Mesh.InnerElements[i]
Dele = np.array(dele + [dele[0]], dtype=int)
Ele = np.array(ele + [ele[0]], dtype=int)
# Get the average state of the considered element
Var1 = np.mean(Solution.Sol[:,dele], axis = 1)
res = np.zeros((1,self.Problem.GoverningEquations.NbVariables))
# Looping on each edge
for f in range(len(ele)):
# Get the physical vertices of the considered edge
L = self.Mesh.ElementsBoundaryWiseWidth[i][f]
Nei = self.Mesh.Neighborhoods[i][f]
Edge = Ele[f:f+2]
# Get the variational information at each Dof lying on the edge
DEdge = Dele[f:f+2]
Flags = Solution.FluidFlag[DEdge]
# Get the average state of the neighboring element
Dele2 = np.array(self.Mesh.ElementsBoundaryDofs[Nei])
Var2 = np.mean(Solution.Sol[:,Dele2], axis=1)
# Compute the input exchange from the neighboring element at the edge midpoint
EdgeMidPoint = np.array(2*[0.5*(np.sum(self.Mesh.CartesianPoints[Edge, :], axis=0))])
MidFluxes = fflx[i,f,:,:]
sprctrad = self.Problem.GoverningEquations.SpectralRadius(np.array([Var1, Var2]).T, Flags, np.array(2*[n]), EdgeMidPoint)
# Compute the pondered contribution
MidFluxes[:,:,1] = -MidFluxes[:,:,1]
avrg = 0.5*np.sum(MidFluxes, axis=2)
FluxN = np.dot(avrg.T,n.T)
# Add the contribution
res[:] += L*(FluxN - np.max(sprctrad)*(Var1[:]-Var2[:]))
# Compute the average of du over the cell and compute its associated
# contribution over the element
dubar = (np.sum(du[:,dele], axis=1)/len(ele))*self.Mesh.InnerElementsVolume[i]
# Distribute it equally to each submit
#Res[:,dele] += (dt*res[:].T+dubar)*self.Mesh.DualAreas[i][:].T/self.Mesh.InnerElementsVolume[i]
Res[:,dele] += (dt*res[:].T+dubar)/3
# Updating each nodal value (here representing the average) by the residual
# pondered by the control volume's area
Res[:, :] = Res[:, :]
return(Res)
#### 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.
Input: SolverSpecs, string, the name of the setting file
Output: None, fills direclty the Parameters data structure
Note: This is only for later plotting purposes.
----------------------------------------------------------------------------------- """
# Stupid reconstruction at the closest value for Dofs #TODO better
Solution.RLSValues = np.zeros(np.shape(Solution.LSValues))
Solution.RSol = np.zeros(np.shape(Solution.Sol))
Solution.RFluidFlag = np.zeros(np.shape(Solution.FluidFlag))
for i in range(self.Mesh.NbMeshPoints):
xy = self.Mesh.CartesianPoints[i,:]
xyDofs = self.Mesh.DofsXY[:,:]
norm = np.linalg.norm(xy-xyDofs,2, axis=1)
index = np.argsort(norm)[0]
Solution.RLSValues[:,i] = Solution.LSValues[:,index]
Solution.RFluidFlag[i] = Solution.FluidFlag[index]
Solution.RSol[:,i] = Solution.Sol[:,index]
"""
################################################################################
# Module that defines a Vertex centred finite volume scheme #
################################################################################
Contains (in scrolling order):
-> User defined data reader tools
-> Structure definition
-> Initialisation routines
Usage:
Note:
-----------------------------------------------------------------------------"""
# ==============================================================================
# Preliminary imports
# ==============================================================================
import numpy as np
# ==============================================================================
# Parameter (and reader) class of scheme properties and runtime informations
# ==============================================================================
# Name of the module to be externally accessible
SchemeName = "Finite Volumes -- Vertex centred with natural flux"
# Mapping to the mesh that should be loaded
class AssociatedMeshType():
""" 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
MeshType, string, the mesh type according to fenics's classification
MeshOrder, integer, the mesh order
Note: This is only a product-type class: no method is available
--------------------------------------------------------------------------------------"""
#### Automatic initialisation routine ####
def __init__(self, VariationalOrder = 0):
""" Automatic initialisation routine defining the type of mesh
to be loaded for this solver."""
self.MeshType = "DG"
self.MeshOrder = 0 # Or dependent on VariationalOrder
self.ControlVolumes = "Primal"
# ==============================================================================
# Class furnishing all the scheme tools
# ==============================================================================
class Scheme():
""" Class furnishing all the schemes tools and registers the scheme's wished
parameters. The main iteration is accessible with the method "Iteration".
--------------------------------------------------------------------------------------"""
#### Automatic initialisation routine ####
def __init__(self, Problem, Mesh, Params = []):
""" Automatic initialisation routine. Just registers the user defined scheme parameters
within the instance"""
self.Problem = Problem
self.Params = Params
self.Mesh = Mesh
#### Emulates a flux computation as done in the iteration routine ####
def ComputeFlux(self, FluidFlag, U, LSValues):
""" Emulates a flux computation as done in the iteration routine in order
to be used in the DeC subtimestep evaluation.
Input:
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)
Output:
fluxes, numpy float (multiD)array, fluxes in the format requires by the Iteration routine
---------------------------------------------------- """
# -------------------------- Initialisations -------------------------------------------------
fluxes = np.zeros((self.Mesh.NbInnerElements, 3, 2, self.Problem.GoverningEquations.NbVariables,1))
# ------------- Flux computation at the control's volume interface ----------------------------
for i in range(self.Mesh.NbInnerElements):
# Retrieving the information on the local degree of freedom
dele = self.Mesh.ElementsInnerDofs[i][0]
ele = self.Mesh.InnerElements[i]
Eele = self.Mesh.InnerElements[i] + [self.Mesh.InnerElements[i][0]]
# Looping on each edge
for f in range(len(ele)):
# Reconstructing the solution's value at the considered point (usually quadrature point)
Flags = np.array([FluidFlag[dele]])
Var = np.array([U[:,dele]])
# Compute the input exchange from the neighboring element at the edge midpoint
EdgeMidPoint = np.array([0.5*np.sum(self.Mesh.CartesianPoints[Eele[f:f+2],:], axis=0)])
fluxes[i,f,:,:,:] = self.Problem.GoverningEquations.Flux(Var, EdgeMidPoint, Flags)
# Returning the full array
return(fluxes)
#### Main routine, definining the scheme ####
def Iteration(self, Solution, fluxes, 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]
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)
---------------------------------------------------------------------"""
# -------- Initialising the mollification vector (residuals) -------------
Resu = np.zeros(np.shape(Solution.Sol))
# ------- Getting the residuals for each element -------------------------
for i in range(self.Mesh.NbInnerElements):
# Retrieving the information on the local degree of freedom
dele = self.Mesh.ElementsInnerDofs[i][0]
ele = self.Mesh.InnerElements[i]
# Looping on each edge
for f in range(len(ele)):
# Get the physical properties of the considered edge
L = self.Mesh.ElementsBoundaryWiseWidth[i][f]
n = self.Mesh.ElementsBoundaryWiseNormals[i][f,:]
# Retrieving the flux computed through DeC
MidFluxes = fluxes[i,f,:,:,0]
# Compute the pondered contribution
FluxN = np.dot(MidFluxes.T,n)
quant = L*FluxN
# Distribute it equally to each submit
Resu[:,dele] += dt*quant
#----- Modifying the residual on each Dof by addin the DeC increment --------
# Adapting the residual to Dec
Resu += du[:,:]/self.Mesh.Lumping
# Returning the result
return(Resu)
#### Routine that maps the values from the Dofs to the Physical vertices of the mesh ####
def ReconstructSolutionAtVertices(self, Solution):
""" Routine that maps the values from the Dofs to the Physical vertices of the mesh
of the FluidSpotter, solution and flags.
Input: SolverSpecs, string, the name of the setting file
Output: None, fills direclty the Parameters 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):
# Retrieve the catesian coordinates of the considered point
xy = self.Mesh.CartesianPoints[i,:]
# Stupid reconstruction at the closest value for Dofs #TODO better
norm = np.linalg.norm(xy-xyDofs,2, axis=1)
index = np.argsort(norm)[0]
# Fill the reconstructed solution upon the found index
Solution.RLSValues[:,i] = Solution.LSValues[:,index]
Solution.RFluidFlag[i] = Solution.FluidFlag[index]
Solution.RSol[:,i] = Solution.Sol[:,index]