From f33caca00dd1a14a2d89735d8fab9f95a475b86a Mon Sep 17 00:00:00 2001 From: Tobias Hienzsch Date: Sun, 1 Sep 2024 23:24:13 +0200 Subject: [PATCH] [python] Add pylint config --- .pylintrc | 641 ++++++++++++++++++ .../pffdtd/air_abs/get_air_absorption.py | 1 - src/python/pffdtd/air_abs/modal_filter.py | 58 +- src/python/pffdtd/air_abs/ola_filter.py | 42 +- .../pffdtd/air_abs/test_air_abs_filters.py | 44 +- src/python/pffdtd/air_abs/visco_filter.py | 40 +- src/python/pffdtd/analysis/response.py | 6 +- src/python/pffdtd/analysis/room_modes.py | 6 +- src/python/pffdtd/analysis/t60.py | 8 +- src/python/pffdtd/analysis/waterfall.py | 18 +- src/python/pffdtd/common/box.py | 19 +- src/python/pffdtd/common/check_version.py | 23 +- src/python/pffdtd/common/myfuncs.py | 222 +++--- src/python/pffdtd/diffusor/measurement.py | 2 +- src/python/pffdtd/diffusor/prd.py | 2 +- src/python/pffdtd/localization/tdoa.py | 4 +- src/python/pffdtd/materials/adm_funcs.py | 318 +++++---- src/python/pffdtd/sim2d/fdtd.py | 39 +- src/python/pffdtd/sim2d/report.py | 4 +- src/python/pffdtd/sim3d/engine.py | 109 +-- src/python/pffdtd/sim3d/sim_constants.py | 3 +- src/python/pffdtd/sim3d/sim_setup.py | 76 ++- src/python/pffdtd/sim3d/sim_signals.py | 27 +- src/python/pffdtd/voxelizer/cart_grid.py | 83 ++- src/python/pffdtd/voxelizer/vox_grid.py | 39 +- src/python/pffdtd/voxelizer/vox_grid_base.py | 18 +- src/python/pffdtd/voxelizer/vox_scene.py | 79 +-- 27 files changed, 1283 insertions(+), 648 deletions(-) create mode 100644 .pylintrc diff --git a/.pylintrc b/.pylintrc new file mode 100644 index 0000000..aa23335 --- /dev/null +++ b/.pylintrc @@ -0,0 +1,641 @@ +[MAIN] + +# Analyse import fallback blocks. This can be used to support both Python 2 and +# 3 compatible code, which means that the block might have code that exists +# only in one or another interpreter, leading to false positives when analysed. +analyse-fallback-blocks=no + +# Clear in-memory caches upon conclusion of linting. Useful if running pylint +# in a server-like mode. +clear-cache-post-run=no + +# Load and enable all available extensions. Use --list-extensions to see a list +# all available extensions. +#enable-all-extensions= + +# In error mode, messages with a category besides ERROR or FATAL are +# suppressed, and no reports are done by default. Error mode is compatible with +# disabling specific errors. +#errors-only= + +# Always return a 0 (non-error) status code, even if lint errors are found. +# This is primarily useful in continuous integration scripts. +#exit-zero= + +# A comma-separated list of package or module names from where C extensions may +# be loaded. Extensions are loading into the active Python interpreter and may +# run arbitrary code. +extension-pkg-allow-list= + +# A comma-separated list of package or module names from where C extensions may +# be loaded. Extensions are loading into the active Python interpreter and may +# run arbitrary code. (This is an alternative name to extension-pkg-allow-list +# for backward compatibility.) +extension-pkg-whitelist= + +# Return non-zero exit code if any of these messages/categories are detected, +# even if score is above --fail-under value. Syntax same as enable. Messages +# specified are enabled, while categories only check already-enabled messages. +fail-on= + +# Specify a score threshold under which the program will exit with error. +fail-under=10 + +# Interpret the stdin as a python script, whose filename needs to be passed as +# the module_or_package argument. +#from-stdin= + +# Files or directories to be skipped. They should be base names, not paths. +ignore=CVS + +# Add files or directories matching the regular expressions patterns to the +# ignore-list. The regex matches against paths and can be in Posix or Windows +# format. Because '\\' represents the directory delimiter on Windows systems, +# it can't be used as an escape character. +ignore-paths= + +# Files or directories matching the regular expression patterns are skipped. +# The regex matches against base names, not paths. The default value ignores +# Emacs file locks +ignore-patterns=^\.# + +# List of module names for which member attributes should not be checked +# (useful for modules/projects where namespaces are manipulated during runtime +# and thus existing member attributes cannot be deduced by static analysis). It +# supports qualified module names, as well as Unix pattern matching. +ignored-modules= + +# Python code to execute, usually for sys.path manipulation such as +# pygtk.require(). +#init-hook= +#init-hook='import sys; sys.path.append(".")' + +# Use multiple processes to speed up Pylint. Specifying 0 will auto-detect the +# number of processors available to use, and will cap the count on Windows to +# avoid hangs. +jobs=1 + +# Control the amount of potential inferred values when inferring a single +# object. This can help the performance when dealing with large functions or +# complex, nested conditions. +limit-inference-results=100 + +# List of plugins (as comma separated values of python module names) to load, +# usually to register additional checkers. +load-plugins= + +# Pickle collected data for later comparisons. +persistent=yes + +# Minimum Python version to use for version dependent checks. Will default to +# the version used to run pylint. +py-version=3.10 + +# Discover python modules and packages in the file system subtree. +recursive=no + +# When enabled, pylint would attempt to guess common misconfiguration and emit +# user-friendly hints instead of false-positive error messages. +suggestion-mode=yes + +# Allow loading of arbitrary C extensions. Extensions are imported into the +# active Python interpreter and may run arbitrary code. +unsafe-load-any-extension=no + +# In verbose mode, extra non-checker-related info will be displayed. +#verbose= + + +[BASIC] + +# Naming style matching correct argument names. +argument-naming-style=snake_case + +# Regular expression matching correct argument names. Overrides argument- +# naming-style. If left empty, argument names will be checked with the set +# naming style. +#argument-rgx= + +# Naming style matching correct attribute names. +attr-naming-style=snake_case + +# Regular expression matching correct attribute names. Overrides attr-naming- +# style. If left empty, attribute names will be checked with the set naming +# style. +#attr-rgx= + +# Bad variable names which should always be refused, separated by a comma. +bad-names=foo, + bar, + baz, + toto, + tutu, + tata + +# Bad variable names regexes, separated by a comma. If names match any regex, +# they will always be refused +bad-names-rgxs= + +# Naming style matching correct class attribute names. +class-attribute-naming-style=any + +# Regular expression matching correct class attribute names. Overrides class- +# attribute-naming-style. If left empty, class attribute names will be checked +# with the set naming style. +#class-attribute-rgx= + +# Naming style matching correct class constant names. +class-const-naming-style=UPPER_CASE + +# Regular expression matching correct class constant names. Overrides class- +# const-naming-style. If left empty, class constant names will be checked with +# the set naming style. +#class-const-rgx= + +# Naming style matching correct class names. +class-naming-style=PascalCase + +# Regular expression matching correct class names. Overrides class-naming- +# style. If left empty, class names will be checked with the set naming style. +#class-rgx= + +# Naming style matching correct constant names. +const-naming-style=UPPER_CASE + +# Regular expression matching correct constant names. Overrides const-naming- +# style. If left empty, constant names will be checked with the set naming +# style. +#const-rgx= + +# Minimum line length for functions/classes that require docstrings, shorter +# ones are exempt. +docstring-min-length=-1 + +# Naming style matching correct function names. +function-naming-style=snake_case + +# Regular expression matching correct function names. Overrides function- +# naming-style. If left empty, function names will be checked with the set +# naming style. +#function-rgx= + +# Good variable names which should always be accepted, separated by a comma. +good-names=i, + j, + k, + ex, + Run, + _ + +# Good variable names regexes, separated by a comma. If names match any regex, +# they will always be accepted +good-names-rgxs= + +# Include a hint for the correct naming format with invalid-name. +include-naming-hint=no + +# Naming style matching correct inline iteration names. +inlinevar-naming-style=any + +# Regular expression matching correct inline iteration names. Overrides +# inlinevar-naming-style. If left empty, inline iteration names will be checked +# with the set naming style. +#inlinevar-rgx= + +# Naming style matching correct method names. +method-naming-style=snake_case + +# Regular expression matching correct method names. Overrides method-naming- +# style. If left empty, method names will be checked with the set naming style. +#method-rgx= + +# Naming style matching correct module names. +module-naming-style=snake_case + +# Regular expression matching correct module names. Overrides module-naming- +# style. If left empty, module names will be checked with the set naming style. +#module-rgx= + +# Colon-delimited sets of names that determine each other's naming style when +# the name regexes allow several styles. +name-group= + +# Regular expression which should only match function or class names that do +# not require a docstring. +no-docstring-rgx=^_ + +# List of decorators that produce properties, such as abc.abstractproperty. Add +# to this list to register other decorators that produce valid properties. +# These decorators are taken in consideration only for invalid-name. +property-classes=abc.abstractproperty + +# Regular expression matching correct type variable names. If left empty, type +# variable names will be checked with the set naming style. +#typevar-rgx= + +# Naming style matching correct variable names. +variable-naming-style=snake_case + +# Regular expression matching correct variable names. Overrides variable- +# naming-style. If left empty, variable names will be checked with the set +# naming style. +#variable-rgx= + + +[CLASSES] + +# Warn about protected attribute access inside special methods +check-protected-access-in-special-methods=no + +# List of method names used to declare (i.e. assign) instance attributes. +defining-attr-methods=__init__, + __new__, + setUp, + __post_init__ + +# List of member names, which should be excluded from the protected access +# warning. +exclude-protected=_asdict, + _fields, + _replace, + _source, + _make + +# List of valid names for the first argument in a class method. +valid-classmethod-first-arg=cls + +# List of valid names for the first argument in a metaclass class method. +valid-metaclass-classmethod-first-arg=mcs + + +[DESIGN] + +# List of regular expressions of class ancestor names to ignore when counting +# public methods (see R0903) +exclude-too-few-public-methods= + +# List of qualified class names to ignore when counting class parents (see +# R0901) +ignored-parents= + +# Maximum number of arguments for function / method. +max-args=5 + +# Maximum number of attributes for a class (see R0902). +max-attributes=7 + +# Maximum number of boolean expressions in an if statement (see R0916). +max-bool-expr=5 + +# Maximum number of branch for function / method body. +max-branches=12 + +# Maximum number of locals for function / method body. +max-locals=15 + +# Maximum number of parents for a class (see R0901). +max-parents=7 + +# Maximum number of public methods for a class (see R0904). +max-public-methods=20 + +# Maximum number of return / yield for function / method body. +max-returns=6 + +# Maximum number of statements in function / method body. +max-statements=50 + +# Minimum number of public methods for a class (see R0903). +min-public-methods=2 + + +[EXCEPTIONS] + +# Exceptions that will emit a warning when caught. +overgeneral-exceptions=builtins.BaseException,builtins.Exception + + +[FORMAT] + +# Expected format of line ending, e.g. empty (any line ending), LF or CRLF. +expected-line-ending-format= + +# Regexp for a line that is allowed to be longer than the limit. +ignore-long-lines=^\s*(# )??$ + +# Number of spaces of indent required inside a hanging or continued line. +indent-after-paren=4 + +# String used as indentation unit. This is usually " " (4 spaces) or "\t" (1 +# tab). +indent-string=' ' + +# Maximum number of characters on a single line. +max-line-length=100 + +# Maximum number of lines in a module. +max-module-lines=1000 + +# Allow the body of a class to be on the same line as the declaration if body +# contains single statement. +single-line-class-stmt=no + +# Allow the body of an if to be on the same line as the test if there is no +# else. +single-line-if-stmt=no + + +[IMPORTS] + +# List of modules that can be imported at any level, not just the top level +# one. +allow-any-import-level= + +# Allow explicit reexports by alias from a package __init__. +allow-reexport-from-package=no + +# Allow wildcard imports from modules that define __all__. +allow-wildcard-with-all=no + +# Deprecated modules which should not be used, separated by a comma. +deprecated-modules= + +# Output a graph (.gv or any supported image format) of external dependencies +# to the given file (report RP0402 must not be disabled). +ext-import-graph= + +# Output a graph (.gv or any supported image format) of all (i.e. internal and +# external) dependencies to the given file (report RP0402 must not be +# disabled). +import-graph= + +# Output a graph (.gv or any supported image format) of internal dependencies +# to the given file (report RP0402 must not be disabled). +int-import-graph= + +# Force import order to recognize a module as part of the standard +# compatibility libraries. +known-standard-library= + +# Force import order to recognize a module as part of a third party library. +known-third-party=enchant + +# Couples of modules and preferred modules, separated by a comma. +preferred-modules= + + +[LOGGING] + +# The type of string formatting that logging methods do. `old` means using % +# formatting, `new` is for `{}` formatting. +logging-format-style=old + +# Logging modules to check that the string format arguments are in logging +# function parameter format. +logging-modules=logging + + +[MESSAGES CONTROL] + +# Only show warnings with the listed confidence levels. Leave empty to show +# all. Valid levels: HIGH, CONTROL_FLOW, INFERENCE, INFERENCE_FAILURE, +# UNDEFINED. +confidence=HIGH, + CONTROL_FLOW, + INFERENCE, + INFERENCE_FAILURE, + UNDEFINED + +# Disable the message, report, category or checker with the given id(s). You +# can either give multiple identifiers separated by comma (,) or put this +# option multiple times (only on the command line, not in the configuration +# file where it should appear only once). You can also use "--disable=all" to +# disable everything first and then re-enable specific checks. For example, if +# you want to run only the similarities checker, you can use "--disable=all +# --enable=similarities". If you want to run only the classes checker, but have +# no Warning level messages displayed, use "--disable=all --enable=classes +# --disable=W". +disable=raw-checker-failed, + bad-inline-option, + locally-disabled, + file-ignored, + suppressed-message, + useless-suppression, + deprecated-pragma, + use-symbolic-message-instead, + invalid-name, + missing-class-docstring, + missing-function-docstring, + missing-module-docstring, + import-outside-toplevel, + import-error, + too-many-arguments, + too-many-locals, + too-many-statements, + wrong-import-order, + line-too-long, + no-member, + attribute-defined-outside-init, + not-an-iterable, + unspecified-encoding, + too-many-instance-attributes, + +# Enable the message, report, category or checker with the given id(s). You can +# either give multiple identifier separated by comma (,) or put this option +# multiple time (only on the command line, not in the configuration file where +# it should appear only once). See also the "--disable" option for examples. +enable=c-extension-no-member + + +[METHOD_ARGS] + +# List of qualified names (i.e., library.method) which require a timeout +# parameter e.g. 'requests.api.get,requests.api.post' +timeout-methods=requests.api.delete,requests.api.get,requests.api.head,requests.api.options,requests.api.patch,requests.api.post,requests.api.put,requests.api.request + + +[MISCELLANEOUS] + +# List of note tags to take in consideration, separated by a comma. +notes=FIXME, + XXX, + TODO + +# Regular expression of note tags to take in consideration. +notes-rgx= + + +[REFACTORING] + +# Maximum number of nested blocks for function / method body +max-nested-blocks=5 + +# Complete name of functions that never returns. When checking for +# inconsistent-return-statements if a never returning function is called then +# it will be considered as an explicit return statement and no message will be +# printed. +never-returning-functions=sys.exit,argparse.parse_error + + +[REPORTS] + +# Python expression which should return a score less than or equal to 10. You +# have access to the variables 'fatal', 'error', 'warning', 'refactor', +# 'convention', and 'info' which contain the number of messages in each +# category, as well as 'statement' which is the total number of statements +# analyzed. This score is used by the global evaluation report (RP0004). +evaluation=max(0, 0 if fatal else 10.0 - ((float(5 * error + warning + refactor + convention) / statement) * 10)) + +# Template used to display messages. This is a python new-style format string +# used to format the message information. See doc for all details. +msg-template= + +# Set the output format. Available formats are text, parseable, colorized, json +# and msvs (visual studio). You can also give a reporter class, e.g. +# mypackage.mymodule.MyReporterClass. +#output-format= + +# Tells whether to display a full report or only the messages. +reports=no + +# Activate the evaluation score. +score=yes + + +[SIMILARITIES] + +# Comments are removed from the similarity computation +ignore-comments=yes + +# Docstrings are removed from the similarity computation +ignore-docstrings=yes + +# Imports are removed from the similarity computation +ignore-imports=yes + +# Signatures are removed from the similarity computation +ignore-signatures=yes + +# Minimum lines number of a similarity. +min-similarity-lines=4 + + +[SPELLING] + +# Limits count of emitted suggestions for spelling mistakes. +max-spelling-suggestions=4 + +# Spelling dictionary name. Available dictionaries: none. To make it work, +# install the 'python-enchant' package. +spelling-dict= + +# List of comma separated words that should be considered directives if they +# appear at the beginning of a comment and should not be checked. +spelling-ignore-comment-directives=fmt: on,fmt: off,noqa:,noqa,nosec,isort:skip,mypy: + +# List of comma separated words that should not be checked. +spelling-ignore-words= + +# A path to a file that contains the private dictionary; one word per line. +spelling-private-dict-file= + +# Tells whether to store unknown words to the private dictionary (see the +# --spelling-private-dict-file option) instead of raising a message. +spelling-store-unknown-words=no + + +[STRING] + +# This flag controls whether inconsistent-quotes generates a warning when the +# character used as a quote delimiter is used inconsistently within a module. +check-quote-consistency=no + +# This flag controls whether the implicit-str-concat should generate a warning +# on implicit string concatenation in sequences defined over several lines. +check-str-concat-over-line-jumps=no + + +[TYPECHECK] + +# List of decorators that produce context managers, such as +# contextlib.contextmanager. Add to this list to register other decorators that +# produce valid context managers. +contextmanager-decorators=contextlib.contextmanager + +# List of members which are set dynamically and missed by pylint inference +# system, and so shouldn't trigger E1101 when accessed. Python regular +# expressions are accepted. +generated-members= + +# Tells whether to warn about missing members when the owner of the attribute +# is inferred to be None. +ignore-none=yes + +# This flag controls whether pylint should warn about no-member and similar +# checks whenever an opaque object is returned when inferring. The inference +# can return multiple potential results while evaluating a Python object, but +# some branches might not be evaluated, which results in partial inference. In +# that case, it might be useful to still emit no-member and other checks for +# the rest of the inferred objects. +ignore-on-opaque-inference=yes + +# List of symbolic message names to ignore for Mixin members. +ignored-checks-for-mixins=no-member, + not-async-context-manager, + not-context-manager, + attribute-defined-outside-init + +# List of class names for which member attributes should not be checked (useful +# for classes with dynamically set attributes). This supports the use of +# qualified names. +ignored-classes=optparse.Values,thread._local,_thread._local,argparse.Namespace + +# Show a hint with possible names when a member name was not found. The aspect +# of finding the hint is based on edit distance. +missing-member-hint=yes + +# The minimum edit distance a name should have in order to be considered a +# similar match for a missing member name. +missing-member-hint-distance=1 + +# The total number of similar names that should be taken in consideration when +# showing a hint for a missing member. +missing-member-max-choices=1 + +# Regex pattern to define which classes are considered mixins. +mixin-class-rgx=.*[Mm]ixin + +# List of decorators that change the signature of a decorated function. +signature-mutators= + + +[VARIABLES] + +# List of additional names supposed to be defined in builtins. Remember that +# you should avoid defining new builtins when possible. +additional-builtins= + +# Tells whether unused global variables should be treated as a violation. +allow-global-unused-variables=yes + +# List of names allowed to shadow builtins +allowed-redefined-builtins= + +# List of strings which can identify a callback function by name. A callback +# name must start or end with one of those strings. +callbacks=cb_, + _cb + +# A regular expression matching the name of dummy variables (i.e. expected to +# not be used). +dummy-variables-rgx=_+$|(_[a-zA-Z0-9_]*[a-zA-Z0-9]+?$)|dummy|^ignored_|^unused_ + +# Argument names that match this expression will be ignored. +ignored-argument-names=_.*|^ignored_|^unused_ + +# Tells whether we should check for unused import in __init__ files. +init-import=no + +# List of qualified module names which can have objects that can redefine +# builtins. +redefining-builtins-modules=six.moves,past.builtins,future.builtins,builtins,io diff --git a/src/python/pffdtd/air_abs/get_air_absorption.py b/src/python/pffdtd/air_abs/get_air_absorption.py index b9243b2..7dabde8 100644 --- a/src/python/pffdtd/air_abs/get_air_absorption.py +++ b/src/python/pffdtd/air_abs/get_air_absorption.py @@ -1,7 +1,6 @@ # SPDX-License-Identifier: MIT import numpy as np -from numpy import array as npa from numpy import log10, exp, sqrt, log, pi def get_air_absorption(freq_vec,temperature_celsius,rel_humidity_pnct,pressure_atmospheric_kPa=101.325): diff --git a/src/python/pffdtd/air_abs/modal_filter.py b/src/python/pffdtd/air_abs/modal_filter.py index 263430e..37eead1 100644 --- a/src/python/pffdtd/air_abs/modal_filter.py +++ b/src/python/pffdtd/air_abs/modal_filter.py @@ -1,19 +1,16 @@ # SPDX-License-Identifier: MIT -import numpy as np import numba as nb -from numpy import array as npa -from numpy import exp, sqrt, log, pi, cos -from scipy.fft import dct,idct #default type2 -from pffdtd.air_abs.get_air_absorption import get_air_absorption -from pffdtd.common.myfuncs import iround, iceil +import numpy as np +from numpy import exp, pi, cos +from scipy.fft import dct, idct # default type2 from tqdm import tqdm -#apply filter, x is (Nchannel,Nsamples) array -#see I3DA paper for details -#Tc is temperature in deg Celsius -#rh is relative humidity -def apply_modal_filter(x,Fs,Tc,rh,pad_t=0.0): +from pffdtd.air_abs.get_air_absorption import get_air_absorption +from pffdtd.common.myfuncs import iceil + + +def apply_modal_filter(x, Fs, Tc, rh, pad_t=0.0): """ This is an implementation of an air absorption filter based on a modal approach. It solves a system of 1-d dissipative wave equations tune to air attenuation @@ -23,13 +20,19 @@ def apply_modal_filter(x,Fs,Tc,rh,pad_t=0.0): Hamilton, B. "Adding air attenuation to simulated room impulse responses: A modal approach", to be presented at I3DA 2021 conference in Bologna Italy. """ + + # apply filter, x is (Nchannel,Nsamples) array + # see I3DA paper for details + # Tc is temperature in deg Celsius + # rh is relative humidity + Ts = 1/Fs x = np.atleast_2d(x) Nt0 = x.shape[-1] Nt = iceil(pad_t/Ts)+Nt0 - xp = np.zeros((x.shape[0],Nt)) - xp[:,:Nt0] = x + xp = np.zeros((x.shape[0], Nt)) + xp[:, :Nt0] = x del x y = np.zeros(xp.shape) @@ -38,7 +41,7 @@ def apply_modal_filter(x,Fs,Tc,rh,pad_t=0.0): wqTs = pi*(np.arange(Nx)/Nx) wq = wqTs/Ts - rd = get_air_absorption(wq/2/pi,Tc,rh) + rd = get_air_absorption(wq/2/pi, Tc, rh) alphaq = rd['absfull_Np'] c = rd['c'] @@ -46,8 +49,8 @@ def apply_modal_filter(x,Fs,Tc,rh,pad_t=0.0): P1 = np.zeros(xp.shape) fx = np.zeros(xp.shape) - fx[:,0] = 1 - Fm = dct(fx,type=2,norm='ortho',axis=-1) + fx[:, 0] = 1 + Fm = dct(fx, type=2, norm='ortho', axis=-1) sigqTs = c*alphaq*Ts a1 = 2*exp(-sigqTs)*cos(wqTs) @@ -55,23 +58,24 @@ def apply_modal_filter(x,Fs,Tc,rh,pad_t=0.0): Fmsig1 = Fm*(1+sigqTs/2)/(1+sigqTs) Fmsig2 = Fm*(1-sigqTs/2)/(1+sigqTs) - u = np.zeros((xp.shape[0],Nt+1)) - u[:,1:] = xp[:,::-1] #flip + u = np.zeros((xp.shape[0], Nt+1)) + u[:, 1:] = xp[:, ::-1] # flip - pbar = tqdm(total=Nt,desc=f'modal filter',ascii=True) + pbar = tqdm(total=Nt, desc='modal filter', ascii=True) - @nb.jit(nopython=True,parallel=True) - def _run_step(P0,P1,a1,a2,Fmsig1,Fmsig2,un1,un0): + @nb.jit(nopython=True, parallel=True) + def _run_step(P0, P1, a1, a2, Fmsig1, Fmsig2, un1, un0): P0[:] = a1*P1 + a2*P0 + Fmsig1*un1 - Fmsig2*un0 for n in range(Nt): - #P0 = a1*P1 + a2*P0 + Fmsig1*u[:,n+1] - Fmsig2*u[:,n] + # P0 = a1*P1 + a2*P0 + Fmsig1*u[:,n+1] - Fmsig2*u[:,n] for i in range(P0.shape[0]): - _run_step(P0[i],P1[i],a1,a2,Fmsig1[i],Fmsig2[i],u[i,n+1],u[i,n]) - if n0) + assert n_start > 0 - y[:,:n_start] = x[:,:n_start] + y[:, :n_start] = x[:, :n_start] Tsg2 = 2*Ts*g Tsg2pi = 2*Ts*g*pi gTs = g*Ts dt_fac = 0.1*log(10)*NdB*gTs - pbar = tqdm(total=Nt,desc=f'visco filter',ascii=True) - for n in range(n_start,Nt0): + pbar = tqdm(total=Nt, desc='visco filter', ascii=True) + for n in range(n_start, Nt0): dt = sqrt(dt_fac*n)/Ts dt_int = iceil(dt) - nv = np.arange(n-dt_int,n+dt_int+1) - assert n>=dt_int - y[:,nv] += (Ts/sqrt(n*Tsg2pi))*x[:,n][:,None] * exp(-((n-nv)*Ts)**2/(n*Tsg2))[None,:] + nv = np.arange(n-dt_int, n+dt_int+1) + assert n >= dt_int + y[:, nv] += (Ts/sqrt(n*Tsg2pi))*x[:, n][:, None] * \ + exp(-((n-nv)*Ts)**2/(n*Tsg2))[None, :] pbar.update(1) pbar.close() - return np.squeeze(y) #squeeze to 1d in case + return np.squeeze(y) # squeeze to 1d in case diff --git a/src/python/pffdtd/analysis/response.py b/src/python/pffdtd/analysis/response.py index 8353919..d3ba72e 100644 --- a/src/python/pffdtd/analysis/response.py +++ b/src/python/pffdtd/analysis/response.py @@ -1,5 +1,5 @@ import click -import scipy.io.wavfile as wavfile +from scipy.io import wavfile import matplotlib.pyplot as plt import numpy as np from matplotlib.ticker import ScalarFormatter @@ -48,9 +48,6 @@ def main(filename, fmin, fmax, label_a, label_b, smoothing, target): file_a = filename[0] file_b = filename[1] - label_a = label_a - label_b = label_b - fs_a, buf_a = wavfile.read(file_a) fs_b, buf_b = wavfile.read(file_b) @@ -75,7 +72,6 @@ def main(filename, fmin, fmax, label_a, label_b, smoothing, target): dB_b += 75.0 if smoothing > 0.0: - smoothing = smoothing dB_a = fractional_octave_smoothing(dB_a, fs_a, nfft, smoothing) dB_b = fractional_octave_smoothing(dB_b, fs_b, nfft, smoothing) diff --git a/src/python/pffdtd/analysis/room_modes.py b/src/python/pffdtd/analysis/room_modes.py index d19d548..d86a441 100644 --- a/src/python/pffdtd/analysis/room_modes.py +++ b/src/python/pffdtd/analysis/room_modes.py @@ -6,7 +6,7 @@ import matplotlib.pyplot as plt import numpy as np from scipy.signal import find_peaks -import scipy.io.wavfile as wavfile +from scipy.io import wavfile from pffdtd.common.myfuncs import iceil from pffdtd.common.plot import plot_styles @@ -18,8 +18,8 @@ def find_nearest(array, value): return array[idx] -def collect_wav_paths(dir, pattern="*.wav"): - return list(sorted(glob.glob(os.path.join(dir, pattern)))) +def collect_wav_paths(folder, pattern="*.wav"): + return list(sorted(glob.glob(os.path.join(folder, pattern)))) def hz_to_note(frequency): diff --git a/src/python/pffdtd/analysis/t60.py b/src/python/pffdtd/analysis/t60.py index b3e6c3f..1526d7f 100644 --- a/src/python/pffdtd/analysis/t60.py +++ b/src/python/pffdtd/analysis/t60.py @@ -4,11 +4,11 @@ import click import numpy as np -import scipy.io.wavfile as wavfile -import scipy.signal as signal import matplotlib.pyplot as plt from matplotlib.axes import Axes from matplotlib.ticker import ScalarFormatter +from scipy import signal +from scipy.io import wavfile from pffdtd.common.plot import plot_styles @@ -73,7 +73,7 @@ def run(files, fmin, fmax, show_all=False, show_tolerance=True, target=None): plt.rcParams.update(plot_styles) - fig, axs = plt.subplots(2, 1) + _, axs = plt.subplots(2, 1) formatter = ScalarFormatter() formatter.set_scientific(False) @@ -140,7 +140,7 @@ def run(files, fmin, fmax, show_all=False, show_tolerance=True, target=None): label="EBU Tech 3000" ) - ax.set_title(f"Tolerance") + ax.set_title("Tolerance") ax.set_ylabel("Difference [s]") ax.set_xlabel("Frequency [Hz]") ax.xaxis.set_major_formatter(formatter) diff --git a/src/python/pffdtd/analysis/waterfall.py b/src/python/pffdtd/analysis/waterfall.py index 17c97ad..9cdb57d 100644 --- a/src/python/pffdtd/analysis/waterfall.py +++ b/src/python/pffdtd/analysis/waterfall.py @@ -15,7 +15,7 @@ def main(filename): frequencies, times, Zxx = stft(rir, fs=fs, nperseg=nfft) decay_time = np.zeros(frequencies.shape) - for i, freq in enumerate(frequencies): + for i, _ in enumerate(frequencies): magnitude = np.abs(Zxx[i, :]) magnitude_db = 20 * np.log10(magnitude / np.max(magnitude)) try: @@ -39,15 +39,15 @@ def main(filename): )]) fig.update_layout( title='Decay Times', - scene=dict( - xaxis_title='X: Time [s]', - yaxis_title='Y: Frequency [Hz]', - zaxis_title='Z: Amplitude [dB]', + scene={ + "xaxis_title": 'X: Time [s]', + "yaxis_title": 'Y: Frequency [Hz]', + "zaxis_title": 'Z: Amplitude [dB]', - xaxis=dict(range=[0, 1],), - yaxis=dict(type='log',), - # zaxis=dict(range=[-60, 0],), - ), + "xaxis": {"range": [0, 1]}, + "yaxis": {"type": "log"}, + # "zaxis": {"range": [-60, 0]}, + } ) Z = 20 * np.log10(np.abs(Zxx/nfft)) diff --git a/src/python/pffdtd/common/box.py b/src/python/pffdtd/common/box.py index 6b38321..7785aa0 100644 --- a/src/python/pffdtd/common/box.py +++ b/src/python/pffdtd/common/box.py @@ -1,24 +1,13 @@ -############################################################################## -# This file is a part of PFFDTD. -# -# PFFTD is released under the MIT License. -# For details see the LICENSE file. -# -# Copyright 2021 Brian Hamilton. -# -# File name: box.py -# -# Description: This is a class for a box. Used in a few places (vox_grid, tri_box). -# -############################################################################## +# SPDX-License-Identifier: MIT import numpy as np from numpy import array as npa import numpy.random as npr -import numpy.linalg as npl from pffdtd.common.myfuncs import mydefault,rotmatrix_ax_ang class Box: + """This is a class for a box. Used in a few places (vox_grid, tri_box). + """ def __init__(self,Lx=None,Ly=None,Lz=None,Rax=None,Rang=None,shift=None,centered=True): #defaults Lx = mydefault(Lx,1.0) @@ -133,7 +122,7 @@ def main(): Lx = npr.random() Ly = npr.random() Lz = npr.random() - print('Lx = %.2f, Ly = %.2f, Lz = %.2f' % (Lx,Ly,Lz)) + print(f'Lx = {Lx:.2f}, Ly = {Ly:.2f}, Lz = {Lz:.2f}') box = Box(Lx,Ly,Lz) print('box 0 ... bmin = %.2f,%.2f,%.2f' % (box.bmin[0],box.bmin[1],box.bmin[2])) diff --git a/src/python/pffdtd/common/check_version.py b/src/python/pffdtd/common/check_version.py index b6c7756..597e1f3 100644 --- a/src/python/pffdtd/common/check_version.py +++ b/src/python/pffdtd/common/check_version.py @@ -1,20 +1,9 @@ -############################################################################## -# This file is a part of PFFDTD. -# -# PFFTD is released under the MIT License. -# For details see the LICENSE file. -# -# Copyright 2021 Brian Hamilton. -# -# File name: check_version.py -# -# Description: checks Python version -# -############################################################################## +# SPDX-License-Identifier: MIT import sys -ATLEASTVERSION39 = (sys.version_info >= (3,9)) -ATLEASTVERSION38 = (sys.version_info >= (3,8)) -ATLEASTVERSION37 = (sys.version_info >= (3,7)) -ATLEASTVERSION36 = (sys.version_info >= (3,6)) + +ATLEASTVERSION39 = sys.version_info >= (3, 9) +ATLEASTVERSION38 = sys.version_info >= (3, 8) +ATLEASTVERSION37 = sys.version_info >= (3, 7) +ATLEASTVERSION36 = sys.version_info >= (3, 6) assert ATLEASTVERSION39 diff --git a/src/python/pffdtd/common/myfuncs.py b/src/python/pffdtd/common/myfuncs.py index e2f18c2..ab6a9c9 100644 --- a/src/python/pffdtd/common/myfuncs.py +++ b/src/python/pffdtd/common/myfuncs.py @@ -1,33 +1,25 @@ -############################################################################## -# This file is a part of PFFDTD. -# -# PFFTD is released under the MIT License. -# For details see the LICENSE file. -# -# Copyright 2021 Brian Hamilton. -# -# File name: myfuncs.py -# -# Description: Miscellaneous python/numpy functions. Not all used or useful. -# -############################################################################## - -import numpy as np -from numpy import array as npa -import numpy.linalg as npl -from numpy import cos,sin,pi -from typing import Any,Optional -from pffdtd.common.myasserts import assert_np_array_float, assert_np_array_complex -from pathlib import Path -import struct -import os +# SPDX-License-Identifier: MIT +"""Miscellaneous python/numpy functions. Not all used or useful. +""" import argparse import hashlib import multiprocessing as mp +import os +from pathlib import Path +from typing import Any +import struct + + +import numpy as np +import numpy.linalg as npl +from numpy import cos, sin, pi import scipy.io.wavfile +from pffdtd.common.myasserts import assert_np_array_complex + EPS = np.finfo(np.float64).eps + def rotmatrix_ax_ang(Rax:Any, Rang:float): assert isinstance(Rax,np.ndarray) assert Rax.shape == (3,) @@ -43,53 +35,56 @@ def rotmatrix_ax_ang(Rax:Any, Rang:float): assert npl.norm(npl.inv(R)-R.T) < 1e-8 return R -def rotate_xyz_deg(thx_d,thy_d,thz_d): - #R applies Rz then Ry then Rx (opposite to wikipedia) - #rotations about x,y,z axes, right hand rule + +def rotate_xyz_deg(thx_d, thy_d, thz_d): + # R applies Rz then Ry then Rx (opposite to wikipedia) + # rotations about x,y,z axes, right hand rule thx = np.deg2rad(thx_d) thy = np.deg2rad(thy_d) thz = np.deg2rad(thz_d) - Rx = npa([[1,0,0],\ - [0,np.cos(thx),-np.sin(thx)],\ - [0,np.sin(thx),np.cos(thx)]]) + Rx = np.array([[1, 0, 0], + [0, np.cos(thx), -np.sin(thx)], + [0, np.sin(thx), np.cos(thx)]]) - Ry = npa([[np.cos(thy),0,np.sin(thy)],\ - [0,1,0],\ - [-np.sin(thy),0,np.cos(thy)]]) + Ry = np.array([[np.cos(thy), 0, np.sin(thy)], + [0, 1, 0], + [-np.sin(thy), 0, np.cos(thy)]]) - Rz = npa([[np.cos(thz),-np.sin(thz),0],\ - [np.sin(thz),np.cos(thz),0],\ - [0,0,1]]) + Rz = np.array([[np.cos(thz), -np.sin(thz), 0], + [np.sin(thz), np.cos(thz), 0], + [0, 0, 1]]) R = Rx @ Ry @ Rz - return R,Rx,Ry,Rz + return R, Rx, Ry, Rz -def rotate_az_el_deg(az_d,el_d): - #R applies Rel then Raz - #Rel is rotation about negative y-axis, right hand rule - #Raz is rotation z-axis, right hand rule - #this uses matlab convention - #note, this isn't the most general rotation - _,_,Ry,Rz = rotate_xyz_deg(0,-el_d,az_d) +def rotate_az_el_deg(az_d, el_d): + # R applies Rel then Raz + # Rel is rotation about negative y-axis, right hand rule + # Raz is rotation z-axis, right hand rule + # this uses matlab convention + # note, this isn't the most general rotation + _, _, Ry, Rz = rotate_xyz_deg(0, -el_d, az_d) Rel = Ry Raz = Rz R = Raz @ Rel - return R,Raz,Rel + return R, Raz, Rel + -#only used in box -def mydefault(a:Any, default:Any): +def mydefault(a: Any, default: Any): + # only used in box if a is None: a = default return a -#seconds to day/hour/min/s -def s2dhms(time:int): - #give time in seconds + +def s2dhms(time: int): + # seconds to day/hour/min/s + # give time in seconds day = time // (24 * 3600) time = time % (24 * 3600) hour = time // 3600 @@ -99,76 +94,94 @@ def s2dhms(time:int): seconds = time return (day, hour, minutes, seconds) -#misc numpy functions -def clamp(x:Any, xmin:Any, xmax:Any): - assert isinstance(x,(np.ndarray,int,np.integer,float)) - assert isinstance(xmin,(np.ndarray,int,np.integer,float)) - assert isinstance(xmax,(np.ndarray,int,np.integer,float)) - return np.where(xxmax,xmax,x)) + +def clamp(x: Any, xmin: Any, xmax: Any): + assert isinstance(x, (np.ndarray, int, np.integer, float)) + assert isinstance(xmin, (np.ndarray, int, np.integer, float)) + assert isinstance(xmax, (np.ndarray, int, np.integer, float)) + return np.where(x < xmin, xmin, np.where(x > xmax, xmax, x)) + def vclamp(x, xmin, xmax): - return np.where(xxmax,xmax,x)) #only this works for vectorized xmin,xmax + # only this works for vectorized xmin,xmax + return np.where(x < xmin, xmin, np.where(x > xmax, xmax, x)) + def dot2(v): - return dotv(v,v) + return dotv(v, v) + + +def dotv(v1, v2): + return np.sum(v1*v2, axis=-1) -def dotv(v1,v2): - return np.sum(v1*v2,axis=-1) def vecnorm(v1): return np.sqrt(dot2(v1)) + def vecnorm2(v1): - return np.sum(v1**2,axis=-1) + return np.sum(v1**2, axis=-1) + -def normalise(v1,eps=EPS): +def normalise(v1, eps=EPS): return (v1.T/(vecnorm(v1)+eps)).T -def ceilint(x:float): - assert isinstance(x,float) + +def ceilint(x: float): + assert isinstance(x, float) return np.int_(np.ceil(x)) -def roundint(x:float): - assert isinstance(x,float) + +def roundint(x: float): + assert isinstance(x, float) return np.int_(np.round(x)) -def floorint(x:float): - assert isinstance(x,float) + +def floorint(x: float): + assert isinstance(x, float) return np.int_(np.floor(x)) -def maxabs(x:Any): + +def maxabs(x: Any): assert_np_array_complex(x) - if x.ndim==2: - return np.amax(np.abs(x),axis=-1)[:,None] + if x.ndim == 2: + return np.amax(np.abs(x), axis=-1)[:, None] else: return np.amax(np.abs(x)) -def to_col_2d(x:Any): - assert x.ndim<3 - x = x.reshape(x.shape[-1],-1) - assert x.shape[0]>x.shape[1] + +def to_col_2d(x: Any): + assert x.ndim < 3 + x = x.reshape(x.shape[-1], -1) + assert x.shape[0] > x.shape[1] return x -def to_row_2d(x:Any): - assert x.ndim<3 - x = x.reshape(-1,x.shape[-1]) - assert x.shape[1]>x.shape[0] + +def to_row_2d(x: Any): + assert x.ndim < 3 + x = x.reshape(-1, x.shape[-1]) + assert x.shape[1] > x.shape[0] return x -def ind2sub3d(ii,Nx,Ny,Nz): + +def ind2sub3d(ii, Nx, Ny, Nz): iz = ii % Nz iy = (ii - iz)//Nz % Ny ix = ((ii - iz)//Nz-iy)//Ny - return ix,iy,iz + return ix, iy, iz -def rel_diff(x0,x1): #relative difference at machine epsilon level + +def rel_diff(x0, x1): + # relative difference at machine epsilon level return (x0-x1)/(2.0**np.floor(np.log2(x0))) + def get_default_nprocs(): - return max(1,int(0.8*mp.cpu_count())) + return max(1, int(0.8*mp.cpu_count())) + -#taken from SO, for argparse bool problem -def str2bool(v:str): +def str2bool(v: str): + # taken from SO, for argparse bool problem assert isinstance(v, str) if v.lower() in ('true'): return True @@ -177,8 +190,9 @@ def str2bool(v:str): else: raise argparse.ArgumentTypeError('Boolean value expected.') + def clear_dat_folder(dat_folder_str=None): - #clear dat folder + # clear dat folder dat_path = Path(dat_folder_str) if not dat_path.exists(): dat_path.mkdir(parents=True) @@ -199,44 +213,49 @@ def clear_folder(folder_str): os.system('rm -rf ' + str(Path(folder_str)) + '/*') print(f'{folder_str} cleared') + def double_to_hex(d): return hex(struct.unpack('=0 - assert nabs<=1 + assert nabs >= 0 + assert nabs <= 1 return np.sqrt(1.0-nabs) -#specific admittance to reflection coefficient + def convert_Yn_to_R(Yn): - assert np.all(Yn>0.0) + """specific admittance to reflection coefficient + """ + assert np.all(Yn > 0.0) R = (1.0-Yn)/(1.0+Yn) return R -#reflection coefficient to specific admittance + def convert_R_to_Yn(R): - assert np.all(R<1.0) + """reflection coefficient to specific admittance + """ + assert np.all(R < 1.0) Yn = (1.0-R)/(1.0+R) return Yn -#reflection coefficient to specific impedance + def convert_R_to_Zn(R): - assert np.all(R<1.0) + """reflection coefficient to specific impedance + """ + assert np.all(R < 1.0) Zn = 1.0/convert_R_to_Yn(R) return Zn -#sabine absorption to specific admittance, Paris formula inversion with Newton method -def convert_Sabs_to_Yn(Sabs,max_iter=100): + +def convert_Sabs_to_Yn(Sabs, max_iter=100): + """sabine absorption to specific admittance, Paris formula inversion with Newton method + """ def _print(fstring): print(f'--MATS: {fstring}') - if Sabs>0.9512: + + def fg(g): + return 8.0*g*(1+g/(1+g)-2*g*np.log((g+1)/g)) + + def fgd(g): + return -8.0*(-4*g**2 - 6*g + 4 * (1+g)**2*g*np.log((g+1)/g)-1)/(1+g)**2 + + if Sabs > 0.9512: _print('warning, Sabs>0.9512 -- not possible for locally-reactive model') Sabs = 0.9512 - fg = lambda g: 8.0*g*(1+g/(1+g)-2*g*np.log((g+1)/g)) - fgd = lambda g: -8.0*(-4*g**2 - 6*g +4*(1+g)**2*g*np.log((g+1)/g)-1)/(1+g)**2 - if Sabs==0: - g=0 - else: #newton solve - x_old = Sabs/8.0 #starting guess + + if Sabs == 0: + g = 0 + else: # newton solve + x_old = Sabs/8.0 # starting guess x_new = 0 err = np.inf - niter=0 - while (niter 1e-6): + niter = 0 + while (niter < max_iter) and (err > 1e-6): x_new = x_old - (fg(x_old)-Sabs)/fgd(x_old) niter += 1 err = np.abs(1-x_new/x_old) - #_print(f'Iteration {niter}: {x_new=}, {err=}') - x_old = x_new; + # _print(f'Iteration {niter}: {x_new=}, {err=}') + x_old = x_new g = x_new return g -#write HDF5 mat file from a specific impedance (frequency-independent) -def write_freq_ind_mat_from_Zn(Zn,filename): + +def write_freq_ind_mat_from_Zn(Zn, filename): + """write HDF5 mat file from a specific impedance (frequency-independent) + """ def _print(fstring): print(f'--MATS: {fstring}') assert ~np.isnan(Zn) - assert ~np.isinf(Zn) #rigid should be specified in scene (no material) - assert Zn>=0 + assert ~np.isinf(Zn) # rigid should be specified in scene (no material) + assert Zn >= 0 filename = Path(filename) - h5f = h5py.File(filename,'w') - DEF = np.array([0,Zn,0]) - assert np.all(np.sum(DEF>0,axis=-1)) #at least one non-zero + h5f = h5py.File(filename, 'w') + DEF = np.array([0, Zn, 0]) + assert np.all(np.sum(DEF > 0, axis=-1)) # at least one non-zero _print(f'{DEF=}') - h5f.create_dataset(f'DEF', data=np.atleast_2d(DEF)) + h5f.create_dataset('DEF', data=np.atleast_2d(DEF)) h5f.close() -#write HDF5 mat file from a specific impedance (frequency-independent) -def write_freq_ind_mat_from_Yn(Yn,filename): + +def write_freq_ind_mat_from_Yn(Yn, filename): + """write HDF5 mat file from a specific impedance (frequency-independent) + """ def _print(fstring): print(f'--MATS: {fstring}') assert ~np.isnan(Yn) assert ~np.isinf(Yn) - assert Yn>0 #rigid should be specified in scene (no material) - write_freq_ind_mat_from_Zn(1/Yn,filename) + assert Yn > 0 # rigid should be specified in scene (no material) + write_freq_ind_mat_from_Zn(1/Yn, filename) + -#write HDF5 mat file from frequency-independent triplet (D=F=0) -def write_freq_dep_mat(DEF,filename): +def write_freq_dep_mat(DEF, filename): + """write HDF5 mat file from frequency-independent triplet (D=F=0) + """ def _print(fstring): print(f'--MATS: {fstring}') _print(f'{DEF=}') DEF = np.atleast_2d(DEF) - assert ~np.any(np.isinf(DEF)) #rigid should be specified in scene (no material) + # rigid should be specified in scene (no material) + assert ~np.any(np.isinf(DEF)) assert ~np.any(np.isnan(DEF)) - assert np.all(DEF>=0) #all non-neg - assert np.all(np.sum(DEF>0,axis=-1)) #at least one non-zero - assert DEF.shape[1]==3 + assert np.all(DEF >= 0) # all non-neg + assert np.all(np.sum(DEF > 0, axis=-1)) # at least one non-zero + assert DEF.shape[1] == 3 filename = Path(filename) - h5f = h5py.File(filename,'w') + h5f = h5py.File(filename, 'w') _print(f'{DEF=}') - h5f.create_dataset(f'DEF', data=DEF) + h5f.create_dataset('DEF', data=DEF) h5f.close() -#write HDF5 mat file from frequency-independent triplet (D=F=0) + def read_mat_DEF(filename): - h5f = h5py.File(Path(filename),'r') + """write HDF5 mat file from frequency-independent triplet (D=F=0) + """ + h5f = h5py.File(Path(filename), 'r') DEF = h5f['DEF'][()] h5f.close() return DEF -#plot some admittance based on RLC triplets (not specific admittance) -def plot_RLC_admittance(RLC,Y0): - plot_DEF_admittance(RLC*Y0) -#plot some admittance based on DEF triplets (specific admittance) -def plot_DEF_admittance(fv,DEF,model_Rf=None,show=True): +def plot_DEF_admittance(fv, DEF, model_Rf=None, show=True): + """plot some admittance based on DEF triplets (specific admittance) + """ DEF = np.atleast_2d(DEF) assert DEF.shape[1] == 3 - D,E,F = DEF.T + D, E, F = DEF.T jw = 1j*fv*2*np.pi - Rf,Yn,Zn_br,Rf_br = compute_Rf_from_DEF(jw,D,E,F) + Rf, Yn, Zn_br, Rf_br = compute_Rf_from_DEF(jw, D, E, F) if model_Rf is not None: model_Yn = convert_R_to_Yn(model_Rf) fig = plt.figure() - #reflection coefficient + # reflection coefficient ax = fig.add_subplot(2, 1, 1) - ax.plot(fv,np.abs(Rf),linestyle='-') - ax.plot(fv,np.abs(Rf_br),linestyle=':') + ax.plot(fv, np.abs(Rf), linestyle='-') + ax.plot(fv, np.abs(Rf_br), linestyle=':') if model_Rf is not None: - ax.plot(fv,np.abs(model_Rf),linestyle='--') + ax.plot(fv, np.abs(model_Rf), linestyle='--') ax.set_xscale('log') ax.set_xlabel('frequency (Hz)') ax.set_ylabel(r'$|R|$') @@ -150,10 +163,10 @@ def plot_DEF_admittance(fv,DEF,model_Rf=None,show=True): ax.grid(which='both', axis='both') ax = fig.add_subplot(2, 1, 2) - ax.plot(fv,np.angle(Rf),linestyle='-') - ax.plot(fv,np.angle(Rf_br),linestyle=':') + ax.plot(fv, np.angle(Rf), linestyle='-') + ax.plot(fv, np.angle(Rf_br), linestyle=':') if model_Rf is not None: - ax.plot(fv,np.angle(model_Rf),linestyle='--') + ax.plot(fv, np.angle(model_Rf), linestyle='--') ax.set_xscale('log') ax.set_xlabel('frequency (Hz)') ax.set_ylabel(r'$\angle R$') @@ -161,12 +174,12 @@ def plot_DEF_admittance(fv,DEF,model_Rf=None,show=True): ax.grid(which='both', axis='both') fig = plt.figure() - #specific (normalised) admittance (abs/ang) + # specific (normalised) admittance (abs/ang) ax = fig.add_subplot(2, 1, 1) - ax.plot(fv,np.abs(Yn),linestyle='-') - ax.plot(fv,np.abs(1.0/Zn_br),linestyle=':') + ax.plot(fv, np.abs(Yn), linestyle='-') + ax.plot(fv, np.abs(1.0/Zn_br), linestyle=':') if model_Rf is not None: - ax.plot(fv,np.abs(model_Yn),linestyle='--') + ax.plot(fv, np.abs(model_Yn), linestyle='--') ax.set_xscale('log') ax.set_xlabel('frequency (Hz)') ax.set_ylabel(r'$|Y|$') @@ -174,10 +187,10 @@ def plot_DEF_admittance(fv,DEF,model_Rf=None,show=True): ax.grid(which='both', axis='both') ax = fig.add_subplot(2, 1, 2) - ax.plot(fv,np.angle(Yn),linestyle='-') - ax.plot(fv,np.angle(1.0/Zn_br),linestyle=':') + ax.plot(fv, np.angle(Yn), linestyle='-') + ax.plot(fv, np.angle(1.0/Zn_br), linestyle=':') if model_Rf is not None: - ax.plot(fv,np.angle(model_Yn),linestyle='--') + ax.plot(fv, np.angle(model_Yn), linestyle='--') ax.set_xscale('log') ax.set_xlabel('frequency (Hz)') ax.set_ylabel(r'$\angle Y$') @@ -185,12 +198,12 @@ def plot_DEF_admittance(fv,DEF,model_Rf=None,show=True): ax.grid(which='both', axis='both') fig = plt.figure() - #specific (normalised) admittance (real/imag) + # specific (normalised) admittance (real/imag) ax = fig.add_subplot(2, 1, 1) - ax.plot(fv,np.real(Yn),linestyle='-') - ax.plot(fv,np.real(1.0/Zn_br),linestyle=':') + ax.plot(fv, np.real(Yn), linestyle='-') + ax.plot(fv, np.real(1.0/Zn_br), linestyle=':') if model_Rf is not None: - ax.plot(fv,np.real(model_Yn),linestyle='--') + ax.plot(fv, np.real(model_Yn), linestyle='--') ax.set_xscale('log') ax.set_xlabel('frequency (Hz)') ax.set_ylabel(r'$\Re(Y)$') @@ -198,10 +211,10 @@ def plot_DEF_admittance(fv,DEF,model_Rf=None,show=True): ax.grid(which='both', axis='both') ax = fig.add_subplot(2, 1, 2) - ax.plot(fv,np.imag(Yn),linestyle='-') - ax.plot(fv,np.imag(1.0/Zn_br),linestyle=':') + ax.plot(fv, np.imag(Yn), linestyle='-') + ax.plot(fv, np.imag(1.0/Zn_br), linestyle=':') if model_Rf is not None: - ax.plot(fv,np.imag(model_Yn),linestyle='--') + ax.plot(fv, np.imag(model_Yn), linestyle='--') ax.set_xscale('log') ax.set_xlabel('frequency (Hz)') ax.set_ylabel(r'$\Im(Y)$') @@ -211,37 +224,43 @@ def plot_DEF_admittance(fv,DEF,model_Rf=None,show=True): if show: plt.show() -#plot compute reflection coefficient from DEF triplets (specific admittance) -def compute_Rf_from_DEF(jw,D,E,F): - Zn_br = jw[:,None]*D[None,:] + E + F[None,:]/jw[:,None] - Yn = np.sum(1.0/Zn_br,axis=-1) + +def compute_Rf_from_DEF(jw, D, E, F): + """plot compute reflection coefficient from DEF triplets (specific admittance) + """ + Zn_br = jw[:, None]*D[None, :] + E + F[None, :]/jw[:, None] + Yn = np.sum(1.0/Zn_br, axis=-1) Rf = (1.0-Yn)/(1.0+Yn) Rf_br = (Zn_br-1.0)/(Zn_br+1.0) - return Rf,Yn,Zn_br,Rf_br + return Rf, Yn, Zn_br, Rf_br -#Ynm is max abs normalised admittance -#dw is half-power bandwidth in radians/s -#w0 is resonant frequency in radians/s -def _to_DEF(Ynm,dw,w0): + +def _to_DEF(Ynm, dw, w0): + # Ynm is max abs normalised admittance + # dw is half-power bandwidth in radians/s + # w0 is resonant frequency in radians/s D = 1.0/Ynm/dw E = 1.0/Ynm F = w0**2/Ynm/dw - return D,E,F + return D, E, F + -def _from_DEF(D,E,F): +def _from_DEF(D, E, F): Ynm = 1.0/E dw = E/D w0 = np.sqrt(F/D) - return Ynm,dw,w0 + return Ynm, dw, w0 + -#fit Yn from 11 Sabine octave-band coefficients (16Hz to 16KHz) -#this is a simple fitting routine that is a starting point for octave-band input data -def fit_to_Sabs_oct_11(Sabs,filename,plot=False): +def fit_to_Sabs_oct_11(Sabs, filename, plot=False): + # fit Yn from 11 Sabine octave-band coefficients (16Hz to 16KHz) + # this is a simple fitting routine that is a starting point for octave-band input data assert Sabs.size == 11 Noct = Sabs.size - fv = np.logspace(np.log10(10),np.log10(20e3),1000) #frequency vector to fit over + # frequency vector to fit over + fv = np.logspace(np.log10(10), np.log10(20e3), 1000) jw = 1j*fv*2*np.pi - fcv = 1000*(2.0**np.arange(-6,5)) + fcv = 1000*(2.0**np.arange(-6, 5)) ymv = np.zeros(Noct) dwv = np.zeros(Noct) w0v = np.zeros(Noct) @@ -250,69 +269,74 @@ def fit_to_Sabs_oct_11(Sabs,filename,plot=False): fc = fcv[j] sa = Sabs[j] Ynm = convert_Sabs_to_Yn(sa) - if j==0: - i1=0 + if j == 0: + i1 = 0 else: - i1 = np.flatnonzero(fv>=fc/np.sqrt(2))[0] - if j==Noct-1: - i2=fv.size + i1 = np.flatnonzero(fv >= fc/np.sqrt(2))[0] + if j == Noct-1: + i2 = fv.size else: - i2 = np.flatnonzero(fv>=fc*np.sqrt(2))[0] + i2 = np.flatnonzero(fv >= fc*np.sqrt(2))[0] Y_target[i1:i2] = Ynm w0 = 2*np.pi*fc - dw = w0/np.sqrt(2) #using resonance with half-octave bandwidth + dw = w0/np.sqrt(2) # using resonance with half-octave bandwidth ymv[j] = Ynm dwv[j] = dw w0v[j] = w0 R_target = (1.0-Y_target)/(1.0+Y_target) + def cost_function_3(x0): - if np.any(x0<0): - return np.finfo(np.float64).max #something big to disallow neg values - else: - x0 = x0.reshape((-1,3)) - ym = x0[:,0] - dw = x0[:,1] - w0 = x0[:,2] - - D,E,F = _to_DEF(ym,dw,w0) - - Rf_opt,Yn_opt,_,_ = compute_Rf_from_DEF(jw,D,E,F) - assert Yn_opt.ndim==1 - assert Rf_opt.ndim==1 - abs_opt = 1-np.abs(Rf_opt)**2 - abs_target = 1-np.abs(R_target)**2 - cost = np.sum(np.abs(abs_opt-abs_target)) #simply fit to absorption coefficients - #other possible cost functions - #cost = np.sum(np.abs(np.abs(Y_target) - np.abs(Yn_opt))) - #cost = np.sum(np.abs(1-np.abs(Y_target)/np.abs(Yn_opt))) - #cost = np.sum(np.abs(Y_target - Yn_opt)) - #cost = np.sum(np.abs(R_target - Rf_opt)) - return cost - - #limit optimisation to just Yabs (keep bandwidths fix and w0 fixed) - cost_function = lambda x0: cost_function_3(np.c_[x0,dwv,w0v].flat[:]) - - #run optimisation + if np.any(x0 < 0): + # something big to disallow neg values + return np.finfo(np.float64).max + + x0 = x0.reshape((-1, 3)) + ym = x0[:, 0] + dw = x0[:, 1] + w0 = x0[:, 2] + + D, E, F = _to_DEF(ym, dw, w0) + + Rf_opt, Yn_opt, _, _ = compute_Rf_from_DEF(jw, D, E, F) + assert Yn_opt.ndim == 1 + assert Rf_opt.ndim == 1 + abs_opt = 1-np.abs(Rf_opt)**2 + abs_target = 1-np.abs(R_target)**2 + # simply fit to absorption coefficients + cost = np.sum(np.abs(abs_opt-abs_target)) + # other possible cost functions + # cost = np.sum(np.abs(np.abs(Y_target) - np.abs(Yn_opt))) + # cost = np.sum(np.abs(1-np.abs(Y_target)/np.abs(Yn_opt))) + # cost = np.sum(np.abs(Y_target - Yn_opt)) + # cost = np.sum(np.abs(R_target - Rf_opt)) + return cost + + # limit optimisation to just Yabs (keep bandwidths fix and w0 fixed) + def cost_function(x0): + return cost_function_3(np.c_[x0, dwv, w0v].flat[:]) + + # run optimisation x0 = ymv initial_cost = cost_function(x0) - res = scpo.minimize(cost_function,x0,method='Nelder-Mead') #'Powell' works well too + # 'Powell' works well too + res = scpo.minimize(cost_function, x0, method='Nelder-Mead') final_cost = cost_function(res.x) - assert (final_cost <= initial_cost) + assert final_cost <= initial_cost ymv_opt = res.x dwv_opt = dwv w0v_opt = w0v - D,E,F = _to_DEF(ymv_opt,dwv_opt,w0v_opt) - DEF = np.c_[D,E,F] + D, E, F = _to_DEF(ymv_opt, dwv_opt, w0v_opt) + DEF = np.c_[D, E, F] - #now save - h5f = h5py.File(filename,'w') - assert np.all(np.sum(DEF>0,axis=-1)) #at least one non-zero + # now save + h5f = h5py.File(filename, 'w') + assert np.all(np.sum(DEF > 0, axis=-1)) # at least one non-zero print(f'{DEF=}') - h5f.create_dataset(f'DEF', data=np.atleast_2d(DEF)) + h5f.create_dataset('DEF', data=np.atleast_2d(DEF)) h5f.close() if plot: - plot_DEF_admittance(fv,np.c_[D,E,F],model_Rf = R_target) + plot_DEF_admittance(fv, np.c_[D, E, F], model_Rf=R_target) diff --git a/src/python/pffdtd/sim2d/fdtd.py b/src/python/pffdtd/sim2d/fdtd.py index 0046b20..d5ca8a3 100644 --- a/src/python/pffdtd/sim2d/fdtd.py +++ b/src/python/pffdtd/sim2d/fdtd.py @@ -46,8 +46,7 @@ def quantize(i): i_q = low if err_low < err_high else high if scale: return i_q*dx, err - else: - return i_q, err + return i_q, err x, y = pos xq, yq = quantize(x), quantize(y) @@ -89,7 +88,7 @@ def add_diffusor(prime, well_width, max_depth, room, in_mask, X, Y, dx, c, verbo print(f" error_d={derr_q/depth*100:.2f}%") print('--SIM-SETUP: Locate diffusor') - depths, g = primitive_root_diffuser(prime, g=None, depth=depth_q) + depths, _ = primitive_root_diffuser(prime, g=None, depth=depth_q) depths = quadratic_residue_diffuser(prime, depth_q) prime = depths.shape[0] for w in range(n): @@ -132,7 +131,7 @@ def stencil_air(u0, u1, u2, mask): @nb.njit(parallel=True) def stencil_boundary_rigid(u0, u1, u2, bn_ixy, adj_bn): - Nx, Ny = u1.shape + _, Ny = u1.shape Nb = bn_ixy.size for i in nb.prange(Nb): ib = bn_ixy[i] @@ -173,7 +172,6 @@ def main(): parser.add_argument('--save', action='store_true') parser.add_argument('--video', action='store_true') parser.add_argument('--verbose', action='store_true') - parser.add_argument('--apply_rigid', action=bool_action, default=True) parser.add_argument('--apply_loss', action=bool_action, default=True) args = parser.parse_args() @@ -192,10 +190,7 @@ def main(): duration = args.duration # seconds refl_coeff = 0.99 # reflection coefficient - apply_rigid = args.apply_rigid apply_loss = args.apply_loss - if apply_loss: - assert apply_rigid Lx, Ly = room[0], room[1] # box dims (with lower corner at origin) x_in, y_in = Lx*0.5, Ly*0.5 # source input position @@ -233,17 +228,16 @@ def main(): arc_center = (x_in, y_in-arc_radius) out_ixy, _ = make_receiver_arc(180, arc_center, arc_radius, dx, Nx, Ny) - if apply_rigid: - print('--SIM-SETUP: Create node ABCs') - # Calculate number of interior neighbours (for interior points only) - K_map = np.zeros((Nx, Ny), dtype=int) - K_map[1:-1, 1:-1] += in_mask[2:, 1:-1] - K_map[1:-1, 1:-1] += in_mask[:-2, 1:-1] - K_map[1:-1, 1:-1] += in_mask[1:-1, 2:] - K_map[1:-1, 1:-1] += in_mask[1:-1, :-2] - K_map[~in_mask] = 0 - bn_ixy = np.where((K_map.flat > 0) & (K_map.flat < 4))[0] - adj_bn = K_map.flat[bn_ixy] + print('--SIM-SETUP: Create node ABCs') + # Calculate number of interior neighbours (for interior points only) + K_map = np.zeros((Nx, Ny), dtype=int) + K_map[1:-1, 1:-1] += in_mask[2:, 1:-1] + K_map[1:-1, 1:-1] += in_mask[:-2, 1:-1] + K_map[1:-1, 1:-1] += in_mask[1:-1, 2:] + K_map[1:-1, 1:-1] += in_mask[1:-1, :-2] + K_map[~in_mask] = 0 + bn_ixy = np.where((K_map.flat > 0) & (K_map.flat < 4))[0] + adj_bn = K_map.flat[bn_ixy] # Grid forcing points inx = int(np.round(x_in / dx + 0.5) + 1) @@ -324,10 +318,9 @@ def main(): if args.video: for nt in tqdm(range(Nt)): stencil_air(u0, u1, u2, in_mask) - if apply_rigid: - stencil_boundary_rigid(u0, u1, u2, bn_ixy, adj_bn) - if apply_loss: - stencil_boundary_loss(u0, u2, bn_ixy, adj_bn, loss_factor) + stencil_boundary_rigid(u0, u1, u2, bn_ixy, adj_bn) + if apply_loss: + stencil_boundary_loss(u0, u2, bn_ixy, adj_bn, loss_factor) u0[inx, iny] = u0[inx, iny] + src_sig[nt] diff --git a/src/python/pffdtd/sim2d/report.py b/src/python/pffdtd/sim2d/report.py index dd43e75..c50ecdf 100644 --- a/src/python/pffdtd/sim2d/report.py +++ b/src/python/pffdtd/sim2d/report.py @@ -3,7 +3,7 @@ import h5py import matplotlib.pyplot as plt import numpy as np -import scipy.signal as signal +from scipy import signal from pffdtd.diffusor.measurement import polar_response @@ -72,7 +72,7 @@ def main(): constrained_layout=True, subplot_kw={'projection': 'polar'} ) - fig.suptitle(f"Diffusion") + fig.suptitle("Diffusion") ax[0][0].plot(np.deg2rad(mic_angles), rms_values[0][0]) ax[0][0].set_title(rms_values[0][1]) diff --git a/src/python/pffdtd/sim3d/engine.py b/src/python/pffdtd/sim3d/engine.py index 0b1b251..9748917 100644 --- a/src/python/pffdtd/sim3d/engine.py +++ b/src/python/pffdtd/sim3d/engine.py @@ -21,16 +21,16 @@ # - Plots simulations (mayavi is best, matplotlib is fallback) # ############################################################################## +import json +import time +from pathlib import Path +import h5py import numpy as np -from numpy import array as npa import numba as nb -from pathlib import Path -from pffdtd.common.timerdict import TimerDict from tqdm import tqdm -import time -import h5py -import json as json + +from pffdtd.common.timerdict import TimerDict from pffdtd.common.myfuncs import ind2sub3d,rel_diff,get_default_nprocs MMb = 12 #max allowed number of branches @@ -75,7 +75,7 @@ def load_h5_data(self): saf_bn = h5f['saf_bn'][...] h5f.close() - ii = (mat_bn>-1) + ii = mat_bn>-1 self.saf_bnl = saf_bn[ii] self.mat_bnl = mat_bn[ii] self.bnl_ixyz = self.bn_ixyz[ii] @@ -110,19 +110,19 @@ def load_h5_data(self): self.print(f'Nr={self.Nr} Ns={self.Ns} Nt={self.Nt}') if self.fcc: - assert (self.Nx%2)==0 - assert (self.Ny%2)==0 - assert (self.Nz%2)==0 - self.print(f'On FCC subgrid') - assert (self.adj_bn.shape[1] == 12) + assert self.Nx%2 ==0 + assert self.Ny%2 ==0 + assert self.Nz%2 ==0 + self.print('On FCC subgrid') + assert self.adj_bn.shape[1] == 12 if self.fcc: self.ssaf_bnl = self.saf_bnl*0.5/np.sqrt(2.0) #rescaled by S*h/V - assert (self.l<=1.0) - assert (self.l2<=1.0) + assert self.l<=1.0 + assert self.l2<=1.0 else: self.ssaf_bnl = self.saf_bnl - assert (self.l<=np.sqrt(1/3)) - assert (self.l2<=1/3) + assert self.l<=np.sqrt(1/3) + assert self.l2<=1/3 h5f = h5py.File(Path(data_dir / Path('materials.h5')),'r') @@ -160,7 +160,7 @@ def _load_abc(self): self.bna_ixyz = bna_ixyz self.Nba = Nba if self.energy_on: - self.V_bna = (2.0**-Q_bna) + self.V_bna = 2.0**-Q_bna def allocate_mem(self): self.print('allocating mem..') @@ -169,7 +169,6 @@ def allocate_mem(self): Nx = self.Nx Ny = self.Ny Nz = self.Nz - Nm = self.Nm u0 = np.zeros((Nx,Ny,Nz),dtype=np.float64) u1 = np.zeros((Nx,Ny,Nz),dtype=np.float64) @@ -214,7 +213,6 @@ def setup_mask(self): self.bn_mask = bn_mask def set_coeffs(self): - l = self.l l2 = self.l2 Ts = self.Ts mat_bnl = self.mat_bnl @@ -223,12 +221,12 @@ def set_coeffs(self): Nm = self.Nm if self.fcc: - assert(l2<=1.0) - a1 = (2.0-l2*3.0) #0.25*12 + assert l2<=1.0 + a1 = 2.0-l2*3.0 #0.25*12 a2 = 0.25*l2 else: - assert(l2<=1/3) - a1 = (2.0-l2*6.0) + assert l2<=1/3 + a1 = 2.0-l2*6.0 a2 = l2 av = np.array([a1,a2],dtype=np.float64) @@ -251,7 +249,7 @@ def set_coeffs(self): Fh = F*Ts b = 1.0/(2.0*Dh+Eh+0.5*Fh) - d = (2.0*Dh-Eh-0.5*Fh) + d = 2.0*Dh-Eh-0.5*Fh assert ~np.any(np.isnan(np.r_[b,d])) assert ~np.any(np.isinf(np.r_[b,d])) @@ -306,8 +304,8 @@ def run_all(self,nsteps=1): timer = TimerDict() pbar = {} - pbar['vox']= tqdm(total=Nt*Npts,desc=f'FDTD run',unit='vox',unit_scale=True,ascii=True,leave=False,position=0,dynamic_ncols=True) - pbar['samples'] = tqdm(total=Nt,desc=f'FDTD run',unit='samples',unit_scale=True,ascii=True,leave=False,position=1,ncols=0) + pbar['vox']= tqdm(total=Nt*Npts,desc='FDTD run',unit='vox',unit_scale=True,ascii=True,leave=False,position=0,dynamic_ncols=True) + pbar['samples'] = tqdm(total=Nt,desc='FDTD run',unit='samples',unit_scale=True,ascii=True,leave=False,position=1,ncols=0) timer.tic('run') for n in range(0,Nt,nsteps): @@ -330,7 +328,6 @@ def run_plot(self,nsteps=1,draw_backend='mayavi',json_model=None): Ny = self.Ny Nz = self.Nz Nt = self.Nt - Npts = Nx*Ny*Nz bn_mask = self.bn_mask in_ixyz = self.in_ixyz ix,iy,iz = ind2sub3d(in_ixyz,Nx,Ny,Nz) @@ -401,15 +398,15 @@ def run_plot(self,nsteps=1,draw_backend='mayavi',json_model=None): mat_str = list(mats_dict.keys()) #for section cuts - edges_xy_0 = npa([[],[]]).T - edges_xy_1 = npa([[],[]]).T - edges_xz_0 = npa([[],[]]).T - edges_xz_1 = npa([[],[]]).T - edges_yz_0 = npa([[],[]]).T - edges_yz_1 = npa([[],[]]).T + edges_xy_0 = np.array([[],[]]).T + edges_xy_1 = np.array([[],[]]).T + edges_xz_0 = np.array([[],[]]).T + edges_xz_1 = np.array([[],[]]).T + edges_yz_0 = np.array([[],[]]).T + edges_yz_1 = np.array([[],[]]).T for mat in mat_str: - pts = npa(mats_dict[mat]['pts'],dtype=np.float64) #no rotation - tris = npa(mats_dict[mat]['tris'],dtype=np.int64) + pts = np.array(mats_dict[mat]['pts'],dtype=np.float64) #no rotation + tris = np.array(mats_dict[mat]['tris'],dtype=np.int64) for j in range(3): ii = np.nonzero((pts[tris[:,j%3],2]>zv[iz_in]) ^ (pts[tris[:,(j+1)%3],2]>zv[iz_in]))[0] edges_xy_0 = np.r_[edges_xy_0,pts[tris[ii,j%3]][:,[0,1]]] @@ -426,18 +423,18 @@ def run_plot(self,nsteps=1,draw_backend='mayavi',json_model=None): edges_yz_1 = np.r_[edges_yz_1,pts[tris[ii,(j+1)%3]][:,[1,2]]] del ii - u = edges_xy_1[:,0]-edges_xy_0[:,0]; - v = edges_xy_1[:,1]-edges_xy_0[:,1]; + u = edges_xy_1[:,0]-edges_xy_0[:,0] + v = edges_xy_1[:,1]-edges_xy_0[:,1] w=np.zeros(u.shape) mlab.quiver3d(edges_xy_0[:,0],edges_xy_0[:,1],np.full(u.shape,zv[iz_in]),u,v,w,mode='2ddash',scale_factor=1.,color=(0,0,0)) - u = edges_xz_1[:,0]-edges_xz_0[:,0]; + u = edges_xz_1[:,0]-edges_xz_0[:,0] v=np.zeros(u.shape) - w = edges_xz_1[:,1]-edges_xz_0[:,1]; + w = edges_xz_1[:,1]-edges_xz_0[:,1] mlab.quiver3d(edges_xz_0[:,0],np.full(u.shape,yv[iy_in]),edges_xz_0[:,1],u,v,w,mode='2ddash',scale_factor=1.,color=(0,0,0)) - v = edges_yz_1[:,0]-edges_yz_0[:,0]; - w = edges_yz_1[:,1]-edges_yz_0[:,1]; + v = edges_yz_1[:,0]-edges_yz_0[:,0] + w = edges_yz_1[:,1]-edges_yz_0[:,1] u=np.zeros(w.shape) mlab.quiver3d(np.full(u.shape,xv[ix_in]),edges_yz_0[:,0],edges_yz_0[:,1],u,v,w,mode='2ddash',scale_factor=1.,color=(0,0,0)) @@ -464,14 +461,24 @@ def run_plot(self,nsteps=1,draw_backend='mayavi',json_model=None): mlab.orientation_axes() - fake_verts = np.c_[npa([xv[0],yv[0],zv[0]]),npa([xv[-1],yv[-1],zv[-1]])] + fake_verts = np.c_[np.array([xv[0],yv[0],zv[0]]),np.array([xv[-1],yv[-1],zv[-1]])] mlab.plot3d(*fake_verts,transparent=True,opacity=0) mlab.axes(xlabel='x', ylabel='y', zlabel='z', color=(0., 0., 0.)) #fix z-up fig.scene.interactor.interactor_style = tvtk.InteractorStyleTerrain() - pbar = tqdm(total=Nt,desc=f'FDTD run',unit='samples',unit_scale=True,ascii=True,leave=False,position=1,ncols=0) + pbar = tqdm( + total=Nt, + desc='FDTD run', + unit='samples', + unit_scale=True, + ascii=True, + leave=False, + position=1, + ncols=0, + ) + for n in range(0,Nt,nsteps): nrun = min(nsteps,Nt-n) self.run_steps(n,nrun) @@ -502,7 +509,7 @@ def run_plot(self,nsteps=1,draw_backend='mayavi',json_model=None): cmax_xy = np.max(np.abs(uxy.flat[:])) cmax_xz = np.max(np.abs(uxz.flat[:])) cmax_yz = np.max(np.abs(uyz.flat[:])) - cmax = np.max(npa([cmax_xy,cmax_yz,cmax_xz])) + cmax = np.max(np.array([cmax_xy,cmax_yz,cmax_xz])) hh_xy.mlab_source.scalars = ldraw(uxy.T) hh_xy.mlab_source.update() @@ -549,7 +556,6 @@ def run_steps(self,nstart,nsteps): in_ixyz = self.in_ixyz out_ixyz = self.out_ixyz energy_on = self.energy_on - av = self.av l = self.l l2 = self.l2 u2b = self.u2b @@ -694,11 +700,10 @@ def print_last_energy(self,Np): def save_outputs(self): data_dir = self.data_dir u_out = self.u_out - Ts = self.Ts out_reorder = self.out_reorder #just raw outputs, recombine elsewhere h5f = h5py.File(data_dir / Path('sim_outs.h5'),'w') - h5f.create_dataset(f'u_out', data=u_out[out_reorder,:]) + h5f.create_dataset('u_out', data=u_out[out_reorder,:]) h5f.close() self.print('saved outputs in {data_dir}') @@ -740,7 +745,7 @@ def nb_stencil_air_fcc(Lu1,u1,bn_mask): @nb.jit(nopython=True,parallel=True) def nb_stencil_bn_fcc(Lu1,u1,bn_ixyz,adj_bn): - Nx,Ny,Nz = u1.shape + _,Ny,Nz = u1.shape Nb = bn_ixyz.size for i in nb.prange(Nb): K = np.sum(adj_bn[i,:]) @@ -762,7 +767,7 @@ def nb_stencil_bn_fcc(Lu1,u1,bn_ixyz,adj_bn): @nb.jit(nopython=True,parallel=True) def nb_stencil_bn_cart(Lu1,u1,bn_ixyz,adj_bn): - Nx,Ny,Nz = u1.shape + _,Ny,Nz = u1.shape Nb = bn_ixyz.size for i in nb.prange(Nb): K = np.sum(adj_bn[i,:]) @@ -879,11 +884,11 @@ def nb_get_abc_ib(bna_ixyz,Q_bna,Nx,Ny,Nz,fcc): if fcc and (ix+iy+iz)%2==1: continue Q = 0 - if ((ix==1) or (ix==Nx-2)): + if ix in (1, Nx-2): Q+=1 - if ((iy==1) or (iy==Ny-2)): + if iy in (1, Ny-2): Q+=1 - if ((iz==1) or (iz==Nz-2)): + if iz in (1, Nz-2): Q+=1 if Q>0: bna_ixyz[ii] = ix*Ny*Nz + iy*Nz + iz diff --git a/src/python/pffdtd/sim3d/sim_constants.py b/src/python/pffdtd/sim3d/sim_constants.py index 5cfdc07..e19577b 100644 --- a/src/python/pffdtd/sim3d/sim_constants.py +++ b/src/python/pffdtd/sim3d/sim_constants.py @@ -1,9 +1,8 @@ # SPDX-License-Identifier: MIT +from pathlib import Path import numpy as np -from numpy import array as npa import h5py -from pathlib import Path class SimConstants: diff --git a/src/python/pffdtd/sim3d/sim_setup.py b/src/python/pffdtd/sim3d/sim_setup.py index c18bfa9..5865a17 100644 --- a/src/python/pffdtd/sim3d/sim_setup.py +++ b/src/python/pffdtd/sim3d/sim_setup.py @@ -12,9 +12,9 @@ # ############################################################################## -import numpy as np -import json as json from pathlib import Path + +import numpy as np from numpy import array as npa from pffdtd.common.room_geo import RoomGeo @@ -69,74 +69,76 @@ def sim_setup( assert duration is not None if (bmin is not None) and (bmax is not None): - #custom bmin/bmax (for open scenes) - bmin = npa(bmin,dtype=np.float64) - bmax = npa(bmax,dtype=np.float64) + # custom bmin/bmax (for open scenes) + bmin = npa(bmin, dtype=np.float64) + bmax = npa(bmax, dtype=np.float64) - #set up room geometry (reads in JSON export, rotates scene) - room_geo = RoomGeo(model_json_file,az_el=rot_az_el,bmin=bmin,bmax=bmax) + # set up room geometry (reads in JSON export, rotates scene) + room_geo = RoomGeo(model_json_file, az_el=rot_az_el, bmin=bmin, bmax=bmax) room_geo.print_stats() - #sources have to be specified in advance (edit JSON if necessary) - Sxyz = room_geo.Sxyz[source_num-1] #one source (one-based indexing) - Rxyz = room_geo.Rxyz #many receivers + # sources have to be specified in advance (edit JSON if necessary) + Sxyz = room_geo.Sxyz[source_num-1] # one source (one-based indexing) + Rxyz = room_geo.Rxyz # many receivers - #some constants for the simulation, in one place - constants = SimConstants(Tc=Tc,rh=rh,fmax=fmax,PPW=PPW,fcc=fcc_flag) + # some constants for the simulation, in one place + constants = SimConstants(Tc=Tc, rh=rh, fmax=fmax, PPW=PPW, fcc=fcc_flag) constants.save(save_folder) - #link up the wall materials to impedance datasets + # link up the wall materials to impedance datasets materials = SimMaterials(save_folder=save_folder) - materials.package(mat_files_dict=mat_files_dict,mat_list=room_geo.mat_str,read_folder=mat_folder) + materials.package(mat_files_dict=mat_files_dict, + mat_list=room_geo.mat_str, read_folder=mat_folder) - #set the cartesian grid (also for FCC) - cart_grid = CartGrid(h=constants.h,offset=3.5,bmin=room_geo.bmin,bmax=room_geo.bmax,fcc=fcc_flag) + # set the cartesian grid (also for FCC) + cart_grid = CartGrid(h=constants.h, offset=3.5, + bmin=room_geo.bmin, bmax=room_geo.bmax, fcc=fcc_flag) cart_grid.print_stats() cart_grid.save(save_folder) - #set up source/receiver positions and input signals - sim_comms = SimSignals(save_folder=save_folder) #reads from cart_grid + # set up source/receiver positions and input signals + sim_comms = SimSignals(save_folder=save_folder) # reads from cart_grid sim_comms.prepare_source_pts(Sxyz) sim_comms.prepare_receiver_pts(Rxyz) - sim_comms.prepare_source_signals(duration,sig_type=insig_type) + sim_comms.prepare_source_signals(duration, sig_type=insig_type) if diff_source: sim_comms.diff_source() sim_comms.save(compress=compress) - #set up the voxel grid (volume hierarchy for ray-triangle intersections) - vox_grid = VoxGrid(room_geo,cart_grid,Nvox_est=Nvox_est,Nh=Nh) + # set up the voxel grid (volume hierarchy for ray-triangle intersections) + vox_grid = VoxGrid(room_geo, cart_grid, Nvox_est=Nvox_est, Nh=Nh) vox_grid.fill(Nprocs=1) vox_grid.print_stats() - #'voxelize' the scene (calculate FDTD mesh adjacencies and identify/correct boundary surfaces) - vox_scene = VoxScene(room_geo,cart_grid,vox_grid,fcc=fcc_flag) + # 'voxelize' the scene (calculate FDTD mesh adjacencies and identify/correct boundary surfaces) + vox_scene = VoxScene(room_geo, cart_grid, vox_grid, fcc=fcc_flag) vox_scene.calc_adj(Nprocs=1) vox_scene.check_adj_full() - vox_scene.save(save_folder,compress=compress) + vox_scene.save(save_folder, compress=compress) - #check that source/receivers don't intersect with boundaries + # check that source/receivers don't intersect with boundaries sim_comms.check_for_clashes(vox_scene.bn_ixyz) - #make copy for sorting/rotation for gpu + # make copy for sorting/rotation for gpu if save_folder_gpu is not None and Path(save_folder_gpu) != Path(save_folder): - copy_sim_data(save_folder,save_folder_gpu) + copy_sim_data(save_folder, save_folder_gpu) if save_folder_gpu is not None: rotate(save_folder_gpu) if fcc_flag: fold_fcc_sim_data(save_folder_gpu) sort_sim_data(save_folder_gpu) - #extra folder for cpu version (if needed for testing) - #if save_folder_cpu is not None and Path(save_folder_cpu) != Path(save_folder): - #copy_sim_data(save_folder,save_folder_cpu) - #if save_folder_cpu is not None: - #rotate(save_folder_cpu) - #if fcc_flag: - #fold_fcc_sim_data(save_folder_gpu) - #sort_sim_data(save_folder_cpu) + # extra folder for cpu version (if needed for testing) + # if save_folder_cpu is not None and Path(save_folder_cpu) != Path(save_folder): + # copy_sim_data(save_folder,save_folder_cpu) + # if save_folder_cpu is not None: + # rotate(save_folder_cpu) + # if fcc_flag: + # fold_fcc_sim_data(save_folder_gpu) + # sort_sim_data(save_folder_cpu) - #draw the voxelisation (use polyscope for dense grids) + # draw the voxelisation (use polyscope for dense grids) if draw_vox: - room_geo.draw(wireframe=False,backend=draw_backend) + room_geo.draw(wireframe=False, backend=draw_backend) vox_scene.draw(backend=draw_backend) room_geo.show(backend=draw_backend) diff --git a/src/python/pffdtd/sim3d/sim_signals.py b/src/python/pffdtd/sim3d/sim_signals.py index 087523d..8a33089 100644 --- a/src/python/pffdtd/sim3d/sim_signals.py +++ b/src/python/pffdtd/sim3d/sim_signals.py @@ -11,17 +11,16 @@ # Description: Class for source/receiver positions and input signals # ############################################################################## +from pathlib import Path +import h5py import numpy as np from numpy import array as npa -from pffdtd.voxelizer.cart_grid import CartGrid -from pathlib import Path -from pffdtd.common.timerdict import TimerDict -from pffdtd.common.myfuncs import iceil -import h5py -from scipy.signal import lfilter,bilinear_zpk from numpy import pi,cos,sin +from scipy.signal import lfilter +from pffdtd.common.timerdict import TimerDict +from pffdtd.common.myfuncs import iceil class SimSignals: def __init__(self,save_folder): #will read h,xv,yv,zv from h5 data @@ -44,9 +43,9 @@ def __init__(self,save_folder): self.fcc = (self.fcc_flag>0) if self.fcc: - assert (self.xv.size%2)==0 - assert (self.yv.size%2)==0 - assert (self.zv.size%2)==0 + assert self.xv.size%2 == 0 + assert self.yv.size%2 == 0 + assert self.zv.size%2 == 0 self.save_folder = save_folder self._diff = False @@ -81,7 +80,7 @@ def prepare_source_signals(self,duration,sig_type='impulse'): elif sig_type=='dhann30': #symmetric, for viz N = 30 n = np.arange(N) - in_sig[:N] = cos(pi*n/N)*sin(pi*n/N); + in_sig[:N] = cos(pi*n/N)*sin(pi*n/N) elif sig_type=='hann5ms': #for consistency checking N = iceil(5e-3/Ts) n = np.arange(N) @@ -113,7 +112,7 @@ def diff_source(self): b = 2/Ts*npa([1.0,-1.0]) #don't need Ts scaling but simpler to keep for receiver post-processing a = npa([1.0,1.0]) - in_sigs = lfilter(b,a,in_sigs,axis=-1); + in_sigs = lfilter(b,a,in_sigs,axis=-1) self._diff = True self.in_sigs = in_sigs @@ -236,14 +235,14 @@ def _check_for_clashes(_ixyz,bn_ixyz): ixyz = np.unique(_ixyz) #can have duplicates in receivers #could speed this up with a numba routine (don't need to know clashes, just say yes or no) assert np.union1d(ixyz.flat[:],bn_ixyz).size == ixyz.size + bn_ixyz.size - self.print(f'intersection with boundaries: passed') + self.print('intersection with boundaries: passed') timer = TimerDict() timer.tic('check in_xyz') - self.print(f'boundary intersection check with in_ixyz..') + self.print('boundary intersection check with in_ixyz..') _check_for_clashes(self.in_ixyz,bn_ixyz) self.print(timer.ftoc('check in_xyz')) timer.tic('check in_xyz') - self.print(f'boundary intersection check with out_ixyz..') + self.print('boundary intersection check with out_ixyz..') _check_for_clashes(self.out_ixyz,bn_ixyz) self.print(timer.ftoc('check in_xyz')) diff --git a/src/python/pffdtd/voxelizer/cart_grid.py b/src/python/pffdtd/voxelizer/cart_grid.py index 0b559ea..aada9ac 100644 --- a/src/python/pffdtd/voxelizer/cart_grid.py +++ b/src/python/pffdtd/voxelizer/cart_grid.py @@ -1,55 +1,48 @@ -############################################################################## -# This file is a part of PFFDTD. -# -# PFFTD is released under the MIT License. -# For details see the LICENSE file. -# -# Copyright 2021 Brian Hamilton. -# -# File name: cart_grid.py -# -# Description: This is a class for a Cartesian grid with grid spacing 'h' and -# some bounds. Saves minimal data in HDF5 file if needed outside. -# -############################################################################## +# SPDX-License-Identifier: MIT import numpy as np from numpy import array as npa from pathlib import Path import h5py + + class CartGrid(): - def __init__(self,h=None,offset=None,bmin=None,bmax=None,fcc=False): + """This is a class for a Cartesian grid with grid spacing 'h' and some bounds. + Saves minimal data in HDF5 file if needed outside. + """ + + def __init__(self, h=None, offset=None, bmin=None, bmax=None, fcc=False): assert h is not None assert offset is not None - assert bmin is not None #bbox min - assert bmax is not None #bbox max + assert bmin is not None # bbox min + assert bmax is not None # bbox max - #ensure three-layer halo, for ABCs, etc. - assert offset>2.0 #distance offset + # ensure three-layer halo, for ABCs, etc. + assert offset > 2.0 # distance offset xyzmin0 = bmin - offset*h xyzmax0 = bmax + offset*h - #64-bit ints signed (for simplicity and numpy indexing) - Nx,Ny,Nz = np.int_(np.ceil((xyzmax0 - xyzmin0)/h))+1 - if fcc: #make all dims even (so we can rotate any and halve) - Nx += (Nx%2) - Ny += (Ny%2) - Nz += (Nz%2) + # 64-bit ints signed (for simplicity and numpy indexing) + Nx, Ny, Nz = np.int_(np.ceil((xyzmax0 - xyzmin0)/h))+1 + if fcc: # make all dims even (so we can rotate any and halve) + Nx += (Nx % 2) + Ny += (Ny % 2) + Nz += (Nz % 2) self.print(f'To use FCC subgrid') - #grid vectors - xv,yv,zv = np.ogrid[0.:Nx,0.:Ny,0.:Nz] + # grid vectors + xv, yv, zv = np.ogrid[0.:Nx, 0.:Ny, 0.:Nz] xv = xv.ravel()*h + xyzmin0[0] yv = yv.ravel()*h + xyzmin0[1] zv = zv.ravel()*h + xyzmin0[2] - xyzmin = np.array([xv[0],yv[0],zv[0]]) - xyzmax = np.array([xv[-1],yv[-1],zv[-1]]) + xyzmin = np.array([xv[0], yv[0], zv[0]]) + xyzmax = np.array([xv[-1], yv[-1], zv[-1]]) - assert np.all(xyzmin==xyzmin0) - assert np.all(xyzmax>=xyzmax0) + assert np.all(xyzmin == xyzmin0) + assert np.all(xyzmax >= xyzmax0) self.print(f'{Nx=}, {Ny=}, {Nz=}, N_grid_pts={Nx*Ny*Nz:g}') if fcc: @@ -67,28 +60,34 @@ def __init__(self,h=None,offset=None,bmin=None,bmax=None,fcc=False): self.Nx = Nx self.Ny = Ny self.Nz = Nz - self.Nxyz = npa([Nx,Ny,Nz]) + self.Nxyz = npa([Nx, Ny, Nz]) self.Npts = np.prod(self.Nxyz) self.xyzmin = xyzmin self.xyzmax = xyzmax - def print(self,fstring): + def print(self, fstring): print(f'--CART_GRID: {fstring}') - #don't use this unless grid is small - def draw_gridpoints(self,backend='mayavi'): + # don't use this unless grid is small + def draw_gridpoints(self, backend='mayavi'): xv = self.xv yv = self.yv zv = self.zv h = self.h - X,Y,Z = np.meshgrid(xv,yv,zv,indexing='ij') - if backend=='mayavi': + X, Y, Z = np.meshgrid(xv, yv, zv, indexing='ij') + if backend == 'mayavi': from mayavi import mlab - mlab.points3d(X.flat[:],Y.flat[:],Z.flat[:],color=(0,0,0),mode='sphere',resolution=4,scale_factor=h/4) + mlab.points3d( + X.flat[:], Y.flat[:], Z.flat[:], + color=(0, 0, 0), + mode='sphere', + resolution=4, + scale_factor=h/4 + ) mlab.draw() else: - raise #todo polyscope + raise # todo polyscope def print_stats(self): cg = self @@ -98,8 +97,8 @@ def print_stats(self): self.print(f'{cg.xyzmin=}') self.print(f'{cg.xyzmax=}') - #save in HDF5 file - def save(self,save_folder): + # save in HDF5 file + def save(self, save_folder): xv = self.xv yv = self.yv zv = self.zv @@ -113,7 +112,7 @@ def save(self,save_folder): assert save_folder.is_dir() kw = {'compression': "gzip", 'compression_opts': 9} - h5f = h5py.File(save_folder / Path('cart_grid.h5'),'w') + h5f = h5py.File(save_folder / Path('cart_grid.h5'), 'w') h5f.create_dataset('xv', data=xv, **kw) h5f.create_dataset('yv', data=yv, **kw) h5f.create_dataset('zv', data=zv, **kw) diff --git a/src/python/pffdtd/voxelizer/vox_grid.py b/src/python/pffdtd/voxelizer/vox_grid.py index cb295a4..5df7f47 100644 --- a/src/python/pffdtd/voxelizer/vox_grid.py +++ b/src/python/pffdtd/voxelizer/vox_grid.py @@ -1,32 +1,23 @@ -############################################################################## -# This file is a part of PFFDTD. -# -# PFFTD is released under the MIT License. -# For details see the LICENSE file. -# -# Copyright 2021 Brian Hamilton. -# -# File name: vox_grid.py -# -# Description: VoxGrid class for voxelizer -# Uses multiprocessing. -# -# Notes: -# - Performance will depend on geometry, grid spacing, voxel size and # processes -# - Simple *heuristic* default auto-tuning provided -# - Not designed for more than 10^6 voxels (initialisation becomes bottleneck) -# - Expected to run on a powerful CPU (+4 cores with SMT). -# -############################################################################## +# SPDX-License-Identifier: MIT + +""" +VoxGrid class for voxelizer. Uses multiprocessing. + +Notes: + - Performance will depend on geometry, grid spacing, voxel size and # processes + - Simple *heuristic* default auto-tuning provided + - Not designed for more than 10^6 voxels (initialisation becomes bottleneck) + - Expected to run on a powerful CPU (+4 cores with SMT). +""" import numpy as np from numpy import array as npa -from pffdtd.common.room_geo import RoomGeo -from pffdtd.voxelizer.vox_grid_base import VoxGridBase,VoxBase -from pffdtd.voxelizer.cart_grid import CartGrid -from pffdtd.common.timerdict import TimerDict from tqdm import tqdm + from pffdtd.common.myfuncs import get_default_nprocs,iceil +from pffdtd.common.room_geo import RoomGeo +from pffdtd.voxelizer.cart_grid import CartGrid +from pffdtd.voxelizer.vox_grid_base import VoxGridBase,VoxBase class Voxel(VoxBase): #using cubic voxels for simplicity diff --git a/src/python/pffdtd/voxelizer/vox_grid_base.py b/src/python/pffdtd/voxelizer/vox_grid_base.py index bbd6724..a8c4371 100644 --- a/src/python/pffdtd/voxelizer/vox_grid_base.py +++ b/src/python/pffdtd/voxelizer/vox_grid_base.py @@ -12,23 +12,15 @@ # Uses multiprocessing # ############################################################################## +import multiprocessing as mp +from multiprocessing import shared_memory import numpy as np -from numpy import array as npa -from pffdtd.common.timerdict import TimerDict -from pffdtd.common.tri_box_intersection import tri_box_intersection_vec -from pffdtd.common.room_geo import RoomGeo -from pffdtd.common.tris_precompute import tris_precompute -import multiprocessing as mp -from pffdtd.common.myfuncs import clear_dat_folder,clear_console -from pffdtd.common.myfuncs import get_default_nprocs from tqdm import tqdm -import pffdtd.common.check_version as cv -import time -assert cv.ATLEASTVERSION38 #for shared memory (but project needs 3.9 anyway) - -from multiprocessing import shared_memory +from pffdtd.common.timerdict import TimerDict +from pffdtd.common.tri_box_intersection import tri_box_intersection_vec +from pffdtd.common.myfuncs import clear_dat_folder, get_default_nprocs #base class for a voxel class VoxBase: diff --git a/src/python/pffdtd/voxelizer/vox_scene.py b/src/python/pffdtd/voxelizer/vox_scene.py index 86b5b80..ea83bc9 100644 --- a/src/python/pffdtd/voxelizer/vox_scene.py +++ b/src/python/pffdtd/voxelizer/vox_scene.py @@ -1,60 +1,49 @@ -############################################################################## -# This file is a part of PFFDTD. -# -# PFFTD is released under the MIT License. -# For details see the LICENSE file. -# -# Copyright 2021 Brian Hamilton. -# -# File name: vox_scene.py -# -# Description: a scene voxelizer for FDTD -# -# Sets up primary input for FDTD (adjacencies, materials, boundary nodes) -# CPU-based and uses multiprocessing. -# Tested up to ~10^10 points. Meant for single-node use (no MPI) -# -# Notes: -# - Performance will depend on geometry, grid spacing, voxel size and # processes -# - Simple *heuristic* default auto-tuning (Nvox_est) provided (manual choice usually better) -# - This temporarily uses disk space for multiprocessing, will warn if not enough space -# - In single precision or with FCC, Disk usage relative to simulation size goes up -# -# About voxelisation: -# - despite the use of term 'voxelizer' this is not necessarily a solid or surface voxelizer -# - this computes a lower-level adjacency graph for points on a grid (a grid of voxels) -# - this allows for non-watertight scenes -# - expects sides of materials properly oriented, boundary conditions depend on this (also FDTD performance) -# - has surface-area corrections to mitigate staircasing errors -# - points too near boundary (within some eps) are set not adjacent to neighbours -# - this is meant for interpretation of FDTD grid as FVTD mesh of voxels/cells -# - this exports data just for boundary nodes (anything with non-adjacency to a neighbour) -# -############################################################################## +# SPDX-License-Identifier: MIT + +"""A scene voxelizer for FDTD + +Sets up primary input for FDTD (adjacencies, materials, boundary nodes) +CPU-based and uses multiprocessing. +Tested up to ~10^10 points. Meant for single-node use (no MPI) + +Notes: + - Performance will depend on geometry, grid spacing, voxel size and # processes + - Simple *heuristic* default auto-tuning (Nvox_est) provided (manual choice usually better) + - This temporarily uses disk space for multiprocessing, will warn if not enough space + - In single precision or with FCC, Disk usage relative to simulation size goes up + +About voxelisation: + - despite the use of term 'voxelizer' this is not necessarily a solid or surface voxelizer + - this computes a lower-level adjacency graph for points on a grid (a grid of voxels) + - this allows for non-watertight scenes + - expects sides of materials properly oriented, boundary conditions depend on this (also FDTD performance) + - has surface-area corrections to mitigate staircasing errors + - points too near boundary (within some eps) are set not adjacent to neighbours + - this is meant for interpretation of FDTD grid as FVTD mesh of voxels/cells + - this exports data just for boundary nodes (anything with non-adjacency to a neighbour) +""" + +from pathlib import Path import numpy as np from numpy import array as npa +import numba as nb +import h5py +import psutil + +import multiprocessing as mp +from multiprocessing import shared_memory +from tqdm import tqdm + from pffdtd.common.room_geo import RoomGeo from pffdtd.common.timerdict import TimerDict from pffdtd.common.tri_ray_intersection import tri_ray_intersection_vec -from pffdtd.common.tri_box_intersection import tri_box_intersection_vec from pffdtd.voxelizer.cart_grid import CartGrid from pffdtd.voxelizer.vox_grid import VoxGrid from pffdtd.common.myfuncs import clear_dat_folder from pffdtd.common.myfuncs import yes_or_no,ind2sub3d from pffdtd.common.myfuncs import dotv from pffdtd.common.myfuncs import get_default_nprocs -from pathlib import Path -import numba as nb -import sys -import h5py - -import psutil - -import multiprocessing as mp -from multiprocessing import shared_memory -from tqdm import tqdm -from memory_profiler import profile as memory_profile F_EPS = np.finfo(np.float64).eps R_EPS = 1e-6 #relative eps (to grid spacing) for near hits