diff --git a/Makefile b/Makefile index b6baaf76a3..4da67a4831 100644 --- a/Makefile +++ b/Makefile @@ -276,6 +276,7 @@ docs: @echo "************************************" -${PROTEUS_ENV} pip3 install sphinx sphinx_rtd_theme breathe exhale + @${PROTEUS_ENV} sh -c 'command -v dot' >/dev/null || { echo "graphviz 'dot' not found; inheritance diagrams need it (conda install graphviz)"; exit 1; } cd docs && ${PROTEUS_ENV} PROTEUS=${PWD} make html @echo "**********************************" diff --git a/docs/Makefile b/docs/Makefile index e2edd87f2a..7c4da2574d 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -51,13 +51,23 @@ clean: rm -rf source/api rm -rf source/capi -html: $(BUILDDIR)/capi/html/files.html source/api/modules.rst +html: $(BUILDDIR)/capi/html/files.html source/api/proteus.rst $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR) @echo @echo "Build finished. The HTML pages are in $(BUILDDIR)." -source/api/modules.rst: - sphinx-apidoc --output-dir=source/api ../proteus "../**/*.so" "../proteus/tests" "../proteus/test_utils" "../proteus/config" --maxdepth=1 --separate --module-first --templatedir=source/_templates/apidoc +source/api/proteus.rst: + # flcbdfWrappers.pyx: extension is commented out in setup.py (unbuilt, + # dead code) -- no compiled .so exists for autodoc to import. + # SWFlow/utils/SWEs_so.py: a parun-driven split-operator driver script, + # not an importable API module -- its top-level code sniffs sys.argv + # for a case file to __import__() by name, which only makes sense when + # invoked by parun itself and raises ModuleNotFoundError otherwise. + # --no-toc: apicapi.rst links directly to api/proteus.html (the + # top-level package page), never to the modules.rst index apidoc + # generates by default -- that page would be orphaned (not reachable + # from any toctree) and Sphinx warns accordingly, so skip generating it. + sphinx-apidoc --output-dir=source/api ../proteus "../**/*.so" "../proteus/tests" "../proteus/test_utils" "../proteus/config" "../proteus/flcbdfWrappers.pyx" "../proteus/SWFlow/utils/SWEs_so.py" --maxdepth=1 --separate --module-first --no-toc --templatedir=source/_templates/apidoc $(BUILDDIR)/capi/html/files.html: mkdir -p ./build diff --git a/docs/source/_templates/apidoc/package.rst_t b/docs/source/_templates/apidoc/package.rst_t index 2274b9fe4f..e084d4336a 100644 --- a/docs/source/_templates/apidoc/package.rst_t +++ b/docs/source/_templates/apidoc/package.rst_t @@ -13,6 +13,10 @@ {%- endfor %} {%- endmacro %} +{%- if pkgname == "proteus" %} +:orphan: + +{% endif %} {%- if is_namespace %} {{- [pkgname, "namespace"] | join(" ") | e | heading }} {% else %} diff --git a/docs/source/conf.py b/docs/source/conf.py index 0ce412a7f8..18d0973ef9 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -166,7 +166,6 @@ def linkcode_resolve(domain, info): 'canonical_url': '', 'analytics_id': '', 'logo_only': True, - 'display_version': True, 'prev_next_buttons_location': 'bottom', 'style_external_links': False, # 'vcs_pageview_mode': '', @@ -290,8 +289,27 @@ def linkcode_resolve(domain, info): # not chapters. #latex_use_parts = False +# LaTeX macros usable in docstring math, for both the HTML (MathJax) and +# LaTeX builds: name -> (definition, number of arguments) +math_macros = { + 'grad': (r'\nabla', 0), + 'deld': (r'\nabla\cdot', 0), + 'outer': (r'\otimes', 0), + 'ds': (r'\,\mathrm{d}s', 0), + 'ten': (r'\mathbf{#1}', 1), + 'mat': (r'\mathbf{#1}', 1), + 'pd': (r'\frac{\partial #1}{\partial #2}', 2), + 'od': (r'\frac{\mathrm{d} #1}{\mathrm{d} #2}', 2), +} +mathjax3_config = { + 'tex': {'macros': {name: [body, nargs] if nargs else body + for name, (body, nargs) in math_macros.items()}}, +} + # Additional stuff for the LaTeX preamble. -latex_elements['preamble'] = '\input{macros}' +latex_elements['preamble'] = '\n'.join( + r'\providecommand{\%s}[%d]{%s}' % (name, nargs, body) + for name, (body, nargs) in math_macros.items()) # If true, show page references after internal links. #latex_show_pagerefs = False @@ -416,8 +434,6 @@ def linkcode_resolve(domain, info): graphviz_output_format = 'svg' -mathjax_pth = "https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-MML-AM_CHTM" - breathe_projects = { "Proteus C++": "../build/capi/xml" } @@ -450,4 +466,4 @@ def linkcode_resolve(domain, info): # "exhaleExecutesDoxygen": True, # #"exhaleUseDoxyfile": True, # "exhaleDoxygenStdin": exhaleDoxygenStdin -#} \ No newline at end of file +#} diff --git a/docs/source/tools/wave_tools.rst b/docs/source/tools/wave_tools.rst index bbd0f8a059..c15601361d 100644 --- a/docs/source/tools/wave_tools.rst +++ b/docs/source/tools/wave_tools.rst @@ -104,7 +104,7 @@ If the relaxation zone method is used, then the class should be passed through t tank.setGenerationZones(x_n=True, waves=wave, dragAlpha=dragAlpha, smoothing = smoothing) -Guidance for using the ``setUnsteadyTwoPhaseVelocityInlet`` and ``setGenerationZones`` functions are given in the :doc:`BoundaryConditions <./boundary_conditions>` and :doc:`Spatial Tools <./spatial_tools>` sections of the documentation +Guidance for using the ``setUnsteadyTwoPhaseVelocityInlet`` and ``setGenerationZones`` functions is given in the ``BoundaryConditions`` API reference and the :doc:`Spatial Tools <./spatial_tools>` section of the documentation Simple examples of usage within the context of a 2D numerical tank can be found in `air-water-vv `_ diff --git a/environment-dev-up.yml b/environment-dev-up.yml index 5d75650a06..7eb4114310 100644 --- a/environment-dev-up.yml +++ b/environment-dev-up.yml @@ -62,3 +62,5 @@ dependencies: - cffi >=1.17.1 - jupyterlab >=4.6 - llvm-openmp >=22.1 + # provides `dot` for the docs inheritance diagrams + - graphviz diff --git a/proteus/DiagUtils.py b/proteus/DiagUtils.py index fabd1f15bd..90e09734a1 100644 --- a/proteus/DiagUtils.py +++ b/proteus/DiagUtils.py @@ -1,8 +1,5 @@ """ Module for diagnostic utilities - -.. inheritance-diagram:: proteus.DiagUtils - :parts: 1 """ from .EGeometry import * from .MeshTools import * diff --git a/proteus/EGeometry.py b/proteus/EGeometry.py index c90f65ff71..3ff3fc4e24 100644 --- a/proteus/EGeometry.py +++ b/proteus/EGeometry.py @@ -1,8 +1,5 @@ """ A simple implementation of basic 3D Euclidean Geometry - -.. inheritance-diagram:: proteus.EGeometry - :parts: 1 """ import math from math import * diff --git a/proteus/FemTools.py b/proteus/FemTools.py index e40169de13..736978f68d 100644 --- a/proteus/FemTools.py +++ b/proteus/FemTools.py @@ -567,54 +567,55 @@ class BernsteinOnCube(LocalFunctionSpace): Bernstein polynomials on the unit nd-cube Bernstein basis functions on the reference nd-cube (nd <=3) with - coordinates xi[0],xi[1],and xi[2]. - The basis functions are numbered as follows: - ############## - # *** 2D *** # - ############## - ########### - # Order=1 # - ########### - 3--2 - | | - 0--1 - ########### - # Order=2 # - ########### - 3-6-2 - | | | - 7-8-5 - | | | - 0-4-1 - ########### - # Order=3 # - ########### - 3--9--8--2 - | | | | - 10-14-15-7 - | | | | - 11-12-13-6 - | | | | - 0--4--5--1 - ########### - # Order=4 # - ########### - 3--12-11-10-2 - | | | | | - 13-22-23-24-9 - | | | | | - 14-19-20-21-8 - | | | | | - 15-16-17-18-7 - | | | | | - 0--4--5--6--1 - - NOTE (mql): *** The polynomials are defined for arbitrary order in 1D. - This definition is carried to multi-D via tensor products. - *** The numbering is defined for arbitrary order (via funMap) in 2D. - *** TODO: Generalize numbering for 1D and 3D. - *** TODO: The quad rule Lobatto edge alt must be generalized. - *** TODO: define hessians and higher order deriavtives + coordinates xi[0],xi[1],and xi[2]. + The basis functions are numbered as follows:: + + ############## + # *** 2D *** # + ############## + ########### + # Order=1 # + ########### + 3--2 + | | + 0--1 + ########### + # Order=2 # + ########### + 3-6-2 + | | | + 7-8-5 + | | | + 0-4-1 + ########### + # Order=3 # + ########### + 3--9--8--2 + | | | | + 10-14-15-7 + | | | | + 11-12-13-6 + | | | | + 0--4--5--1 + ########### + # Order=4 # + ########### + 3--12-11-10-2 + | | | | | + 13-22-23-24-9 + | | | | | + 14-19-20-21-8 + | | | | | + 15-16-17-18-7 + | | | | | + 0--4--5--6--1 + + NOTE (mql): *** The polynomials are defined for arbitrary order in 1D. + This definition is carried to multi-D via tensor products. + *** The numbering is defined for arbitrary order (via funMap) in 2D. + *** TODO: Generalize numbering for 1D and 3D. + *** TODO: The quad rule Lobatto edge alt must be generalized. + *** TODO: define hessians and higher order deriavtives """ from math import factorial diff --git a/proteus/LinearAlgebraTools.py b/proteus/LinearAlgebraTools.py index 5752e86f36..678870db97 100644 --- a/proteus/LinearAlgebraTools.py +++ b/proteus/LinearAlgebraTools.py @@ -106,14 +106,14 @@ def petsc_load_IS(filename): def csr_2_petsc(size,csr): """ Create an petsc4py matrix from size and CSR information. - Parameters: + Parameters ---------- size : tuple A 2-tuple with the number of matrix rows and columns. csr : tuple A 3-tuple with the sparse matrix csr information. - Returns: + Returns -------- matrix : PETSc4py aij matrix """ @@ -365,14 +365,14 @@ def dense_numpy_2_petsc4py(dense_numpy, eps = 1.e-12): def csr_2_petsc_mpiaij(size,csr): """ Create an MPIaij petsc4py matrix from size and CSR information. - Parameters: + Parameters ---------- size : tuple Two entires: (num_rows, num_cols) csr : tuple (row_idx, col_idx, vals) - Returns: + Returns -------- matrix : PETSc4py MPIaij matrix """ @@ -389,11 +389,11 @@ def split_PETSc_Mat(mat): """ Decompose a PETSc matrix into a symmetric and skew-symmetric matrix - Parameters: + Parameters ---------- mat : :class: `PETSc4py Matrix` - Returns: + Returns -------- H : :class: `PETSc4py Matrix` Symmetric (or Hermitian) component of mat @@ -552,6 +552,9 @@ def print_info(cls): class ParMat_petsc4py(p4pyPETSc.Mat): """ Parallel matrix based on petsc4py's wrappers for PETSc. + + Parameters + ---------- ghosted_csr_mat : :class:`proteus.superluWrappers.SparseMatrix` Primary CSR information for the ParMat. par_bs : int diff --git a/proteus/LinearSolvers.py b/proteus/LinearSolvers.py index e68145d8a1..f0f31e8603 100644 --- a/proteus/LinearSolvers.py +++ b/proteus/LinearSolvers.py @@ -1292,6 +1292,7 @@ class DofOrderInfo(object): """Base class for managing DOF ordering information Parameters + ---------- dof_order_type : str This describes the type of dof ordering that will be constructed. Currently supports: 'blocked' @@ -4386,7 +4387,7 @@ def apply_ns(self, Arguments --------- - par_b : :vec:`petsc4py_vec` + par_b : petsc4py vector """ # Check whether a global null space vector for a constant # pressure has been created. If not, create one. @@ -4454,7 +4455,7 @@ def apply_ns(self, Arguments --------- - par_b : :vec:`petsc4py_vec` + par_b : petsc4py vector """ # Check whether a global null space vector for a constant # has been created. If not, create one. diff --git a/proteus/NonlinearSolvers.py b/proteus/NonlinearSolvers.py index 82af700d9e..f14e5ead10 100644 --- a/proteus/NonlinearSolvers.py +++ b/proteus/NonlinearSolvers.py @@ -809,7 +809,7 @@ def solve(self,u,r=None,b=None,par_u=None,par_r=None): class ExplicitLumpedMassMatrix(Newton): """ - This is a fake solver meant to be used with optimized code + This is a fake solver meant to be used with optimized code A simple iterative solver that is Newton's method if you give it the right Jacobian """ @@ -851,7 +851,7 @@ def no_solve(self,u,r=None,b=None,par_u=None,par_r=None): class ExplicitLumpedMassMatrixForRichards(Newton): """ - This is a fake solver meant to be used with optimized code + This is a fake solver meant to be used with optimized code A simple iterative solver that is Newton's method if you give it the right Jacobian """ @@ -877,7 +877,7 @@ def solve(self,u,r=None,b=None,par_u=None,par_r=None): class ExplicitConsistentMassMatrixWithRedistancing(Newton): """ - This is a fake solver meant to be used with optimized code + This is a fake solver meant to be used with optimized code A simple iterative solver that is Newton's method if you give it the right Jacobian """ @@ -982,7 +982,7 @@ def solve(self,u,r=None,b=None,par_u=None,par_r=None): class ExplicitConsistentMassMatrixForVOF(Newton): """ - This is a fake solver meant to be used with optimized code + This is a fake solver meant to be used with optimized code A simple iterative solver that is Newton's method if you give it the right Jacobian """ @@ -1018,7 +1018,7 @@ def solve(self,u,r=None,b=None,par_u=None,par_r=None): class ExplicitConsistentMassMatrixForRichards(Newton): """ - This is a fake solver meant to be used with optimized code + This is a fake solver meant to be used with optimized code A simple iterative solver that is Newton's method if you give it the right Jacobian """ @@ -1054,7 +1054,7 @@ def solve(self,u,r=None,b=None,par_u=None,par_r=None): class NewtonWithL2ProjectionForMassCorrection(Newton): """ - This is a fake solver meant to be used with optimized code + This is a fake solver meant to be used with optimized code A simple iterative solver that is Newton's method if you give it the right Jacobian """ diff --git a/proteus/Norms.py b/proteus/Norms.py index d2c0180a3f..c935425955 100644 --- a/proteus/Norms.py +++ b/proteus/Norms.py @@ -1,8 +1,5 @@ """ Tools for calculating norms on function spaces. - -.. inheritance-diagram:: proteus.Norms - :parts: 1 """ from .FemTools import * from .Quadrature import * diff --git a/proteus/NumericalFlux.py b/proteus/NumericalFlux.py index 767b7432b3..8782b4eb81 100644 --- a/proteus/NumericalFlux.py +++ b/proteus/NumericalFlux.py @@ -4511,7 +4511,7 @@ class RusanovNumericalFlux(RusanovNumericalFlux_Diagonal): .. math:: - \bar{\lamda}= max |\vec f^{\prime} \cdot n\| + \bar{\lambda}= max |\vec f^{\prime} \cdot n\| Then we just apply apply numerical flux diff --git a/proteus/PostProcessingTools.py b/proteus/PostProcessingTools.py index 23bc166521..0fc0464d90 100644 --- a/proteus/PostProcessingTools.py +++ b/proteus/PostProcessingTools.py @@ -519,7 +519,7 @@ def evaluateElementVelocityField(self,x,ci): class VPP_PWL_RT0(VelocityPostProcessingAlgorithmBase): - """Local Larson-Niklasson method with RT_0 representation for element + r"""Local Larson-Niklasson method with RT_0 representation for element velocities This is not the correct postprocessing for optimized NS codes @@ -531,7 +531,7 @@ class VPP_PWL_RT0(VelocityPostProcessingAlgorithmBase): where :math:`p_i` is the vertex across from face i, :math:`|E|` is the volume of the element, and d is the space dimension the degrees of freedom - are :math:`V^i = \\int_{e_i}\vec v\\dot n_{i}\\ds` + are :math:`V^i = \int_{e_i}\vec v\dot n_{i}\ds` """ # TODO: @@ -1263,7 +1263,7 @@ def evaluateElementVelocityField(self,x,ci): class VPP_PWL_BDM(VPP_PWL_RT0): - """Local Larson-Niklasson method with BDM1 representation for element + r"""Local Larson-Niklasson method with BDM1 representation for element velocities This is not the correct postprocessing for optimized NS codes @@ -1274,7 +1274,7 @@ class VPP_PWL_BDM(VPP_PWL_RT0): .. math:: - \vec N_i = \\lambda_{i/d}\vec e_{i%d} + \vec N_i = \lambda_{i/d}\vec e_{i \bmod d} That is the dofs are locally (say in 2d) :math:`[v^x_0,v^y_0,v^x_1,v^y_1,v^x_2,v^y_2]` @@ -1980,7 +1980,7 @@ def evaluateElementVelocityField(self,x,ci): return vx class VPP_PWL_BDM_OPT(VPP_PWL_RT0_OPT): - """Local Larson-Niklasson method with BDM1 representation for element + r"""Local Larson-Niklasson method with BDM1 representation for element velocities for optimized parallel codes Only difference from VPP_PWL_RT0_OPT should be steps for local @@ -1990,7 +1990,7 @@ class VPP_PWL_BDM_OPT(VPP_PWL_RT0_OPT): .. math:: - \vec N_i = \\lambda_{i%(d+1)} \vec e_{i/(d+1)} + \vec N_i = \lambda_{i \bmod (d+1)} \vec e_{i/(d+1)} Have to use BDM projection to get degrees of freedom @@ -2109,7 +2109,7 @@ def evaluateElementVelocityField(self,x,ci): class VPP_PWC_RT0(VelocityPostProcessingAlgorithmBase): - """Global piecewise constant Larson-Niklasson method with :math:`RT_0` + r"""Global piecewise constant Larson-Niklasson method with :math:`RT_0` representation for element velocities uses RT_0 representation in terms of local basis @@ -2121,7 +2121,7 @@ class VPP_PWC_RT0(VelocityPostProcessingAlgorithmBase): where :math:`p_i` is the vertex across from face :math:`i`, :math:`|E|` is the volume of the element, and :math:`d` is the space dimension the degrees of freedom are :math:`V^i = - \\int_{e_i}\vec v\\dot n_{i}\\ds` + \int_{e_i}\vec v\dot n_{i}\ds` """ # TOOD: @@ -2327,7 +2327,7 @@ def evaluateElementVelocityField(self,x,ci): class VPP_SUN_RT0(VelocityPostProcessingAlgorithmBase): - """Global piecewise constant Sun-Wheeler method with :math:`RT_0` + r"""Global piecewise constant Sun-Wheeler method with :math:`RT_0` representation for element velocities uses RT_0 representation in terms of local basis @@ -2338,7 +2338,7 @@ class VPP_SUN_RT0(VelocityPostProcessingAlgorithmBase): where :math:`p_i` is the vertex across from face i, :math:`|E|` is the volume of the element, and d is the space dimension the - degrees of freedom are :math:`V^i = \\int_{e_i}\vec v\\dot n_{i}\\ds` + degrees of freedom are :math:`V^i = \int_{e_i}\vec v\dot n_{i}\ds` """ # TOOD: @@ -2588,7 +2588,7 @@ def evaluateElementVelocityField(self,x,ci): class VPP_SUN_GS_RT0(VelocityPostProcessingAlgorithmBase): - """Local Sun-Wheeler "Gauss-Seidel" method with :math:`RT_0` representation + r"""Local Sun-Wheeler "Gauss-Seidel" method with :math:`RT_0` representation for element velocities uses :math:`RT_0` representation in terms of local basis @@ -2599,7 +2599,7 @@ class VPP_SUN_GS_RT0(VelocityPostProcessingAlgorithmBase): where p_i is the vertex across from face i, :math:`|E|` is the volume of the element, and d is the space dimension the degrees of - freedom are :math:`V^i = \\int_{e_i}\vec v\\dot n_{i}\\ds` + freedom are :math:`V^i = \int_{e_i}\vec v\dot n_{i}\ds` """ def __init__(self,vectorTransport=None,vtComponents=[0]): @@ -2853,14 +2853,14 @@ def evaluateElementVelocityField(self,x,ci): return vx class VPP_DG_BDM(VPP_DG_RT0): - """Project DG solution onto BDM1 local representation which is + r"""Project DG solution onto BDM1 local representation which is :math:`[P^1]^d` locally with continuous linear fluxes on each face use standard basis .. math:: - \vec N_i = \\lambda_{i%(d+1)} \vec e_{i/(d+1)} + \vec N_i = \lambda_{i \bmod (d+1)} \vec e_{i/(d+1)} Have to use BDM projection to get degrees of freedom @@ -3018,12 +3018,12 @@ def __init__(self,vectorTransport=None,vtComponents=[0]): vtComponents=vtComponents) def postprocess_component(self,ci,verbose=0): - """compute velocity field by just using point evaluation of trial + r"""compute velocity field by just using point evaluation of trial solution and coefficients .. math:: - v = -\ten{a}_h\\grad \\phi_h + \vec f_h + v = -\ten{a}_h\grad \phi_h + \vec f_h """ # TODO @@ -4996,13 +4996,13 @@ def getConservationJacobianPWL_opt(self,ci): self.vt.ebq_global['n'], self.nodeStarFactors[ci]) def postprocessPointEval(self,ci,verbose=0): - """ + r""" compute velocity field by just using point evaluation of trial solution and coefficients .. math:: - v = -\ten{a}_h\\grad \\phi_h + \vec f_h + v = -\ten{a}_h\grad \phi_h + \vec f_h """ self.q[('velocity'),ci][:]=0.0 diff --git a/proteus/Profiling.py b/proteus/Profiling.py index 8c576efafc..9bf2821802 100644 --- a/proteus/Profiling.py +++ b/proteus/Profiling.py @@ -193,7 +193,7 @@ def __call__(self, func, func_args, func_kwargs, profile_name=''): def profile_function(self, func, args, kwargs, profile_name): """ - Profile a Proteus function using the call: func(*args, **kwargs) + Profile a Proteus function using the call: ``func(*args, **kwargs)`` Returns the output of the function call. """ diff --git a/proteus/RefUtils.py b/proteus/RefUtils.py index 5253363a1d..b61ca22256 100644 --- a/proteus/RefUtils.py +++ b/proteus/RefUtils.py @@ -2,9 +2,6 @@ """ Classes, functions, and some global data that are useful for doing FEM calculations on reference elements, etc - -.. inheritance-diagram:: proteus.RefUtils - :parts: 1 """ from .EGeometry import * from .Quadrature import * diff --git a/proteus/SimTools.py b/proteus/SimTools.py index c6d5b8c16b..849ec69e01 100644 --- a/proteus/SimTools.py +++ b/proteus/SimTools.py @@ -1402,11 +1402,13 @@ def getTensorExteriorElementBoundaryStorageKeys(self,mlvt,tsim): return tensorElementStorageKeys def stepPlotElementQuantitiesEnsight(self,mlvt,tsim): """ - sort through desired quantities in quadrature dictionaries like m, dm, to plot + sort through desired quantities in quadrature dictionaries like m, dm, to plot:: + p --- problem definition n --- numerics definition mlvt --- multilevel vector transport that holds the quantities to measure tsim --- simulation time + assumes this is the correct time to plot and plotOffSet is set correctly """ @@ -1464,11 +1466,13 @@ def stepStoreQuantities(self,mlvt,tsim): if multilevelModel is set, stores basically everything - Input + Input:: + p --- problem definition n --- numerics definition mlvt --- multilevel vector transport that holds the quantities to measure tsim --- simulation time + assumes this is the correct time to store TODO: add option for storage directory @@ -1586,11 +1590,13 @@ def computeNodalQuadratureInfo(self,mlvt,t): # def stepPlotEnsight(self,mlvt,tsim): """ - plot solution and solution 'velocity' for ensight/paraview + plot solution and solution 'velocity' for ensight/paraview:: + p --- problem definition n --- numerics definition mlvt --- multilevel vector transport that holds the quantities to measure tsim --- simulation time + assumes this is the correct time to plot and plotOffSet is set correctly @@ -1696,17 +1702,18 @@ def writeVectorElementFunctionHeaderEnsight(self,ckey,filename,append=False,firs def plotScalarElementQuantityEnsight(self,ckey,mlvt,tsim): """ - Ensight plotting routine to look at scalar quantity stored in element quad dictionary q + Ensight plotting routine to look at scalar quantity stored in element quad dictionary q:: + ckey --- what should be plotted p --- problem definition n --- numerics definition mlvt --- multilevel vector transport that holds the quantities to measure tsim --- simulation time + assumes this is the correct time to plot and plotOffSet is set correctly assumes that initial case,sos, geo files set somewhere else - """ p = self.pFile; n = self.nFile @@ -1731,12 +1738,14 @@ def plotScalarElementQuantityEnsight(self,ckey,mlvt,tsim): self.stepPlotCalled['ensightElementQuantities']= True def plotVectorElementQuantityEnsight(self,ckey,mlvt,tsim,scaleOutput=None): """ - Ensight plotting routine to look at vector quantity stored in element quad dictionary q + Ensight plotting routine to look at vector quantity stored in element quad dictionary q:: + ckey --- what should be plotted p --- problem definition n --- numerics definition mlvt --- multilevel vector transport that holds the quantities to measure tsim --- simulation time + assumes this is the correct time to plot and plotOffSet is set correctly @@ -1773,17 +1782,18 @@ def plotVectorElementQuantityEnsight(self,ckey,mlvt,tsim,scaleOutput=None): def plotScalarGlobalElementBoundaryQuantityEnsight(self,ckey,mlvt,tsim): """ - Ensight plotting routine to look at scalar quantity stored in element quad dictionary ebq_global + Ensight plotting routine to look at scalar quantity stored in element quad dictionary ebq_global:: + ckey --- what should be plotted p --- problem definition n --- numerics definition mlvt --- multilevel vector transport that holds the quantities to measure tsim --- simulation time + assumes this is the correct time to plot and plotOffSet is set correctly assumes that initial case,sos, geo files set somewhere else - """ p = self.pFile; n = self.nFile @@ -1810,12 +1820,14 @@ def plotScalarGlobalElementBoundaryQuantityEnsight(self,ckey,mlvt,tsim): self.stepPlotCalled['ensightElementQuantities']= True def plotVectorGlobalElementBoundaryQuantityEnsight(self,ckey,mlvt,tsim,scaleOutput=None): """ - Ensight plotting routine to look at vector quantity stored in element quad dictionary ebq_global + Ensight plotting routine to look at vector quantity stored in element quad dictionary ebq_global:: + ckey --- what should be plotted p --- problem definition n --- numerics definition mlvt --- multilevel vector transport that holds the quantities to measure tsim --- simulation time + assumes this is the correct time to plot and plotOffSet is set correctly diff --git a/proteus/StepControl.py b/proteus/StepControl.py index 01c7d37cf0..9f4ee4cad3 100644 --- a/proteus/StepControl.py +++ b/proteus/StepControl.py @@ -854,24 +854,23 @@ def retryStep_errorFailure(self): return retry class HeuristicNL_dt_controller(SC_base): - """ + r""" Classical Heuristic step controller that picks time step based on threshholds in nonlinear - solver iterations - - if nnl < nonlinearIterationsFloor: - dt *= dtNLgrowFactor - else if nnl > nonlinearIterationsCeil: - dt *= dtNLreduceFactor - end + solver iterations:: - if the nonlinear solver fails, the time step is modified using + if nnl < nonlinearIterationsFloor: + dt *= dtNLgrowFactor + else if nnl > nonlinearIterationsCeil: + dt *= dtNLreduceFactor + end - dt *= dtNLfailureReduceFactor + if the nonlinear solver fails, the time step is modified using:: - Also includes simple linear predictor for initial guess + dt *= dtNLfailureReduceFactor - y^{n+1,p} = y^{n} + (y^{n}-y^{n-1})/(\\Delta t^{n})(t - t^{n-1}) + Also includes simple linear predictor for initial guess:: + y^{n+1,p} = y^{n} + (y^{n}-y^{n-1})/(\Delta t^{n})(t - t^{n-1}) TODO: Implementation: @@ -977,39 +976,40 @@ def errorFailure(self): return False class GustafssonFullNewton_dt_controller(SC_base): - """ + r""" Try version of basic Gustafsson and Soederlind 97 time step selection strategy that accounts for nonlinear solver performance assuming a full Newton nonlinear solver Also includes error control based on classical "deadbeat" control - Right now, decisions based on finest level solve + Right now, decisions based on finest level solve:: - input: dt_prev - output: dt + input: dt_prev + output: dt - get time step estimate based on temporal error --> dt_e - get convergence rate estimate from nonlinear solver --> a - get number of iterations from nonlinear solver --> nnl + get time step estimate based on temporal error --> dt_e + get convergence rate estimate from nonlinear solver --> a + get number of iterations from nonlinear solver --> nnl - if nonlinear solver converges + if nonlinear solver converges - r_a = phi(a_ref/a) - r_a = min(r_a_max,max(r_a,r_a_min)) + r_a = phi(a_ref/a) + r_a = min(r_a_max,max(r_a,r_a_min)) - else + else - if a_ref < a #convergence rate ok but took too many iterations anyway - r_a = phi(nnl_ref/nnl) - else - r_a = phi(a_ref/a) - # - r_a = min(r_a_max,max(r_a,r_a_min)) - # + if a_ref < a #convergence rate ok but took too many iterations anyway + r_a = phi(nnl_ref/nnl) + else + r_a = phi(a_ref/a) + # + r_a = min(r_a_max,max(r_a,r_a_min)) + # - dt = min(dt_e,r_a dt_prev) + dt = min(dt_e,r_a dt_prev) Here, + a_ref -- target convergence rate nnl_ref -- target number of nonlinear iterations r_a_max -- max growth rate @@ -1017,10 +1017,9 @@ class GustafssonFullNewton_dt_controller(SC_base): phi -- limiter function, defaut is phi(x) = x - Also includes simple linear predictor for initial guess - - y^{n+1,p} = y^{n} + (y^{n}-y^{n-1})/(\\Delta t^{n})(t - t^{n-1}) + Also includes simple linear predictor for initial guess:: + y^{n+1,p} = y^{n} + (y^{n}-y^{n-1})/(\Delta t^{n})(t - t^{n-1}) TODO: Implementation: diff --git a/proteus/SubgridError.py b/proteus/SubgridError.py index e24c80a528..b758606539 100644 --- a/proteus/SubgridError.py +++ b/proteus/SubgridError.py @@ -1679,16 +1679,18 @@ def calculateSubgridError(self,q): class AdvectionDiffusionReactionTransientSubscales_ASGS(AdvectionDiffusionReaction_ASGS): - """ + r""" track subgrid scales in time with Backward Euler - \\delta u^{n+1} = -\tau_t\tilde{R}_h + .. math:: + + \delta u^{n+1} &= -\tau_t\tilde{R}_h - \tilde{R}_h = R_h - m^{\\prime,k}\frac{\\delta u^{n}}{\\Delta t^{n+1}} + \tilde{R}_h &= R_h - m^{\prime,k}\frac{\delta u^{n}}{\Delta t^{n+1}} - \tau_t = \frac{\\Delta t^{n+1}\tau_s}{m^{prime,n+1}\tau_s + \\Delta t^{n+1}} + \tau_t &= \frac{\Delta t^{n+1}\tau_s}{m^{prime,n+1}\tau_s + \Delta t^{n+1}} - \tau_s = normal spatial tau, supposed to have \tau_s \approx \\mathcal{L}^{-1}_{s} + \tau_s &= \text{normal spatial tau, supposed to have } \tau_s \approx \mathcal{L}^{-1}_{s} for now m^{prime} evaluated at k=n for subgrid error but not sure if this is right or not diff --git a/proteus/SubsurfaceTransportCoefficients.py b/proteus/SubsurfaceTransportCoefficients.py index 08eb778c9e..6253a8c2ac 100644 --- a/proteus/SubsurfaceTransportCoefficients.py +++ b/proteus/SubsurfaceTransportCoefficients.py @@ -1209,12 +1209,12 @@ def setMaterialTypes(self, # class TwophaseDarcy_fc(TwophaseDarcyFlow_base): - """continuity equation for each phase + r"""continuity equation for each phase .. math:: - \\pd{m_w}{t} - \\deld (\ten{a}_w \\grad \\phi_w) + r_w = 0 - \\pd{m_n}{t} - \\deld (\ten{a}_n \\grad \\phi_n) + r_n = 0 + \pd{m_w}{t} - \deld (\ten{a}_w \grad \phi_w) + r_w = 0 + \pd{m_n}{t} - \deld (\ten{a}_n \grad \phi_n) + r_n = 0 """ # (normalized) mass for each phase @@ -1587,12 +1587,12 @@ def __init__(self, class TwophaseDarcy_fc_pp(TwophaseDarcyFlow_base): - """continuity equation for each phase + r"""continuity equation for each phase .. math:: - \\pd{m_w}{t} - \\deld (\ten{a}_w \\grad \\phi_w) + r_w = 0 - \\pd{m_n}{t} - \\deld (\ten{a}_n \\grad \\phi_n) + r_n = 0 + \pd{m_w}{t} - \deld (\ten{a}_w \grad \phi_w) + r_w = 0 + \pd{m_n}{t} - \deld (\ten{a}_n \grad \phi_n) + r_n = 0 """ diff --git a/proteus/TimeIntegration.py b/proteus/TimeIntegration.py index 44604f6031..2010a03195 100644 --- a/proteus/TimeIntegration.py +++ b/proteus/TimeIntegration.py @@ -11,16 +11,16 @@ from .Comm import globalMax class TI_base(object): - """ + r""" The base class for time discretizations of tightly coupled systems of transport equations. The purpose of time integration objects is to construct - approximations: + approximations:: - m_t, dm_t -- the time derivatives and their derivatives w.r.t. u - u_*, du_* -- the solution variables and their derivatives w.r.t. u - f_*, df_*,... --the coefficients and their derivatives w.r.t. u + m_t, dm_t -- the time derivatives and their derivatives w.r.t. u + u_*, du_* -- the solution variables and their derivatives w.r.t. u + f_*, df_*,... --the coefficients and their derivatives w.r.t. u The objects provide a set_dt function for setting the time step as well as functions for SUGGESTING the time step. The StepController @@ -37,7 +37,9 @@ class TI_base(object): method and results in the steady state scalar transport equation: - \\deld ( f - a \\grad phi) + r = 0 + .. math:: + + \deld ( f - a \grad \phi ) + r = 0 NoIntegration is an alias for this class. """ @@ -1513,7 +1515,7 @@ def postAdaptUpdate(self,oldTime): self.dt_history[:] = oldTime.dt_history[:] def computeErrorEstimate(self): - """calculate :math:`\vec{e}^{n+1}` + r"""calculate :math:`\vec{e}^{n+1}` depending on order of approximation To be consistent, this must be called after step taken but before update time history @@ -1521,17 +1523,16 @@ def computeErrorEstimate(self): Initially, use .. math:: - + (\vec y^{n+1}-\vec y^{n+1,p})/2 for first order and - .. math: - + .. math:: + r/(1+r)(\vec y^{n+1}-(1+r)\vec y^{n}+r\vec y^{n-1}) for second order. - """ #mwf debug #import pdb @@ -1612,38 +1613,38 @@ def setFromOptions(self,nOptions): #mwf hack for cfl calculations from . import cfemIntegrals class ExplicitRK_base(TI_base): - """ + r""" base class for explicit RK schemes which can be written in the Shu Osher 88 (Cockburn Shu 89) style generic formula u^0 = u^n - u^l = \\sum_{k=0}^{l-1}\alpha_{l,k}u^k + \beta_{l,k} \\Delta t L(u^{k},t^n + d_{k} \\Delta t^n) + u^l = \sum_{k=0}^{l-1}\alpha_{l,k}u^k + \beta_{l,k} \Delta t L(u^{k},t^n + d_{k} \Delta t^n) u^{n+1} = u^{m} We will try to include these schemes in a Rothe framework rather than MOL so in calculate element coefficients, at stage number l, u_t will be replaced with - u_t --> (u - \\sum_{k=0}^{l-1} \alpha_{l,k} u^k)/\\Delta t + u_t --> (u - \sum_{k=0}^{l-1} \alpha_{l,k} u^k)/\Delta t or - u_t --> (u - \\sum_{k=0}^{l-1} \alpha_{l,k} u^k) + u_t --> (u - \sum_{k=0}^{l-1} \alpha_{l,k} u^k) - depending on whether \\Delta t is assumed to be in the mass matrix or not + depending on whether \Delta t is assumed to be in the mass matrix or not The spatial residual at stage l will be replaced with (in calculateElementSpatialResidual) - r --> \\sum_{k=0,l-1} \beta_{l,k} r_k + r --> \sum_{k=0,l-1} \beta_{l,k} r_k or - r --> \\Delta t \\sum_{k=0,l-1} \beta_{l,k} r_k + r --> \Delta t \sum_{k=0,l-1} \beta_{l,k} r_k - depending on whether \\Delta t is assumed to be in the mass matrix or not + depending on whether \Delta t is assumed to be in the mass matrix or not Here r_k corresponds to - L(u^{k},t^n + d_{k} \\Delta t^n) + L(u^{k},t^n + d_{k} \Delta t^n) Substeps in NumericalSolution framework correspond to stage values @@ -1779,12 +1780,12 @@ def resetOrder(self,order,stages): #end ci # def setCoefficients(self): - """set alpha,beta, d for + r"""set alpha,beta, d for .. math:: u^0 = u^n - u^l = \\sum_{k=0}^{l-1}\alpha_{l,k}u^k + \beta_{l,k} \\Delta t L(u^{k},t^n + d_{k} \\Delta t^n) + u^l = \sum_{k=0}^{l-1}\alpha_{l,k}u^k + \beta_{l,k} \Delta t L(u^{k},t^n + d_{k} \Delta t^n) u^{n+1} = u^{m} must be implemented in derived class @@ -1797,7 +1798,7 @@ def setCoefficients(self): beta[l,k] --> beta_{l+1,k} dcoefs[k] --> d_{k+1} - because of the whole caching\\delayed eval deal second index + because of the whole caching\delayed eval deal second index (k) is actual level """ @@ -1816,10 +1817,10 @@ def updateStage(self): if self.lstage < 0 or self.lstage > self.nStages: logEvent("""stage= %d out of allowed bounds [0,%d]""" % (self.lstage,self.nStages)) def calculateElementCoefficients(self,q): - """ + r""" set m_t as - m_t --> (m - \\sum_{k=0}^{l-1} \alpha_{l,k} m^k)/\\Delta t + m_t --> (m - \sum_{k=0}^{l-1} \alpha_{l,k} m^k)/\Delta t """ if self.dt <= 1.0e-24: @@ -1853,10 +1854,10 @@ def calculateElementCoefficients(self,q): #end if q has key #end ci def calculateElementCoefficientsNoDtInMass(self,q): - """ + r""" set m_t as - m_t --> (m - \\sum_{k=0}^{l-1} \alpha_{l,k} m^k) + m_t --> (m - \sum_{k=0}^{l-1} \alpha_{l,k} m^k) """ assert self.dt > 1.0e-24,"dt = %d too small in calculateElementCoefficients" % self.dt @@ -1878,10 +1879,10 @@ def calculateElementCoefficientsNoDtInMass(self,q): #end ci def calculateElementSpatialResidual(self,elementResidual): - """ + r""" The spatial residual at stage l will be replaced with (in calculateElementSpatialResidual) - r --> \\sum_{k=0,l-1} \beta_{l,k} r_k + r --> \sum_{k=0,l-1} \beta_{l,k} r_k """ @@ -1898,10 +1899,10 @@ def calculateElementSpatialResidual(self,elementResidual): #end ci #end calculateElementSpatialResidual def calculateElementSpatialResidualNoDtInMass(self,elementResidual): - """ + r""" The spatial residual at stage l will be replaced with (in calculateElementSpatialResidual) - r --> \\Delta t \\sum_{k=0,l-1} \beta_{l,k} r_k + r --> \Delta t \sum_{k=0,l-1} \beta_{l,k} r_k """ @@ -2110,13 +2111,13 @@ def setFromOptions(self,nOptions): class LinearSSPRKintegration(ExplicitRK_base): - """ + r""" implementation of Explicit RK schemes for ODE that are SSP for linear operators (aka) spatial discretizations in .. math:: - \\od{\vec y}{t} = \\mat{L}\vec y + \od{\vec y}{t} = \mat{L}\vec y See Gottlieb, Shu, Tadmor siam review article and notes @@ -2125,22 +2126,21 @@ class LinearSSPRKintegration(ExplicitRK_base): .. math:: u^0 = u^n - u^l = u^{l-1} + \\Delta t L(u^{l-1}) l = 1,...,m-1 - u^m = \\sum_{k=0}^{l-1}\alpha_{m,k}u^k + - \alpha_{m,m-1} \\Delta t L(u^{m-1}) + u^l = u^{l-1} + \Delta t L(u^{l-1}) l = 1,...,m-1 + u^m = \sum_{k=0}^{l-1}\alpha_{m,k}u^k + + \alpha_{m,m-1} \Delta t L(u^{m-1}) u^{n+1} = u^m so \alpha_{l,l-1} = 1, \beta_{l,l-1} = 1.0 l < m, - \beta_{m,m-1} = \alph_{m,m-1} + \beta_{m,m-1} = \alpha_{m,m-1} Apparently, .. math:: - d_{l} = l \\Delta t l < m, + d_{l} = l \Delta t l < m, d_m = 1.0 - - which doesn't make a lot of sense for time dependent problems + which doesn't make a lot of sense for time dependent problems """ def __init__(self,transport,order=1,runCFL=0.1,usingSSPRKNewton=False): ExplicitRK_base.__init__(self,transport,timeOrder=order,nStages=order,runCFL=runCFL, @@ -2148,26 +2148,29 @@ def __init__(self,transport,order=1,runCFL=0.1,usingSSPRKNewton=False): def setCoefficients(self): - """ - alpha matrix in Gottlieb, Shu, Tadmor review paper + r""" + alpha matrix in Gottlieb, Shu, Tadmor review paper:: - u^{l} = \\sum_{k=0}^{l-1}\alpha_{l,k}u^k + u^{l} = \sum_{k=0}^{l-1}\alpha_{l,k}u^k so \alpha_{l,l-1} = 1 for l < m - beta matrix multiplies - \beta_{l,k} \\Delta t L(u^{k},t^n + d_{k} \\Delta t^n) + beta matrix multiplies:: + + \beta_{l,k} \Delta t L(u^{k},t^n + d_{k} \Delta t^n) + + so:: - so \beta_{l,l-1} = 1.0 for l < m \beta_{m,m-1} = \alpha_{m,m-1} - Note that + Note that:: + alpha[l,k] --> alpha_{l+1,k} beta[l,k] --> beta_{l+1,k} - dcoefs[k] --> d_{k+1} because of the whole caching\\ delayed eval deal - second index (k) is actual level + dcoefs[k] --> d_{k+1} because of the whole caching\ delayed eval deal + second index (k) is actual level """ order = self.timeOrder; stages = self.nStages assert order == stages, "LinearSSPRK problem order (%s) != stages (%s) " % (order,stages) @@ -2212,46 +2215,52 @@ def setFromOptions(self,nOptions): class SSPRKPIintegration(ExplicitRK_base): - """So called SSP RK integration with limiting step at end of each + r"""So called SSP RK integration with limiting step at end of each phase This one holds the original TVD RK schemes up to 3rd order 1st order is Forward Euler 2nd order .. math:: - u^1 = u^n + \\Delta t L(u^n) - u^2 = 0.5(u^n + u^1) + 0.5\\Delta t L(u^1,t^n + \\Delta t^n), same stage values etc linear SSPRK + + u^1 = u^n + \Delta t L(u^n) + + u^2 = 0.5(u^n + u^1) + 0.5\Delta t L(u^1,t^n + \Delta t^n), same stage values etc linear SSPRK 3rd order .. math:: - u^1 = u^n + \\Delta t L(u^n,t^n) - u^2 = 0.75 u^n + 0.25u^1 + 0.25\\Delta t L(u^1,t^n + \\Delta t^n), - u^3 = 1/3 u^n + 2/3 u^2 + 2/3 \\Delta t L(u^2,t^n + 0.5\\Delta t^n) + + u^1 = u^n + \Delta t L(u^n,t^n) + + u^2 = 0.75 u^n + 0.25u^1 + 0.25\Delta t L(u^1,t^n + \Delta t^n), + + u^3 = 1/3 u^n + 2/3 u^2 + 2/3 \Delta t L(u^2,t^n + 0.5\Delta t^n) generic formula is .. math:: u^0 = u^n - u^l = \\sum_{k=0}^{l-1}\alpha_{l,k}u^k + \beta_l \\Delta t L(u^{l-1},t^n + d_{l-1} \\Delta t^n) + + u^l = \sum_{k=0}^{l-1}\alpha_{l,k}u^k + \beta_l \Delta t L(u^{l-1},t^n + d_{l-1} \Delta t^n) + u^{n+1} = u^{m} - so :math:`beta_l --> beta_{l,k-1}` in our basic framework + so :math:`beta_l \to beta_{l,k-1}` in our basic framework Try to put this into Rothe paradigm with linear implicit mass matrix evaluation by saying at stage l - + .. math:: - m_t \approx \frac{u^l - \\sum_{k=0}^{l-1}\alpha_{l,k}u^k}{\\Delta t} + m_t \approx \frac{u^l - \sum_{k=0}^{l-1}\alpha_{l,k}u^k}{\Delta t} spatial residual --> spatial residual * beta_l and solve system for :math:`u^l` The limiting right now assumes same space for each component - """ def __init__(self,transport,order=1,runCFL=0.1,limiterType=None,usingSSPRKNewton=False): ExplicitRK_base.__init__(self,transport,timeOrder=order,nStages=order,runCFL=runCFL, diff --git a/proteus/TransportCoefficients.py b/proteus/TransportCoefficients.py index d8ff6d40d6..991daa88a6 100644 --- a/proteus/TransportCoefficients.py +++ b/proteus/TransportCoefficients.py @@ -6602,8 +6602,8 @@ class kEpsilon(TC_base): .. math:: \pd{\bar{\vec v}}{t} + \deld \left(\bar{\vec v} \outer \bar{\vec v}\right) - -\nu \deld \ten \bar{D} + \frac{1}{\rho}\grad \bar p - - \frac{1}{rho}\deld \ten{R} = 0 + -\nu \deld \ten{\bar{D}} + \frac{1}{\rho}\grad \bar p + - \frac{1}{\rho}\deld \ten{R} = 0 Reynolds stress term @@ -6611,7 +6611,7 @@ class kEpsilon(TC_base): \ten R = -\rho <\vec v^{'}\outer \vec v^{'}> \frac{1}{\rho}\ten{R} = 2 \nu_t \bar{D} - \frac{2}{3}k\ten{I} - D_{ij}(\vec v) = \frac{1}{2} \left( \pd{v_i}{x_j} + \pd{v_j}{x_i}) + D_{ij}(\vec v) = \frac{1}{2} \left( \pd{v_i}{x_j} + \pd{v_j}{x_i}\right) \ten D \bar{\ten D} = D(<\vec v>), \ten D^{'} = \ten D(\vec v^{'}) k-epsilon tranport equations @@ -6936,8 +6936,8 @@ class kEpsilon_k(TC_base): \deld \bar{\vec v} = 0 \pd{\bar{\vec v}}{t} + \deld \left(\bar{\vec v} \outer \bar{\vec v}\right) - -\nu \deld \ten \bar{D} + \frac{1}{\rho}\grad \bar p - - \frac{1}{rho}\deld \ten{R} = 0 + -\nu \deld \ten{\bar{D}} + \frac{1}{\rho}\grad \bar p + - \frac{1}{\rho}\deld \ten{R} = 0 Reynolds stress term @@ -6945,7 +6945,7 @@ class kEpsilon_k(TC_base): \ten R = -\rho <\vec v^{'}\outer \vec v^{'}> \frac{1}{\rho}\ten{R} = 2 \nu_t \bar{D} - \frac{2}{3}k\ten{I} - D_{ij}(\vec v) = \frac{1}{2} \left( \pd{v_i}{x_j} + \pd{v_j}{x_i}) + D_{ij}(\vec v) = \frac{1}{2} \left( \pd{v_i}{x_j} + \pd{v_j}{x_i}\right) \ten D \bar{\ten D} = D(<\vec v>), \ten D^{'} = \ten D(\vec v^{'}) k-epsilon tranport equations @@ -7207,8 +7207,8 @@ class kEpsilon_epsilon(TC_base): \deld \bar{\vec v} = 0 \pd{\bar{\vec v}}{t} + \deld \left(\bar{\vec v} \outer \bar{\vec v}\right) - -\nu \deld \ten \bar{D} + \frac{1}{\rho}\grad \bar p - - \frac{1}{rho}\deld \ten{R} = 0 + -\nu \deld \ten{\bar{D}} + \frac{1}{\rho}\grad \bar p + - \frac{1}{\rho}\deld \ten{R} = 0 Reynolds stress term @@ -7216,7 +7216,7 @@ class kEpsilon_epsilon(TC_base): \ten R = -\rho <\vec v^{'}\outer \vec v^{'}> \frac{1}{\rho}\ten{R} = 2 \nu_t \bar{D} - \frac{2}{3}k\ten{I} - D_{ij}(\vec v) = \frac{1}{2} \left( \pd{v_i}{x_j} + \pd{v_j}{x_i}) + D_{ij}(\vec v) = \frac{1}{2} \left( \pd{v_i}{x_j} + \pd{v_j}{x_i}\right) \ten D \bar{\ten D} = D(<\vec v>), \ten D^{'} = \ten D(\vec v^{'}) k-epsilon tranport equations @@ -7990,16 +7990,16 @@ def evaluate(self,t,c): #ci from proteus.ctransportCoefficients import evaluateBuckleyLeverettLiuExample class BuckleyLeverettLiuExample(TC_base): - """ + r""" 5 spot example from Liu 93 Siam paper. - S_t + \\deld (\vec a f(S)) = 0 + S_t + \deld (\vec a f(S)) = 0 f(S) = S^2 / (0.2 - 0.4 S + 1.2 S^2) - \vec a = (\\pd{\\phi}{x},\\pd{phi}{y}) + \vec a = (\pd{\phi}{x},\pd{phi}{y}) - \\phi = 0.01 \\log(\\sqrt{x^2 + y^2}) + \phi = 0.01 \log(\sqrt{x^2 + y^2}) However, note that his example is for S=0 injecting into S=1 background saturation using his definition of f diff --git a/proteus/Viewers.py b/proteus/Viewers.py index b95bab4723..24169f0afe 100644 --- a/proteus/Viewers.py +++ b/proteus/Viewers.py @@ -555,11 +555,13 @@ def stepPlotExact(self,mlvt,tsim): #end def def stepPlotElementQuantities(self,mlvt,tsim): """ - sort through desired quantities in quadrature dictionaries like m, dm, to plot + sort through desired quantities in quadrature dictionaries like m, dm, to plot:: + p --- problem definition n --- numerics definition mlvt --- multilevel vector transport that holds the quantities to measure tsim --- simulation time + assumes this is the correct time to plot and plotOffSet is set correctly """ @@ -929,15 +931,17 @@ def plotVectorGlobalElementBoundaryQuantity(self,ckey,mlvt,tsim): def plotVectorElementQuantity(self,ckey,mlvt,tsim,nVectorPlotPointsPerElement=1): """ plotting routine to look at vector quantity stored in global element quad dictionary q - input : + + input:: + ckey --- what should be plotted p --- problem definition n --- numerics definition mlvt --- multilevel vector transport that holds the quantities to measure tsim --- simulation time + assumes this is the correct time to plot and plotOffSet is set correctly - """ from proteusGraphical import vtkViewers p = self.p; n = self.n @@ -1182,18 +1186,16 @@ def viewGlobalScalarPointData(self,cmdFile,nSpace,q,ckey,name=None, storeMeshData=True,figureNumber=1,title=None): """ input scalar variable and coordinates stored in dictionary - q['x'], q[ckey] + q['x'], q[ckey] respectively, generate global continuous interpolant should work for q, ebq_global, and ebqe quadrature dictionaries uses delaunay triangulation in 2d and 3d - scalar data is stored in - name_q + scalar data is stored in name_q - if storeMeshData = True, writes out - name_x_q -- point data - tri_name_q -- Delaunay representation (2d,3d) + if storeMeshData = True, writes out name_x_q (point data) and + tri_name_q (Delaunay representation, 2d/3d) returns number of figures actually plotted """ @@ -1300,18 +1302,16 @@ def viewGlobalVectorPointData(self,cmdFile,nSpace,q,ckey,name=None, storeMeshData=True,figureNumber=1,title=None): """ input vector variable and coordinates stored in dictionary - q['x'], q[ckey] + q['x'], q[ckey] respectively, generate global continuous interpolant should work for q, ebq_global, and ebqe quadrature dictionaries uses delaunay triangulation in 2d and 3d - vector data is stored in - name_q + vector data is stored in name_q - if storeMeshData = True, writes out - name_x_q -- point data - tri_name_q -- Delaunay representation (2d,3d) + if storeMeshData = True, writes out name_x_q (point data) and + tri_name_q (Delaunay representation, 2d/3d) returns number of figures actually generated """ @@ -1427,18 +1427,16 @@ def viewLocalScalarPointData(self,cmdFile,nSpace,q,ckey,name=None, storeMeshData=True,figureNumber=1,title=None): """ input scalar variable and coordinates stored in dictionary - q['x'], q[ckey] + q['x'], q[ckey] respectively, generate local interpolant that is element-wise continuous should work for q, ebq_global, and ebqe quadrature dictionaries uses delaunay triangulation in 2d and 3d on each element (stored in cell array) - scalar data is stored in - name_q + scalar data is stored in name_q - if storeMeshData = True, writes out - name_x_q -- point data - tri_name_q -- Delaunay representation (2d,3d) + if storeMeshData = True, writes out name_x_q (point data) and + tri_name_q (Delaunay representation, 2d/3d) returns number of figures actually plotted """ diff --git a/proteus/WaveTools.py b/proteus/WaveTools.py index fb3f00dd12..51f1c6405d 100644 --- a/proteus/WaveTools.py +++ b/proteus/WaveTools.py @@ -438,18 +438,19 @@ def vel_mode(x, t, kDir, kAbs, omega, phi, amplitude, mwl, depth, vDir, gA def sigma(omega,omega0): """Calculates sigma function for JONSWAP spectrum - See http://www.wikiwaves.org/Ocean-Wave_Sectra + See https://www.wikiwaves.org/index.php/Ocean-Wave_Spectra + Parameters ---------- omega : numpy.ndarray Angular frequency array omega0 : numpy.ndarray Peak angular frequency + Returns -------- numpy.ndarray 1D Numpy array of simga function with respect to f - """ sigmaReturn = np.where(omega > omega0,0.09,0.07) return sigmaReturn @@ -902,19 +903,20 @@ class MonochromaticWaves(object): Fenton: uses BCoeffs/YCoeffs provided by user autoFenton: bool autoFenton=True: uses waveheight, period, depth, and g to - calculate coeffs + calculate coeffs autoFenton=False: uses BCoeffs/YCoeffs provided by user autoFentonOpts: dict options for autoFenton. The dictionary must contain the following - entries (here the default values if autoFentonOpts is None): - autoFentonOpts = {'mode': 'Period', - 'current_criterion': 1, - 'height_steps': 1, - 'niter': 40, - 'conv_crit': 1e-05, - 'points_freesurface': 50, - 'points_velocity': 16, - 'points_vertical': 20} + entries (here the default values if autoFentonOpts is None):: + + autoFentonOpts = {'mode': 'Period', + 'current_criterion': 1, + 'height_steps': 1, + 'niter': 40, + 'conv_crit': 1e-05, + 'points_freesurface': 50, + 'points_velocity': 16, + 'points_vertical': 20} Ycoeff : numpy.ndarray Fenton Fourier coefficients for free-surface elevation Bcoeff : numpy.ndarray diff --git a/proteus/deim_utils.py b/proteus/deim_utils.py index 4dc46002ed..ab7b282725 100644 --- a/proteus/deim_utils.py +++ b/proteus/deim_utils.py @@ -58,7 +58,7 @@ def generate_svd_decomposition(archive,nsnap,val_name,outbase): return U,s,V def calculate_deim_indices(Uin): - """ + r""" input: Uin n x m array of basis vectors for nonlinear function snapshots output: rho, m vector of indices \rho_i for extracting $\vec F$ values diff --git a/proteus/mbd/CouplingFSI.pyx b/proteus/mbd/CouplingFSI.pyx index 7e411ca25c..580a73ecbc 100644 --- a/proteus/mbd/CouplingFSI.pyx +++ b/proteus/mbd/CouplingFSI.pyx @@ -1179,13 +1179,15 @@ cdef class ProtChBody: Angular velocity of body (default: False. Set to True to record). ang_acc: bool Angular acceleration of body (default: False. Set to True to record). + Notes ----- - To add another value manually, add to dictionary self.record_dict: - key: header of the column in .csv - value: list of length 2: [variable name, index within variable] - (if no index, use None) - e.g. self.record_dict['m']['mass', None] + To add another value manually, add to dictionary self.record_dict:: + + key: header of the column in .csv + value: list of length 2: [variable name, index within variable] + (if no index, use None) + e.g. self.record_dict['m']['mass', None] """ if all_values is True: pos = rot = F = M = acc = vel = ang_acc = ang_vel = h_predict = True @@ -1921,7 +1923,8 @@ cdef class ProtChSystem: def log_bodies_text(self, d_time, l_logging_info): """ Logs the chrono information into a text file at each timestep - Parameters + + Parameters ---------- self: object ProtChSystem being referenced @@ -1931,6 +1934,7 @@ cdef class ProtChSystem: Contains the information to be logged. Structure is [body number, type] with multiple entries being included as a 2d list. Valid types are 'position', 'rotation', 'force', and 'torque'. + Returns ------- None. Data is logged to a text file @@ -1984,6 +1988,7 @@ cdef class ProtChSystem: def log_bodies_h5(self, l_logging_info): """ Logs the chrono information into a h5 file at each timestep + Parameters ---------- self: object @@ -1992,6 +1997,7 @@ cdef class ProtChSystem: Contains the information to be logged. Structure is [body number, type] with multiple entries being included as a 2d list. Valid types are 'position', 'rotation', 'force', and 'torque'. + Returns ------- None. Data is saved to an h5 file. diff --git a/proteus/mprans/BodyDynamics.py b/proteus/mprans/BodyDynamics.py index 3d28052189..ea326ccf6d 100644 --- a/proteus/mprans/BodyDynamics.py +++ b/proteus/mprans/BodyDynamics.py @@ -13,6 +13,7 @@ class RigidBody(AuxiliaryVariables.AV_base, object): class instance acting as a rigid body. To set a shape as a rigid body, use shape.setRigidBody(). The class instance is created automatically when shape.setRigidBody() has been called and after calling assembleDomain(). + Parameters ---------- shape: proteus.mprans.SpatialTools.Shape_RANS @@ -342,6 +343,7 @@ def getAngularDisplacement(self, dt): def setSprings(self, springs, Kx, Ky, Krot, Cx, Cy, Crot, Kz=0.0, Cz=0.0): """ Sets a system of uniform springs to model soil's reactions (for moving bodies) + Parameters ---------- spring: string @@ -425,7 +427,6 @@ def inputMotion(self, InputMotion=False, pivot=None, def imposeSinusoidalMotion(self): """ Motion is imposed rather than calculated. - Parameters """ t = self.model.stepController.t_model_last Tra = np.array([0., 0., 0.]) @@ -453,6 +454,7 @@ def imposeSinusoidalMotion(self): def step(self, dt): """ Step for rigid body calculations in Python + Parameters ---------- dt: float @@ -493,6 +495,7 @@ def step(self, dt): def setConstraints(self, free_x, free_r): """ Sets constraints on the Shape (for moving bodies) + Parameters ---------- free_x: array_like @@ -506,6 +509,7 @@ def setConstraints(self, free_x, free_r): def setMass(self, mass): """ Set mass of the shape. + Parameters ---------- mass: float @@ -516,6 +520,7 @@ def setMass(self, mass): def setInertiaTensor(self, It): """ Set the inertia tensor of the shape + Parameters ---------- It: array_like, float @@ -537,6 +542,7 @@ def setInertiaTensor(self, It): def getInertia(self, vec=(0., 0., 1.), pivot=None): """ Gives the inertia of the shape from an axis and a pivot + Parameters ---------- vec: array_like @@ -586,6 +592,7 @@ def setRecordValues(self, filename=None, all_values=False, pos=False, """ Sets the rigid body attributes that are to be recorded in a csv file during the simulation. + Parameters ---------- filename: Optional[string] @@ -616,13 +623,15 @@ def setRecordValues(self, filename=None, all_values=False, pos=False, Angular velocity of body (default: False. Set to True to record). ang_acc: bool Angular acceleration of body (default: False. Set to True to record). + Notes ----- - To add another value manually, add to dictionary self.record_dict: - key: header of the column in .csv - value: list of length 2: [variable name, index within variable] - (if no index, use None) - e.g. self.record_dict['m']['mass', None] + To add another value manually, add to dictionary self.record_dict:: + + key: header of the column in .csv + value: list of length 2: [variable name, index within variable] + (if no index, use None) + e.g. self.record_dict['m']['mass', None] """ if all_values is True: pos = rot = F = M = acc = vel = ang_acc = ang_vel = True @@ -821,6 +830,7 @@ def _store_last_values(self): def step(self, dt, substeps=20): """ Step for rigid body calculations in Python + Parameters ---------- dt: float @@ -873,6 +883,7 @@ def step(self, dt, substeps=20): def getInertia(self, vec=(0., 0., 1.), pivot=None): """ Gives the inertia of the shape from an axis and a pivot + Parameters ---------- vec: array_like @@ -924,6 +935,7 @@ def getInertia(self, vec=(0., 0., 1.), pivot=None): def setFriction(self, friction, m_static, m_dynamic, tolerance, grainSize): """ Sets material properties for sliding and overturning modules + Parameters ---------- friction: string @@ -949,6 +961,7 @@ def setFriction(self, friction, m_static, m_dynamic, tolerance, grainSize): def setOverturning(self, overturning): """ Sets material properties for sliding and overturning modules + Parameters ---------- overturning: string @@ -959,6 +972,7 @@ def setOverturning(self, overturning): def setSprings(self, springs, Kx, Ky, Krot, Cx, Cy, Crot, Kz=0.0, Cz=0.0): """ Sets a system of uniform springs to model soil's reactions (for moving bodies) + Parameters ---------- spring: string @@ -989,6 +1003,7 @@ def setSprings(self, springs, Kx, Ky, Krot, Cx, Cy, Crot, Kz=0.0, Cz=0.0): def setNumericalScheme(self, scheme): """ Sets the numerical scheme used to solve motion. + Parameters ---------- scheme: string @@ -1000,6 +1015,7 @@ def setNumericalScheme(self, scheme): def friction_module(self, dt): """ Calculate sliding motion modelling frictional force. + Parameters ---------- dt : Time step. @@ -1027,6 +1043,7 @@ def friction_module(self, dt): def dynamic_case(self, sign, Fx, Fv, mass, m): """ Set a dynamic friction. + Parameters ---------- sign : It's function of horizontal force. @@ -1144,6 +1161,7 @@ def dynamic_case(self, sign, Fx, Fv, mass, m): def overturning_module(self, dt): """ Calculate overturning motion modelling soil foundation reactions. + Parameters ---------- dt : Time step. @@ -1156,6 +1174,7 @@ def overturning_module(self, dt): def calculate_rotation(self, floating, h): """ Calculate rotation. + Parameters ---------- floating : enable floating body calculation with NO restrictions offered by friction or springs. @@ -1255,6 +1274,7 @@ def setRecordValues(self, filename=None, all_values=False, pos=False, """ Sets the rigid body attributes that are to be recorded in a csv file during the simulation. + Parameters ---------- filename: Optional[string] @@ -1285,13 +1305,15 @@ def setRecordValues(self, filename=None, all_values=False, pos=False, Angular velocity of body (default: False. Set to True to record). ang_acc: bool Angular acceleration of body (default: False. Set to True to record). + Notes ----- - To add another value manually, add to dictionary self.record_dict: - key: header of the column in .csv - value: list of length 2: [variable name, index within variable] - (if no index, use None) - e.g. self.record_dict['m']['mass', None] + To add another value manually, add to dictionary self.record_dict:: + + key: header of the column in .csv + value: list of length 2: [variable name, index within variable] + (if no index, use None) + e.g. self.record_dict['m']['mass', None] """ if all_values is True: pos = rot = F = M = acc = vel = ang_acc = ang_vel = elasticPlastic = True @@ -1407,6 +1429,7 @@ def _store_last_values(self): def step(self, dt, substeps=20): """ Step for rigid body calculations in Python + Parameters ---------- dt: float @@ -1492,12 +1515,13 @@ def inputMotion(self, InputMotion=False, pivot=None, def imposeSinusoidalMotion(self, tt = None): """ Motion is imposed rather than calculated. - ---------- + + Parameters + ---------- tt: None or float - tt should be None for other than testing purposes. If not None, + tt should be None for other than testing purposes. If not None, time is externally set. - - """ + """ if tt is not None: t= tt else: @@ -1560,6 +1584,7 @@ def getEulerAngles(coord_system): def runge_kutta(u0, v0, a0, dt, substeps, F, K, C, m, velCheck): """ Function that applies Runge Kutta's scheme for motion calculation. + Parameters ---------- u0 : translational or rotational displacement. diff --git a/proteus/mprans/BoundaryConditions.py b/proteus/mprans/BoundaryConditions.py index 55f52834d5..35652bedd7 100644 --- a/proteus/mprans/BoundaryConditions.py +++ b/proteus/mprans/BoundaryConditions.py @@ -322,6 +322,7 @@ def setAtmosphere(self, orientation=None, vof_air=1.,kInflow=None,dInflow=None): """ Sets atmosphere boundary conditions (water can come out) (!) pressure dirichlet set to 0 for this BC + Parameters ---------- orientation: Optional[array_like] @@ -454,10 +455,10 @@ def setWallFunction(self, wall, shearStress=False): At the moment version with shearStress=False is the only one that returns good results. Keep it False at the moment! + - When True, the wall function prescribes diffusive boundaries for velocity and kappa. It's like imposing the shear stress. - If False, the wall function prescribes dirichlet conditions. - """ wf = wall self.reset() @@ -643,6 +644,7 @@ def setTwoPhaseVelocityInlet(self, U, waterLevel, smoothing, Uwind=None, """ Imposes a velocity profile lower than the sea level and an open boundary for higher than the sealevel. + Parameters ---------- U: list. @@ -653,7 +655,7 @@ def setTwoPhaseVelocityInlet(self, U, waterLevel, smoothing, Uwind=None, water level at global coordinate system. smoothing: float. range within smoothing function is valid. - [3.0 times mesh element size can be a good value] + [3.0 times mesh element size can be a good value] vert_axis: optional. index of vertical in position vector, must always be aligned with gravity, by default set to 1]. @@ -669,13 +671,16 @@ def setTwoPhaseVelocityInlet(self, U, waterLevel, smoothing, Uwind=None, Air K inflow value for turbulent model imposed at the boundary. dissipationInflowAir: float (optional). Air dissipation inflow value for turbulent model imposed at the boundary. + + Notes + ----- Below the seawater level, the condition returns the _dirichlet and p_advective condition according to the inflow velocity. Above the sea water level, the condition returns the gravity as zero, and sets _dirichlet condition to zero, only if there is a zero inflow velocity component. (!) This condition is best used for boundaries and gravity aligned with - one of the main axes. + one of the main axes. """ self.reset() self.BC_type = 'TwoPhaseVelocityInlet' @@ -777,7 +782,7 @@ def setHydrostaticPressureOutletWithDepth(self, seaLevel, rhoUp, rhoDown, g, If the boundary is aligned with one of the main axes, sets the tangential velocity components to zero as well. (!) This condition is best used for boundaries and gravity aligned with - one of the main axes. + one of the main axes. Parameters ---------- @@ -1225,7 +1230,7 @@ def __init__(self, turbModel, kWall, Y, Yplus, U0, nu=1.004e-6, Cmu=0.09, K=0.41 impose the wall shear stress. - k is assumed to be constant in the fully turbulent region close to the wall, - in this way kv = kp. + in this way kv = kp. - dissipation is calculated. Parameters diff --git a/proteus/mprans/MoveMeshMonitor.py b/proteus/mprans/MoveMeshMonitor.py index 3d731051f8..04729b7df2 100644 --- a/proteus/mprans/MoveMeshMonitor.py +++ b/proteus/mprans/MoveMeshMonitor.py @@ -50,7 +50,7 @@ class Coefficients(TransportCoefficients.PoissonEquationCoefficients): grading: double grading of mesh with distance function. 1.-> no grading, 1.1->10% increase (!) if using grading, func must be a distance function - otherwise func is an area function + otherwise func is an area function grading_type: int type of grading: 0 -> no grading diff --git a/proteus/mprans/Pres.py b/proteus/mprans/Pres.py index 7e4a9bb71b..79972f555d 100644 --- a/proteus/mprans/Pres.py +++ b/proteus/mprans/Pres.py @@ -30,14 +30,14 @@ def __init__(self, class Coefficients(TC_base): - """ + r""" The coefficients for pressure solution Update is given by .. math:: - p^{k+1} - p^{k} - phi^{k+1} + \nabla\\cdot(\\mu \\mathbf{u}^{k+1}) = 0 + p^{k+1} - p^{k} - phi^{k+1} + \nabla\cdot(\mu \mathbf{u}^{k+1}) = 0 """ def __init__(self, diff --git a/proteus/mprans/PresInc.py b/proteus/mprans/PresInc.py index d6699db4d8..0d28623332 100644 --- a/proteus/mprans/PresInc.py +++ b/proteus/mprans/PresInc.py @@ -31,14 +31,14 @@ def __init__(self, class Coefficients(TC_base): - """ + r""" The coefficients for pressure increment solution Update is given by .. math:: - \nabla\\cdot( -a \nabla \\phi^{k+1} - \\mathbf{q^t}^{k+1} ) = 0 + \nabla\cdot( -a \nabla \phi^{k+1} - \mathbf{q^t}^{k+1} ) = 0 a = \frac{\tau (1-\theta_s)}{\rho_f} + \frac{\tau \theta_s}{\rho_s} q^t = (1-\theta_s) v_f + \theta_s v_s """ diff --git a/setup.py b/setup.py index c1aa8f7c3f..96c7e22a2a 100644 --- a/setup.py +++ b/setup.py @@ -414,7 +414,17 @@ def _host_supports_avx(): language='c++', extra_compile_args=PROTEUS_OPT+['-std=c++20'], ), - + Extension( + # proteus/richards/ADR.py does `from .cADR import *`. This entry was + # dropped in merge 9e2aca55 (2026-08-12); see erdc/proteus#1285. + 'richards.cADR', + sources=['proteus/richards/cADR.cpp'], + depends=['proteus/richards/ADR.h', 'proteus/mprans/ArgumentsDict.h', 'proteus/ModelFactory.h', 'proteus/CompKernel.h'], + include_dirs=get_xtensor_include(), + language='c++', + extra_compile_args=PROTEUS_OPT+['-std=c++20'], + ), + Extension( 'm_comp_co2.cm_comp_co2', sources=['proteus/m_comp_co2/cm_comp_co2.cpp'],