diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9e4d58a3..3edc10ff 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -75,7 +75,7 @@ jobs: - name: Install dependencies run: | brew update && (brew list cmake || brew install cmake) - brew install openblas eigen libomp lapack + brew install openblas eigen@3 libomp lapack sudo ln -sf /opt/homebrew/bin/gfortran-14 /opt/homebrew/bin/gfortran - name: Export Environment Variables @@ -92,13 +92,14 @@ jobs: export CPPFLAGS+=" -I/opt/homebrew/opt/libomp/include" export LDFLAGS+=" -L/opt/homebrew/opt/lapack/lib" export CPPFLAGS+=" -I/opt/homebrew/opt/lapack/include" + export CPPFLAGS+=" -I/opt/homebrew/opt/eigen@3/include" export PKG_CONFIG_PATH+=" /opt/homebrew/opt/lapack/lib/pkgconfig" - name: Create build directory run: mkdir -p build - name: Run CMake - run: cmake -DCMAKE_Fortran_COMPILER="/opt/homebrew/bin/gfortran-14" -S . -B build + run: cmake -DEigen3_DIR=/opt/homebrew/Cellar/eigen@3/3.4.1/share/eigen3/cmake -DCMAKE_Fortran_COMPILER="/opt/homebrew/bin/gfortran-14" -S . -B build - name: Build MOLE library run: cmake --build build diff --git a/CMakeLists.txt b/CMakeLists.txt index 06ab105c..50b6bb30 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,23 +14,27 @@ if (CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang") set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -L/usr/local/opt/libomp/lib -L/opt/homebrew/opt/libomp/lib -lomp") message(STATUS "Using AppleClang-specific flags.") include_directories("/usr/local/opt/libomp/include" "/opt/homebrew/opt/libomp/include") + message(STATUS "Adding Eigen3 include directory for AppleClang.") + set(EIGEN3_INCLUDE_DIR "/opt/homebrew/include/eigen3") + include_directories(${EIGEN3_INCLUDE_DIR}) + elseif (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM") - set(CMAKE_CXX_FLAGS "-O3 -qopenmp -DARMA_DONT_USE_WRAPPER -DARMA_USE_SUPERLU -diag-disable=10430") - set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") - message(STATUS "Using non-Clang compiler flags.") - # Get MKLROOT from environment or fallback to default - if(DEFINED ENV{MKLROOT}) - set(MKLROOT $ENV{MKLROOT}) - else() - set(MKLROOT "/opt/intel/oneapi/mkl/latest") - endif() - - # Automatically set Armadillo to link against MKL - set(BLAS_LIBRARIES "${MKLROOT}/lib/intel64/libmkl_rt.so" CACHE STRING "BLAS library path for MKL") - set(LAPACK_LIBRARIES "${MKLROOT}/lib/intel64/libmkl_rt.so" CACHE STRING "LAPACK library path for MKL") - set(ARMA_USE_WRAPPER OFF CACHE BOOL "Disable Armadillo wrapper to directly use MKL") - - message(STATUS "Using MKL from: ${MKLROOT}") + set(CMAKE_CXX_FLAGS "-O3 -qopenmp -DARMA_DONT_USE_WRAPPER -DARMA_USE_SUPERLU -diag-disable=10430") + set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") + message(STATUS "Using non-Clang compiler flags.") + # Get MKLROOT from environment or fallback to default + if(DEFINED ENV{MKLROOT}) + set(MKLROOT $ENV{MKLROOT}) + else() + set(MKLROOT "/opt/intel/oneapi/mkl/latest") + endif() + # Automatically set Armadillo to link against MKL + set(BLAS_LIBRARIES "${MKLROOT}/lib/intel64/libmkl_rt.so" CACHE STRING "BLAS library path for MKL") + set(LAPACK_LIBRARIES "${MKLROOT}/lib/intel64/libmkl_rt.so" CACHE STRING "LAPACK library path for MKL") + set(ARMA_USE_WRAPPER OFF CACHE BOOL "Disable Armadillo wrapper to directly use MKL") + + message(STATUS "Using MKL from: ${MKLROOT}") + else() set(CMAKE_CXX_FLAGS "-O3 -fopenmp -DARMA_DONT_USE_WRAPPER -DARMA_USE_SUPERLU") set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") diff --git a/doc/sphinx/Makefile b/doc/sphinx/Makefile index 76d9c969..d00af8b6 100644 --- a/doc/sphinx/Makefile +++ b/doc/sphinx/Makefile @@ -129,15 +129,24 @@ $(SPHINX_TARGETS): .PHONY: copy-images copy-images: + @echo "[IMAGE_COPY_DEBUG] ===== Starting image copy process =====" + @echo "[IMAGE_COPY_DEBUG] Copying images from ../assets/img/ to source/_images/" mkdir -p source/_images + @ls -la ../assets/img/*.png 2>/dev/null | grep -E "MOLE_pillars|MOLE_OSE_circles" || echo "[IMAGE_COPY_DEBUG] Source images not found" cp -r ../assets/img/* source/_images/ || true + @echo "[IMAGE_COPY_DEBUG] Images copied to source/_images/:" + @ls -la source/_images/*.png 2>/dev/null | grep -E "MOLE_pillars|MOLE_OSE_circles" || echo "[IMAGE_COPY_DEBUG] Images not in source/_images/" # Copy README images to the _images directory for proper reference + @echo "[IMAGE_COPY_DEBUG] Copying images to source/intros/doc/assets/img/" mkdir -p source/intros/doc/assets/img cp -r ../assets/img/* source/intros/doc/assets/img/ || true + @echo "[IMAGE_COPY_DEBUG] Images copied to source/intros/doc/assets/img/:" + @ls -la source/intros/doc/assets/img/*.png 2>/dev/null | grep -E "MOLE_pillars|MOLE_OSE_circles" || echo "[IMAGE_COPY_DEBUG] Images not in source/intros/doc/assets/img/" # Create figure directories for SVG files mkdir -p source/api/examples/md/figures mkdir -p source/api/examples-m/md/figures mkdir -p source/math_functions/figures + @echo "[IMAGE_COPY_DEBUG] ===== Image copy process complete =====" .PHONY: clean clean: diff --git a/doc/sphinx/source/_ext/image_path_logger.py b/doc/sphinx/source/_ext/image_path_logger.py new file mode 100644 index 00000000..d203e1f4 --- /dev/null +++ b/doc/sphinx/source/_ext/image_path_logger.py @@ -0,0 +1,147 @@ +""" +Sphinx extension to log image path resolution for debugging image inclusion issues. + +This extension hooks into Sphinx's document reading and image processing to log: +1. When documents with includes are processed +2. Image references found in markdown files +3. Path resolution attempts +4. Final image paths used by Sphinx +""" +import os +import logging +from pathlib import Path +from docutils import nodes +from sphinx.util import logging as sphinx_logging + +logger = sphinx_logging.getLogger(__name__) + +def log_document_read(app, doctree): + """Log when a document is read and what images it contains.""" + docname = app.env.docname + + # Get the source file path + source_file = app.env.doc2path(docname, base=None) + + # Check if this is an OSE organization wrapper document + if 'ose_organization' in docname: + logger.info("=" * 80) + logger.info(f"[IMAGE_DEBUG] Processing document: {docname}") + logger.info(f"[IMAGE_DEBUG] Source file: {source_file}") + logger.info(f"[IMAGE_DEBUG] Source directory: {os.path.dirname(source_file)}") + + # Check for image nodes + for node in doctree.traverse(nodes.image): + uri = node.get('uri', 'NO_URI') + logger.info(f"[IMAGE_DEBUG] Found image node with URI: {uri}") + + # Try to resolve the full path + if hasattr(app.env, 'relfn2path'): + try: + rel_fn, abs_fn = app.env.relfn2path(uri, docname) + logger.info(f"[IMAGE_DEBUG] Resolved relative path: {rel_fn}") + logger.info(f"[IMAGE_DEBUG] Resolved absolute path: {abs_fn}") + logger.info(f"[IMAGE_DEBUG] File exists: {os.path.exists(abs_fn)}") + except Exception as e: + logger.warning(f"[IMAGE_DEBUG] Failed to resolve path: {e}") + + logger.info("=" * 80) + + +def log_missing_reference(app, env, node, contnode): + """Log missing reference attempts (including images).""" + if isinstance(node, nodes.image): + uri = node.get('uri', 'UNKNOWN') + logger.warning(f"[IMAGE_DEBUG] Missing image reference: {uri}") + logger.warning(f"[IMAGE_DEBUG] Current docname: {env.docname}") + logger.warning(f"[IMAGE_DEBUG] Current doc path: {env.doc2path(env.docname)}") + return None + + +def log_image_copying(app, exception): + """Log image copying at the end of the build.""" + if exception is not None: + return + + logger.info("=" * 80) + logger.info("[IMAGE_DEBUG] Build finished - checking copied images") + + # Check what images ended up in _images + build_images_dir = Path(app.outdir) / "_images" + if build_images_dir.exists(): + logger.info(f"[IMAGE_DEBUG] Images directory: {build_images_dir}") + logger.info(f"[IMAGE_DEBUG] Images in _images directory:") + for img in sorted(build_images_dir.glob("*.png")): + logger.info(f"[IMAGE_DEBUG] - {img.name} ({img.stat().st_size} bytes)") + else: + logger.warning(f"[IMAGE_DEBUG] Images directory does not exist: {build_images_dir}") + + logger.info("=" * 80) + + +def trace_myst_include_processing(app, docname, source): + """ + Hook into source-read event to log MyST include processing. + This runs BEFORE MyST parses the document. + """ + if 'ose_organization' in docname: + logger.info("=" * 80) + logger.info(f"[IMAGE_DEBUG] Source-read event for: {docname}") + + # Check if this file has an include directive + if '{include}' in source[0]: + logger.info(f"[IMAGE_DEBUG] Document contains include directive") + + # Extract the included file path + import re + include_match = re.search(r'\{include\}\s+([^\s\n]+)', source[0]) + if include_match: + included_path = include_match.group(1) + logger.info(f"[IMAGE_DEBUG] Including file: {included_path}") + + # Resolve the full path + source_dir = Path(app.env.doc2path(docname, base=None)).parent + full_included_path = (source_dir / included_path).resolve() + logger.info(f"[IMAGE_DEBUG] Source document dir: {source_dir}") + logger.info(f"[IMAGE_DEBUG] Resolved include path: {full_included_path}") + logger.info(f"[IMAGE_DEBUG] Include file exists: {full_included_path.exists()}") + + if full_included_path.exists(): + # Read the included file and look for image references + with open(full_included_path, 'r') as f: + included_content = f.read() + + # Find markdown image references + image_matches = re.findall(r'!\[([^\]]*)\]\(([^)]+)\)', included_content) + if image_matches: + logger.info(f"[IMAGE_DEBUG] Found {len(image_matches)} image(s) in included file:") + for alt_text, img_path in image_matches: + logger.info(f"[IMAGE_DEBUG] - Alt: '{alt_text}', Path: '{img_path}'") + + # Check if image path is relative to included file or source file + # Path relative to included file's location + rel_to_included = (full_included_path.parent / img_path).resolve() + logger.info(f"[IMAGE_DEBUG] Relative to included file: {rel_to_included}") + logger.info(f"[IMAGE_DEBUG] Exists (rel to included): {rel_to_included.exists()}") + + # Path relative to source document + rel_to_source = (source_dir / img_path).resolve() + logger.info(f"[IMAGE_DEBUG] Relative to source doc: {rel_to_source}") + logger.info(f"[IMAGE_DEBUG] Exists (rel to source): {rel_to_source.exists()}") + + logger.info("=" * 80) + + +def setup(app): + """Register the extension.""" + # Connect to various events + app.connect('doctree-read', log_document_read) + app.connect('missing-reference', log_missing_reference) + app.connect('build-finished', log_image_copying) + app.connect('source-read', trace_myst_include_processing) + + return { + 'version': '0.1', + 'parallel_read_safe': True, + 'parallel_write_safe': True, + } + diff --git a/doc/sphinx/source/conf.py b/doc/sphinx/source/conf.py index 6e8e6e25..c05552d4 100644 --- a/doc/sphinx/source/conf.py +++ b/doc/sphinx/source/conf.py @@ -371,10 +371,58 @@ def fix_math_environments(app, docname, source): source[0] = src +def copy_governance_images(app, config): + """ + Copy governance and organization images before Sphinx build starts. + + This ensures images referenced in OSE_ORGANIZATION.md are available + when Sphinx processes the document, regardless of whether the build + is run via Makefile or directly via sphinx-build (e.g., on ReadTheDocs). + + The images are copied to two locations: + 1. source/_images/ - Standard Sphinx image directory + 2. source/intros/doc/assets/img/ - Where MyST resolves relative paths + from the including file (ose_organization_wrapper.md) + + See: GitHub Issue #222 + """ + import shutil + from pathlib import Path + + # Determine paths relative to conf.py location + conf_dir = Path(app.confdir) + + # Source: doc/assets/img/ (relative to repo root) + img_source = conf_dir.parent.parent.parent / "doc" / "assets" / "img" + + # Destinations + img_dest_1 = conf_dir / "_images" + img_dest_2 = conf_dir / "intros" / "doc" / "assets" / "img" + + if not img_source.exists(): + print(f"Warning: Image source directory not found: {img_source}") + return + + # Copy to both locations + for dest in [img_dest_1, img_dest_2]: + dest.mkdir(parents=True, exist_ok=True) + + # Copy all image files + for pattern in ["*.png", "*.jpg", "*.jpeg", "*.gif", "*.svg"]: + for img in img_source.glob(pattern): + dest_file = dest / img.name + try: + shutil.copy2(img, dest_file) + except Exception as e: + print(f"Warning: Could not copy {img.name}: {e}") + def setup(app): """Setup function for Sphinx extension.""" app.add_js_file('mathconf.js') + # Copy governance images before build starts (Fix for Issue #222) + app.connect('config-inited', copy_governance_images) + # Add capability to replace problematic math environments app.connect('source-read', fix_math_environments) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Homogeneous-Dirichlet.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Homogeneous-Dirichlet.md index 97fd4a69..dc28f08a 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Homogeneous-Dirichlet.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Homogeneous-Dirichlet.md @@ -41,10 +41,11 @@ $$ This example is implemented in: - [MATLAB/ OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/elliptic1DHomogeneousDirichlet.cpp) Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Non-Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Dirichlet, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightRobin.m) - [Left Neumann, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) - [Left Neumann, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Dirichlet-Right-Robin.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Dirichlet-Right-Robin.md index 53a8a4c6..8c8a9952 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Dirichlet-Right-Robin.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Dirichlet-Right-Robin.md @@ -51,7 +51,7 @@ This example is implemented in: Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) - [Non-Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Neumann, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) - [Left Neumann, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightRobin.m) - [Left Robin, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftRobinRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Dirichlet-Right-Neumann.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Dirichlet.md similarity index 85% rename from doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Dirichlet-Right-Neumann.md rename to doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Dirichlet.md index 6c95118b..04d78de4 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Dirichlet-Right-Neumann.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Dirichlet.md @@ -1,6 +1,6 @@ -### Elliptic1D Left Dirichlet and Right Neumann Boundary Conditions +### Elliptic1D Left Neumann and Right Dirichlet Boundary Conditions -Solves the 1D Poisson equation with left Dirichlet and right Neumann boundary conditions. +Solves the 1D Poisson equation with left Neumann and right Dirichlet boundary conditions. $$ -\nabla^2 u(x) = 1 @@ -43,7 +43,8 @@ $$ --- This example is implemented in: -- [MATLAB/ OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [MATLAB/ OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/elliptic1DLeftNeumannRightDirichlet.cpp) Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Neumann.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Neumann.md index 1aaa61b3..b4c5dabd 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Neumann.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Neumann.md @@ -47,11 +47,12 @@ $$ This example is implemented in: - [MATLAB/ OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/elliptic1DLeftNeumannRightNeumann.cpp) Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) - [Non-Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Dirichlet, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightRobin.m) - [Left Neumann, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightRobin.m) - [Left Robin, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftRobinRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Robin.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Robin.md index 414559a7..944286fb 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Robin.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Neumann-Right-Robin.md @@ -54,7 +54,7 @@ This example is implemented in: Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) - [Non-Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Dirichlet, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightRobin.m) - [Left Neumann, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) - [Left Robin, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftRobinRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Robin-Right-Robin.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Robin-Right-Robin.md index be844749..a4290e8d 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Robin-Right-Robin.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Left-Robin-Right-Robin.md @@ -53,7 +53,7 @@ This example is implemented in: Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) - [Non-Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Dirichlet, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightRobin.m) - [Left Neumann, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) - [Left Neumann, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Non-Homogeneous-Dirichlet.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Non-Homogeneous-Dirichlet.md index 911a380e..def7b829 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Non-Homogeneous-Dirichlet.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Non-Homogeneous-Dirichlet.md @@ -41,10 +41,11 @@ $$ This example is implemented in: - [MATLAB/ OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/elliptic1DNonHomogeneousDirichlet.cpp) Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Dirichlet, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightRobin.m) - [Left Neumann, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) - [Left Neumann, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Non-Periodic-BC.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Non-Periodic-BC.md index 555ffd28..6f9567b1 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Non-Periodic-BC.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Non-Periodic-BC.md @@ -48,7 +48,7 @@ This example is implemented in: Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) - [Non-Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Dirichlet, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightRobin.m) - [Left Neumann, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) - [Left Neumann, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Periodic-BC.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Periodic-BC.md index a672363f..267d4601 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Periodic-BC.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-Periodic-BC.md @@ -37,7 +37,7 @@ This example is implemented in: Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) - [Non-Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Dirichlet, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightRobin.m) - [Left Neumann, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) - [Left Neumann, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-add-Scalar-BC.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-add-Scalar-BC.md index d1bba6cb..208956f3 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-add-Scalar-BC.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D-add-Scalar-BC.md @@ -46,7 +46,7 @@ This example is implemented in: Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Non-Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Dirichlet, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightRobin.m) - [Left Neumann, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) - [Left Neumann, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D.md b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D.md index 82c40a2c..a524b0f2 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D.md +++ b/doc/sphinx/source/examples/Elliptic/1D/Elliptic1D.md @@ -33,7 +33,7 @@ This example is implemented in: Additional MATLAB/ OCTAVE variants of this example with different boundary conditions: - [Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DHomogeneousDirichlet.m) - [Non-Homogeneous Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DNonHomogeneousDirichlet.m) -- [Left Dirichlet, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m) +- [Left Neumann, Right Dirichlet](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m) - [Left Dirichlet, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftDirichletRightRobin.m) - [Left Neumann, Right Neumann](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightNeumann.m) - [Left Neumann, Right Robin](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/elliptic1DLeftNeumannRightRobin.m) diff --git a/doc/sphinx/source/examples/Elliptic/1D/index.md b/doc/sphinx/source/examples/Elliptic/1D/index.md index ee1af391..94e76e15 100644 --- a/doc/sphinx/source/examples/Elliptic/1D/index.md +++ b/doc/sphinx/source/examples/Elliptic/1D/index.md @@ -8,7 +8,7 @@ These examples demonstrate the solution of elliptic equations in one dimension, Elliptic1D Elliptic1D-add-Scalar-BC Elliptic1D-Homogeneous-Dirichlet -Elliptic1D-Left-Dirichlet-Right-Neumann +Elliptic1D-Left-Neumann-Right-Dirichlet Elliptic1D-Left-Dirichlet-Right-Robin Elliptic1D-Left-Neumann-Right-Neumann Elliptic1D-Left-Neumann-Right-Robin diff --git a/doc/sphinx/source/examples/Sturm-Liouville/Bessel.md b/doc/sphinx/source/examples/Sturm-Liouville/Bessel.md new file mode 100644 index 00000000..02452029 --- /dev/null +++ b/doc/sphinx/source/examples/Sturm-Liouville/Bessel.md @@ -0,0 +1,54 @@ +# Bessel Sturm-Liouville Problem + +This example solves the Bessel differential equation, which is a classic Sturm-Liouville problem: + +$$ x^2 u'' + x u' + (x^2 - \nu^2) u = 0, \quad 0 < x 1 $$ + +with Dirichlet boundary conditions: +$$ u(0) = 0, \quad u(1) = J_\nu(1) $$ + +The exact solution to this problem is the Bessel function of the first kind of order $\nu$, denoted as $J_\nu(x)$. For $\nu = 3$, the solution is $ J_3(x) = \frac{1}{\pi}\int^{\pi}_0 \cos(3\tau) - x\sin(\tau)\,d\tau $. + +## Mathematical Background + +Bessel's differential equation is a special case of the Sturm-Liouville problem, which has the general form: + +$$\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) + q(x)u + \lambda r(x)u = 0$$ + +For Bessel's equation, we have: + +- $p(x) = x^2$ +- $q(x) = x$ +- $r(x) = \frac{1}{x}$ +- $\lambda = -\nu^2$ + +## Discretization + +The equation is discretized using mimetic finite difference operators. The spatial derivative operators are constructed with a specified order of accuracy $k$. + +The discrete system is: + +$$A u = b$$ + +where: + +- $A = x^2 L + x I_{FC} G + (x^2 - \nu^2) I$ +- $L$ is the mimetic Laplacian +- $G$ is the mimetic gradient +- $I_{FC}$ is the interpolation operator from faces to centers +- $I$ is the identity matrix + +Boundary conditions are applied using `RobinBC`. + +--- + +This example is implemented in: + +- [MATLAB/OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/sturmLiouvilleBessel.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/sturmLiouvilleBessel.cpp) + +## Results + +The numerical solution closely matches the exact solution, which is the Bessel function of the first kind of order $\nu$: + +$J_\nu(x) = \frac{1}{\pi}\int^{\pi}_0 \cos(\nu\tau) - x\sin(\tau)\,d\tau$. diff --git a/doc/sphinx/source/examples/Sturm-Liouville/Chebyshev.md b/doc/sphinx/source/examples/Sturm-Liouville/Chebyshev.md index 226f6354..10edecb3 100644 --- a/doc/sphinx/source/examples/Sturm-Liouville/Chebyshev.md +++ b/doc/sphinx/source/examples/Sturm-Liouville/Chebyshev.md @@ -16,6 +16,7 @@ Chebyshev's differential equation is a special case of the Sturm-Liouville probl $$\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) + q(x)u + \lambda r(x)u = 0$$ For Chebyshev's equation, we have: + - $p(x) = 1-x^2$ - $q(x) = 0$ - $r(x) = 1$ @@ -30,23 +31,28 @@ The discrete system is: $$A u = b$$ where: -- $A = (1-x^2) L - x I G + n^2 I$ + +- $A = (1-x^2) L - x I_{FC} G + n^2 I$ - $L$ is the mimetic Laplacian - $G$ is the mimetic gradient -- $I$ is the interpolation operator from faces to centers +- $I_{FC}$ is the interpolation operator from faces to centers +- $I$ is the identity matrix Boundary conditions are applied using the `addScalarBC1D` function. --- This example is implemented in: + - [MATLAB/OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/sturmLiouvilleChebyshev.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/sturmLiouvilleChebyshev.cpp) ## Results -The numerical solution closely matches the exact solution, which is the Chebyshev polynomial $T_2(x) = 2x^2 - 1$. +The numerical solution closely matches the exact solution, which is the Chebyshev polynomial $T_2(x) = 2x^2 - 1$. Chebyshev polynomials are important in numerical analysis and approximation theory because they: + 1. Minimize the maximum error in polynomial approximation 2. Have roots that are optimal interpolation points (Chebyshev nodes) -3. Are closely related to the Fourier cosine series +3. Are closely related to the Fourier cosine series diff --git a/doc/sphinx/source/examples/Sturm-Liouville/Helmholtz.md b/doc/sphinx/source/examples/Sturm-Liouville/Helmholtz.md new file mode 100644 index 00000000..b6778ce1 --- /dev/null +++ b/doc/sphinx/source/examples/Sturm-Liouville/Helmholtz.md @@ -0,0 +1,56 @@ +# Helmholtz Sturm-Liouville Problem + +This example solves the Helmholtz differential equation, which is a classic Sturm-Liouville problem: + +$$ u'' + u = 0, \quad 0 < x < 3 $$ +or +$$ u'' + \mu^2 u = 0, \quad 0 < x < 1 $$ + +with boundary conditions: +$$ u(0) = 0, \quad u(3) = \sin(3) $$ +or +$$ u'(0) = 0, \quad u(1) + u'(1) = \cos(\mu) - \mu \sin(\mu) $$ + +The exact solution to this problem is $\sin(x)$ or $\cos(\mu x)$. + +## Mathematical Background + +Helmholtz's differential equation is a special case of the Sturm-Liouville problem, which has the general form: + +$$\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) + q(x)u + \lambda r(x)u = 0$$ + +For Helmholtz's equation, we have: + +- $p(x) = 1$ +- $q(x) = 0$ +- $r(x) = 1$ +- $\lambda = \mu^2$ + +## Discretization + +The equation is discretized using mimetic finite difference operators. The spatial derivative operators are constructed with a specified order of accuracy $k$. + +The discrete system is: + +$$A u = b$$ + +where: + +- $A = L + I$ +- $L$ is the mimetic Laplacian +- $I$ is the identity matrix + +Boundary conditions are applied using `RobinBC` or `MixedBC`. + +--- + +These examples are implemented in: + +- [MATLAB/OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/sturmLiouvilleHelmholtzDirichletDirichlet.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/sturmLiouvilleHelmholtzDirichletDirichlet.cpp) +- [MATLAB/OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/sturmLiouvilleHelmholtzDirichletRobin.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/sturmLiouvilleHelmholtzDirichletRobin.cpp) + +## Results + +The numerical solutions closely match the exact solutions, which are $\sin(x)$ or $\cos(\mu x)$. diff --git a/doc/sphinx/source/examples/Sturm-Liouville/Hermite.md b/doc/sphinx/source/examples/Sturm-Liouville/Hermite.md new file mode 100644 index 00000000..93162ecd --- /dev/null +++ b/doc/sphinx/source/examples/Sturm-Liouville/Hermite.md @@ -0,0 +1,52 @@ +# Hermite Sturm-Liouville Problem + +This example solves the Hermite differential equation, which is a classic Sturm-Liouville problem: + +$$ u'' - 2 u' + 2 m u = 0, \quad -1 < x < 1 $$ + +with Dirichlet boundary conditions: +$$ u(-1) = H_m(-1), \quad u(1) = H_m(1) $$ + +The exact solution to this problem is the Hermite polynomial of degree $m$, denoted as $H_m(x)$. For $m = 4$, the solution is $H_4(x) = e^{x^2} \frac{d^4}{dx^4}e^{-x^2}$. + +## Mathematical Background + +Hermite's differential equation is a special case of the Sturm-Liouville problem, which has the general form: + +$$\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) + q(x)u + \lambda r(x)u = 0$$ + +For Hermite's equation, we have: + +- $p(x) = e^{-x^2}$ +- $q(x) = 0$ +- $r(x) = e^{-x^2}$ +- $\lambda = m^2$ + +## Discretization + +The equation is discretized using mimetic finite difference operators. The spatial derivative operatores are constructed with a specified order of accuracy $k$. + +The discrete system is: + +$$A u = b$$ + +where + +- $A = L - 2 x I_{FC} G + 2 m I$ +- $L$ is the mimetic Laplacian +- $G$ is the mimetic gradient +- $I_{FC}$ is the interpolation operator from faces to centers +- $I$ is the identity matrix + +Boundary conditions are applied using `RobinBC`. + +--- + +This example is implemented in: + +- [MATLAB/OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/sturmLiouvilleHermite.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/sturmLiouvilleHermite.cpp) + +## Results + +The numerical solution closely matches the exact solution, which is the Hermite polynomial $H_4(x) = e^{x^2} \frac{d^4}{dx^4}e^{-x^2}$. diff --git a/doc/sphinx/source/examples/Sturm-Liouville/Laguerre.md b/doc/sphinx/source/examples/Sturm-Liouville/Laguerre.md new file mode 100644 index 00000000..26674a25 --- /dev/null +++ b/doc/sphinx/source/examples/Sturm-Liouville/Laguerre.md @@ -0,0 +1,52 @@ +# Laguerre Sturm-Liouville Problem + +This example solves the Laguerre differential equation, which is a classic Sturm-Liouville problem: + +$$x u'' + (1 - x) u' + n u = 0, \quad 0 < x < 2$$ + +with Dirichlet boundary conditions: +$$ u(0) = L_n(2), \quad u(2) = L_n(2) $$ + +The exact solution to this problem is the Laguerre polynomial of degree $n$, denoted as $L_n(x)$. For $n = 4$, the solution is $L_4(x) = \frac{e^x}{4!} \frac{d^4}{dx^4} (e^{-x}x^4) $. + +## Mathematical Background + +Laguerre's differential equation is a special case of the Sturm-Liouville problem, which has the general form: + +$$\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) + q(x)u + \lambda r(x)u = 0$$ + +For Laguerre's equation, we have: + +- $p(x) = xe^{-x}$ +- $q(x) = 0$ +- $r(x) = e^{-x}$ +- $\lambda = n$ + +## Discretization + +The equation is discretized using mimetic finite difference operators. The spatial derivative operators are constructed with a specified order of accuracy $k$. + +The discrete system is: + +$$A u = b$$ + +where: + +- $A = x L + (1 - x) I_{FC} G + n I$ +- $L$ is the mimetic Laplacian +- $G$ is the mimetic gradient +- $I_{FC}$ is the interpolation operator from faces to centers +- $I$ is the identity matrix + +Boundary conditions are applied using `RobinBC`. + +--- + +This example is implemented in: + +- [MATLAB/OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/sturmLiouvilleLaguerre.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/sturmLiouvilleLaguerre.cpp) + +## Results + +The numerical solution closely matches the exact solution, which is the Laguerre polynomial $L_4(x) = \frac{e^x}{4!} \frac{d^4}{dx^4} (e^{-x}x^4) $. diff --git a/doc/sphinx/source/examples/Sturm-Liouville/Legendre.md b/doc/sphinx/source/examples/Sturm-Liouville/Legendre.md new file mode 100644 index 00000000..f697b5fa --- /dev/null +++ b/doc/sphinx/source/examples/Sturm-Liouville/Legendre.md @@ -0,0 +1,52 @@ +# Legendre Sturm-Liouville Problem + +This example solves the Legendre differential equation, which is a classic Sturm-Liouville problem: + +$$ (1 - x^2) u'' - 2 x u' + n (n + 1) u = 0, \quad -1 < x < 1 $$ + +with Dirichlet boundary conditions: +$$ u(-1) = -1, \quad u(1) = 1 $$ + +The exact solutoin to this problem is the Legendre polynomial of degree $n$, denoted as $L_n(x)$. For $n = 4$, the solution is $L_4(x) = \frac{1}{8}(35x^4 -30x^2 + 3)$. + +## Mathematical Background + +Legendre's differential equation is a special case of the Sturm-Liouville problem, which has the general form: + +$$\frac{d}{dx}\left(p(x)\frac{du}{dx}\right) + q(x)u + \lambda r(x)u = 0$$ + +For Legendre's equation, we have: + +- $p(x) = 1-x^2$ +- $q(x) = 0$ +- $r(x) = 1$ +- $\lambda = n(n+1)$ + +## Discretization + +The equation is discretized using mimetic finite difference operators. The spatial derivative operators are constructed with a specified order of accuracy $k$. + +The discrete system is: + +$$A u = b$$ + +where: + +- $A = (1 - x^2) L - 2 x I_{FC} G + n(n+1) I$ +- $L$ is the mimetic Laplacian +- $G$ is the mimetic gradient +- $I_{FC}$ is the interpolation operator from faces to centers +- $I$ is the identity matrix + +Boundary conditions are applied using `RobinBC`. + +--- + +This example is implemented in: + +- [MATLAB/OCTAVE](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/sturmLiouvilleLegendre.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/sturmLiouvilleLegendre.cpp) + +## Results + +The numerical solution closely matches the exact solution, which is the Legendre polynomial $L_4(x) = \frac{1}{8}(35x^4 -30x^2 + 3)$. diff --git a/doc/sphinx/source/examples/Sturm-Liouville/index.md b/doc/sphinx/source/examples/Sturm-Liouville/index.md index b5db197a..9545b9b4 100644 --- a/doc/sphinx/source/examples/Sturm-Liouville/index.md +++ b/doc/sphinx/source/examples/Sturm-Liouville/index.md @@ -6,5 +6,10 @@ Sturm-Liouville theory deals with a class of second-order linear differential eq :maxdepth: 1 :caption: Contents +Bessel Chebyshev -``` +Helmholtz +Hermite +Laguerre +Legendre +``` diff --git a/doc/sphinx/source/examples/Time-Integrators/backward_euler.md b/doc/sphinx/source/examples/Time-Integrators/backward_euler.md new file mode 100644 index 00000000..7e383451 --- /dev/null +++ b/doc/sphinx/source/examples/Time-Integrators/backward_euler.md @@ -0,0 +1,43 @@ +### Backward Euler + +This example solves the first-order ordinary differential equation using the backward Euler method. + +$$ +\frac{dy}{dt} = \sin^2(t) \cdot y +$$ + +with initial condition $y(0) = 0$ over the time interval $[0,5]$. + +#### Mathematical Background + +The backward Euler is an implicit method for solving first-order differential equations where we are evaluating the function at the next point in time $y_{i+1}$. + +$$ +y_{i+1} = y_i + h \cdot f(t_{i+1}, y_{i+1}) +$$ + +where: +- $h$ is the step-size +- $t_i$ is the current time +- $y_i$ is the solution at time $t_i$ +- $y_{i+1}$ is the solution at the next time step $t_{i+1}$ +- $f(t_{i+1}, y_{i+1})$ is the function at time step $t_{i+1}$ and $y_{i+1}$ + +#### Implementation Details + +Note that $f(t_{i+1}, y_{i+1})$ is not known due to $y_{i+1}$ being on both the left and right hand side. Therefore, we can solve for $y_{i+1}$ by finding the root: +$$ +y_{i+1} - y_i + h \cdot f(t_{i+1}, y_{i+1}) = 0 +$$ + +In the C++ example, we used the fixed-point iteration method for rootfinding due to it's simpler nature. + +#### Results + +![Backward Euler gnuplot](../../_images/backward_euler.png) + +--- + +This example is implemented in: +- [MATLAB/Octave](https://github.com/csrc-sdsu/mole/blob/main/examples/matlab_octave/backward_euler.m) +- [C++](https://github.com/csrc-sdsu/mole/blob/main/examples/cpp/backward_euler.cpp) diff --git a/doc/sphinx/source/examples/Tutorials/Two_Way_Wave_Equation/Two_Way_Wave_Equation.md b/doc/sphinx/source/examples/Tutorials/Two_Way_Wave_Equation/Two_Way_Wave_Equation.md new file mode 100644 index 00000000..bf6868cf --- /dev/null +++ b/doc/sphinx/source/examples/Tutorials/Two_Way_Wave_Equation/Two_Way_Wave_Equation.md @@ -0,0 +1,243 @@ +# Mimetic Difference Operators Vs Finite Differences +### Solving the 1D Two-Way Wave Equation + +This tutorial demonstrates how to solve the **1D two-way wave equation** using both **standard finite differences (FD)** and **mimetic finite differences (MD)**. +We will use three main functions: + +1. [`finite_diff_two_way_wave_eq`](#finite-diff-two-way-wave-eq) +2. [`mimetic_diff_two_way_wave_eq`](#mimetic-diff-two-way-wave-eq) +3. [`comparison_two_way_wave_md_vs_fd`](#comparison-two-way-wave-md-vs-fd) + +The goal is to understand the **differences in accuracy, stability, and numerical properties** between traditional finite and mimetic difference operators. The two numerical functions have been created to be as similar as possible, so that only the changes necessary for mimetic difference operators are visible. + +While each code can be run independanlty, the easiest way to compare the two methods is with `comparison_two_way_wave_md_vs_fd` + +--- + +## Overview + +The **1D two-way wave equation** is given by: + +$$ +\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}, \quad x \in [-2, 2] \qquad \qquad \qquad (1) +$$ + +Both solvers use a **centered-in-time, centered-in-space (CTCS)** scheme, which is **second-order accurate in both time and space**. + +We assume no boundary effects (large domain), allowing a clean comparison of the spatial discretization methods. Here is a little bit more information about the two seperate functions. + +--- +## `finite diff two way wave eq` + +**Purpose:** This function solves the the wave equation using **standard finite differences (FD)**. + +**Time-stepping scheme (CTCS):** +If we discretize the equation with standard centered second order finite differences, Equation(1) turns into + +$$ +\frac{U^{k-1}_i - 2U^{k}_i + U^{k+1}_i}{\Delta t^2} = c^2\frac{U_{i+1}^k - 2U_{i}^k + U_{i+1}^k}{\Delta x^2} +$$ + +Rearranging the equation for the unknown value $U^{k+1}_i$ we get: + +$$ +U^{k+1}_i = 2U^k_i - U^{k-1}_i + \frac{c^2\Delta t^2}{\Delta x^2}\Big(U^k_{i-1} - 2U^k_{i} + U^k_{i+1} \Big) +$$ + +Which can be expressed as a matrix $D_{fd}$ times the vector ${U^k}$ + +$$ +U^{k+1} = 2U^k - U^{k-1} + D_{fd}U^k +$$ + + +The matrix is written out in the code using sparse matrix operations and is equivalent to: + +$$ +D_{fd} = \frac{c^2\Delta t^2}{\Delta x^2} +\begin{bmatrix} +-2 & 1 & 0 & \cdots \\ +1 & -2 & 1 & \cdots \\ +0 & 1 & -2 & \cdots \\ +\vdots & \vdots & \vdots & \ddots +\end{bmatrix} +$$ + +The new values are updated each time step with values from the previous two time steps (k, k-1). + +$$ +U^{k+1} = 2U^k - U^{k-1} + D_{fd}U^k +$$ + +### Function Outputs +| Variable | Description | +|-----------|--------------| +| `U2_fd` | Final solution at the last time step | +| `error_fd` | Norm of error vs. analytic/reference solution | +| `walltime_fd` | Wall-clock time (seconds) | +| `flops_fd` | Estimated floating-point operation count | + +### Notes +We are using an explicit two-step scheme (leapfrog). The CFL condition must be satisfied: $ c \Delta t^2 / \Delta x^2 \le 1 $. If you change the stepping, or increase the number of cells, be sure to change this value. No boundary conditions are applied, the domain is large enough to avoid reflection effects. This is to test just the spatial scheme without any boundary considerations. + +--- +## `mimetic diff two way wave eq` + +**Purpose:** Solve the wave equation using **mimetic difference operators (MD)**. + +The mimetic Laplacian operator `L` replaces the finite-difference stencil, automatically enforcing discrete analogs of conservation and symmetry properties. + +**Time-stepping scheme:** +The mimetic operator directly takes the place fo the second derivative (Laplacian), therefore we only need to discretize the time derivative. Doing so to Equation (1) yields: + +$$ +\frac{U^{k-1}_i - 2U^{k}_i + U^{k+1}_i}{\Delta t^2} = c^2 L\,(U^k) +$$ + +Here, `L` is constructed from mimetic gradient (`G`) and divergence (`D`) operators such that: + +$$ +L = D \, G +$$ + +where + +$$ +\mathbf{D} = \frac{1}{\Delta x} +\begin{bmatrix} + -\tfrac{1}{2} & -\tfrac{1}{2} & 0 & \cdots & 0\\ + 1 & -1 & 0 & \cdots & 0\\ + 0 & 1 & -1 & \cdots & 0\\ + \vdots & & \ddots & \ddots & \vdots\\ + 0 & \cdots & 0 & \tfrac{1}{2} & -\tfrac{1}{2} +\end{bmatrix}. +$$ + +and + +$$ +\mathbf{G} = \frac{1}{\Delta x} +\begin{bmatrix} + -\tfrac{1}{2} & 1 & 0 & 0 & \cdots & 0 & 0\\[4pt] + -\tfrac{1}{2} & -1 & 1 & 0 & \cdots & 0 & 0\\[4pt] + 0 & 0 & -1 & 1 & \cdots & 0 & 0\\[4pt] + \vdots & & \ddots & \ddots & \ddots & & \vdots\\[4pt] + 0 & \cdots & 0 & -1 & 1 & 0 & 0\\[4pt] + 0 & \cdots & 0 & 0 & -1 & -1 & 1\\[4pt] + 0 & \cdots & 0 & 0 & 0 & -\tfrac{1}{2} & \tfrac{1}{2} +\end{bmatrix}. +$$ + +ensuring energy and flux consistency at the discrete level. The mimetic library Laplacian operator is just the composition $DG$ under the hood. + +Rearranging our equation for the unknown value $U^{k+1}$ leads us to: + +$$ +U^{k+1}_i = 2U^k_i - U^{k-1}_i + c^2\Delta t^2 L(U^k) +$$ + +To save computations in the code, `L` is premultiplied by $C^2\Delta t^2$ before the computation loop. + +### Outputs +| Variable | Description | +|-----------|--------------| +| `U2_md` | Final solution at the last time step | +| `error_md` | Norm of error vs. analytic/reference solution | +| `walltime_md` | Wall-clock time (seconds) | +| `flops_md` | Estimated floating-point operation count | + +### Notes +This uses the same CTCS time integration scheme as FD. The **order of accuracy** can be increased by adjusting the `k` parameter (e.g. `k=2`, `k=4`, etc.). Mimetic operators preserve **discrete conservation laws**, often improving numerical stability and physical fidelity. If you change the number of cells, be sure to change the time step, to conform to the CFL condition. + +--- +## `comparison two way wave md vs fd` + +**Purpose:** Compare **mimetic** and **finite difference** results for the same PDE setup. + +This script runs both solvers and measures differences in: +- Numerical accuracy +- Computational cost +- Wall-clock time +- Error convergence rate + +Assorted plots are generated showing differences between the two methods. + + +## Grid Comparison + +The key difference between FD and MD methods lies in the spatial grid. Here is a slightly more involved explanation to help the user understand the difference in grids. + +### Finite Difference Grid (Uniform) + +Each cell is the same width Δx between **A** and **B**, with points at cell boundaries (0,1,2,...N-1,N,N+1): +``` +A B + <-----dx-----> <-----dx-----> ... <-----dx-----> <-----dx-----> +|----cell 1----|----cell 2----| ... |----cell N-1--|----cell N----| +0 1 2 ... N-1 N N+1 +``` + +For example, from 0 to 1 with 5 cells has six values (x0,x1,x2,x3,x4,x5) each 0.2 away from each other: +``` +0.0 0.2 0.4 0.6 0.8 1.0 + o-----o-----o-----o-----o-----o + x0 x1 x2 x3 x4 x5 +``` +### Mimetic Difference Grid (Staggered) + +Mimetic operators use a staggered grid—half-step offsets at boundaries improve conservation and symmetry. Here, note that **cell 1** and **cell N** are not the sam size as the internal cells, and that there are now **N+2** points. +``` +A B + <--dx/2--> <-----dx-----> ... <-----dx-----> <--dx/2--> +|--cell 1--|----cell 2----| ... |----cell N-1--|--cell N--| +0 1 2 ... N N+1 N+2 +``` + +For the same interval (0–1) with 5 cells, a staggered grid will have smaller dt/s at the beginning and end, and also one mroe point (x0,x1,x2,x3,x4,x5,x6): +``` +0.0 0.1 0.3 0.5 0.7 0.9 1.0 + o-----o-----o-----o-----o-----o-----o + x0 x1 x2 x3 x4 x5 x6 +``` + +This staggered layout provides **better discrete analogs** of differential operators, leading to improved conservation and often reduced numerical dispersion. + +--- + +## Convergence and Accuracy + +Both FD and MD methods use a **second-order accurate** CTCS scheme in time. + +| Scheme | Spatial Order | Time Order | Conservation | Grid Type | +|---------|----------------|-------------|---------------|------------| +| Finite Difference | 2 | 2 | Approximate | Uniform | +| Mimetic Difference | 2 (or higher) | 2 | Exact (discrete) | Staggered | + +--- + +## Practical Notes + +- This setup isolates **spatial discretization effects**—no boundary reflections or external forcing. +- Increasing `N` refines the grid and reduces spatial error. +- For fair comparison, both solvers should use the same `dt`, `dx`, and runtime length. +- The mimetic operator’s flexibility allows easy order-of-accuracy experiments by simply changing `k` in `mimetic_diff_two_way_wave_eq.m`. + +--- + +## Summary + +| Feature | Finite Difference | Mimetic Difference | +|----------|------------------|--------------------| +| Grid | Uniform | Staggered | +| Operator | Explicit stencil | Discrete mimetic operator | +| Conservation | Approximate | Exact (discrete) | +| Flexibility | Fixed 2nd order | Adjustable order (k) | +| Accuracy | 2nd order | ≥ 2nd order | +| Ease of implementation | Simple | Slightly more complex | +| Best use case | Quick prototyping | Physically consistent PDE solvers | + +--- + +## References +- Castillo, J. E., *Mimetic Methods for Partial Differential Equations*. Cambridge University Press, 2008. +- LeVeque, R. J., *Finite Difference Methods for Ordinary and Partial Differential Equations*. SIAM, 2007. \ No newline at end of file diff --git a/doc/sphinx/source/examples/Tutorials/index.md b/doc/sphinx/source/examples/Tutorials/index.md new file mode 100644 index 00000000..5be70266 --- /dev/null +++ b/doc/sphinx/source/examples/Tutorials/index.md @@ -0,0 +1,9 @@ +# Tutorials + +Tutorials for using the mimetic operators including comparison with other standard ways of solving systems of PDEs. These examples are heavily commented, and show how standard workflows could be modified to use mimetic operators. + +```{toctree} +:maxdepth: 1 +:caption: Contents +Two_Way_Wave_Equation/Two_Way_Wave_Equation +``` diff --git a/doc/sphinx/source/examples/index.md b/doc/sphinx/source/examples/index.md index d0107ba5..92d7a3e2 100644 --- a/doc/sphinx/source/examples/index.md +++ b/doc/sphinx/source/examples/index.md @@ -17,6 +17,7 @@ Navier-Stokes/index Mixed/index Integration/index Time-Integrators/index +Tutorials/index ``` More examples will be added in the future. diff --git a/examples/cpp/backward_euler.cpp b/examples/cpp/backward_euler.cpp new file mode 100644 index 00000000..4dee7c96 --- /dev/null +++ b/examples/cpp/backward_euler.cpp @@ -0,0 +1,66 @@ +/** + * @file backward_euler.cpp + * @brief Solves a first-order ODE with backward Euler with an initial value condition. + * Equation: $\frac{dy}{dt} = \sin^2(t) \cdot y$ + * Domain: N/A + * Boundary Conditions: N/A + */ + +#include "mole.h" +#include + +// Function declaration for f(t, y) +double f(double t, double y); + +int main() { + // Problem parameters + constexpr double h = 0.01; // Step-size h + arma::vec t = arma::regspace(0, h, 5); // Calculates up to y(5) at step-size h + arma::mat y = arma::mat(1, t.size(), arma::fill::zeros); + constexpr double tol = 1e-6; // tolerance for fixed-point iteration + + // Initial conditions + y(0) = 2.0; + + // Backward Euler + for (int i = 0; i < t.size()-1; i++) { + double y_old = y(i); + for (int k = 0; k < 100; k++) { // fixed-point iteration for rootfinding + y(i+1) = y_old + h*f(t(i+1), y_old); // Backward Euler + if (std::abs(y(i+1) - y_old) < tol) { // Stopping criteria for approximate relative error + break; + } + } + } + + // Create a GNUplot script file + std::ofstream plot_script("plot.gnu"); + if (!plot_script) { + std::cerr << "Error: Failed to create GNUplot script.\n"; + return 1; + } + plot_script << "set title 'Approximation to y(t) using Backward Euler'\n"; + plot_script << "set xlabel 't'\n"; + plot_script << "set ylabel 'y'\n"; + plot_script << "plot '-' using 1:2 with lines\n"; + + // Print the time and solution values + for (int i = 0; i < t.size(); ++i) { + // output to stdout + std::cout << t(i) << " " << y(i) << "\n"; + // AND output to plot_script (plot.gnu) + plot_script << t(i) << " " << y(i) << "\n"; + } + plot_script.close(); + + // Execute GNUplot using the script + if (system("gnuplot -persist plot.gnu") != 0) { + std::cerr << "Error: Failed to execute GNUplot.\n"; + return 1; + } +} + +// Function definition for f(t, y) +double f(double t, double y) { + return std::pow(std::sin(t), 2) * y; +} \ No newline at end of file diff --git a/examples/cpp/elliptic1DHomogeneousDirichlet.cpp b/examples/cpp/elliptic1DHomogeneousDirichlet.cpp new file mode 100644 index 00000000..8ffe7965 --- /dev/null +++ b/examples/cpp/elliptic1DHomogeneousDirichlet.cpp @@ -0,0 +1,87 @@ +/** + * @file elliptic1DHomogeneousDirichlet.cpp + * @brief Solves the 1D linear equation -u'' = 1 + * + * ## Spatial Domains: + * - The spatial domain is [0, 1] + * - Interior points are spaced by dx = (b - a) / m. + * + * ## Boundary Conditions: + * u(0) = 0, u(1) = 0 + * + * The solution is computed numerically, and the result is compared with the exact solution: + * u_exact(x) = x(1 - x) / 2 + * + * The results are saved to a file "plot.gnu" and visualized using GNUplot. +*/ + +#include "mole.h" +#include + +int main() { + + const int k = 2; + const int m = 2 * k + 1; + const Real dx = 1.0 / m; + const Real a = 0; // left boundary + const Real b = 1; // right boundary + + // Mimetic operators + Laplacian L(k, m, dx); + Real d = 1; // Dirichlet coeff + Real n = 0; // Neumann coeff + RobinBC BC(k, m, dx, d, n); + L = L + BC; + + // 1D grid + arma::vec grid(m + 2); + grid(0) = a; + grid(1) = grid(0) + dx / 2.0; + for (int i = 2; i <= m; i++) { + grid(i) = grid(i - 1) + dx; + } + grid(m + 1) = b; + + // RHS + arma::vec rhs(m+2); + rhs.fill(-1.0); + rhs(0) = 0.0; + rhs(m+1) = 0.0; + + // Solve the system + #ifdef EIGEN + // Use eigen if available + arma::vec sol = Utils:: spsolve_eigen(L, rhs); + #else + arma::vec sol = spsolve(L, rhs); + #endif + + // Create a GNUplot script file + std::ofstream plot_script("plot.gnu"); + if (!plot_script) { + std::cerr << "Error: Failed to create GNUplot script.\n"; + return 1; + } + plot_script << "set title \"-u'' = 1, u(0) = 0, u(1) = 0\"\n"; + plot_script << "set xlabel 'x'\n"; + plot_script << "set ylabel 'y'\n"; + plot_script << "plot '-' using 1:2 with lines title \"Estimated Solution\", " + << "x*(1-x)/2 with lines title \"Exact Solution\"\n"; + + + for (int i = 0; i <= m + 1; ++i) { + plot_script << grid(i) << " " << sol(i) << "\n"; + } + + plot_script.close(); + + // Execute GNUplot using the script + if (system("gnuplot -persist plot.gnu") != 0) { + std::cerr << "Error: Failed to execute GNUplot.\n"; + return 1; + } + + cout << sol; + + return 0; +} \ No newline at end of file diff --git a/examples/cpp/elliptic1DLeftNeumannRightDirichlet.cpp b/examples/cpp/elliptic1DLeftNeumannRightDirichlet.cpp new file mode 100644 index 00000000..d1f0a42a --- /dev/null +++ b/examples/cpp/elliptic1DLeftNeumannRightDirichlet.cpp @@ -0,0 +1,86 @@ +/** + * @file elliptic1DLeftNeumannRightDirichlet.cpp + * @brief Solves the 1D linear equation -u'' = 1 + * + * ## Spatial Domains: + * - The spatial domain is [0, 1] + * - Interior points are spaced by dx = (b - a) / m. + * + * ## Boundary Conditions: + * u'(0) = 0, u(1) = 0 + * + * The solution is computed numerically, and the result is compared with the exact solution: + * u_exact(x) = (1 - x^2)/2 + * + * The results are saved to a file "plot.gnu" and visualized using GNUplot. +*/ + +#include "mole.h" +#include + +int main() { + + const int k = 2; + const int m = 2 * k + 1; + const Real dx = 1.0 / m; + const Real a = 0; // left boundary + const Real b = 1; // right boundary + + // Mimetic operators + Laplacian L(k, m, dx); + MixedBC BC(k, m, dx, + "Neumann", std::vector{1.0}, + "Dirichlet", std::vector{1.0}); + L += BC; + + // 1D grid + arma::vec grid(m + 2); + grid(0) = a; + grid(1) = grid(0) + dx / 2.0; + for (int i = 2; i <= m; i++) { + grid(i) = grid(i - 1) + dx; + } + grid(m + 1) = b; + + // RHS + arma::vec rhs(m+2); + rhs.fill(-1.0); + rhs(0) = 0.0; + rhs(m+1) = 0.0; + + // Solve the system + #ifdef EIGEN + // Use eigen if available + arma::vec sol = Utils:: spsolve_eigen(L, rhs); + #else + arma::vec sol = spsolve(L, rhs); + #endif + + // Create a GNUplot script file + std::ofstream plot_script("plot.gnu"); + if (!plot_script) { + std::cerr << "Error: Failed to create GNUplot script.\n"; + return 1; + } + plot_script << "set title \"-u'' = 1, u(0) = 0, u'(1) = 0\"\n"; + plot_script << "set xlabel 't'\n"; + plot_script << "set ylabel 'y'\n"; + plot_script << "plot '-' using 1:2 with lines title \"Estimated Solution\", " + << "(1 - x**2)/2 with lines title \"Exact Solution\"\n"; + + for (int i = 0; i <= m + 1; ++i) { + plot_script << grid(i) << " " << sol(i) << "\n"; + } + plot_script.close(); + + // Execute GNUplot using the script + if (system("gnuplot -persist plot.gnu") != 0) { + std::cerr << "Error: Failed to execute GNUplot.\n"; + return 1; + } + + + cout << sol; + + return 0; +} \ No newline at end of file diff --git a/examples/cpp/elliptic1DLeftNeumannRightNeumann.cpp b/examples/cpp/elliptic1DLeftNeumannRightNeumann.cpp new file mode 100644 index 00000000..6261866c --- /dev/null +++ b/examples/cpp/elliptic1DLeftNeumannRightNeumann.cpp @@ -0,0 +1,93 @@ +/** + * @file elliptic1DLeftNeumannRightNeumann.cpp + * @brief Solves the 1D linear equation -u'' = x - 1/2 + * + * ## Spatial Domains: + * - The spatial domain is [0, 1] + * - Interior points are spaced by dx = (b - a) / m. + * + * ## Boundary Conditions: + * u'(0) = 0, u'(1) = 0 + * + * The solution is computed numerically, and the result is compared with the exact solution: + * u_exact(x) = x^2/4 - x^3/6 + * + * The results are saved to a file "plot.gnu" and visualized using GNUplot. +*/ + +#include "mole.h" +#include + +int main() { + + const int k = 2; + const int m = 10; + const Real dx = 1.0 / m; + const Real a = 0; // left boundary + const Real b = 1; // right boundary + + // Mimetic operators + Laplacian L(k, m, dx); + const Real d = 0; // Dirichlet coeff + const Real n = 1; // Neumann coeff + RobinBC BC(k, m, dx, d, n); + L += BC; + + // 1D grid + arma::vec grid(m + 2); + grid(0) = a; + grid(1) = grid(0) + dx / 2.0; + for (int i = 2; i <= m; i++) { + grid(i) = grid(i - 1) + dx; + } + grid(m + 1) = b; + + // RHS + arma::vec rhs(m+2); rhs.zeros(); + for (int i=1; i<=m; ++i) { + rhs(i) = -(grid(i) - 0.5); + } + rhs(0) = 0.0; + rhs(m+1) = 0.0; + + // Solve the system + #ifdef EIGEN + // Use eigen if available + arma::vec sol = Utils:: spsolve_eigen(L, rhs); + #else + arma::vec sol = spsolve(L, rhs); + #endif + + // shift + const Real val = sol(0); + for (int j = 0; j < m+2; ++j) { + sol(j) -= val; + } + + // Create a GNUplot script file + std::ofstream plot_script("plot.gnu"); + if (!plot_script) { + std::cerr << "Error: Failed to create GNUplot script.\n"; + return 1; + } + plot_script << "set title \"-u'' = x - 1/2, u'(0) = 0, u'(1) = 0\"\n"; + plot_script << "set xlabel 't'\n"; + plot_script << "set ylabel 'y'\n"; + plot_script << "plot '-' using 1:2 with lines title \"Estimated Solution\", " + << "x**2 / 4 - x**3 / 6 with lines title \"Exact Solution\"\n"; + + for (int i = 0; i <= m + 1; ++i) { + plot_script << grid(i) << " " << sol(i) << "\n"; + } + plot_script.close(); + + // Execute GNUplot using the script + if (system("gnuplot -persist plot.gnu") != 0) { + std::cerr << "Error: Failed to execute GNUplot.\n"; + return 1; + } + + cout << sol; + + return 0; +} \ No newline at end of file diff --git a/examples/cpp/elliptic1DNonHomogeneousDirichlet.cpp b/examples/cpp/elliptic1DNonHomogeneousDirichlet.cpp new file mode 100644 index 00000000..3365d4d0 --- /dev/null +++ b/examples/cpp/elliptic1DNonHomogeneousDirichlet.cpp @@ -0,0 +1,87 @@ + /** + * @file elliptic1DNonHomogeneousDirichlet.cpp + * @brief Solves the 1D linear equation -u'' = 1 + * + * ## Spatial Domains: + * - The spatial domain is [0, 1] + * - Interior points are spaced by dx = (b - a) / m. + * + * ## Boundary Conditions: + * u(0) = 1/2, u(1) = 1/2 + * + * The solution is computed numerically, and the result is compared with the exact solution: + * u_exact(x) = (-x^2 + x + 1)/2 + * + * The results are saved to a file "plot.gnu" and visualized using GNUplot. +*/ + +#include "mole.h" +#include + +int main() { + + const int k = 2; + const int m = 2 * k + 1; + const Real dx = 1.0 / m; + const Real a = 0; // left boundary + const Real b = 1; // right boundary + const Real left_dirichlet = 0.5; + const Real right_dirichlet = 0.5; + + // Mimetic operators + Laplacian L(k, m, dx); + const Real d = 1; // Dirichlet coeff + const Real n = 0; // Neumann coeff + RobinBC BC(k, m, dx, d, n); + L += BC; + + // 1D grid + arma::vec grid(m + 2); + grid(0) = a; + grid(1) = grid(0) + dx / 2.0; + for (int i = 2; i <= m; i++) { + grid(i) = grid(i - 1) + dx; + } + grid(m + 1) = b; + + // RHS + arma::vec rhs(m+2); + rhs.fill(-1.0); + rhs(0) = left_dirichlet; + rhs(m+1) = right_dirichlet; + + // Solve the system + #ifdef EIGEN + // Use eigen if available + arma::vec sol = Utils:: spsolve_eigen(L, rhs); + #else + arma::vec sol = spsolve(L, rhs); + #endif + + // Create a GNUplot script file + std::ofstream plot_script("plot.gnu"); + if (!plot_script) { + std::cerr << "Error: Failed to create GNUplot script.\n"; + return 1; + } + plot_script << "set title \"-u'' = 1, u(0) = 1/2, u(1) = 1/2\"\n"; + plot_script << "set xlabel 't'\n"; + plot_script << "set ylabel 'y'\n"; + plot_script << "plot '-' using 1:2 with lines title \"Estimated Solution\", " + << "(-x**2 + x + 1)/2 with lines title \"Exact Solution\"\n"; + + for (int i = 0; i <= m + 1; ++i) { + plot_script << grid(i) << " " << sol(i) << "\n"; + } + plot_script.close(); + + // Execute GNUplot using the script + if (system("gnuplot -persist plot.gnu") != 0) { + std::cerr << "Error: Failed to execute GNUplot.\n"; + return 1; + } + + cout << sol; + + return 0; +} \ No newline at end of file diff --git a/examples/cpp/hyperbolic1D_upwind.cpp b/examples/cpp/hyperbolic1D_upwind.cpp index 476a309e..4cda8af8 100644 --- a/examples/cpp/hyperbolic1D_upwind.cpp +++ b/examples/cpp/hyperbolic1D_upwind.cpp @@ -124,6 +124,9 @@ int main() scriptFile.close(); // Run GNUplot with the script - system("gnuplot -persistent gp_script"); + if (system("gnuplot -persistent gp_script") != 0) { + cerr << "Error: Failed to execute GNUplot.\n"; + return EXIT_FAILURE; + } return EXIT_SUCCESS; } diff --git a/examples/cpp/sturmLiouvilleBessel.cpp b/examples/cpp/sturmLiouvilleBessel.cpp new file mode 100644 index 00000000..c1cde3ac --- /dev/null +++ b/examples/cpp/sturmLiouvilleBessel.cpp @@ -0,0 +1,76 @@ +/** + * @file sturmLiouvilleBessel.cpp + * @brief Solves the 1D Bessel function of the first kind and third order in Sturm-Liouville form + * + * The equation being solved is: + * $$ x^2 * u'' + x u' + (x^2 - \nu^2) * u = 0 $$ + * + * ## Spatial Domain: + * - The spatial domain is $x \in [0, 1]$ + * - The grid spacing is $dx = (2 k + 1)^{-1}$ + * + * ## Boundary Conditions: + * - $u(0) = 0$ + * - $u(1) = J_3(1)$ + * + * The solution is computed using mimetic finite difference operators and is compared with the exact result: + * $$ u(x) = J_3(x) $$ + * + * The norms of each solution and the error are printed + */ + + +#include "mole.h" +#include + +int main() +{ + // Parameters + const int k = 2; // Order of accuracy + const int m = 2 * k + 1; // Number of cells + const Real dx = 1.0 / m; // Grid spacing + const Real nu = 3.0; + + // Build grid of cell centers + arma::vec xc(m+2); + xc(0) = 0.0; + xc(1) = dx / 2.0; + for (int i = 2; i <= m; i++) + { + xc(i) = xc(i - 1) + dx; + } + xc(m + 1) = 1.0; + + // Exact solution -- J_3(xc) + arma::vec ue("0.0 0.000020820315755 0.000559343047749 0.002563729994587 0.006929654826751 0.014434028475866 0.019563353982668"); + + // Mimetic Operators + Interpol I(false, m, 0.5); // Interpolates from faces to centers + Gradient G(k, m, dx); + Laplacian L(k, m, dx); + RobinBC robin(k, m, dx, 1, 0); // Dirichlet BC + + // Set up system of equations + arma::sp_mat xc_mat = arma::sp_mat( arma::diagmat(xc) ); // x + arma::sp_mat xc_mat_sq = arma::sp_mat( arma::diagmat( arma::pow(xc, 2) ) ); // x^2 + arma::sp_mat xc_mat_sq_sub = arma::sp_mat( arma::diagmat( arma::pow(xc, 2) - nu*nu ) ); // x^2 - nu^2 + + arma::sp_mat A = xc_mat_sq * (arma::sp_mat)L + xc_mat * (arma::sp_mat)I * (arma::sp_mat)G + xc_mat_sq_sub * arma::speye(m+2, m+2); + + // Apply BC + A.row(0).zeros(); + A.row(A.n_rows - 1).zeros(); + A = A + (arma::sp_mat)robin; + + arma::vec b(m+2); + b(m+1) = 0.019563353982668; // J_3(1) + + // Solve + arma::vec sol = arma::spsolve(A, b); + + arma::vec diff = sol - ue; + std::cout << "norm(u_numerical) = " << arma::norm(sol) << std::endl; + std::cout << "norm(u_exact) = " << arma::norm(ue) << std::endl; + std::cout << "norm(u_numerical - u_exact) = " << arma::norm(diff) << std::endl; + +} \ No newline at end of file diff --git a/examples/cpp/sturmLiouvilleChebyshev.cpp b/examples/cpp/sturmLiouvilleChebyshev.cpp new file mode 100644 index 00000000..0f4bbe67 --- /dev/null +++ b/examples/cpp/sturmLiouvilleChebyshev.cpp @@ -0,0 +1,74 @@ +/** + * @file sturmLiouvilleChebyshev.cpp + * @brief Solves the Chebyshev equation in Sturm-Liouville form + * + * The equation being solved is: + * $$ (1 - x^2) * u'' - x * u' + p^2 * u = 0 $$ + * + * ## Spatial Domain: + * - The spatial domain is $x \in [-1, 1]$ + * - The grid spacing is $dx = 2 / m$ + * + * ## Boundary Conditions: + * - $u(-1) = u(1) = 1$ + * + * The solution is computed using mimetic finite difference operators and is compared with the exact result: + * $$ u(x) = T_2(x) $$ (Chebyshev Polynomial of degree 2) + * + * The norms of each solution and the error are printed + */ + + +#include "mole.h" +#include + +int main() +{ + // Parameters + const int k = 2; // Order of accuracy + const int m = 2 * k + 1; // Number of cells + const Real dx = 2.0 / m; // Grid spacing + const Real p = 2.0; + + // Build grid of cell centers + arma::vec xc(m+2); + xc(0) = -1.0; + xc(1) = -1.0 + dx / 2.0; + for (int i = 2; i <= m; i++) + { + xc(i) = xc(i - 1) + dx; + } + xc(m + 1) = 1.0; + + // Exact solution -- T_2(xc) + arma::vec ue("1.0 0.28 -0.68 -1.0 -0.68 0.28 1.0"); + + // Mimetic Operators + Interpol I(false, m, 0.5); // Interpolates from faces to centers + Gradient G(k, m, dx); + Laplacian L(k, m, dx); + RobinBC robin(k, m, dx, 1, 0); // Dirichlet BC + + // Set up system of equations + arma::sp_mat xc_mat = arma::sp_mat( arma::diagmat(xc) ); // x + arma::sp_mat xc_mat_sq = arma::sp_mat( arma::diagmat( 1.0 - arma::pow(xc, 2) ) ); // 1 - x^2 + arma::sp_mat A = xc_mat_sq * (arma::sp_mat)L - xc_mat * (arma::sp_mat)I * (arma::sp_mat)G + p * p * arma::speye(m+2, m+2); + + // Apply BC + A.row(0).zeros(); + A.row(A.n_rows - 1).zeros(); + A = A + (arma::sp_mat)robin; + + arma::vec b(m+2); + b(0) = 1.0; + b(m+1) = 1.0; + + // Solve + arma::vec sol = arma::spsolve(A, b); + + arma::vec diff = sol - ue; + std::cout << "norm(u_numerical) = " << arma::norm(sol) << std::endl; + std::cout << "norm(u_exact) = " << arma::norm(ue) << std::endl; + std::cout << "norm(u_numerical - u_exact) = " << arma::norm(diff) << std::endl; + +} \ No newline at end of file diff --git a/examples/cpp/sturmLiouvilleHelmholtzDirichletDirichlet.cpp b/examples/cpp/sturmLiouvilleHelmholtzDirichletDirichlet.cpp new file mode 100644 index 00000000..9383c6dd --- /dev/null +++ b/examples/cpp/sturmLiouvilleHelmholtzDirichletDirichlet.cpp @@ -0,0 +1,70 @@ +/** + * @file sturmLiouvilleHelmholtzDirichletDirichlet.cpp + * @brief Solves the 1D Helmholtz equation in Sturm-Liouville form + * + * The equation being solved is: + * $$ u'' + u = 0 $$ + * + * ## Spatial Domain: + * - The spatial domain is $x \in [0, 3]$ + * - The grid spacing is $dx = 3 / m$ + * + * ## Boundary Conditions: + * - $u(0) = 0$ + * - $u(3) = \sin(3)$ + * + * The solution is computed using mimetic finite difference operators and is compared with the exact result: + * $$ u(x) = \sin(x) $$ + * + * The norms of each solution and the error are printed + */ + + +#include "mole.h" +#include + +int main() +{ + // Parameters + const int k = 2; // Order of accuracy + const int m = 40; // Number of cells + const Real dx = 3.0 / m; // Grid spacing + + // Build grid of cell centers + arma::vec xc(m+2); + xc(0) = 0.0; + xc(1) = dx / 2.0; + for (int i = 2; i <= m; i++) + { + xc(i) = xc(i - 1) + dx; + } + xc(m + 1) = 3.0; + + // Exact solution -- sin(x) + arma::vec ue = sin(xc); + + // Mimetic Operators + Laplacian L(k, m, dx); + RobinBC robin(k, m, dx, 1, 0); // Dirichlet BC + + // Set up system of equations + arma::sp_mat A = (arma::sp_mat)L + arma::speye(m+2, m+2); + + // Apply BC + A.row(0).zeros(); + A.row(A.n_rows - 1).zeros(); + A = A + (arma::sp_mat)robin; + + arma::vec b(m+2); + b(0) = 0.0; + b(m+1) = sin(3.0); + + // Solve + arma::vec sol = arma::spsolve(A, b); + + arma::vec diff = sol - ue; + std::cout << "norm(u_numerical) = " << arma::norm(sol) << std::endl; + std::cout << "norm(u_exact) = " << arma::norm(ue) << std::endl; + std::cout << "norm(u_numerical - u_exact) = " << arma::norm(diff) << std::endl; + +} \ No newline at end of file diff --git a/examples/cpp/sturmLiouvilleHelmholtzDirichletRobin.cpp b/examples/cpp/sturmLiouvilleHelmholtzDirichletRobin.cpp new file mode 100644 index 00000000..2a5799c0 --- /dev/null +++ b/examples/cpp/sturmLiouvilleHelmholtzDirichletRobin.cpp @@ -0,0 +1,78 @@ +/** + * @file sturmLiouvilleHelmholtzDirichletRobin.cpp + * @brief Solves the 1D Helmholtz equation in Sturm-Liouville form + * + * The equation being solved is: + * $$ u'' + \mu^2 u = 0 $$ + * + * ## Spatial Domain: + * - The spatial domain is $x \in [0, 1] + * - The grid spacing is $dx = 1 / m$ + * + * ## Boundary Conditions: + * - $u'(0) = 0$ + * - $u(1) + u'(1) = \cos(\mu) - mu * \sin(\mu)$ + * + * The solution is computed using mimetic finite difference operators and is compared with the exact result: + * $$ u(x) = \cos(\mu * x) $$ + * + * The norms of each solution and the error are printed + */ + + +#include "mole.h" +#include + +int main() +{ + // Parameters + const int k = 2; // Order of accuracy + const int m = 150; // Number of cells + const Real dx = 1.0 / m; // Grid spacing + const Real mu = 0.86; + + // Build grid of cell centers + arma::vec xc(m+2); + xc(0) = 0.0; + xc(1) = dx / 2.0; + for (int i = 2; i <= m; i++) + { + xc(i) = xc(i - 1) + dx; + } + xc(m + 1) = 1.0; + + // Exact solution -- cos(mu * x) + arma::vec ue = cos(mu * xc); + + // Mimetic Operators + Laplacian L(k, m, dx); + const std::string left = "Neumann"; + const std::string right = "Robin"; + std::vector lbc(2); + std::vector rbc(2); + lbc[0] = 1.0; + rbc[0] = 1.0; + rbc[1] = 1.0; + MixedBC mbc(k, m, dx, "Neumann", lbc, "Robin", rbc); // Left side Neumann, right side Robin DC + + // Set up system of equations + arma::sp_mat A = (arma::sp_mat)L + mu * mu * arma::speye(m+2, m+2); + + // Apply BC + A.row(0).zeros(); + A.row(A.n_rows - 1).zeros(); + A = A + (sp_mat)mbc; + + arma::vec b(m+2); + b(0) = 0.0; + b(m+1) = cos(mu) - mu * sin(mu); + + // Solve + arma::vec sol = arma::spsolve(A, b); + + arma::vec diff = sol - ue; + std::cout << "norm(u_numerical) = " << arma::norm(sol) << std::endl; + std::cout << "norm(u_exact) = " << arma::norm(ue) << std::endl; + std::cout << "norm(u_numerical - u_exact) = " << arma::norm(diff) << std::endl; + +} \ No newline at end of file diff --git a/examples/cpp/sturmLiouvilleHermite.cpp b/examples/cpp/sturmLiouvilleHermite.cpp new file mode 100644 index 00000000..7786a60f --- /dev/null +++ b/examples/cpp/sturmLiouvilleHermite.cpp @@ -0,0 +1,76 @@ +/** + * @file sturmLiouvilleHermite.cpp + * @brief Solves the 1D Hermite equation in Sturm-Liouvulle form + * + * The equation being solved is: + * $$ u'' - 2 * u' + 2 * m * u = 0 $$ + * + * ## Spatial Domain: + * - The spatial domain is $x \in [-1, 1] + * - The grid spacing is $dx = 2 / m$ + * + * ## Boundary Conditions: + * - $u(-1) = H_4(-1)$ + * - $u(1) = H_4(1)$ + * + * The solution is computed using mimetic finite difference operators and is compared with the exact result: + * $$ H_4(x) $$ (Hermite Polynomial of order 4) + * + * The norms of each solution and the error are printed + */ + + +#include "mole.h" +#include + +int main() +{ + // Parameters + const int k = 2; // Order of accuracy + const int m = 20; // Number of cells + const Real dx = 2.0 / m; // Grid spacing + const Real Hm = 4.0; + + // Build grid of cell centers + arma::vec xc(m+2); + xc(0) = -1.0; + xc(1) = -1.0 + dx / 2.0; + for (int i = 2; i <= m; i++) + { + xc(i) = xc(i - 1) + dx; + } + xc(m + 1) = 1.0; + + // Exact solution -- H_4(xc) + arma::vec ue("-20.0 -18.2879 -14.3279 -9.9375 -5.4239 -1.0559 2.9361 6.3601 " + "9.0625 10.9281 11.8801 11.8801 10.9281 9.0625 6.3601 2.9361 " + "-1.0559 -5.4239 -9.9375 -14.3279 -18.2879 -20.0"); + + // Mimetic Operators + Interpol I(false, m, 0.5); // Interpolates from faces to centers + Gradient G(k, m, dx); + Laplacian L(k, m, dx); + RobinBC robin(k, m, dx, 1, 0); // Dirichlet BC + + // Set up system of equations + arma::sp_mat xc_mat = arma::sp_mat( arma::diagmat(xc) ); + arma::sp_mat A = (arma::sp_mat)L - 2 * xc_mat * (arma::sp_mat)I * (arma::sp_mat)G + 2.0 * Hm * arma::speye(m+2, m+2); + + // Apply BC + A.row(0).zeros(); + A.row(A.n_rows - 1).zeros(); + A = A + (arma::sp_mat)robin; + + arma::vec b(m+2); + b(0) = -20.0; + b(m+1) = -20.0; + + // Solve + arma::vec sol = arma::spsolve(A, b); + + arma::vec diff = sol - ue; + std::cout << "norm(u_numerical) = " << arma::norm(sol) << std::endl; + std::cout << "norm(u_exact) = " << arma::norm(ue) << std::endl; + std::cout << "norm(u_numerical - u_exact) = " << arma::norm(diff) << std::endl; + +} \ No newline at end of file diff --git a/examples/cpp/sturmLiouvilleLaguerre.cpp b/examples/cpp/sturmLiouvilleLaguerre.cpp new file mode 100644 index 00000000..9dc9c200 --- /dev/null +++ b/examples/cpp/sturmLiouvilleLaguerre.cpp @@ -0,0 +1,84 @@ +/** + * @file sturmLiouvilleLaguerre.cpp + * @brief Solves the 1D Laguerre equation in Sturm-Liouville form + * + * The equation being solved is: + * $$ x * u'' + (1 - x) * u' + n * u = 0 $$ + * + * ## Spatial Domain: + * - The spatial domain is $x \in [0, 2] + * - The grid spacing is $dx = 2 / m$ + * + * ## Boundary Conditions: + * - $u(0) = L_4(0)$ + * - $u(2) = L_4(2)$ + * + * The solution is computed using mimetic finite difference operators and is compared with the exact result: + * $$ L_4(x) $$ (Laguerre Polynomial of order 4) + * + * The norms of each solution and the error are printed + */ + + +#include "mole.h" +#include + +int main() +{ + // Parameters + const int k = 2; // Order of accuracy + const int m = 30; // Number of cells + const Real dx = 2.0 / m; // Grid spacing + const Real n = 4.0; + + // Build gird of cell centers + arma::vec xc(m+2); + xc(0) = 0.0; + xc(1) = dx / 2.0; + for (int i = 2; i <= m; i++) + { + xc(i) = xc(i - 1) + dx; + } + xc(m + 1) = 2.0; + + // Exact solution -- L_4(xc) + arma::vec ue("1.0 0.869975360082304 0.6293375 0.413612397119342 " + "0.221654372427984 0.0523375 -0.095444393004115 " + "-0.222777726337449 -0.330729166666667 -0.420345627572016 " + "-0.492654269547325 -0.5486625 -0.589357973251029 " + "-0.615708590534979 -0.6286625 -0.629148096707819 " + "-0.618074022633745 -0.596329166666667 -0.564782664609053 " + "-0.524283899176955 -0.4756625 -0.419728343621399 " + "-0.357271553497942 -0.2890625 -0.215851800411523 " + "-0.138370318930041 -0.057329166666667 0.026580298353909 " + "0.112686471193416 0.2003375 0.288901286008230 0.333333333333333"); + + // Mimetic Operators + Interpol I(false, m, 0.5); // Interpolates from faces to centers + Gradient G(k, m, dx); + Laplacian L(k, m, dx); + RobinBC robin(k, m, dx, 1, 0); // Dirichlet BC + + // Set up system of equations + arma::sp_mat xc_mat = arma::sp_mat( arma::diagmat(xc) ); // x + arma::sp_mat xc_mat_sub = arma::sp_mat( arma::diagmat(1 - xc) ); // 1 - x + arma::sp_mat A = xc_mat * (arma::sp_mat)L + xc_mat_sub * (arma::sp_mat)I * (arma::sp_mat)G + n * arma::speye(m+2, m+2); + + // Apply BC + A.row(0).zeros(); + A.row(A.n_rows - 1).zeros(); + A = A + (arma::sp_mat)robin; + + arma::vec b(m+2); + b(0) = 1.0; + b(m+1) = 0.333333333333333; + + // Solve + arma::vec sol = arma::spsolve(A, b); + + arma::vec diff = sol - ue; + std::cout << "norm(u_numerical) = " << arma::norm(sol) << std::endl; + std::cout << "norm(u_exact) = " << arma::norm(ue) << std::endl; + std::cout << "norm(u_numerical - u_exact) = " << arma::norm(diff) << std::endl; + +} \ No newline at end of file diff --git a/examples/cpp/sturmLiouvilleLegendre.cpp b/examples/cpp/sturmLiouvilleLegendre.cpp new file mode 100644 index 00000000..b9d83a1a --- /dev/null +++ b/examples/cpp/sturmLiouvilleLegendre.cpp @@ -0,0 +1,79 @@ +/** + * @file sturmLiouvilleLegendre.cpp + * @brief Solves the 1D Legendre equation in Sturm-Liouville form + * + * The equation being solved is: + * $$ (1 - x^2) * u'' - 2 * x * u' + n * (n + 1) * u = 0 $$ + * + * ## Spatial Domain: + * - The spatial domain is $x \in [-1, 1]$ + * - The grid spacing is $dx = 2 / m$ + * + * ## Boundary Conditions: + * - $u(-1) = -1$ + * - $u(1) = 1$ + * + * The solution is computed using mimetic finite difference operators and is compared with the exact result: + * $$ L_3(x) $$ (Legendre Polynomial of order 3) + * + * The norms of each solution and the error are printed + */ + + +#include "mole.h" +#include + +int main() +{ + // Parameters + const int k = 2; // Order of accuracy + const int m = 20; // Number of cells + const Real dx = 2.0 / m; // Grid spacing + const Real n = 3.0; + + // Build grid of cell centers + arma::vec xc(m+2); + xc(0) = -1.0; + xc(1) = -1.0 + dx / 2.0; + for (int i = 2; i <= m; i++) + { + xc(i) = xc(i - 1) + dx; + } + xc(m + 1) = 1.0; + + // Exact solution -- L_3(xc) + arma::vec ue("-1.0 -0.7184375 -0.2603125 0.0703125 " + "0.2884375 0.4090625 0.4471875 0.4178125 " + "0.3359375 0.2165625 0.0746875 -0.0746875 " + "-0.2165625 -0.3359375 -0.4178125 -0.4471875 " + "-0.4090625 -0.2884375 -0.0703125 0.2603125 0.7184375 1.0"); + + // Mimetic operators + Interpol I(false, m, 0.5); // Interpolates from faces to centers + Gradient G(k, m, dx); + Laplacian L(k, m, dx); + RobinBC robin(k, m, dx, 1, 0); // Dirichlet BC + + // Set up system of equations + arma::sp_mat xc_mat = arma::sp_mat( arma::diagmat(xc) ); // x + arma::sp_mat xc_mat_sq = arma::sp_mat( arma::diagmat(1 - arma::pow(xc, 2)) ); // 1 - x^2 + arma::sp_mat A = xc_mat_sq * (arma::sp_mat)L - 2 * xc_mat * (arma::sp_mat)I * (arma::sp_mat)G + n * (n-1) * arma::speye(m+2, m+2); + + // Apply BC + A.row(0).zeros(); + A.row(A.n_rows - 1).zeros(); + A = A + (arma::sp_mat)robin; + + arma::vec b(m+2); + b(0) = -1.0; + b(m+1) = 1.0; + + // Solve + vec sol = arma::spsolve(A, b); + + vec diff = sol - ue; + std::cout << "norm(u_numerical) = " << arma::norm(sol) << std::endl; + std::cout << "norm(u_exact) = " << arma::norm(ue) << std::endl; + std::cout << "norm(u_numerical - u_exact) = " << arma::norm(diff) << std::endl; + +} \ No newline at end of file diff --git a/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m b/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m similarity index 93% rename from examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m rename to examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m index 53daff0f..27e7275d 100755 --- a/examples/matlab_octave/elliptic1DLeftDirichletRightNeumann.m +++ b/examples/matlab_octave/elliptic1DLeftNeumannRightDirichlet.m @@ -1,5 +1,5 @@ % ====================== Test 3 ===================== -% 1D Poisson BVP: Dirichlet, Neumann Homogeneous BC +% 1D Poisson BVP: Neumann, Dirichlet Homogeneous BC % - u'' = 1, 0 < x < 1, u'(0) = 0, u(1) = 0 % exact solution: u(x) = (1 - x^2)/2 % =================================================== diff --git a/examples/matlab_octave/tutorials/two_way_wave_equation/comparison_two_way_wave_md_vs_fd.m b/examples/matlab_octave/tutorials/two_way_wave_equation/comparison_two_way_wave_md_vs_fd.m new file mode 100644 index 00000000..8facb992 --- /dev/null +++ b/examples/matlab_octave/tutorials/two_way_wave_equation/comparison_two_way_wave_md_vs_fd.m @@ -0,0 +1,160 @@ +%% COMPARISON_TWO_WAY_WAVE_MD_VS_FD +% Compare 1D two-way wave equation solutions using Mimetic Differences (MD) +% and standard Finite Differences (FD). +% +% PDE: +% d^2u/dt^2 = c^2 * d^2u/dx^2, -2 <= x <= 2 +% +% METHODS: +% 1. mimetic_diff_two_way_wave_eq.m +% - Uses mimetic finite differences for the spatial Laplacian. +% - Preserves discrete conservation properties. +% +% 2. finite_diff_two_way_wave_eq.m +% - Uses standard finite differences (tridiagonal Laplacian). +% +% PURPOSE: +% - Compare numerical accuracy, computational cost, and runtime of MD vs FD. +% - Demonstrate differences in error, stability, and sparsity handling. +% +% DOMAIN & DISCRETIZATION: +% - Spatial domain: [-2, 2] +% - N spatial points / dx spacing (user-defined) +% - Time integration: centered-in-time, centered-in-space (CTCS) +% +% OUTPUTS (printed or plotted by this script): +% - U2_md : Final solution using mimetic differences +% - U2_fd : Final solution using finite differences +% - error_md : Error norm for MD +% - error_fd : Error norm for FD +% - walltime_md : Computation time for MD +% - walltime_fd : Computation time for FD +% - flops_md : Estimated FLOPs for MD +% - flops_fd : Estimated FLOPs for FD +% +% NOTES: +% - Both solvers use a two-step leapfrog (centered-time) scheme. +% - CFL condition must be satisfied: c*dt/dx <= 1 +% - Visualization and comparison are performed at the final time step. +% +% SEE ALSO: +% mimetic_diff_two_way_wave_eq, finite_diff_two_way_wave_eq +% +% EXPLANATION: +% We set the time step to be constant (dt=0.001) and vary the grid spacing to +% analyze the error of the spacial numerical scheme. The simulation runs such +% that the wave traverses 1 unit. The domain is set to be 4 units wide. +% This removes any boundary, so we can compare just the methods, without +% having to do anything for the boundary. A straightforward comparison of +% just the methods. +% +% The main difference between the methods is how the grid is discretized. +% Finite differences use a standard grid, where the domain from A to B is +% broken into N equal sized cells, with n+1 grid points demarking the +% boundaries. +% +% A B +% <-----dx-----> <-----dx-----> <-----dx-----> <-----dx-----> -- space +% |----cell 1----|----cell 2----| ... |----cell N-1--|----cell N----| -- cells +% 0 1 2 ...N-1 N N+1 -- index +% +% Where the mimetic difference uses a staggered grid, with a half step at each +% boundary. +% A B +% <--dx/2--> <-----dx-----> <-----dx-----> <--dx/2--> -- space +% |--cell 1--|----cell 2----| ... |----cell N-1--|--cell N--| -- cells +% 0 1 2 ... N N+1 N+2 -- index +% +% In numbers, from 0 to 1, with 5 cells, an FD grid and MD grid would be +% FD +% 0.0 0.2 0.4 0.6 0.8 1.0 +% o----------o----------o----------o----------o----------o +% x0 x1 x2 x3 x4 x5 +% +% MD +% 0.0 0.1 0.3 0.5 0.7 0.9 1.0 +% o-----o---------o----------o----------o----------o-----o +% x0 x1 x2 x3 x4 x5 x6 +% +% +% The CTCS has a second order accurate time scheme, so it will never get >2 +% convergence. The user can change the order of the mimetic scheme very easily +% in the mimetic_diff_two_way_wave_eq.m file (k=2), which shows how easy it is +% to switch orders of the experiment, for very little cost. +% +% ---------------------------------------------------------------------------- +% SPDX-License-Identifier: GPL-3.0-or-later +% © 2008-2024 San Diego State University Research Foundation (SDSURF). +% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +% ---------------------------------------------------------------------------- +% + +clc +close all +clear + +addpath('../../../../src/matlab_octave'); + +%% Problem definition - These values aremimicked in the FD and MD functions +k = 2; % Mimetic Order of Accuracy, can change to 4,6,8 +c = 0.1; % Velocity, 1 makes FD sceme exact +west = -2.0; % Domain's leftmost limits +east = 2.0; % Domain's rightmost limit + +% Number of cells to try, points is cells+1 +num_cells = [ 10, 20, 40, 80, 160 ]; + +%% Run each of the methods over the different grids +[ U2_fd, error_fd, walltime_fd, flops_fd ] = finite_diff_two_way_wave_eq(); +[ U2_md, error_md, walltime_md, flops_md ] = mimetic_diff_two_way_wave_eq(); + +%% Error analysis + +% We kept the time step constant, so any error, or better, any improvement +% measured is due to the spacial scheme. here we fit a regression plot to the +% log error of the lines. Orders of accuracy (k) are O^k, so we need to extract +% the exponent, hence the log of the errors and spacial steps. +% Regression fit line to error plot, should be close to 2. +p_fd = polyfit(log(1 ./ num_cells), log(error_fd), 1); +p_md = polyfit(log(1 ./ num_cells), log(error_md), 1); +fprintf("Finite Difference Error Convergence Slope: %3.4f\n", p_fd(1) ); +fprintf("Mimetic Difference Error Convergence Slope: %3.4f\n", p_md(1) ); + +% Plot the comparison of the errors of each scheme - on a log-log plot, +% so we can capture the slope of the lines == the order of accuracy of the +% spacial scheme. +figure(1) +loglog(1 ./ num_cells, error_fd, 'LineWidth', 2); hold on; +loglog(1 ./ num_cells, error_md, 'LineWidth', 2) +xlim([1 / num_cells(end) 1 / num_cells(1)]); +xlabel('dx'); ylabel('error'); +grid on; +title(['Error: FD slope~', num2str(p_fd(1)), ', MD slope~', num2str(p_md(1))]); +legend('FD Error','MD error') +set(gca, "linewidth", 2, "fontsize", 16) + +% Walltime, there is some spooling MATLAB does, so may not be +% super accurate. Given here for reference. +figure(2) +plot(num_cells, walltime_fd, 'LineWidth', 2); hold on; +plot(num_cells, walltime_md, 'LineWidth', 2); +xlim([num_cells(1) num_cells(end)]); +ylim([min( min(walltime_fd), min(walltime_md) ) max( max(walltime_fd), max(walltime_md) )]); +xlabel('num points'); ylabel('walltime [s]'); +grid on; +title(['Walltime, MD is order ', num2str(k)]); +legend('FD time', 'MD time'); +set(gca, "linewidth", 2, "fontsize", 16) + +% Floating point operations needed for each method. The mimetic difference +% methods use at least one more point, and have boundary calculations built in. +figure(3) +loglog(num_cells, flops_fd, 'LineWidth', 2); hold on; +loglog(num_cells, flops_md, 'LineWidth', 2); +xlim([ num_cells(1) num_cells(end) ]); ylim([ 1e5 1e8 ]); +xlabel('num cells'); ylabel('FLOPs'); +grid on; +title(['FLOPs for each method, mimetic order:', num2str(k)]); +legend('FD FLOPs', 'MD FLOPs'); +set(gca, "linewidth", 2, "fontsize", 16) + diff --git a/examples/matlab_octave/tutorials/two_way_wave_equation/finite_diff_two_way_wave_eq.m b/examples/matlab_octave/tutorials/two_way_wave_equation/finite_diff_two_way_wave_eq.m new file mode 100644 index 00000000..4533a985 --- /dev/null +++ b/examples/matlab_octave/tutorials/two_way_wave_equation/finite_diff_two_way_wave_eq.m @@ -0,0 +1,132 @@ +function [ U2_fd, error_fd, walltime_fd, flops_fd ] = finite_diff_two_way_wave_eq() +%% FINITE_DIFF_TWO_WAY_WAVE_EQ +% Solve the 1D two-way wave equation using standard finite differences. +% +% PDE: d^2u/dt^2 = c^2 * d^2u/dx^2 on domain [-2, 2] +% Scheme: Centered-in-Time, Centered-in-Space (CTCS) +% Explicit second-order accurate in time and space +% +% Time-stepping update: +% U^{n+1} = 2*U^{n} - U^{n-1} + (c^2 * dt^2 * D_fd) * U^{n} +% where D_fd is the standard second-derivative finite-difference matrix: +% +% D_fd = (1/dx^2) * tridiag([1, -2, 1]) +% +% OUTPUTS: +% U2_fd - Final solution vector at last time step +% error_fd - Norm of error vs. analytic or reference solution +% walltime_fd - Wall-clock time (seconds) for the simulation +% flops_fd - Estimated floating-point operation count +% +% SEE ALSO: +% mimetic_diff_two_way_wave_eq, comparison_two_way_wave_md_vs_fd +% +% NOTES: +% We use a larger domain for graphing wave movement without boundary issues. +% This removes any boundary, so we can compare just the methods, without +% having to do anything for the boundary. A straightforward comparison of +% just the methods. +% +% The discretization has a second order accurate time scheme, so it will never +% get >2 convergence. This file is commented for first time users, to explain +% and show differences with the mimetic difference method. +% +% ---------------------------------------------------------------------------- +% SPDX-License-Identifier: GPL-3.0-or-later +% © 2008-2024 San Diego State University Research Foundation (SDSURF). +% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +% ---------------------------------------------------------------------------- +% + +%% Problem definition +c = 0.1; % Velocity, 1 makes FD scheme exact +west = -2.0; % Domain's leftmost limits +east = 2.0; % Domain's rightmost limit + +%% Number of cells to try, grid points is cells+1 +num_cells = [ 10, 20, 40, 80, 160 ]; + +% generic holders for metrics within the loops +error_fd = zeros(size(num_cells)); +walltime_fd = zeros(size(num_cells)); +flops_fd = zeros(size(num_cells)); + +% Initial Condition Function +f = @(x) ( (x > -0.5) & (x < 0.5) ) .* (cos(pi * x).^2); + +% Wave solution using d'Alembert +u = @(x,t) 0.5 * ( f(x - c * t) + f(x + c * t) ); + +for cell_index = 1 : numel(num_cells) + + %% Setup the domain + m = num_cells(cell_index); % Number of cells, mimetics uses 'cells' + nx = m+1; % number of grid points + + dx = (east - west) / m; % spacial discretization + dt = 0.001; % Time step constant for error analysis + + r2_fd = c^2 * (dt^2 / dx^2); % c in the equation + + t = ceil( 1/(c * dt) ); % first step euler, so t is one less + + % FD grid + grid_fd = west : dx : east; + + %% Initial Displacement is cos curve + U0_fd = u(grid_fd', 0); + + % The analytic solution, so we can check the error + analytic_fd = u(grid_fd', (t+1) * dt); + + %% Centered Scheme Matrix L + % L is (I + cD), a tridiagonal sparse matrix: + % + % L = + % [2-2c c 0 0 0 + % c 2-2c c 0 0 + % 0 α 2-2c c 0 + % 0 0 c 2-2c c + % 0 0 0 c 2-2c] + % + % where c = r2_fd = c^2 * (dt^2 / dx^2) + + % For step 1 + % Make a vector of ones, this will be the -1,+1 diagonal of matrix LTCS + B = ( r2_fd ) * ones(nx,1); + + % The next vector of ones will be the main diagonal of LTCS + A = (2 - (2 * r2_fd)) * ones(nx,1); + + % Build sparse matrix consisting of diagonals A and D + LTCS = spdiags([B A B], [-1,0,1], nx, nx); + + % Boundary Conditions, not needed with time<2, but here for reference + LTCS(1,1) = 1; LTCS(1,2) = 0; + LTCS(end,end) = 1; LTCS(end,end-1) = 0; + + % Step one is different, precomputed using Euler method + U1_fd = 0.5 * LTCS * U0_fd; + + U2_fd = U1_fd; % Increment time step + + %% Loop over time values + nnz_fd = nnz(LTCS); + tic + for i = 1 : t + + U2_fd = (LTCS * U1_fd) - U0_fd; + + % Shift everyone back for leapfrog scheme + U0_fd = U1_fd; + U1_fd = U2_fd; + + end + walltime_fd(cell_index) = toc; + + % Number of flops + flops_fd(cell_index) = (2 * nnz_fd + length(U0_fd)) * t; + error_fd(cell_index) = max(U2_fd-analytic_fd); + +end + diff --git a/examples/matlab_octave/tutorials/two_way_wave_equation/mimetic_diff_two_way_wave_eq.m b/examples/matlab_octave/tutorials/two_way_wave_equation/mimetic_diff_two_way_wave_eq.m new file mode 100644 index 00000000..045ac22b --- /dev/null +++ b/examples/matlab_octave/tutorials/two_way_wave_equation/mimetic_diff_two_way_wave_eq.m @@ -0,0 +1,113 @@ +function [ U2_md, error_md, walltime_md, flops_md ] = mimetic_diff_two_way_wave_eq() +%% MIMETIC_DIFF_TWO_WAY_WAVE_EQ +% Solve the 1D two-way wave equation using mimetic finite differences. +% +% PDE: d^2u/dt^2 = c^2 * d^2u/dx^2 on domain [-2, 2] +% Scheme: Centered-in-Time, Centered-in-Space (CTCS) +% Explicit second-order accurate in time and space +% +% Time-stepping update: +% U^{n+1} = 2*U^{n} - U^{n-1} + (c^2 * dt^2 * L) * U^{n} +% where L is the mimetic discrete Laplacian operator. Note there is no spacial +% discretization, the L takes care of that. +% +% OUTPUTS: +% U2_md - Final solution vector at last time step +% error_md - Norm of error vs. analytic or reference solution +% walltime_md - Wall-clock time (seconds) for the simulation +% flops_md - Estimated floating-point operation count +% +% SEE ALSO: +% finite_diff_two_way_wave_eq, comparison_two_way_wave_md_vs_fd +% +% NOTES: +% We use a larger domain for graphing wave movement without bondary issues. +% This removes any boundary, so we can compare just the methods, without +% having to do anything for the boundary. A straightforward comparison of +% just the methods. +% +% The discretization has a second order accurate time scheme, so it will never +% get >2 convergence. This file is commented for first time users, to explain +% and show differences with the mimetic difference method. +% +% The order of accuracy of the mimetic difference scheme can be changed with the +% k parameter. +% +% [ Compare with finite_diff_two_way_wave_eq.m ] +% +% ---------------------------------------------------------------------------- +% SPDX-License-Identifier: GPL-3.0-or-later +% © 2008-2024 San Diego State University Research Foundation (SDSURF). +% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +% ---------------------------------------------------------------------------- +% + +%% Problem definition +k = 2; % Mimetic Order of Accuracy, can change to 4,6,8 +c = 0.1; % Velocity, 1 makes FD sceme exact +west = -2.0; % Domain's leftmost limits +east = 2.0; % Domain's rightmost limit + +%% Number of cells to try, points is cells+1 +num_cells = [ 10, 20, 40, 80, 160 ]; + +% generic holders for loop info +error_md = zeros(size(num_cells)); +walltime_md = zeros(size(num_cells)); +flops_md = zeros(size(num_cells)); + +% Initial Condition Function +f = @(x) ( (x > -0.5) & (x < 0.5) ) .* (cos(pi * x).^2); + +% Wave solution using d'Almbert +u = @(x,t) 0.5 * ( f(x - c * t) + f(x + c * t) ); + +for cell_index = 1 : numel(num_cells) + + %% Setup the domain + m = num_cells(cell_index); % Number of cells, mimetics uses 'cells' + nx = m+1; % number of grid points + + dx = (east - west) / m; % spacial discretization + dt = 0.001; % Time step constant for error analysis + + r2_md = c^2 * (dt^2); % c in the equation, dx is built into the + % mimetic operator Laplacian (L) + + t = ceil( 1/(c * dt) ); % first step euler, so t is one less + + % MD grid, note the extra staggered (dx/2) step near the boundaries + grid_md = [ west west+dx/2 : dx : east-dx/2 east ]; + + %% Initial Displacement is cos curve + U0_md = u(grid_md', 0); + + % The analytic solution, so we can check the error + analytic_md = u(grid_md', (t+1) * dt); + + %% Mimetic Scheme Matrix L + L = lap(k,m,dx); % mimetic Laplacian operator + L = r2_md * L; % premultiply by dt^2 * c^2, knowing dx is already in L + + U1_md = U0_md + (0.5 * L * U0_md); % First step, Euler + U2_md = U1_md; % Increment time step + + %% Loop over time values + nnz_md = nnz(L); + tic + for i = 1 : t + + U2_md = (2.0 * U1_md) + (L * U1_md) - U0_md; + + % Shift everyone back for leapfrog scheme + U0_md = U1_md; + U1_md = U2_md; + + end + walltime_md(cell_index) = toc; + + % Number of flops + flops_md(cell_index) = (2 * nnz_md + (3 * length(U0_md)) ) * t; + error_md(cell_index) = max(U2_md-analytic_md); + +end diff --git a/src/cpp/mole.h b/src/cpp/mole.h index 3fe0a3e1..fae7777c 100644 --- a/src/cpp/mole.h +++ b/src/cpp/mole.h @@ -24,5 +24,7 @@ #include "operators.h" #include "robinbc.h" #include "utils.h" +#include "weightsQ.h" +#include "weightsP.h" #endif // MOLE_H diff --git a/src/cpp/utils.cpp b/src/cpp/utils.cpp index b6f02d6d..4642cc86 100644 --- a/src/cpp/utils.cpp +++ b/src/cpp/utils.cpp @@ -18,8 +18,9 @@ #include "utils.h" #include -#ifdef EIGEN +//#ifdef EIGEN #include +#include vec Utils::spsolve_eigen(const sp_mat &A, const vec &b) { Eigen::SparseMatrix eigen_A(A.n_rows, A.n_cols); @@ -47,7 +48,34 @@ vec Utils::spsolve_eigen(const sp_mat &A, const vec &b) { return vec(eigen_x.data(), eigen_x.size()); } -#endif + +vec Utils::spsolve_eigenQR(const sp_mat &A, const vec &b) { + Eigen::SparseMatrix eigen_A(A.n_rows, A.n_cols); + std::vector> triplets; + Eigen::SparseQR, Eigen::COLAMDOrdering> solver; + + Eigen::VectorXd eigen_x(A.n_rows); + triplets.reserve(5 * A.n_rows); + + auto it = A.begin(); + while (it != A.end()) { + triplets.push_back(Eigen::Triplet(it.row(), it.col(), *it)); + ++it; + } + + eigen_A.setFromTriplets(triplets.begin(), triplets.end()); + triplets.clear(); + + auto b_ = conv_to>::from(b); + Eigen::Map eigen_b(b_.data(), b_.size()); + + solver.analyzePattern(eigen_A); + solver.factorize(eigen_A); + eigen_x = solver.solve(eigen_b); + + return vec(eigen_x.data(), eigen_x.size()); +} +//#endif // Basic implementation of Kronecker product /* diff --git a/src/cpp/utils.h b/src/cpp/utils.h index dc860a85..e3fb27ac 100644 --- a/src/cpp/utils.h +++ b/src/cpp/utils.h @@ -68,6 +68,8 @@ class Utils { */ static vec spsolve_eigen(const sp_mat &A, const vec &b); + static vec spsolve_eigenQR(const sp_mat &A, const vec &b); + /** * @brief An analog to the MATLAB/Octave 2D meshgrid operation * diff --git a/src/cpp/weightsP.cpp b/src/cpp/weightsP.cpp new file mode 100644 index 00000000..9203d699 --- /dev/null +++ b/src/cpp/weightsP.cpp @@ -0,0 +1,34 @@ +/* +* SPDX-License-Identifier: GPL-3.0-or-later +* © 2008-2024 San Diego State University Research Foundation (SDSURF). +* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +*/ + +/* + * @file weights.cpp + * + * @brief Mimetic Differences Q Weights + * + * @date 2025/10/14 + * + */ + +#include "utils.h" +#include "divergence.h" +#include "gradient.h" +#include "weightsP.h" + + +WeightsP::WeightsP(u16 k, u32 m, Real dx) +{ + Gradient G(k, m, dx); + + vec b(m+2); + b.at(0) = -1.0; + b.at(m+1) = 1.0; + + + sp_mat Gtranspose = G.t(); + *this = Utils::spsolve_eigenQR(Gtranspose,b); + +} \ No newline at end of file diff --git a/src/cpp/weightsP.h b/src/cpp/weightsP.h new file mode 100644 index 00000000..ea490978 --- /dev/null +++ b/src/cpp/weightsP.h @@ -0,0 +1,42 @@ +/* +* SPDX-License-Identifier: GPL-3.0-or-later +* © 2008-2024 San Diego State University Research Foundation (SDSURF). +* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +*/ + + /* + * @file weights.h + * + * @brief Mimetic Differences Weights + * + * @date 2025/10/14 + * + */ + +#ifndef WEIGHTSP_H +#define WEIGHTSP_H + +#include "utils.h" +#include + +/** + * @brief Generate P Weights + * + */ +class WeightsP : public sp_vec { + +public: + using sp_vec::operator=; + + /** + * @brief P Weights Constructor + * + * @param k Order of accuracy + * @param m Number of cells + * @param dx Spacing between cells + */ + WeightsP(u16 k, u32 m, Real dx); + +}; + +#endif // WEIGHTSP_H diff --git a/src/cpp/weightsQ.cpp b/src/cpp/weightsQ.cpp new file mode 100644 index 00000000..23bef9dc --- /dev/null +++ b/src/cpp/weightsQ.cpp @@ -0,0 +1,43 @@ +/* +* SPDX-License-Identifier: GPL-3.0-or-later +* © 2008-2024 San Diego State University Research Foundation (SDSURF). +* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +*/ + +/* + * @file weights.cpp + * + * @brief Mimetic Differences Q Weights + * + * @date 2025/10/14 + * + */ + +#include "utils.h" +#include "divergence.h" +#include "gradient.h" +#include "weightsQ.h" + +WeightsQ::WeightsQ(u16 k, u32 m, Real dx): sp_vec(m + 1) +{ + Divergence D(k, m, dx); + D.shed_row(0); + D.shed_row(m); + + vec b(m+1); + b.at(0) = -1.0; + b.at(m) = 1.0; + + sp_mat Dtranspose = D.t(); + sp_vec Q_prime = Utils::spsolve_eigenQR(Dtranspose, b); + sp_vec Q(m+2); + Q.at(0) = 1.0; + Q.at(m+1) = 1.0; + for (int i = 1; i < m+1; ++i) { + Q.at(i) = Q_prime.at(i-1); + } + + *this = Q; + +} + diff --git a/src/cpp/weightsQ.h b/src/cpp/weightsQ.h new file mode 100644 index 00000000..37d6b09b --- /dev/null +++ b/src/cpp/weightsQ.h @@ -0,0 +1,43 @@ +/* +* SPDX-License-Identifier: GPL-3.0-or-later +* © 2008-2024 San Diego State University Research Foundation (SDSURF). +* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +*/ + + /* + * @file weights.h + * + * @brief Mimetic Differences Weights + * + * @date 2025/10/14 + * + */ + +#ifndef WEIGHTSQ_H +#define WEIGHTSQ_H + +#include "utils.h" +#include + +/** + * @brief Generate Q Weights + * + */ +class WeightsQ : public sp_vec { + +public: + using sp_vec::operator=; + + /** + * @brief Q Weights Constructor + * + * @param k Order of accuracy + * @param m Number of cells + * @param dx Spacing between cells + */ + WeightsQ(u16 k, u32 m, Real dx); + +}; + + +#endif // WEIGHTSQ_H diff --git a/src/dat/pweights.csv b/src/dat/pweights.csv deleted file mode 100644 index f19ba6bc..00000000 --- a/src/dat/pweights.csv +++ /dev/null @@ -1,78 +0,0 @@ -2,5,0.375000000,1.125000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,6,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,7,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,8,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,9,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,10,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,11,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,12,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,13,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,14,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,15,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,16,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,17,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,18,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,19,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,20,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,21,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,22,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,23,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,24,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,25,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,26,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,27,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,28,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,29,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -2,30,0.375000000,1.125000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.125000000,0.375000000 -4,9,0.354134518,1.228459404,0.898117099,1.018547095,1.000741884,1.000741884,1.018547095,0.898117099,1.228459404,0.354134518 -4,10,0.354134518,1.228459403,0.898117060,1.018546076,1.000715427,1.000055033,1.000715427,1.018546076,0.898117060,1.228459403,0.354134518 -4,11,0.354134518,1.228459403,0.898117058,1.018546036,1.000714408,1.000028576,1.000028576,1.000714408,1.018546036,0.898117058,1.228459403,0.354134518 -4,12,0.354134518,1.228459403,0.898117058,1.018546035,1.000714369,1.000027557,1.000002120,1.000027557,1.000714369,1.018546035,0.898117058,1.228459403,0.354134518 -4,13,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027518,1.000001101,1.000001101,1.000027518,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,14,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001061,1.000000082,1.000001061,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,15,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000042,1.000000042,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,16,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000003,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,17,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,18,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,19,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,20,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,21,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,22,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,23,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,24,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,25,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,26,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,27,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,28,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,29,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,30,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,31,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,32,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,33,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -4,34,0.354134518,1.228459403,0.898117058,1.018546035,1.000714367,1.000027516,1.000001060,1.000000041,1.000000002,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000002,1.000000041,1.000001060,1.000027516,1.000714367,1.018546035,0.898117058,1.228459403,0.354134518 -6,13,0.315722926,1.390677390,0.629532984,1.234237262,0.919144328,1.009804636,1.000880474,1.000880474,1.009804636,0.919144328,1.234237262,0.629532984,1.390677390,0.315722926 -6,14,0.315722926,1.390677390,0.629532983,1.234237260,0.919144539,1.009807801,1.000871112,1.000011976,1.000871112,1.009807801,0.919144539,1.234237260,0.629532983,1.390677390,0.315722926 -6,15,0.315722926,1.390677390,0.629532983,1.234237259,0.919144537,1.009808013,1.000874278,1.000002614,1.000002614,1.000874278,1.009808013,0.919144537,1.234237259,0.629532983,1.390677390,0.315722926 -6,16,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808011,1.000874489,1.000005780,0.999993252,1.000005780,1.000874489,1.009808011,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,17,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874487,1.000005991,0.999996418,0.999996418,1.000005991,1.000874487,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,18,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005989,0.999996629,0.999999583,0.999996629,1.000005989,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,19,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996627,0.999999794,0.999999794,0.999996627,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,20,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999793,1.000000006,0.999999793,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,21,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000004,1.000000004,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,22,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999791,1.000000003,1.000000002,1.000000003,0.999999791,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,23,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,24,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,25,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,26,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,27,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,28,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,29,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,30,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,31,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,32,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,33,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,34,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,35,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,36,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,37,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 -6,38,0.315722926,1.390677390,0.629532983,1.234237259,0.919144536,1.009808010,1.000874486,1.000005988,0.999996626,0.999999792,1.000000003,1.000000001,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000001,1.000000003,0.999999792,0.999996626,1.000005988,1.000874486,1.009808010,0.919144536,1.234237259,0.629532983,1.390677390,0.315722926 diff --git a/src/dat/qweights.csv b/src/dat/qweights.csv deleted file mode 100644 index fed7e11b..00000000 --- a/src/dat/qweights.csv +++ /dev/null @@ -1,78 +0,0 @@ -2,5,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,6,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,7,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,8,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,9,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,10,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,11,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,12,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,13,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,14,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,15,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,16,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,17,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,18,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,19,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,20,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,21,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,22,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,23,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,24,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,25,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,26,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,27,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,28,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,29,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -2,30,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 -4,9,1.000000000,1.125064293,0.751414447,1.162097097,0.962852897,0.997142531,0.962852897,1.162097097,0.751414447,1.125064293,1.000000000 -4,10,1.000000000,1.125064297,0.751414525,1.162099135,0.962905810,0.998516232,0.998516232,0.962905810,1.162099135,0.751414525,1.125064297,1.000000000 -4,11,1.000000000,1.125064297,0.751414528,1.162099214,0.962907849,0.998569145,0.999889934,0.998569145,0.962907849,1.162099214,0.751414528,1.125064297,1.000000000 -4,12,1.000000000,1.125064297,0.751414528,1.162099217,0.962907927,0.998571184,0.999942847,0.999942847,0.998571184,0.962907927,1.162099217,0.751414528,1.125064297,1.000000000 -4,13,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571262,0.999944885,0.999995760,0.999944885,0.998571262,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,14,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944964,0.999997799,0.999997799,0.999944964,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,15,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997877,0.999999837,0.999997877,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,16,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999915,0.999999915,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,17,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999994,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,18,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,19,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,20,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,21,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,22,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,23,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,24,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,25,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,26,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,27,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,28,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,29,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,30,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,31,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,32,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,33,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -4,34,1.000000000,1.125064297,0.751414528,1.162099217,0.962907930,0.998571265,0.999944967,0.999997880,0.999999918,0.999999997,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999997,0.999999918,0.999997880,0.999944967,0.998571265,0.962907930,1.162099217,0.751414528,1.125064297,1.000000000 -6,13,1.000000000,1.188528786,0.464036031,1.670433241,0.529401810,1.170898848,0.978586058,0.996230453,0.978586058,1.170898848,0.529401810,1.670433241,0.464036031,1.188528786,1.000000000 -6,14,1.000000000,1.188528786,0.464036033,1.670433245,0.529401355,1.170891964,0.978605184,0.998103433,0.998103433,0.978605184,1.170891964,0.529401355,1.670433245,0.464036033,1.188528786,1.000000000 -6,15,1.000000000,1.188528786,0.464036034,1.670433247,0.529401360,1.170891509,0.978598300,0.998122559,0.999976412,0.998122559,0.978598300,1.170891509,0.529401360,1.670433247,0.464036034,1.188528786,1.000000000 -6,16,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891513,0.978597845,0.998115674,0.999995538,0.999995538,0.998115674,0.978597845,1.170891513,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,17,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891515,0.978597849,0.998115220,0.999988654,1.000014664,0.999988654,0.998115220,0.978597849,1.170891515,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,18,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115224,0.999988200,1.000007780,1.000007780,0.999988200,0.998115224,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,19,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988204,1.000007326,1.000000896,1.000007326,0.999988204,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,20,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007330,1.000000441,1.000000441,1.000007330,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,21,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000446,0.999999987,1.000000446,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,22,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999991,0.999999991,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,23,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999995,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,24,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,25,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,26,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,27,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,28,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,29,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,30,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,31,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,32,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,33,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,34,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,35,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,36,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,37,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 -6,38,1.000000000,1.188528786,0.464036034,1.670433247,0.529401362,1.170891516,0.978597852,0.998115226,0.999988206,1.000007332,1.000000448,0.999999993,0.999999998,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,0.999999998,0.999999993,1.000000448,1.000007332,0.999988206,0.998115226,0.978597852,1.170891516,0.529401362,1.670433247,0.464036034,1.188528786,1.000000000 diff --git a/tests/cpp/test6.cpp b/tests/cpp/test6.cpp new file mode 100644 index 00000000..c4d3a776 --- /dev/null +++ b/tests/cpp/test6.cpp @@ -0,0 +1,54 @@ +#include "mole.h" +#include + +vec Q_baseline[3]; +vec P_baseline[3]; +Real tolerance = 1.0e-08; + +void run_Qtest(int k) { + Real dx = 1.0; + int m = k*2 + 1; + int j = k/2 - 1; + WeightsQ Q(k,m,dx); + Real total_error = 0.0; + + for (int i = 0; i < m+2; ++i) { + total_error += abs(Q[i] - Q_baseline[j][i]); + } + ASSERT_LE(total_error, tolerance); +} + +void run_Ptest(int k) { + Real dx = 1.0; + int m = k*2 + 1; + int j = k/2 - 1; + Real total_error = 0.0; + + WeightsP P(k,m,dx); + for (int i = 0; i < m; ++i) { + total_error += abs(P[i] - P_baseline[j][i]); + } + ASSERT_LE(total_error, tolerance); + +} + +TEST(WeightTests, Accuracy) { + + // Baseline weights - Generated from MatLab/Octave + Q_baseline[0] = { 1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000,1.000000000}; + Q_baseline[1] = { 1.000000000,1.125064293,0.751414447,1.162097097,0.962852897,0.997142531,0.962852897, + 1.162097097,0.751414447,1.125064293,1.000000000 }; + Q_baseline[2] = { 1.000000000,1.188528786,0.464036031,1.670433241,0.529401810,1.170898848,0.978586058, + 0.996230453,0.978586058,1.170898848,0.529401810,1.670433241,0.464036031,1.188528786,1.000000000 }; + + P_baseline[0] = { 0.375000000,1.125000000,1.000000000,1.000000000,1.125000000,0.375000000 }; + P_baseline[1] = { 0.354134518,1.228459404,0.898117099,1.018547095,1.000741884,1.000741884, + 1.018547095,0.898117099,1.228459404,0.354134518 }; + P_baseline[2] = { 0.315722926,1.390677390,0.629532984,1.234237262,0.919144328,1.009804636, + 1.000880474,1.000880474,1.009804636,0.919144328,1.234237262,0.629532984, + 1.390677390,0.315722926 }; + for (int k : {6}) { + run_Qtest(k); + run_Ptest(k); + } +}