Source code for wolfhece.validation.sensitivity_mesh
"""
Module for Sensitivity Analyses.
This module provides tools and functions to perform sensitivity analyses
on various simulation outputs.
Author: Utashi Ciraane Docile
"""
# Libraries ( modules)
# *************************************************************************
import numpy as np
import logging
import scipy.optimize as opt
import matplotlib.pyplot as plt
from typing import Literal
from tqdm import tqdm
from wolfhece.wolf_array import WolfArray
# from fishydrocodes.tools import refine_WolfArray
# Constants
# *************************************************************************
# Objects
# *************************************************************************
[docs]
class RepCell:
"""
Class to compute the representative cell length
based on mesh characteristics.
Currently, only a number of cell or WolfArray are totally supported.
"""
def __init__(self,
cells: WolfArray|int,
dimensions: Literal['1D','2D','3D']= None,
work_in_relative: bool= False):
# Defining attributes
# *********************************************************
# Setting and checking inputs
# *********************************************************
self.work_in_relative = work_in_relative
if isinstance(cells, int):
self._nb_cells = cells
assert dimensions in ['1D','2D','3D'], "Dimensions must be either '1D', '2D' or '3D'."
self._dimensions = dimensions
elif isinstance(cells, WolfArray):
self._get_info_from_WolfArray(cells)
else:
raise TypeError("cells must be either an integer or a WolfArray object.")
@ property
[docs]
def array(self) -> WolfArray:
"""
Return the provided mesh.
:return: The WolfArray object.
:rtype: WolfArray
"""
return self._array
@ property
[docs]
def dimensions(self) -> int:
"""
Return the number of dimensions of the mesh.
:return: The number of dimensions.
:rtype: int
"""
return self._dimensions
@ property
[docs]
def nb_cells(self) -> int:
"""
Return the number of cells in the mesh.
In case the mesh is a WolfArray,
it returns the number of non-null cells.
:return: The number of cells.
:rtype: int
"""
return self._nb_cells
@ property
[docs]
def work_in_relative(self) -> bool:
"""
Return whether the computation of
the representative cell length is done in relative
terms (True) or absolute terms (False).
Relative computation assumes a unit domain size,
meaning the total length/area/volume is 1.
This formula is mostly used when the mesh is uniform
and square/cubic.
One can consult https://turbmodels.larc.nasa.gov/bump_sa.html,
for more information.
:return: Whether to compute the representative cell length in relative terms.
:rtype: bool
"""
return self._work_in_relative
@ work_in_relative.setter
def work_in_relative(self, value: bool):
"""
Set whether the computation of
the representative cell length is done in relative
terms (True) or absolute terms (False).
Relative computation assumes a unit domain size,
meaning the total length/area/volume is 1.
This formula is mostly used when the mesh is uniform
and square/cubic.
One can consult https://turbmodels.larc.nasa.gov/bump_sa.html,
for more information.
:param value: Whether to compute the representative cell length in relative terms.
:type value: bool
:return: None
"""
assert isinstance(value, bool), "work_in_relative must be a boolean."
self._work_in_relative = value
@ property
[docs]
def length(self) -> float:
"""
Compute the representative cell length based on
the number of cells and dimensions of the mesh.
FIXME Only uniform square/cubic meshes are currently supported
for wolfarrays.
:return: The representative cell length.
:rtype: float
"""
assert isinstance(self.nb_cells, (int, np.integer)),\
"nb_cells must be an integer not {}".format(type(self.nb_cells))
if self.work_in_relative or self.array is None:
if self.dimensions == '1D':
return 1/self.nb_cells
elif self.dimensions == '2D':
return 1.0/np.sqrt(self.nb_cells)
elif self.dimensions == '3D':
return 1.0/np.cbrt(self.nb_cells)
else:
raise ValueError("Dimensions must be either '1D', '2D', or '3D'.")
else:
if isinstance(self.array, WolfArray):
if self.dimensions == '1D':
return (1/self.nb_cells)*(self.nb_cells*self.array.dx)
elif self.dimensions == '2D':
return np.sqrt((1/self.nb_cells)*(self.nb_cells*self.array.dx*self.array.dy))
elif self.dimensions == '3D':
return np.cbrt((1/self.nb_cells)*(self.nb_cells*self.array.dx*self.array.dy*self.array.dz))
area = self.array.dx * self.array.dy
return np.sqrt(area)
[docs]
def _get_info_from_WolfArray(self,
wolfarray: WolfArray):
"""
Get the necessary information from a WolfArray object
and populate the attributes accordingly.
:param cells: The WolfArray object containing the mesh data.
:type cells: WolfArray
"""
self._array = wolfarray
self._dimensions = f"{self._array.nbdims}D"
self._nb_cells = int(self._array.count())
[docs]
def compute_representative_cell_length_2D_unstructured_grid(self,
nb_cells: int,
sum_of_root_of_average_areas: float|int) -> float:
"""
Compute the representative cell size of a 2D unstructured mesh.
Formula: h = (sum of root of average areas)/N or $h = 1/N * ΣA_p^{1/2}$
FIXME: Kept for unstructured grids or multiblock arrays when needed.
:param nb_cells: The number of cells in the mesh.
:type nb_cells: int
:param sum_of_root_of_average_areas: The sum of the root of average areas of the cells in the mesh.
:type sum_of_root_of_average_areas: float|int
:return: The representative cell size.
:rtype: float
"""
assert isinstance(nb_cells, (int, np.integer)), "nb_cells must be an integer not {}".format(type(nb_cells))
assert isinstance(sum_of_root_of_average_areas, (float, int, np.floating, np.integer)),\
"sum_of_root_of_average_areas must be a float or integer not {}".format(type(sum_of_root_of_average_areas))
return sum_of_root_of_average_areas/nb_cells
[docs]
def compute_representative_cell_length_3D_unstructured_grid(self,
nb_cells: int,
sum_of_cubic_root_of_average_volumes: float|int) -> float:
"""
Compute the representative cell size of a 3D unstructured mesh.
Formula: h = (sum of cubic root of average volumes)/N or $h = 1/N * ΣV_p^{1/3}$
FIXME: Kept for unstructured grids or multiblock arrays when needed.
:param nb_cells: The number of cells in the mesh.
:type nb_cells: int
:param sum_of_cubic_root_of_average_volumes: The sum of the cubic root of average volumes of the cells in the mesh.
:type sum_of_cubic_root_of_average_volumes: float|int
:return: The representative cell size.
:rtype: float
"""
assert isinstance(nb_cells, (int, np.integer)), "nb_cells must be an integer not {}".format(type(nb_cells))
assert isinstance(sum_of_cubic_root_of_average_volumes, (float, int, np.floating, np.integer)),\
"sum_of_cubic_root_of_average_volumes must be a float or integer not {}".format(type(sum_of_cubic_root_of_average_volumes))
return sum_of_cubic_root_of_average_volumes/nb_cells
[docs]
class RefinementRatio:
"""
Class to compute the refinement ratio between two meshes
based on the number of cells and dimensions of the meshes.
Only for uniform square/cubic meshes.
"""
def __init__(self,
cell_length_fine: int|float,
cell_length_coarse: int|float,
):
# Defining attributes
# *********************************************************
self.cells_length_fine = cell_length_fine
self.cells_length_coarse = cell_length_coarse
@ property
[docs]
def cells_length_fine(self) -> int|float:
"""
Return the cell length of the fine mesh.
:return: The cell length of the fine mesh.
:rtype: int|float
"""
return self._cells_length_fine
@ cells_length_fine.setter
def cells_length_fine(self, value: int|float):
"""
Set the cell length of the fine mesh.
:param value: The cell length of the fine mesh.
:type value: int|float
:return: None
"""
if self._cells_length_coarse is None:
self._cells_length_fine = value
else:
self._set_cell_lengths(value,
self._cells_length_coarse)
@ property
[docs]
def cells_length_coarse(self) -> int|float:
"""
Return the cell length of the coarse mesh.
:return: The cell length of the coarse mesh.
:rtype: int|float
"""
return self._cells_length_coarse
@ cells_length_coarse.setter
def cells_length_coarse(self, value: int|float):
"""
Set the cell length of the coarse mesh.
:param value: The cell length of the coarse mesh.
:type value: int|float
:return: None
"""
if self._cells_length_fine is None:
self._cells_length_coarse = value
else:
self._set_cell_lengths(self._cells_length_fine,
value)
@ property
[docs]
def ratio(self) -> float:
"""
Compute the refinement ratio between two meshes.
:return: The refinement ratio.
:rtype: float
"""
return self.cells_length_coarse / self.cells_length_fine
@ property
[docs]
def ratio_checked(self)-> float:
"""
Check if the cell size ratios are good for sensitivity analysis (Richardson extrapolation).
A good ratio is typically greater than 1.3 see Celik et al. (2008).
Celik I, Ghia U, Roache P and Freitas C, (2008) ’
Procedure for Estimation and Reporting of Uncertainty due to
Discretisation in CFD Applications’,
Journal of Fluids Engineering, 130(7), 078001. DOI: 10.1115/1.2960953
:return: The cell size ratio if it is good.
:rtype: float
"""
assert self.ratio > 1.3,\
"A good ratio must be greater than 1.3, now it is {}".format(self.ratio)
return self.ratio
[docs]
def _set_cell_lengths(self,
cell_length_fine: int|float,
cell_length_coarse: int|float):
"""
Set the cell lengths ensuring that the fine mesh
has a smaller cell length than the coarse mesh.
"""
assert isinstance(cell_length_fine, (int, float, np.integer, np.floating)),\
"cell_length_fine must be an integer or float not {}".format(type(cell_length_fine))
assert isinstance(cell_length_coarse, (int, float, np.integer, np.floating)),\
"cell_length_coarse must be an integer or float not {}".format(type(cell_length_coarse))
if cell_length_fine <= cell_length_coarse:
self._cells_length_fine = cell_length_fine
self._cells_length_coarse = cell_length_coarse
else:
self._cells_length_fine = cell_length_coarse
self._cells_length_coarse = cell_length_fine
[docs]
class OrderOfConvergence:
"""
Class to compute the order of convergence p
used in Richardson extrapolation and
Grid Convergence Index (GCI).
"""
def __init__(self,
r21: float,
r32: float,
s: int,
epsilon21: float,
epsilon32: float):
assert self.convergence_check, \
f"Grid convergence not achieved {self.convergence_ratio}\
cannot compute order of convergence p."
@ property
[docs]
def r21(self) -> float:
"""
Return the refinement ratio
between fine and medium meshes (grids)."""
return self._r21
@ property
@ property
[docs]
def s(self) -> int:
"""
Return the sign of grid convergence
- (1) for Monotonic convergence (same sign of errors),
- (-1) for Oscillatory convergence (different sign of errors).
"""
return self._s
@ property
[docs]
def epsilon21(self) -> float:
"""
Return the absolute difference (error)
between fine and medium mesh solutions.
"""
return self._epsilon21
@ property
[docs]
def epsilon32(self) -> float:
"""
Return the absolute difference (error)
between medium and coarse mesh solutions.
"""
return self._epsilon32
@ property
[docs]
def convergence_ratio(self) -> float:
"""
Compute the convergence ratio.
$C_R = \epsilon_{21} / \epsilon_{32}$
:return: The convergence ratio.
:rtype: float
"""
return self.epsilon21/ self.epsilon32
@ property
[docs]
def convergence_type(self) -> str:
if self.convergence_ratio >= 1:
return MONOTONIC_DIVERGENCE
if 0 <= self.convergence_ratio < 1:
return MONOTONIC_CONVERGENCE
if -1 <= self.convergence_ratio < 0:
return OSCILLATORY_CONVERGENCE
if self.convergence_ratio < -1:
return OSCILLATORY_DIVERGENCE
@ property
[docs]
def convergence_check(self) -> bool:
"""
Check if the grid convergence is achieved.
:return: True if grid convergence is achieved, False otherwise.
:rtype: bool
"""
if self.convergence_type in [MONOTONIC_CONVERGENCE, OSCILLATORY_CONVERGENCE]:
return True
else:
return False
[docs]
def compute_q(self,
p: float) -> float:
"""
Compute the q value based on the order of convergence and refinement ratios.
:param p: The order of convergence.
:type p: float
:return: The computed q value.
:rtype: float
"""
return np.log((np.power(self.r21, p) - self.s) / (np.power(self.r32, p) - self.s))
[docs]
def compute_residual_of_p(self,
p: float) -> float:
"""
Compute the residual of the order of convergence.
:param p: The order of convergence.
:type p: float
:return: The computed residual.
:rtype: float
"""
q = self.compute_q(p)
return (1.0/np.log(self.r21)) * np.abs((np.log(np.abs(self.epsilon32/self.epsilon21))) + q) - p
[docs]
def bisection(self):
"""
Compute the order of convergence using the bisection method.
:return: The computed order of convergence.
:rtype: float
"""
initial_guess_a = 1e-11
initial_guess_b = 2.
check = 0
if isinstance(self.epsilon21, (WolfArray)) and isinstance(self.epsilon32, (WolfArray)):
initial_guess_a = np.full_like(self.epsilon21, initial_guess_a)
initial_guess_b = np.full_like(self.epsilon21, initial_guess_b)
tolerance = 1e-11
max_iterations = 1000000
residual_initial_guess_a = self.compute_residual_of_p(initial_guess_a)
residual_initial_guess_b = self.compute_residual_of_p(initial_guess_b)
if np.all(residual_initial_guess_a * residual_initial_guess_b >= check):
raise ValueError(f"The root is not bounded in the probable interval - sign of residuals\
but got [{residual_initial_guess_a}, {residual_initial_guess_b}]")
return opt.bisect(self.compute_residual_of_p,
initial_guess_a,
initial_guess_b,
xtol=tolerance,
maxiter=max_iterations)
[docs]
def newton_raphson(self) -> float:
"""
Compute the order of convergence using the Newton-Raphson method.
:param initial_guess: The initial guess for the order of convergence.
:type initial_guess: float
:param tolerance: The tolerance for convergence.
:type tolerance: float
:param max_iterations: The maximum number of iterations.
:type max_iterations: int
:return: The computed order of convergence.
:rtype: float
"""
initial_guess = 1
max_iterations = 1000000
tolerance = 1e-11
return opt.newton(self.compute_residual_of_p,
initial_guess,
maxiter=max_iterations,
tol=tolerance)
[docs]
class RichardsonExtrapolation:
"""
Class to perform Richardson extrapolation and related calculations.
"""
def __init__(self,
array_fine: WolfArray|int,
array_medium: WolfArray|int,
array_coarse: WolfArray|int,
solution_fine: WolfArray|float|int,
solution_medium: WolfArray|float|int,
solution_coarse: WolfArray|float|int,
dimensions: Literal['2D','3D']= None,
computation_method: Literal["Bisection", "Newton-Raphson"]= None,
work_in_relative: bool= False):
"""
Initialize the RichardsonExtrapolation class with mesh arrays and solutions.
:param array_fine: The fine mesh array (WolfArray or number of cells).
:type array_fine: WolfArray|int
:param array_medium: The medium mesh array (WolfArray or number of cells).
:type array_medium: WolfArray|int
:param array_coarse: The coarse mesh array (WolfArray or number of cells).
:type array_coarse: WolfArray|int
:param solution_fine: The solution on the fine mesh (WolfArray or float).
:type solution_fine: WolfArray|float|int
:param solution_medium: The solution on the medium mesh (WolfArray or float).
:type solution_medium: WolfArray|float|int
:param solution_coarse: The solution on the coarse mesh (WolfArray or float).
:type solution_coarse: WolfArray|float|int
:param dimensions: The number of dimensions ('2D' or '3D'). Default is '2D'.
:type dimensions: Literal['2D','3D'], optional
:param computation_method: The method to compute the order of convergence
('Bisection' or 'Newton-Raphson'). Default is 'Newton-Raphson'.
:type computation_method: Literal["Bisection", "Newton-Raphson"], optional
"""
self.computation_method = computation_method
if dimensions is None:
self._dimensions = '2D'
else:
assert dimensions in ['2D','3D'], "Dimensions must be either '2D' or '3D'."
self._dimensions = dimensions
self.set_arrays_and_solutions(array_fine=array_fine,
array_medium=array_medium,
array_coarse=array_coarse,
solution_fine=solution_fine,
solution_medium=solution_medium,
solution_coarse=solution_coarse)
@ property
[docs]
def computation_method(self) -> str:
"""
Return the computation method for the order of convergence.
:return: The computation method.
:rtype: str
"""
return self._computation_method
@ computation_method.setter
def computation_method(self, method: Literal["Bisection", "Newton-Raphson"] = None):
"""
Set the computation method for the order of convergence.
:param method: The computation method ('Bisection' or 'Newton-Raphson').
:type method: Literal["Bisection", "Newton-Raphson"], optional
:return: None
"""
if method is None:
self._computation_method = NEWTON_RAPHSON
else:
assert method in ["Bisection", "Newton-Raphson"],\
"Computation method must be either 'Bisection' or 'Newton-Raphson'."
self._computation_method = method
@ property
[docs]
def work_in_relative(self) -> bool:
"""
Return whether the computation of
the representative cell length is done in relative
terms (True) or absolute terms (False).
Relative computation assumes a unit domain size,
meaning the total length/area/volume is 1.
This formula is mostly used when the mesh is uniform
and square/cubic.
One can consult https://turbmodels.larc.nasa.gov/bump_sa.html,
for more information.
:return: Whether to compute the representative cell length in relative terms.
:rtype: bool
"""
return self._work_in_relative
@ work_in_relative.setter
def work_in_relative(self, value: bool):
"""
Set whether the computation of
the representative cell length is done in relative
terms (True) or absolute terms (False).
Relative computation assumes a unit domain size,
meaning the total length/area/volume is 1.
This formula is mostly used when the mesh is uniform
and square/cubic.
One can consult https://turbmodels.larc.nasa.gov/bump_sa.html,
for more information.
:param value: Whether to compute the representative cell length in relative terms.
:type value: bool
:return: None
"""
assert isinstance(value, bool), "work_in_relative must be a boolean."
self._work_in_relative = value
self._array_fine.work_in_relative = self.work_in_relative
self._array_medium.work_in_relative = self.work_in_relative
self._array_coarse.work_in_relative = self.work_in_relative
@ property
[docs]
def e21_diff_fine_medium(self) -> WolfArray|float:
"""
Return the difference between fine and medium solutions.
:return: The difference between fine and medium solutions.
:rtype: WolfArray|float """
return self.difference_between_solutions(coarse= self._solution_medium,
fine= self._solution_fine)
@ property
[docs]
def e32_diff_medium_coarse(self) -> WolfArray|float:
"""
Return the difference between medium and coarse solutions.
:return: The difference between medium and coarse solutions.
:rtype: WolfArray|float
"""
return self.difference_between_solutions(coarse= self._solution_coarse,
fine= self._solution_medium)
@ property
[docs]
def s_sign_of_convergence(self) -> WolfArray|int:
"""
Return the sign of the convergence based on the differences
between solutions.
:return: The sign of convergence.
:rtype: WolfArray|int
"""
return self.sign_of_convergence(epsilon21= self.e21_diff_fine_medium,
epsilon32= self.e32_diff_medium_coarse)
@ property
[docs]
def h1_representative_cell_length_fine(self) -> float:
"""
Return the representative cell length of the fine mesh.
:return: The representative cell length of the fine mesh.
:rtype: float"""
return self._array_fine.length
@ property
[docs]
def h2_representative_cell_length_medium(self) -> float:
"""
Return the representative cell length of the medium mesh.
:return: The representative cell length of the medium mesh.
:rtype: float
"""
return self._array_medium.length
@ property
[docs]
def h3_representative_cell_length_coarse(self) -> float:
"""
Return the representative cell length of the coarse mesh.
:return: The representative cell length of the coarse mesh.
:rtype: float
"""
return self._array_coarse.length
@ property
[docs]
def r21_refinement_ratio(self) -> float:
"""
Return the refinement ratio between fine and medium meshes.
:return: The refinement ratio between fine and medium meshes.
:rtype: float
"""
return RefinementRatio(self.h1_representative_cell_length_fine,
self.h2_representative_cell_length_medium).ratio_checked
@ property
[docs]
def r32_refinement_ratio(self) -> float:
"""
Return the refinement ratio between medium and coarse meshes.
:return: The refinement ratio between medium and coarse meshes.
:rtype: float
"""
return RefinementRatio(self.h2_representative_cell_length_medium,
self.h3_representative_cell_length_coarse).ratio_checked
[docs]
def set_arrays_and_solutions(self,
array_fine: WolfArray|int,
array_medium: WolfArray|int,
array_coarse: WolfArray|int,
solution_fine: WolfArray|float|int,
solution_medium: WolfArray|float|int,
solution_coarse: WolfArray|float|int):
"""
Set the mesh arrays and mesh solutions (characteristics)
for Richardson extrapolation.
:param array_fine: The fine mesh array (WolfArray or number of cells).
:type array_fine: WolfArray|int
:param array_medium: The medium mesh array (WolfArray or number of cells).
:type array_medium: WolfArray|int
:param array_coarse: The coarse mesh array (WolfArray or number of cells).
:type array_coarse: WolfArray|int
:param solution_fine: The solution on the fine mesh (WolfArray or float).
:type solution_fine: WolfArray|float|int
:param solution_medium: The solution on the medium mesh (WolfArray or float).
:type solution_medium: WolfArray|float|int
:param solution_coarse: The solution on the coarse mesh (WolfArray or float).
:type solution_coarse: WolfArray|float|int
"""
self._array_fine = RepCell(array_fine, self._dimensions, self._work_in_relative)
self._array_medium = RepCell(array_medium, self._dimensions, self._work_in_relative)
self._array_coarse = RepCell(array_coarse, self._dimensions, self._work_in_relative)
if self._array_fine.array is None:
assert isinstance(solution_fine, (float, int, np.floating, np.integer)),\
"solution_fine must be a float or integer not {}".format(type(solution_fine))
self._solution_fine = solution_fine
else:
self._solution_fine = solution_fine
if self._array_medium.array is None:
assert isinstance(solution_medium, (float, int, np.floating, np.integer)),\
"solution_medium must be a float or integer not {}".format(type(solution_medium))
self._solution_medium = solution_medium
else:
self._solution_medium = solution_medium
if self._array_coarse.array is None:
assert isinstance(solution_coarse, (float, int, np.floating, np.integer)),\
"solution_coarse must be a float or integer not {}".format(type(solution_coarse))
self._solution_coarse = solution_coarse
else:
self._solution_coarse = solution_coarse
# Refininining solutions size to allow the computation of differences
if isinstance(self._solution_fine, WolfArray) and isinstance(self._solution_medium, WolfArray) and isinstance(self._solution_coarse, WolfArray):
# refined_solution_medium = refine_WolfArray(self._solution_medium,
# self.factor_21())
# refined_solution_coarse = refine_WolfArray(self._solution_coarse,
# self.factor_31())
self._solution_medium.rebin(factor= 1/self.factor_21())
self._solution_coarse.rebin(factor= 1/self.factor_31())
bounds = self._solution_fine.get_bounds()
# self._solution_medium = refined_solution_medium.crop_array(bounds)
# self._solution_coarse = refined_solution_coarse.crop_array(bounds)
self._solution_medium = self._solution_medium.crop_array(bounds)
self._solution_coarse = self._solution_coarse.crop_array(bounds)
assert self._solution_fine.shape == self._solution_medium.shape == self._solution_coarse.shape,\
"Refined solutions must have the same shape."
[docs]
def difference_between_solutions(self,
coarse: WolfArray|float,
fine: WolfArray|float) -> WolfArray|float:
"""
Compute the difference between two solutions.
:param phi1: The first solution value.
:type phi1: float
:param phi2: The second solution value.
:type phi2: float
:return: The difference between the two solutions.
:rtype: float
"""
if isinstance(coarse, WolfArray) and isinstance(fine, WolfArray):
assert coarse.shape == fine.shape,\
"Coarse and fine solutions must have the same shape."
new_array = WolfArray(mold= coarse)
new_array.array = coarse.array - fine.array
return new_array
return coarse - fine
[docs]
def sign_of_convergence(self,
epsilon21: float|WolfArray,
epsilon32: float|WolfArray) -> WolfArray|int:
"""
Determine the sign of convergence based on the ratio of errors.
:param epsilon21: Error between fine and medium mesh solutions.
:type epsilon21: float
:param epsilon32: Error between medium and coarse mesh solutions.
:type epsilon32: float
:return: 1 if the ratio is positive (indicating monotonic convergence),
-1 otherwise (indicating oscillatory convergence).
:rtype: int
"""
if isinstance(epsilon21, WolfArray) and isinstance(epsilon32, WolfArray):
assert epsilon21.shape == epsilon32.shape,\
"epsilon21 and epsilon32 must have the same shape."
new_array = WolfArray(mold= epsilon21)
new_array.array = np.sign(epsilon32.array / epsilon21.array)
return new_array
return np.sign(epsilon32 / epsilon21)
[docs]
def factor_21(self) -> float:
"""
Compute the refinement factor between fine and medium meshes.
:return: The refinement factor between fine and medium meshes.
:rtype: float
"""
return self._array_medium.length / self._array_fine.length
[docs]
def factor_31(self) -> float:
"""
Compute the refinement factor between fine and coarse meshes.
:return: The refinement factor between fine and coarse meshes.
:rtype: float
"""
return self._array_coarse.length / self._array_fine.length
@ property
[docs]
def order_of_convergence(self,
compute_method: Literal['Bisection','Newton-Raphson']= None) -> WolfArray|float:
"""
Compute the order of convergence p using the specified method.
:param compute_method: The method to compute the order of convergence
('Bisection' or 'Newton-Raphson'). Default is Newton-Raphson to avoid convergence issues.
:return: The order of convergence p.
:rtype: float|WolfArray
"""
if compute_method is None:
compute_method= self.computation_method
if isinstance(self.e21_diff_fine_medium, WolfArray) and isinstance(self.e32_diff_medium_coarse, WolfArray):
e21 = self.e21_diff_fine_medium.array
e32 = self.e32_diff_medium_coarse.array
s = self.s_sign_of_convergence.array
s = np.where(s == 0, 1, s) # Replace zeros with ones to avoid division by zero
s= np.where(np.isnan(s), 1, s) # Replace NaNs with ones to avoid issues
p_wolfarray = WolfArray(mold= self.e21_diff_fine_medium)
else:
e21 = self.e21_diff_fine_medium
e32 = self.e32_diff_medium_coarse
s = self.s_sign_of_convergence
p_order_of_convergence = OrderOfConvergence(self.r21_refinement_ratio,
self.r32_refinement_ratio,
s,
e21,
e32)
if compute_method == BISECTION:
p = p_order_of_convergence.bisection()
elif compute_method == NEWTON_RAPHSON:
p = p_order_of_convergence.newton_raphson()
else:
raise ValueError("Method must be either 'Bisection' or 'Newton-Raphson'.")
if isinstance(self.e21_diff_fine_medium, WolfArray) and isinstance(self.e32_diff_medium_coarse, WolfArray):
p_wolfarray.array = p
return p_wolfarray
return p
@ property
[docs]
def extrapolated_solution(self) -> float|WolfArray:
"""
Compute the extrapolated solution (ideal mesh) using Richardson extrapolation.
:param phi_fine: The solution on the fine mesh.
:type phi_fine: float
:param phi_medium: The solution on the medium mesh.
:type phi_medium: float
:param r21: The grid refinement ratio between the fine and medium meshes.
:type r21: float
:param p: The order of convergence.
:type p: float
:return: The extrapolated solution.
:rtype: float
"""
phi_fine = self._solution_fine
phi_medium = self._solution_medium
r21 = self.r21_refinement_ratio
p = self.order_of_convergence
if isinstance(phi_fine, WolfArray) and isinstance(phi_medium, WolfArray) and isinstance(r21, WolfArray) and isinstance(p, WolfArray):
assert phi_fine.shape == phi_medium.shape == r21.shape == p.shape,\
"phi_fine, phi_medium, r21, and p must have the same shape."
new_array = WolfArray(mold= phi_fine)
new_array.array = ((np.power(r21.array, p.array) * phi_fine.array) - phi_medium.array) / (np.power(r21.array, p.array) - 1)
return new_array
return ((np.power(r21,p)*phi_fine) - phi_medium)/(np.power(r21,p)-1)
@ property
[docs]
def extrapolated_relative_error_21(self) -> float|WolfArray:
"""
Compute the relative error between the fine mesh solution and the extrapolated solution.
:param phi_fine: The solution on the fine mesh.
:type phi_fine: float
:param phi_extrapolated: The extrapolated solution.
:type phi_extrapolated: float
:return: The relative error.
:rtype: float
"""
phi_fine = self._solution_fine
phi_extrapolated = self.extrapolated_solution
if isinstance(phi_fine, WolfArray) and isinstance(phi_extrapolated, WolfArray):
assert phi_fine.shape == phi_extrapolated.shape,\
"phi_fine and phi_extrapolated must have the same shape."
new_array = WolfArray(mold= phi_fine)
new_array.array = np.abs((phi_extrapolated.array - phi_fine.array) / phi_extrapolated.array)
return new_array
return np.abs((phi_extrapolated - phi_fine) / phi_extrapolated)
[docs]
def compute_relative_error(self,
coarse: float|WolfArray,
fine: float|WolfArray) -> float|WolfArray:
"""
Compute the relative error between two solutions.
:param phi1: The first solution value.
:type phi1: float
:param phi2: The second solution value.
:type phi2: float
:return: The relative error between the two solutions.
:rtype: float
"""
if isinstance(coarse, WolfArray) and isinstance(fine, WolfArray):
assert coarse.shape == fine.shape,\
"Coarse and fine solutions must have the same shape."
new_array = WolfArray(mold= coarse)
new_array.array = np.abs((fine.array - coarse.array) / fine.array)
return new_array
return np.abs((fine - coarse) / fine)
@ property
[docs]
def relative_error_21(self) -> float|WolfArray:
"""
Return the relative error between the fine and medium mesh solutions.
:return: The relative error between the fine and medium mesh solutions.
:rtype: float|WolfArray
"""
return self.compute_relative_error(coarse= self._solution_medium,
fine= self._solution_fine)
@ property
[docs]
def relative_error_32(self) -> float|WolfArray:
"""
Return the relative error between the medium and coarse mesh solutions.
:return: The relative error between the medium and coarse mesh solutions.
:rtype: float|WolfArray
"""
return self.compute_relative_error(coarse= self._solution_coarse,
fine= self._solution_medium)
@ property
[docs]
def grid_convergence_index_21(self) -> float|WolfArray:
"""
Compute the Grid Convergence Index (GCI) for the fine and medium mesh solutions.
see Roache (1994) for more details.
Roache P. J. (1994), ’Perspective: A method for uniform reporting of grid refinement studies’, Journal
of Fluids Engineering, 116, pp. 405-413
:param rel_error_21: The relative error between the fine and medium mesh solutions.
:type rel_error_21: float
:param r21: The grid refinement ratio between the fine and medium meshes.
:type r21: float
:param p: The order of convergence.
:type p: float
:return: The Grid Convergence Index (GCI).
:rtype: float
"""
if isinstance(self.relative_error_21, WolfArray) and isinstance(self.r21_refinement_ratio, WolfArray) and isinstance(self.order_of_convergence, WolfArray):
assert self.relative_error_21.shape == self.r21_refinement_ratio.shape == self.order_of_convergence.shape,\
"self.relative_error_21, self.r21_refinement_ratio, and p must have the same shape."
new_array = WolfArray(mold= self.relative_error_21)
new_array.array = 1.25 * self.relative_error_21.array / (np.power(self.r21_refinement_ratio.array, self.order_of_convergence.array) - 1)
return new_array
return 1.25 * self.relative_error_21 / (np.power(self.r21_refinement_ratio, self.order_of_convergence) - 1)
@ property
[docs]
def grid_convergence_index_32(self) -> float|WolfArray:
"""
Compute the Grid Convergence Index (GCI) for the medium and coarse mesh solutions.
see Roache (1994) for more details.
Roache P. J. (1994), ’Perspective: A method for uniform reporting of grid refinement studies’, Journal
of Fluids Engineering, 116, pp. 405-413
:param rel_error_32: The relative error between the medium and coarse mesh solutions.
:type rel_error_32: float
:param r32: The grid refinement ratio between the medium and coarse meshes.
:type r32: float
:param p: The order of convergence.
:type p: float
:return: The Grid Convergence Index (GCI).
:rtype: float
"""
if isinstance(self.relative_error_32, WolfArray) and isinstance(self.r32_refinement_ratio, WolfArray) and isinstance(self.order_of_convergence, WolfArray):
assert self.relative_error_32.shape == self.r32_refinement_ratio.shape == self.order_of_convergence.shape,\
"self.relative_error_32, self.r32_refinement_ratio, and p must have the same shape."
new_array = WolfArray(mold= self.relative_error_32)
new_array.array = 1.25 * self.relative_error_32.array / (np.power(self.r32_refinement_ratio.array, self.order_of_convergence.array) - 1)
return new_array
return 1.25 * self.relative_error_32 / (np.power(self.r32_refinement_ratio, self.order_of_convergence) - 1)
@ property
[docs]
def AsymptoticRange(self) -> bool|WolfArray:
"""
return the investment on computational cost"""
return self.grid_convergence_index_32/((self.relative_error_21**self.order_of_convergence)*self.grid_convergence_index_21)
[docs]
def optimal_refinement_ratio_from_GCI(self,
GCI_target: float) -> float:
"""
Compute the optimal refinement ratio to achieve a target GCI.
:param GCI_target: The target Grid Convergence Index (GCI).
:type GCI_target: float
:return: The optimal refinement ratio.
:rtype: float
"""
return (GCI_target/self.grid_convergence_index_21)**(1/self.order_of_convergence)
[docs]
def plot_extrapolated_solution(self,
ax=None,
figure_size:tuple=(8,6),
color = 'red',
title: str="Extrapolated Solution",
ylabel: str="Solution Value",
ylim: tuple|None=None,
xlim: tuple|None=None,
yfactor: int|float=None,
label: str = 'Solutions',
show: bool=False) -> plt.Axes:
if isinstance(self.extrapolated_solution, WolfArray):
attributes = self.extrapolated_solution.plot_matplotlib(figax = ax)
if show:
plt.tight_layout()
plt.show()
return attributes[1]
# show = False
labels = True # show axes labels and title
if ax is None:
fig, ax = plt.subplots(figsize=figure_size)
labels = False
# show = True
representative_cell_lengths = [
self.h1_representative_cell_length_fine,
self.h2_representative_cell_length_medium,
self.h3_representative_cell_length_coarse
]
solutions = [
self._solution_fine,
self._solution_medium,
self._solution_coarse
]
all_cell_lengths = [0.] + representative_cell_lengths
all_solutions = [self.extrapolated_solution] + solutions
if yfactor is not None:
solutions = [solution / yfactor for solution in solutions]
all_solutions = [solution / yfactor for solution in all_solutions]
points_annotations = ['Infinite'] + ['Fine', 'Medium', 'Coarse']
ax.plot(representative_cell_lengths,
solutions,
marker='o',
ls = '-' ,
lw=2,
color=color,
label=label)
ax.plot(all_cell_lengths,
all_solutions,
marker = 'o',
ls = '--' ,
color=color,
lw=2,
# label='Extrapolated Solution [Ideal Mesh]',
)
for x, y, label in zip(all_cell_lengths, all_solutions, points_annotations):
ax.text(x,
y,
label,
fontsize=10,
ha='left',
va='bottom',
bbox=dict(facecolor='white', alpha=0.1, boxstyle='square,pad=0.2'),
)
if labels:
if xlim is not None:
ax.set_xlim(xlim)
else:
ax.set_xlim(0, max(representative_cell_lengths)*1.1)
if ylim is not None:
ax.set_ylim(ylim)
ax.set_xlabel("h - Representative Cell Length")
ax.set_ylabel(ylabel)
ax.set_title(title)
ax.legend()
ax.grid(lw = 0.3, ls = '--', color = 'grey', alpha=0.7)
if show:
plt.tight_layout()
plt.show()
return ax
[docs]
def get_errors_summary(self, precision = 4) -> dict:
"""
Get a dictionary of errors and related metrics.
in the dictionary there are:
- E_extr %: extrapolated relative error in percentage
- E_rel %: relative error in percentage
- GCI %: grid convergence index in percentage
- p - Order of convergence
- GCI_32: grid convergence index between medium and coarse meshes in percentage
:return: A dictionary containing errors and related metrics.
:rtype: dict
"""
return {
"E_extr %": round(self.extrapolated_relative_error_21 * 100, precision),
"E_rel %": round(self.relative_error_21 * 100, precision),
"GCI %": round(self.grid_convergence_index_21 * 100, precision),
"GCI_32 %": round(self.grid_convergence_index_32 * 100, precision),
"p - Order of converegence": round(self.order_of_convergence, precision),
}
[docs]
def get_results_summary(self,
precision = None,
precision_error = None) -> dict:
"""
in the dictionary there are:
- h fine: representative cell length of the fine mesh
- h medium: representative cell length of the medium mesh
- h coarse: representative cell length of the coarse mesh
- h extrapolated: representative cell length of the extrapolated mesh (infinite cells)
- phi fine: solution value for the fine mesh
- phi medium: solution value for the medium mesh
- phi coarse: solution value for the coarse mesh
- phi extrapolated: extrapolated solution value
:param self: Description
:param precision: Description
:return: Description
:rtype: dict
"""
if precision is None:
precision = 5
if precision_error is None:
precision_error = 2
return{
"h fine": round(self.h1_representative_cell_length_fine, precision),
"h medium": round(self.h2_representative_cell_length_medium, precision),
"h coarse": round(self.h3_representative_cell_length_coarse, precision),
"h extrapolated": 0,
"phi fine": round(self._solution_fine, precision) if not isinstance(self._solution_fine, WolfArray) else self._solution_fine,
"phi medium": round(self._solution_medium, precision) if not isinstance(self._solution_medium, WolfArray) else self._solution_medium,
"phi coarse": round(self._solution_coarse, precision) if not isinstance(self._solution_coarse, WolfArray) else self._solution_coarse,
"phi extrapolated": round(self.extrapolated_solution, precision) if not isinstance(self.extrapolated_solution, WolfArray) else self.extrapolated_solution,
"E_extr %": round(self.extrapolated_relative_error_21 * 100, precision_error),
"E_rel %": round(self.relative_error_21 * 100, precision_error),
"GCI %": round(self.grid_convergence_index_21 * 100, precision_error),
"GCI_32 %": round(self.grid_convergence_index_32 * 100, precision_error),
}
# FIXME: add plot funtions: representative cell length vs solution
# representative cell length vs order of convergence
# representative cell length vs GCI
# representative cell length vs relative error
# representative cell length vs extrapolated solution