kawin.diffusion.mesh
kawin.kawin.diffusion.mesh.MeshBase
def noChangeAtNode(Ds):
Default function for at node diffusivity
def arithmeticMean(Ds):
Arithmetic mean - D_avg = 1/2*(D_i + D_j)
Ds should be in the shape of (m x N x y)
m - number of items to average over
N - number of nodes
y - number of responses
Average is taken along m
def geometricMean(Ds):
Geometric mean - D_avg = sqrt(D_i * D_j)
Ds should be in the shape of (m x N x y)
m - number of items to average over
N - number of nodes
y - number of responses
Average is taken along m
Note: this is only valid when D_i and D_j are the same sign
def logMean(Ds):
Log mean - D_avg = exp(1/2*(log(D_i)+log(D_j)))
Ds should be in the shape of (m x N x y)
m - number of items to average over
N - number of nodes
y - number of responses
Average is taken along m
Note: this is only valid for positive D
def harmonicMean(Ds):
Harmonic mean - D_avg = (1/2 * (D_i^-1 + D_j^-1))^-1
Ds should be in the shape of (m x N x y)
m - number of items to average over
N - number of nodes
y - number of responses
Average is taken along m
Note: in the limit of D_i or D_j -> 0, the mean is 0. In practice, this is undefined, so we
convert all inf/nan to 0
def _getResponseVariables(responses: Union[int, list[str]]) -> tuple[int, list[str]]:
This returns the number of response variables and names for each variable
Parameters
----------
responses: int | list[str]
If int, then this is the number of responses and the names will be R{i}
If list[str], then theses are the response names and the number of responses will be the list length
def _getResponseIndex(responseVar: Union[int, str], responses):
Given the response variable, return the corresponding index
def _formatSpatial(z):
BoundaryCondition(ABC)
def BoundaryCondition.setInitialResponse(self, y):
def BoundaryCondition.adjustFluxes(self, fluxes):
def BoundaryCondition.adjustdXdt(self, dXdt):
ProfileFunction(Protocol)
ProfileBuilder
class ProfileBuilder(self, steps = []):
Stores build steps to construct a response profile in a mesh
Note that all steps are additive to the response profile
TODO: consider removing this since it doesn't really do anything other
than store a list that the user can create themselves
Parameters
----------
steps: list[tuple(callable, int|str | list[int|str])] (optional)
List of build steps to initialize
def ProfileBuilder.addBuildStep(self, func: ProfileFunction, responseVar: int|str|list[int|str] = 0):
Adds a build step to construct the initial profile
Parameters
----------
func: ProfileFunction
Function that takes in a set of coordinates and returns the profile values
at coordinates
responseVar: int or str or list of int or str (optional)
Response variable(s) that func will return
This does not have to match exactly to the mesh variables, but only needs
to be a subset
def ProfileBuilder.clearBuildSteps(self):
Removes all build steps
ConstantProfile(ProfileFunction)
class ConstantProfile(self, value):
Constant value across profile
def ConstantProfile.call(self, z):
DiracDeltaProfile(ProfileFunction)
class DiracDeltaProfile(self, z, value):
Zero at all values except at z
Note: if a mesh does not have a node exactly at z, it will choose the closest node
def DiracDeltaProfile.call(self, z):
GaussianProfile(ProfileFunction)
class GaussianProfile(self, z, sigma, maxValue):
Gaussian function centered around z with standard deviation of sigma
Scaled to maxValue
def GaussianProfile.call(self, z):
BoundedRectangleProfile(ProfileFunction)
class BoundedRectangleProfile(self, lowerZ, upperZ, innerValue, outerValue = 0):
Defines rectangle with lower and upper corners
innerValue corresponds to value inside the rectangle while
outerValue corresponds to value outside the rectangle (defaults to 0)
def BoundedRectangleProfile.call(self, z):
BoundedEllipseProfile(ProfileFunction)
class BoundedEllipseProfile(self, z, r, innerValue, outerValue = 0):
Defines ellipse centered at z with radii r
innerValue corresponds to value inside the rectangle while
outerValue corresponds to value outside the rectangle (defaults to 0)
def BoundedEllipseProfile.call(self, z):
AbstractMesh (ABC)
class AbstractMesh (self, responses):
Abstract mesh class that defines basic methods to be used in a diffusion model
Parameters
----------
responses: int | list[str]
If int, then this is the number of responses and the names will be R{i}
If list[str], then theses are the response names and the number of responses will be the list length
Attributes
----------
numResponses: int
Number of response variables
responses: list[str]
Names of response variables
N: int
Total number of nodes
dims: int
Number of dimensions in z
y: np.ndarray
Initial values of response variables
This is intended not to change in a model. Time evolution
of response variables should be stored in MeshData
z: np.ndarray
Values of spatial coordinates corresponding to y
dz: float
In general, the smallest distance between two nodes
Default assumptions in this mesh:
y has shape [N,e] - e is number of responses
z has shape [N,d] - d is number of dimensiosn
def AbstractMesh .setResponseProfile(self, profileBuilder: ProfileBuilder, boundaryConditions: BoundaryCondition = None):
Creates initial profile on mesh based off a series of build steps
def AbstractMesh .computeFluxes(self, pairs):
def AbstractMesh .computedXdt(self, pairs: list[DiffusionPair]):
Given list of diffusivity and response pairs, compute dX/dt from diffusion fluxes
Parameters
----------
pairs: list[DiffusionPair]
Tuple will contain (diffusivity, response, averaging function, at node function)
Note that the diffusivity and response will be in the form of [N,e] (corresponding
to getDiffusivityCoordinates and getResponseCoordinates), which may not be the
shape of self.y. This function is responsible for accounting for this difference
and the shape of dxdt and y should be the same
Returns
-------
dXdt: np.ndarray
Must be same shape as y
def AbstractMesh .flattenResponse(self, y, numResponses = None):
Converts y [internal shape, ...] to yFlat [N, e, ...]
Parameters
----------
y: np.ndarray
Input y response in internal mesh shape with 1 dimension corresponding to responses
numResponses: int (optional)
Number of responses in y, defaults to total number of responses in mesh
Returns
-------
yFlat: np.ndarray
Shape of [N, e] with N being the total number of nodes and e being the number of responses
def AbstractMesh .flattenSpatial(self, z):
Converts z [internal shape, ...] to zFlat [N, d, ...]
Parameters
----------
z: np.ndarray
Input z response in internal mesh shape with 1 dimension
corresponding to the number of spatial dimensions
Returns
-------
zFlat: np.ndarray
Shape of [N, d] with N being the total number of nodes and d being the number of spatial dimensions
def AbstractMesh .unflattenResponse(self, yFlat, numResponses = None):
Converts yFlat [N, e, ...] to y [internal shape, ...]
Parameters
----------
yFlat: np.ndarray
Shape of [N, e] with N being the total number of nodes and e being the number of responses
numResponses: int (optional)
Number of responses in y, defaults to total number of responses in mesh
Returns
-------
y: np.ndarray
Input y response in internal mesh shape with 1 dimension corresponding to responses
def AbstractMesh .unflattenSpatial(self, zFlat):
Converts zFlat [N, d, ...] to z [internal shape, ...]
Parameters
----------
zFlat: np.ndarray
Shape of [N, d] with N being the total number of nodes and d being the number of spatial dimensions
Returns
-------
z: np.ndarray
Input z response in internal mesh shape with 1 dimension
corresponding to the number of spatial dimensions
def AbstractMesh .getResponseCoordinates(self, y):
Returns y as [N,e] and z as [N,D] arrays to compute response terms
def AbstractMesh .getDiffusivityCoordinates(self, y):
Returns y as [N,e] and z as [N,D] arrays to compute diffusivity terms
We decouple this from the response coordinates to allow for different ways
of computing D^(1/2)
a) By computing D^0 and D^1 and averaging
b) By computing D^(1/2) at (y^1/2, z^1/2)
By default, this will use option a since its an easier implementation, but the user can override this to use option b
def AbstractMesh ._diffusiveFlux(self, D, rHigh, rLow, dz):
Flux = -D (r1 - r0) / dz
FiniteVolumeGrid(AbstractMesh)
class FiniteVolumeGrid(self, responses, zlims, Ns, dims):
Functions for the finite volume method on a rectangular grid
def FiniteVolumeGrid.defineZCoordinates(self):
def FiniteVolumeGrid.flattenResponse(self, y, numResponses = None):
Converts y [n x m x ..., e, ...] to yFlat [N, e, ...]
def FiniteVolumeGrid.flattenSpatial(self, z):
Converts z [n x m x ..., d, ...] to zFlat [N, d, ...]
def FiniteVolumeGrid.unflattenResponse(self, yFlat, numResponses = None):
Converts yFlat [N, e, ...] to y [n x m x ..., e, ...]
def FiniteVolumeGrid.unflattenSpatial(self, zFlat):
Converts zFlat [N, e, ...] to z [n x m x ..., e, ...]
MeshData
class MeshData(self, mesh: AbstractMesh, record: bool | int = False):
Stores time and response variables for a mesh
Parameters
----------
mesh : AbstractMesh
Mesh to take dimensions of response variables from
record : bool | int
Record interval or whether to record
If False, only the current state of the mesh will be stored
If int, then every n iterations will be stored
def MeshData.reset(self):
Resets arrays
def MeshData.record(self, time, y, force = False):
Stores current state of time and response variables
Response variables should be in the flattened state
This is the same format as used in GenericModel
But for the mesh, we must call flattenResponse
def MeshData.setResponseAtN(self, y, N = -1):
Sets response profile at index N
If N is not supplied, then this sets at the current index
def MeshData.finalize(self):
Removes extra padding
def MeshData.y(self, time = None):
Returns reponse variable at time
If recording is disabled, then this will return the current state
kawin.kawin.diffusion.mesh.FVM1D
PeriodicBoundary1D(BoundaryCondition)
This doesn't do anything, but rather is used in the Mesh to compute the fluxes accordingly
Unfortunately, we can't just wrap the fluxes around and instead have to wrap the response variable
around. Which means that we would have to have knowledge on how to compute the flux (which is done in the mesh)
def PeriodicBoundary1D.setInitialResponse(self, y):
def PeriodicBoundary1D.adjustFluxes(self, fluxes):
def PeriodicBoundary1D.adjustdXdt(self, dXdt):
MixedBoundary1D(BoundaryCondition)
class MixedBoundary1D(self, responses: Union[int, list[str]]):
Currently, left and right boundary conditions both are defined here
We'll need a way to define boundary conditions at any edge, which will
especially be the case if we add a 2D mesh
Parameters
----------
responses: int | list[str]
If int, then this is the number of responses and the names will be R{i}
If list[str], then theses are the response names and the number of responses will be the list length
def MixedBoundary1D._setBC(self, responseVar, bcType, value, bcTypeArray, bcValueArray):
Sets boundary condition to given side (left or right)
def MixedBoundary1D.setLBC(self, responseVar, bcType, value):
Left boundary condition
Parameters
----------
responseVar: int | str
Response variable name or index
bcType: str
For flux/Neumann condition, it could be [flux, neumman, MixedBoundary1D.NEUMANN]
For composition/dirichlet condition, it could be [composition, dirichlet, MixedBoundary1D.DIRICHLET]
value: float
Value of boundary condition
def MixedBoundary1D.setRBC(self, responseVar, bcType, value):
Left boundary condition
Parameters
----------
responseVar: int | str
Response variable name or index
bcType: str
For flux/Neumann condition, it could be [flux, neumman, MixedBoundary1D.NEUMANN]
For composition/dirichlet condition, it could be [composition, dirichlet, MixedBoundary1D.DIRICHLET]
value: float
Value of boundary condition
def MixedBoundary1D.setInitialResponse(self, y):
Sets dirichlet boundary conditions onto y
def MixedBoundary1D.adjustFluxes(self, fluxes):
Adjust flux values using boundary conditions
For Neumann boundary conditions, this overrides the flux with the BC value
def MixedBoundary1D.adjustdXdt(self, dXdt):
Adjust dx/dt using boundary conditions
For Dirichlet boundary conditions, this sets dx/dt to 0
StepProfile1D(ProfileFunction)
class StepProfile1D(self, z, leftValue, rightValue):
Step function at z
def StepProfile1D.call(self, z):
LinearProfile1D(ProfileFunction)
class LinearProfile1D(self, leftZ, leftValue, rightZ, rightValue, lowerLeftValue = None, upperRightValue = None):
Linear function from (leftZ, leftValue) to (rightZ, rightValue)
lowerLeftValue and upperRightValue corresponds to the value outside the bounds of
leftZ and rightZ. By default they will use leftValue and rightValue
def LinearProfile1D.call(self, z):
ExperimentalProfile1D(ProfileFunction)
class ExperimentalProfile1D(self, z, values, left=None, right=None):
Given experimental values at z, this will interpolate in between
left and right are the values used if the mesh goes beyond z. By default, they
will use the left-most and right-most values
def ExperimentalProfile1D.call(self, z):
FiniteVolumeMidPointCalculator(ABC)
In the finite volume method, we need to compute diffusivity at the
edges of out nodes, but the response variables are always at the cell centers
This class is intended to provide functionality of how to we approach computing
diffusivities at the edges
def FiniteVolumeMidPointCalculator.getDiffusivityCoordinates(self, *args, **kwargs):
def FiniteVolumeMidPointCalculator.getDMid(self, *args, **kwargs):
FVM1DMidpoint(FiniteVolumeMidPointCalculator)
def FVM1DMidpoint.getDiffusivityCoordinates(y, z, zEdge, isPeriodic = False, *args, **kwargs):
This returns the mid point of y and z (which is zEdge)
Note: first element in ymid is second element of zEdge
For periodic conditions, we also include y value between the first and last node
def FVM1DMidpoint.getDMid(D, isPeriodic = False, atNodeFunc = noChangeAtNode, *args, **kwargs):
For periodic conditions, the last D will correspond to between the first and last node
FVM1DEdge(FiniteVolumeMidPointCalculator)
def FVM1DEdge.getDiffusivityCoordinates(y, z, zEdge, *args, **kwargs):
Return y and z since we compute diffusivity at these coordinates
def FVM1DEdge.getDMid(D, isPeriodic = False, avgFunc = arithmeticMean, *args, **kwargs):
For periodic conditions, we average the first and last D
FiniteVolume1D(FiniteVolumeGrid)
class FiniteVolume1D(self, responses, zlim, N, computeMidpoint = False):
1D finite volume mesh
Cell will be a volume with thickness dz
y - value at node center
z - spatial coordinate at node center
zEdge - spatial coordinate at node edge
Parameters
----------
zlim: list[float]
Left and right boundary position of mesh
N: int
Number of cells
responses: int | list[str]
If int, then this is the number of responses and the names will be R{i}
If list[str], then theses are the response names and the number of responses will be the list length
boundaryConditions: BoundaryCondition1D
Boundary conditions on mesh
computeMidpoint: bool
Whether to compute diffusivity at average midpoint composition (True) or
average diffusivities computed on nearby cell centers (False)
def FiniteVolume1D.defineZCoordinates(self):
def FiniteVolume1D.setResponseProfile(self, profileBuilder, boundaryConditions = None):
def FiniteVolume1D.computeAtMidpoint(self, compute: bool):
Sets midpoint calculator
Parameters
----------
compute: bool
If True, diffusivity is computed at average response at node edge
If False, diffusivity is the average diffusivity computed at the neighboring node centers
def FiniteVolume1D._fluxTodXdt(self, fluxes):
Given fluxes, compute dx/dt. This depends on coordinate system
def FiniteVolume1D.computeFluxes(self, pairs: list[DiffusionPair]):
Compute fluxes from (diffusivity, response) pairs on a 1D FVM mesh
def FiniteVolume1D.computedXdt(self, pairs: list[DiffusionPair]):
Computes fluxes, correct fluxes for boundary conditions, compute dx/dt and correct dx/dt for boundary conditions
def FiniteVolume1D.getDiffusivityCoordinates(self, y):
Return y and z to compute diffusivities at
Cartesian1D(FiniteVolume1D)
1D finite volume mesh in cartesian coordinates
Cell will be a volume with thickness dz
y - value at node center
z - spatial coordinate at node center
zEdge - spatial coordinate at node edge
Parameters
----------
zlim: list[float]
Left and right boundary position of mesh
N: int
Number of cells
responses: int | list[str]
If int, then this is the number of responses and the names will be R{i}
If list[str], then theses are the response names and the number of responses will be the list length
boundaryConditions: BoundaryCondition1D
Boundary conditions on mesh
computeMidpoint: bool
Whether to compute diffusivity at average midpoint composition (True) or
average diffusivities computed on nearby cell centers (False)
def Cartesian1D._fluxTodXdt(self, fluxes):
For cartesian: dx/dt = -dJ/dz
Cylindrical1D(FiniteVolume1D)
class Cylindrical1D(self, responses, rlim, N, computeMidpoint = False):
1D finite volume mesh in cylindrical coordinates
Cell will be a volume with thickness dz
y - value at node center
z - spatial coordinate at node center
zEdge - spatial coordinate at node edge
Parameters
----------
rlim: list[float]
Left and right boundary position of mesh in terms of radius
Both min and max radius must be above 0
N: int
Number of cells
responses: int | list[str]
If int, then this is the number of responses and the names will be R{i}
If list[str], then theses are the response names and the number of responses will be the list length
boundaryConditions: BoundaryCondition1D
Boundary conditions on mesh
computeMidpoint: bool
Whether to compute diffusivity at average midpoint composition (True) or
average diffusivities computed on nearby cell centers (False)
def Cylindrical1D.setResponseProfile(self, profileBuilder, boundaryConditions=None):
def Cylindrical1D._fluxTodXdt(self, fluxes):
For cylindrical: dx/dt = -1/r dJ/dz
Spherical1D(FiniteVolume1D)
class Spherical1D(self, responses, rlim, N, computeMidpoint = False):
1D finite volume mesh in spherical coordinates
Cell will be a volume with thickness dz
y - value at node center
z - spatial coordinate at node center
zEdge - spatial coordinate at node edge
Parameters
----------
rlim: list[float]
Left and right boundary position of mesh in terms of radius
Both min and max radius must be above 0
N: int
Number of cells
responses: int | list[str]
If int, then this is the number of responses and the names will be R{i}
If list[str], then theses are the response names and the number of responses will be the list length
boundaryConditions: BoundaryCondition1D
Boundary conditions on mesh
computeMidpoint: bool
Whether to compute diffusivity at average midpoint composition (True) or
average diffusivities computed on nearby cell centers (False)
def Spherical1D.setResponseProfile(self, profileBuilder, boundaryConditions=None):
def Spherical1D._fluxTodXdt(self, fluxes):
For spherical: dx/dt = -1/r^2 dJ/dz
kawin.kawin.diffusion.mesh.FVM2D
Cartesian2D(FiniteVolumeGrid)
class Cartesian2D(self, responses, zx, Nx, zy, Ny):
2D finite volume mesh
y - value at node center (Nx, Ny, dims)
z - spatial coordinate at node center (Nx, Ny, 2)
zCorner - spatial coordinate at node corners (Nx+1, Ny+1, 2)
dz - thickness of node in both dimensions
TODO: boundary conditions not supported yet, so only no flux conditions
Parameters
----------
responses: int | list[str]
If int, then this is the number of responses and the names will be R{i}
If list[str], then theses are the response names and the number of responses will be the list length
zx: list[float]
Left and right boundary position of mesh
Nx: int
Number of cells along x
zy: list[float]
Top and bottom boundary position of mesh
Ny: int
Number of cells along y
def Cartesian2D.defineZCoordinates(self):
def Cartesian2D.computeFluxes(self, pairs: list[DiffusionPair]):
Compute fluxes from (diffusivity, response) pairs on a 2D FVM mesh
def Cartesian2D.computedXdt(self, pairs: list[DiffusionPair]):
dx/dt = -dJx/dz + -dJy/dz