diff --git a/.github/workflows/test_pvade.yaml b/.github/workflows/test_pvade.yaml new file mode 100644 index 00000000..3a1fd2d7 --- /dev/null +++ b/.github/workflows/test_pvade.yaml @@ -0,0 +1,38 @@ +name: test_pvade + +on: + push: + branches: [ main, dev ] + pull_request: + branches: [ main, dev ] + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + # Job 1 of 2 - run pytest on code + pytest: + name: Pytest + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest] + runs-on: ${{ matrix.os }} + steps: + - uses: actions/checkout@v3 + - name: Build Conda environment + uses: mamba-org/provision-with-micromamba@main + with: + environment-file: environment.yaml + - name: Run pytest + shell: bash -l {0} + run: pytest -sv pvade/tests/ + + # Job 2 of 2 - enforce Black formatting + formatting: + name: Black formatting + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: psf/black@stable diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..64852b99 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +.* +!/.gitignore +!/.readthedocs.yaml +*.pyc +*.egg-info/ + +docs/_build/ diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 00000000..b40a06d2 --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,14 @@ +# .readthedocs.yaml + +version: 2 + +build: + os: "ubuntu-20.04" + tools: + python: "mambaforge-4.10" + +sphinx: + configuration: docs/conf.py + +conda: + environment: environment.yaml diff --git a/README.md b/README.md index fe367098..c2e5a71f 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,15 @@ -# PVOPT -This project will develop an open source fluid-structure interaction model which can be used as part of a modeling chain to provide stressor inputs to mechanical module models to study the physics of failure for degradation mechanisms such as cell cracking, weathering of cracked cells, and glass breakage. +# PVade -https://github.nrel.gov/pages/warsalan/PVOPT/ +`PVade` is an open source fluid-structure interaction model which can be used to study wind loading and stability on solar-tracking PV arrays. `PVade` can be used as part of a larger modeling chain to provide stressor inputs to mechanical module models to study the physics of failure for degradation mechanisms such as cell cracking, weathering of cracked cells, and glass breakage. For more information, visit the [PVade Documentation](TODO). + +[![test_pvade](https://github.com/NREL/PVade/actions/workflows/test_pvade.yaml/badge.svg)](https://github.com/NREL/PVade/actions/workflows/test_pvade.yaml) + +## Getting Started + +New users are encouraged to review the [Getting Started](TODO) guide which describes how to create the Conda environment and run the example simulations. + +## Developer Quick Start -# Getting Started 1. To use this software, begin by creating a Conda environment using the provided `environment.yaml` file: ```bash conda env create -n my_env_name -f environment.yaml diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 00000000..d4bb2cbb --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/background/index.rst b/docs/background/index.rst new file mode 100644 index 00000000..d748404a --- /dev/null +++ b/docs/background/index.rst @@ -0,0 +1,20 @@ +Background +========== + +A collection of resources for further reading on PVade topics and the theory and implementation for solver elements. + +Fluid Solver +------------ + + + +Structural Solver +----------------- + + + +PV Wind Loading and Stability +----------------------------- + +* E. Young, "Predicting wind loading and instability in solar tracking PV arrays", PV Reliability Workshop (2023), https://www.nrel.gov/docs/fy23osti/85567.pdf. + diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 00000000..f946c59e --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,59 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys + +sys.path.insert(0, os.path.abspath("../pvade")) + + +# -- Project information ----------------------------------------------------- + +project = "PVade" +copyright = "2023, Walid Arsalane, Ethan Young" +author = "Walid Arsalane, Ethan Young" + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + "sphinx.ext.autodoc", + "sphinx.ext.napoleon", + "sphinx.ext.autosectionlabel", +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ["_templates"] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +# html_theme = 'alabaster' +html_theme = "sphinx_rtd_theme" + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ["_static"] + +autoclass_content = "both" diff --git a/docs/how_to_guides/getting_started.rst b/docs/how_to_guides/getting_started.rst new file mode 100644 index 00000000..31437782 --- /dev/null +++ b/docs/how_to_guides/getting_started.rst @@ -0,0 +1,23 @@ +Getting Started +=============== + +Building Conda Environment +-------------------------- + +To use this software, begin by creating a Conda environment using the provided ``environment.yaml`` file:: + + conda env create -n my_env_name -f environment.yaml + +where ``my_env_name`` can be replaced with a short name for your Conda environment. When the environment finishes installing, activate it with:: + + conda activate my_env_name + +from within your activate Conda environment, a simulation can be executed with:: + + python main.py --command_line_arg value + + +Running Examples +---------------- + +.. Fill in with walkthrough pointing to an example \ No newline at end of file diff --git a/docs/how_to_guides/hpc_jobs.rst b/docs/how_to_guides/hpc_jobs.rst new file mode 100644 index 00000000..aa8dda01 --- /dev/null +++ b/docs/how_to_guides/hpc_jobs.rst @@ -0,0 +1,4 @@ +HPC Jobs +======== + +.. Fill in with walkthrough pointing to an example \ No newline at end of file diff --git a/docs/how_to_guides/index.rst b/docs/how_to_guides/index.rst new file mode 100644 index 00000000..bc86b34f --- /dev/null +++ b/docs/how_to_guides/index.rst @@ -0,0 +1,8 @@ +How To Guides +============= + +.. toctree:: + :maxdepth: 1 + + getting_started + hpc_jobs diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 00000000..5cf63f81 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,30 @@ +.. PVade documentation master file, created by + sphinx-quickstart on Fri Mar 10 11:27:06 2023. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to PVade's documentation! +================================= + +PVade is an open source fluid-structure interaction model which can be used to study wind loading and stability on solar-tracking PV arrays. PVade can be used as part of a larger modeling chain to provide stressor inputs to mechanical module models to study the physics of failure for degradation mechanisms such as cell cracking, weathering of cracked cells, and glass breakage. + +Organization +------------ + +Documentation is currently organized into three main categories: + +* :ref:`How to Guides`: User guides covering basic topics and use cases for the PVade software +* :ref:`Technical Reference`: Programming details on the PVade API and functions +* :ref:`Background`: Information and research sources for fluid and structural solvers and PV topics + +New users may find it helpful to review the :ref:`Getting Started` materials first. + +Contents +-------- + +.. toctree:: + :maxdepth: 2 + + how_to_guides/index + technical_reference/index + background/index diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 00000000..32bb2452 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/technical_reference/data_stream.rst b/docs/technical_reference/data_stream.rst new file mode 100644 index 00000000..4fd32c77 --- /dev/null +++ b/docs/technical_reference/data_stream.rst @@ -0,0 +1,6 @@ +DataStream +========== + +.. automodule:: pvade.DataStream + :members: + :private-members: \ No newline at end of file diff --git a/docs/technical_reference/domain_creation.rst b/docs/technical_reference/domain_creation.rst new file mode 100644 index 00000000..c0849cb3 --- /dev/null +++ b/docs/technical_reference/domain_creation.rst @@ -0,0 +1,6 @@ +DomainCreation +============== + +.. automodule:: pvade.geometry.template.TemplateDomainCreation + :members: + :private-members: \ No newline at end of file diff --git a/docs/technical_reference/domain_creation_panels.rst b/docs/technical_reference/domain_creation_panels.rst new file mode 100644 index 00000000..8c585df9 --- /dev/null +++ b/docs/technical_reference/domain_creation_panels.rst @@ -0,0 +1,6 @@ +DomainCreationPanels +==================== + +.. automodule:: pvade.geometry.panels.DomainCreation + :members: + :private-members: \ No newline at end of file diff --git a/docs/technical_reference/flow_manager.rst b/docs/technical_reference/flow_manager.rst new file mode 100644 index 00000000..a9d0bb36 --- /dev/null +++ b/docs/technical_reference/flow_manager.rst @@ -0,0 +1,6 @@ +FlowManager +=========== + +.. automodule:: pvade.FlowManager + :members: + :private-members: \ No newline at end of file diff --git a/docs/technical_reference/index.rst b/docs/technical_reference/index.rst new file mode 100644 index 00000000..e42a7a56 --- /dev/null +++ b/docs/technical_reference/index.rst @@ -0,0 +1,14 @@ +Technical Reference +=================== + +.. toctree:: + :maxdepth: 1 + + parameters + data_stream + domain_creation_panels + domain_creation + mesh_manager + mesh_manager_2d + mesh_manager_3d + flow_manager diff --git a/docs/technical_reference/mesh_manager.rst b/docs/technical_reference/mesh_manager.rst new file mode 100644 index 00000000..5ad67595 --- /dev/null +++ b/docs/technical_reference/mesh_manager.rst @@ -0,0 +1,6 @@ +MeshManager +=========== + +.. automodule:: pvade.geometry.MeshManager + :members: + :private-members: \ No newline at end of file diff --git a/docs/technical_reference/mesh_manager_2d.rst b/docs/technical_reference/mesh_manager_2d.rst new file mode 100644 index 00000000..8b6b2ad1 --- /dev/null +++ b/docs/technical_reference/mesh_manager_2d.rst @@ -0,0 +1,6 @@ +MeshManager2D +============= + +.. automodule:: pvade.geometry.MeshManager2d + :members: + :private-members: \ No newline at end of file diff --git a/docs/technical_reference/mesh_manager_3d.rst b/docs/technical_reference/mesh_manager_3d.rst new file mode 100644 index 00000000..4bafc60a --- /dev/null +++ b/docs/technical_reference/mesh_manager_3d.rst @@ -0,0 +1,6 @@ +MeshManager3d +============= + +.. automodule:: pvade.geometry.MeshManager3d + :members: + :private-members: \ No newline at end of file diff --git a/docs/technical_reference/parameters.rst b/docs/technical_reference/parameters.rst new file mode 100644 index 00000000..72e57dd1 --- /dev/null +++ b/docs/technical_reference/parameters.rst @@ -0,0 +1,6 @@ +Parameters +========== + +.. automodule:: pvade.Parameters + :members: + :private-members: \ No newline at end of file diff --git a/environment.yaml b/environment.yaml index 03a0ce3c..dd4adad9 100644 --- a/environment.yaml +++ b/environment.yaml @@ -1,4 +1,4 @@ -name: pvopt-env +name: pvade channels: - conda-forge - defaults @@ -10,13 +10,16 @@ dependencies: - jsonschema - jupyter - matplotlib + - meshio - pandas - pytest - python-gmsh - pyyaml - scipy - snakeviz + - sphinx - tqdm - pip - pip: - - -e . # Installs the pvopt package + - sphinx-rtd-theme + - -e . # Installs the pvade package diff --git a/inputs/2d_cyld.yaml b/inputs/2d_cyld.yaml index d2e60d05..756e70b6 100644 --- a/inputs/2d_cyld.yaml +++ b/inputs/2d_cyld.yaml @@ -40,8 +40,9 @@ fluid: turbulence_model: none periodic: false # c_s: 0.325 - bc_ywall_max: noslip # slip noslip free - bc_ywall_min: noslip # slip noslip free + bc_y_max: noslip # slip noslip free + bc_y_min: noslip # slip noslip free + use_eddy_viscosity: false structure: youngs: 1.0e+09 # Boundary_conditions: diff --git a/inputs/3d_cyld.yaml b/inputs/3d_cyld.yaml index 2ee634c4..baf96a09 100644 --- a/inputs/3d_cyld.yaml +++ b/inputs/3d_cyld.yaml @@ -12,7 +12,7 @@ domain: y_max: 0.41 z_min: 0.0 z_max: 0.41 - l_char: 0.1 #.02 + l_char: 2.0 #.02 cyld_radius: 0.05 pv_array: num_rows: 1 # not used @@ -33,6 +33,9 @@ solver: solver1_pc: jacobi solver2_pc: jacobi solver3_pc: jacobi + save_text_interval: 0.1 + save_xdmf_interval: 0.1 + fluid: u_ref: 0.45 nu: 0.001 # Establish re = 20 with diam = 0.1 and u_bar = u_ref * (4/9) @@ -40,10 +43,11 @@ fluid: turbulence_model: none periodic: false # c_s: 0.325 - bc_ywall_max: noslip # slip noslip free - bc_ywall_min: noslip # slip noslip free - bc_zwall_max: noslip # slip noslip free - bc_zwall_min: noslip # slip noslip free + bc_y_max: noslip # slip noslip free + bc_y_min: noslip # slip noslip free + bc_z_max: noslip # slip noslip free + bc_z_min: noslip # slip noslip free + use_eddy_viscosity: True structure: youngs: 1.0e+09 # Boundary_conditions: diff --git a/inputs/sim_params_alt.yaml b/inputs/sim_params_alt.yaml index e199ee0d..80b389bc 100644 --- a/inputs/sim_params_alt.yaml +++ b/inputs/sim_params_alt.yaml @@ -1,5 +1,5 @@ general: - example: panels + example: panels3d output_dir_mesh: output/panels/mesh output_dir_sol: output/panels/solution create_mesh: True @@ -12,7 +12,7 @@ domain: y_max: 27 z_min: 0 z_max: 20 - l_char: 10 + l_char: 5 pv_array: num_rows: 7 elevation: 1.5 @@ -32,18 +32,19 @@ solver: solver1_pc: jacobi solver2_pc: jacobi solver3_pc: jacobi - save_text_interval: 0.1 - save_xdmf_interval: 0.1 + save_text_interval: 0.01 + save_xdmf_interval: 0.01 fluid: u_ref: 8.0 nu: 0.01 dpdx_val: 0.0 turbulence_model: none periodic: false - bc_ywall_max: slip # slip noslip free - bc_ywall_min: slip # slip noslip free - bc_zwall_max: slip # slip noslip free - bc_zwall_min: noslip # slip noslip free + bc_y_max: slip # slip noslip free + bc_y_min: slip # slip noslip free + bc_z_max: slip # slip noslip free + bc_z_min: noslip # slip noslip free + use_eddy_viscosity: True # c_s: 0.325 structure: youngs: 1.0e+09 diff --git a/inputs/sim_params_alt_2D.yaml b/inputs/sim_params_alt_2D.yaml index 0c3872ed..fffb2070 100644 --- a/inputs/sim_params_alt_2D.yaml +++ b/inputs/sim_params_alt_2D.yaml @@ -38,8 +38,8 @@ fluid: dpdx_val: 0.0 turbulence_model: none periodic: false - bc_ywall_max: slip # slip noslip free - bc_ywall_min: noslip # slip noslip free + bc_y_max: slip # slip noslip free + bc_y_min: noslip # slip noslip free # c_s: 0.325 structure: youngs: 1.0e+09 diff --git a/ns_main.py b/ns_main.py index adf81f66..c76a30fb 100644 --- a/ns_main.py +++ b/ns_main.py @@ -1,84 +1,40 @@ -# <<<<<<< HEAD -# ## from altFlowManager import * -# from pvopt.FlowManager_copy import * -# #from pvopt.IOManager import DataStream -# ======= - -from pvopt.FlowManager import Flow -from pvopt.DataStream import DataStream - -from pvopt.Parameters import SimParams +from pvade.fluid.FlowManager import Flow +from pvade.DataStream import DataStream +from pvade.Parameters import SimParams +from pvade.Utilities import get_input_file, write_metrics +from pvade.geometry.MeshManager import FSIDomain +from dolfinx.common import TimingType, list_timings import cProfile import sys -import numpy as np -from mpi4py import MPI -import time -from dolfinx.common import TimingType, list_timings import tqdm.autonotebook def main(): - # problems solver are: - # cylinder3d - - # cylinder2d - # panels - # problem = "cylinder2d" - - # problem = "panels" - problem = sys.argv[2] - - # use correct inpit file - if problem == "cylinder3d": - params = SimParams("inputs/3d_cyld.yaml") - dim = 3 - elif problem == "cylinder2d": - params = SimParams("inputs/2d_cyld.yaml") - dim = 2 - elif problem == "panels": - params = SimParams("inputs/sim_params_alt.yaml") - dim = 3 - elif problem == "panels2d": - params = SimParams("inputs/sim_params_alt_2D.yaml") - dim = 2 - else: - raise ValueError(f"Problem is not defined, please specify problem type using --problem $problem \n$problem: [cylinder3d,cylinder2d,panels,panels2d] ") - - if MPI.COMM_WORLD.rank == 0: - print("Currently Solving:",problem) - print("command used is: ",str(sys.argv)) - - if dim == 3: - from pvopt.geometry.MeshManager3d import FSIDomain - elif dim == 2: - from pvopt.geometry.MeshManager2d import FSIDomain - else: - print("dimension not defined") - exit() - + # Get the path to the input file from the command line + input_file = get_input_file() + # input_file = "inputs/2d_cyld.yaml" # get_input_file() + # Load the parameters object specified by the input file + params = SimParams(input_file) # Initialize the domain and construct the initial mesh domain = FSIDomain(params) - if params.general.create_mesh == True: - domain.build() - domain.write_mesh_file() + if params.general.create_mesh == True: + domain.build(params) + domain.write_mesh_file(params) elif params.general.read_mesh == True: domain.read() else: raise ValueError(f"Error in creating/loading mesh, please correct inputs") - # domain.test_mesh_functionspace() - if params.general.mesh_only == True: - # dolfinx.cpp.common.TimingType() - list_timings(MPI.COMM_WORLD, [TimingType.wall]) + list_timings(params.comm, [TimingType.wall]) exit() # Check to ensure mesh node matching for periodic simulations # if domain.periodic_simulation: - # domain.check_mesh_periodicity(mpi_info) + # domain.check_mesh_periodicity(params) # Initialize the function spaces for the flow flow = Flow(domain) @@ -91,17 +47,11 @@ def main(): dataIO = DataStream(domain, flow, params) - - # with XDMFFile(domain.comm, results_filename, "w") as xdmf: - # xdmf.write_mesh(domain.msh) - # xdmf.write_function(flow.u_k, 0.0) - # xdmf.write_function(flow.p_k, 0.0) - # # print(np.max(flow.u_k.x.array)) - t_steps = int(params.solver.t_final / params.solver.dt) - # save_vtk_every_n = int(params.solver.save_vtk_every / params.solver.dt) - tic = time.time() if domain.rank == 0: - progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=params.solver.t_steps) + progress = tqdm.autonotebook.tqdm( + desc="Solving PDE", total=params.solver.t_steps + ) + for k in range(params.solver.t_steps): if domain.rank == 0: progress.update(1) @@ -110,7 +60,7 @@ def main(): flow.solve(params) # adjust pressure to avoid dissipation of pressure profile - # flow.adjust_dpdx_for_constant_flux(mpi_info) + # flow.adjust_dpdx_for_constant_flux(params) if (k + 1) % params.solver.save_xdmf_interval_n == 0: if domain.rank == 0: print( @@ -118,61 +68,25 @@ def main(): ) dataIO.save_XDMF_files(flow, (k + 1) * params.solver.dt) - toc = time.time() - if domain.rank == 0: - print(f"Total solve time = {toc-tic:.2f} s.") - print(f"Average time per iteration = {(toc-tic)/t_steps:.2f} s.") - list_timings(MPI.COMM_WORLD, [TimingType.wall]) - -def write_metrics(): - if MPI.COMM_WORLD.Get_rank() == 0: - with open('profiling.txt', 'r') as output_file: - #solver_line = [line for line in output_file if "(solve)" in line] - solver_line = [line for line in output_file if "(solve)" in line] - print(solver_line) - solver_line = solver_line[0].split() - print("solve: total time = ",solver_line[3], " time per call = ",solver_line[4], " num calls = ",solver_line[0]) - - with open('profiling.txt', 'r') as output_file_1: - #solver_line = [line for line in output_file if "(solve)" in line] - solver1_line = [line for line in output_file_1 if "_solver_step_1" in line] - print(solver1_line) - solver1_line = solver1_line[0].split() - print("solver 1: total time = ",solver1_line[3], " time per call = ",solver1_line[4], " num calls = ",solver1_line[0]) - - with open('profiling.txt', 'r') as output_file: - solver2_line = [line for line in output_file if "_solver_step_2" in line] - print(solver2_line) - solver2_line = solver2_line[0].split() - print("solver 2: total time = ",solver2_line[3], " time per call = ",solver2_line[4], " num calls = ",solver2_line[0]) - - with open('profiling.txt', 'r') as output_file: - solver3_line = [line for line in output_file if "_solver_step_3" in line] - print(solver3_line) - solver3_line = solver3_line[0].split() - print("solver 3: total time = ",solver3_line[3], " time per call = ",solver3_line[4], " num calls = ",solver3_line[0]) - - with open('profiling.txt', 'r') as output_file: - meshread_line = [line for line in output_file if "(read)" in line] - if meshread_line: - meshread_line =meshread_line[0].split() - print("mesh read: total time = ",meshread_line[3], " time per call = ",meshread_line[4], " num calls = ",meshread_line[0]) - - with open('profiling.txt', 'r') as output_file: - meshbuild_line = [line for line in output_file if "(build)" in line] - if meshbuild_line: - meshbuild_line = meshbuild_line[0].split() - print("mesh build: total time = ",meshbuild_line[3], " time per call = ",meshbuild_line[4], " num calls = ",meshbuild_line[0]) + + list_timings(params.comm, [TimingType.wall]) + + return params + + # Print profiling results if __name__ == "__main__": profiler = cProfile.Profile() profiler.enable() - main() + params = main() profiler.disable() - profiler.dump_stats("profiling.prof") - with open("profiling.txt", "w") as output_file: - sys.stdout = output_file - profiler.print_stats(sort="cumtime") - sys.stdout = sys.__stdout__ + if params.rank == 0: + profiler.dump_stats("profiling.prof") + + with open("profiling.txt", "w") as output_file: + sys.stdout = output_file + profiler.print_stats(sort="cumtime") + sys.stdout = sys.__stdout__ + write_metrics() diff --git a/pvade/DataStream.py b/pvade/DataStream.py new file mode 100644 index 00000000..326676e2 --- /dev/null +++ b/pvade/DataStream.py @@ -0,0 +1,78 @@ +# from dolfinx import * +import numpy as np +import time +import os +import shutil +from dolfinx.io import XDMFFile, VTKFile +from mpi4py import MPI +from pathlib import Path +import pytest + + +# test actions +class DataStream: + + """Input/Output and file writing class + + This class manages the generation and updating of the saved xdmf files + from each run along with command line outputs and log files. + + Attributes: + comm (MPI Communicator): An MPI communicator used by all PVade objects + ndim (int): The number of dimensions in the problem + num_procs (int): The total number of processors being used to solve the problem + rank (int): A unique ID for each process in the range `[0, 1, ..., num_procs-1]` + results_filename (str): The full path to the directory in which all saved files will be written. + + """ + + def __init__(self, domain, flow, params): + """Initialize the DataStream object + + This initializes an object that manages the I/O for all of PVade. + + Args: + domain (:obj:`pvade.geometry.MeshManager.Domain`): A Domain object + flow (:obj:`pvade.fluid.FlowManager.Flow`): A Flow object + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + + """ + self.comm = params.comm + self.rank = params.rank + self.num_procs = params.num_procs + + self.ndim = domain.fluid.msh.topology.dim + + self.results_filename = f"{params.general.output_dir_sol}/solution.xdmf" + self.log_filename = f"{params.general.output_dir_sol}/log.txt" + + with XDMFFile(self.comm, self.results_filename, "w") as xdmf_file: + tt = 0.0 + xdmf_file.write_mesh(domain.fluid.msh) + xdmf_file.write_function(flow.u_k, 0.0) + xdmf_file.write_function(flow.p_k, 0.0) + + if self.rank == 0: + with open(self.log_filename, "w") as fp: + fp.write("Run Started.\n") + + def save_XDMF_files(self, flow, tt): + """Write additional timestep to XDMF file + + This function appends the state of the flow at time `tt` to an existing XDMF file. + + Args: + flow (:obj:`pvade.fluid.FlowManager.Flow`): A Flow object + tt (float): The time at which the current solution exists + + """ + with XDMFFile(self.comm, self.results_filename, "a") as xdmf_file: + xdmf_file.write_function(flow.u_k, tt) + xdmf_file.write_function(flow.p_k, tt) + + def print_and_log(self, string_to_print): + if self.rank == 0: + print(string_to_print) + + with open(self.log_filename, "a") as fp: + fp.write(f"{string_to_print}\n") diff --git a/pvade/Parameters.py b/pvade/Parameters.py new file mode 100644 index 00000000..3a8db190 --- /dev/null +++ b/pvade/Parameters.py @@ -0,0 +1,409 @@ +import os +import yaml +import argparse + +from mpi4py import MPI +from pandas import json_normalize +from jsonschema import validate + + +class SimParams: + """Manages and validates all simulation parameters + + ``SimParams`` provides a set of methods to obtain user-input values from + both input files and the command line. It also provides input validation + via a yaml schema file. All parameters can later be access from PVade + objects like ``params.parent_group.option`` for example, + ``params.solver.dt``. + + + Attributes: + comm (MPI Communicator): An MPI communicator that will be shared by all PVade objects. + flat_schema_dict (dict): A flattened dictionary version of the yaml schema file used to set default values and provide valid command-line options. + input_dict (dict): The aggregated input dictionary containing all command-line updated values and input-file values. + input_file_dict (dict): The dictionary representing the input file. + num_procs (int): The total number of processors being used to solve the problem + rank (int): A unique ID for each process in the range `[0, 1, ..., num_procs-1]` + schema_dict (dict): The unflattened (nested) representation of the yaml schema dictionary. + """ + + def __init__(self, input_file_path=None): + """Create a SimParams object + + This method begins by reading both the schema file and the + user-specified file into dictionaries. It overrides entries in the + user-specified dictionary with command-line inputs and finally + validates the complete dictionary against the schema file. + + Args: + input_file_path (str, optional): The complete path to the input yaml file defining the problem. + """ + # Get MPI communicators + self.comm = MPI.COMM_WORLD + self.rank = self.comm.Get_rank() + self.num_procs = self.comm.Get_size() + + # Open the schema file for reading and load its contents into a dictionary + pvade_dir = os.path.dirname(os.path.abspath(__file__)) + with open(os.path.join(pvade_dir, "input_schema.yaml"), "r") as fp: + self.schema_dict = yaml.safe_load(fp) + + # Flatten the schema dictionary to make navigation and CLI parsing easier + self._flatten_schema_dict() + + # Initialize the single dict that will hold and combine ALL inputs + self.input_dict = {} + + if input_file_path is not None: + # Assert that the file exists + assert os.path.isfile( + input_file_path + ), f"Could not find file '{input_file_path}', check it exists" + + # Open the input file for reading and load its contents into a dictionary + with open(input_file_path, "r") as fp: + if self.rank == 0: + print(f"Reading problem definition from {input_file_path}") + self.input_file_dict = yaml.safe_load(fp) + + # Set values as specified by the input yaml file + self._set_user_params_from_file() + + # Notes from WALID: + # - Print a warning if the value set from the file isn't in the schema + # - Print the values that don't fit + # - everything you expect to store in parameters should be in the yaml schema + + else: + # Assign all default parameters using the flattened schema + + # We can still support this, but maybe not the official getting + # started strategy(provide a yaml file for demos) + self._initialize_params_to_default() + + # Override any previous value with a command line input + self._set_user_params_from_cli() + + # Check that the complete parameter spec conforms to the schema + self._validate_inputs() + + # Store the nested dictionary as attributes on this object for easy + # access e.g., params.domain.x_max instead of params['domain'] + # ['x_max'] + self._store_dict_as_attrs() + + self._add_derived_quantities() + + def _flatten_schema_dict(self): + """Flatten the schema dictionary + + This method returns a flattened version of the schema dictionary where + nested dictionary key:val pairs are separated with dots (`.`). For + example:: + + general: + name: test_sim + size: + x: 10 + y: 20 + + will be converted to:: + + general.name: test_sim + general.size.x: 10 + general.size.y: 20 + + Args: + None + + """ + flat_schema_raw = json_normalize(self.schema_dict, sep=".").to_dict() + + self.flat_schema_dict = {} + + for key in flat_schema_raw.keys(): + if "default" in key: + root = key.split(".default")[0] + + short_key = root.split(".") + short_key = [sk for sk in short_key if sk != "properties"] + short_key = ".".join(short_key) + + self.flat_schema_dict[short_key] = {} + + for subkey in ["default", "type", "units", "description"]: + try: + val = flat_schema_raw[f"{root}.{subkey}"][0] + except: + val = None + + self.flat_schema_dict[short_key][subkey] = val + + def _set_user_params_from_file(self): + """Store the input-file supplied values + + This method stores the parameters supplied via input yaml file in the + aggregated input dictionary. + + Args: + None + + """ + flat_input_file_dict = json_normalize(self.input_file_dict, sep=".").to_dict() + + for key, val in flat_input_file_dict.items(): + path_to_input = key.split(".") + self._set_nested_dict_value( + self.input_dict, path_to_input, val[0], error_on_missing_key=False + ) + + def _initialize_params_to_default(self): + """Initialize values to default + + This method uses the `default` entry associated with each parameter in + the yaml schema to provide a sane value where none is supplied by the + input yaml file. + + Args: + None + + """ + for key, val in self.flat_schema_dict.items(): + path_to_input = key.split(".") + self._set_nested_dict_value( + self.input_dict, + path_to_input, + val["default"], + error_on_missing_key=False, + ) + + def _set_user_params_from_cli(self): + """Store values from command line + + This method uses the flattened schema dictionary, where nested options + are stored as a one-layer dictionary with keys separated by dots + (`.`), to generate a set of valid command line inputs that can be + used to override options that would otherwise be set by the input + yaml file or defaults. For example:: + + python main.py --solver.dt 0.01 + + will override whatever value of ``dt`` is present in the yaml file + with the value 0.01. Note the double hyphen to set off an argument + name, either with or without the assignment operator is allowed + (e.g., both ``--solver.dt 0.01`` and ``--solver.dt=0.01`` are valid) + + Args: + None + + """ + parser = argparse.ArgumentParser() + + ignore_list = ["input_file"] + + parser.add_argument( + "--input_file", + metavar="", + type=str, + help="The full path to the input file, e.g., 'intputs/my_params.yaml'", + ) + + for key, value in self.flat_schema_dict.items(): + help_message = f"{value['description']} (data type = {value['type']}, Units = {value['units']})" + + if value["type"] == "string": + cli_type = str + elif value["type"] == "number": + cli_type = float + elif value["type"] == "integer": + cli_type = int + else: + cli_type = None + + parser.add_argument( + f"--{key}", metavar="", type=cli_type, help=help_message + ) + + command_line_inputs, unknown = parser.parse_known_args() + + # Find any command line arguments that were used and replace those entries in params + for key, value in vars(command_line_inputs).items(): + if key not in ignore_list and value is not None: + path_to_input = key.split(".") + self._set_nested_dict_value( + self.input_dict, path_to_input, value, error_on_missing_key=False + ) + + if self.rank == 0: + print(f"| Setting {key} = {value} from command line.") + + for key in unknown: + if self.rank == 0: + print(f"| Got unknown option {key}, skipping.") + + def _validate_inputs(self): + """Validate the input dictionary + + This method uses the validate function from the jsonschema package to + validate the complete user-supplied parameter set. This checks that + variables are the expected type, fall within prescribed ranges, and + also ensure that dependencies between variables are satisfied + (e.g., if variable ``a`` is specified, then variable ``b`` must also + be specified.) + + Args: + None + + """ + # This compares the input dictionary against the yaml schema + # to ensure all values are set correctly + validate(self.input_dict, self.schema_dict) + + # self._pprint_dict(self.input_dict) + + def _store_dict_as_attrs(self): + """Converts dictionary notation to nested object attributes + + This method re-structures the parameters dictionary into the slightly + easier-to-navigate ``object.attribute`` format. For example, rather + than having a parameters dictionary and accessing entries like + ``params['general']['name']``, we create a parameters object with + attributes accessed like ``params.general.name``. + + Args: + None + + """ + flat_input_dict = json_normalize(self.input_dict, sep=".").to_dict() + + for key, val in flat_input_dict.items(): + path_to_input = key.split(".") + self._rec_settattr(self, path_to_input, val[0], error_on_missing_key=False) + + def _set_nested_dict_value( + self, parent_obj, path_to_input, value, error_on_missing_key=True + ): + """Set a nested value in a dictionary + + This method uses the flattened path to a dictionary value to create + the nested structure/path and then store the value. + + Args: + parent_obj (dict): The dictionary containing the value to be set + path_to_input (list, str): A list of strings containing the key names to be used at each level + value (int, float, object): The value to store at the final level of nesting + error_on_missing_key (bool, optional): A flag to raise an error if any of the keys or sub-keys is not already present + + """ + key = path_to_input[-1] + path_to_input = path_to_input[:-1] + + for p in path_to_input: + if error_on_missing_key: + assert ( + p in parent_obj + ), f"Could not find option '{p}' in set of valid inputs." + parent_obj = parent_obj.setdefault(p, {}) + + if error_on_missing_key: + assert ( + key in parent_obj + ), f"Could not find option '{key}' in set of valid inputs." + parent_obj[key] = value + + def _rec_settattr( + self, parent_obj, path_to_input, value, error_on_missing_key=True + ): + """Convert nested dictionary to object.attribute format + + This method recursively creates nested objects to hold attributes in the same manner as the nested dictionary. + + Args: + parent_obj (object): The parent object that will contain this value (self when first called) + path_to_input (list, str): A list of strings containing the attribute names to be used at each level + value (int, float, object): The value to store at the final level of nesting + error_on_missing_key (bool, optional): A flag to raise an error if any of the keys or sub-keys is not already present + + """ + + class ParamGroup: + + """Empty class to enable nesting and attribute storage""" + + pass + + if len(path_to_input) == 1: + if error_on_missing_key: + assert hasattr( + parent_obj, path_to_input[0] + ), f"Attribute '{path_to_input[0]}' is invalid." + + setattr(parent_obj, path_to_input[0], value) + + else: + if error_on_missing_key: + child_obj = getattr(parent_obj, path_to_input[0]) + else: + child_obj = getattr(parent_obj, path_to_input[0], ParamGroup()) + + setattr(parent_obj, path_to_input[0], child_obj) + parent_obj = child_obj + self._rec_settattr( + parent_obj, + path_to_input[1:], + value, + error_on_missing_key=error_on_missing_key, + ) + + def _pprint_dict(self, d, indent=0): + """Pretty print a dictionary + + Print a dictionary where nesting is visualized with indentation. + Useful for validation and debugging. + + Args: + d (dict): The dictionary to print + indent (int, optional): The amount of indentation to use for this level, default=0 + + """ + for key, val in d.items(): + for k in range(indent): + print("| ", end="") + + if isinstance(val, dict): + print(f"{key}:") + indent += 1 + self._pprint_dict(val, indent) + indent -= 1 + + else: + print(f"{key}: {val} ({type(val)})") + + def _add_derived_quantities(self): + """Add derived quantities for convenience + + For convenience, it is helpful to store certain derived quantities in + the SimParams object that are not explicitly provided by the user. + For example, ``params.solver.t_steps`` is the integer number of steps + implied by the user-provided entries for ``t_final`` and ``dt``. Such + derived and/or convenience quantities are accumulated here. Note that + no error checking is performed on these values since they do not + appear in the yaml schema file, the expectation being that valid + quantities up till this point will result in valid derived values. + + Args: + None + + """ + self.solver.t_steps = int(self.solver.t_final / self.solver.dt) + self.solver.save_xdmf_interval_n = int( + self.solver.save_xdmf_interval / self.solver.dt + ) + self.solver.save_text_interval_n = int( + self.solver.save_text_interval / self.solver.dt + ) + + if hasattr(self.domain, "z_min"): + self.domain.dim = 3 + else: + self.domain.dim = 2 diff --git a/pvade/Utilities.py b/pvade/Utilities.py new file mode 100644 index 00000000..a21b654a --- /dev/null +++ b/pvade/Utilities.py @@ -0,0 +1,105 @@ +import argparse + + +def get_input_file(): + parser = argparse.ArgumentParser() + + parser.add_argument( + "--input_file", + metavar="", + type=str, + help="The full path to the input file, e.g., 'intputs/my_params.yaml'", + ) + + command_line_inputs, unknown = parser.parse_known_args() + + try: + input_file_path = vars(command_line_inputs)["input_file"] + assert input_file_path is not None + + except: + raise ValueError("No input file specified, quitting.") + + return input_file_path + + +def write_metrics(): + with open("profiling.txt", "r") as output_file: + # solver_line = [line for line in output_file if "(solve)" in line] + solver_line = [line for line in output_file if "(solve)" in line] + print(solver_line) + solver_line = solver_line[0].split() + print( + "solve: total time = ", + solver_line[3], + " time per call = ", + solver_line[4], + " num calls = ", + solver_line[0], + ) + + with open("profiling.txt", "r") as output_file_1: + # solver_line = [line for line in output_file if "(solve)" in line] + solver1_line = [line for line in output_file_1 if "_solver_step_1" in line] + print(solver1_line) + solver1_line = solver1_line[0].split() + print( + "solver 1: total time = ", + solver1_line[3], + " time per call = ", + solver1_line[4], + " num calls = ", + solver1_line[0], + ) + + with open("profiling.txt", "r") as output_file: + solver2_line = [line for line in output_file if "_solver_step_2" in line] + print(solver2_line) + solver2_line = solver2_line[0].split() + print( + "solver 2: total time = ", + solver2_line[3], + " time per call = ", + solver2_line[4], + " num calls = ", + solver2_line[0], + ) + + with open("profiling.txt", "r") as output_file: + solver3_line = [line for line in output_file if "_solver_step_3" in line] + print(solver3_line) + solver3_line = solver3_line[0].split() + print( + "solver 3: total time = ", + solver3_line[3], + " time per call = ", + solver3_line[4], + " num calls = ", + solver3_line[0], + ) + + with open("profiling.txt", "r") as output_file: + meshread_line = [line for line in output_file if "(read)" in line] + if meshread_line: + meshread_line = meshread_line[0].split() + print( + "mesh read: total time = ", + meshread_line[3], + " time per call = ", + meshread_line[4], + " num calls = ", + meshread_line[0], + ) + + with open("profiling.txt", "r") as output_file: + meshbuild_line = [line for line in output_file if "(build)" in line] + if meshbuild_line: + meshbuild_line = meshbuild_line[0].split() + print( + "mesh build: total time = ", + meshbuild_line[3], + " time per call = ", + meshbuild_line[4], + " num calls = ", + meshbuild_line[0], + ) diff --git a/pvopt/__init__.py b/pvade/__init__.py similarity index 64% rename from pvopt/__init__.py rename to pvade/__init__.py index e11e63a8..ebf194a8 100644 --- a/pvopt/__init__.py +++ b/pvade/__init__.py @@ -1,8 +1,8 @@ -"""This is the base package for ``pvopt``. +"""This is the base package for ``pvade``. """ -class PvoptError(Exception): +class PVadeError(Exception): """The base error for the package. All custom exceptions raised will derive from this exception. """ diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py new file mode 100644 index 00000000..deeefe25 --- /dev/null +++ b/pvade/fluid/FlowManager.py @@ -0,0 +1,568 @@ +"""Summary +""" +import dolfinx +import ufl + +from petsc4py import PETSc +from mpi4py import MPI + +import numpy as np +import scipy.interpolate as interp + +import warnings + +from pvade.fluid.boundary_conditions import ( + build_velocity_boundary_conditions, + build_pressure_boundary_conditions, +) + + +class Flow: + """This class solves the CFD problem""" + + def __init__(self, domain): + """Initialize the fluid solver + + This method initialize the Flow object, namely, it creates all the + necessary function spaces on the mesh, initializes key counting and + boolean variables and records certain characteristic quantities like + the minimum cell size and the number of degrees of freedom attributed + to both the pressure and velocity function spaces. + + Args: + domain (:obj:`pvade.geometry.MeshManager.Domain`): A Domain object + + """ + # Store the comm and mpi info for convenience + self.comm = domain.comm + self.rank = domain.rank + self.num_procs = domain.num_procs + + # Pressure (Scalar) + P1 = ufl.FiniteElement("Lagrange", domain.fluid.msh.ufl_cell(), 1) + self.Q = dolfinx.fem.FunctionSpace(domain.fluid.msh, P1) + + # Velocity (Vector) + P2 = ufl.VectorElement("Lagrange", domain.fluid.msh.ufl_cell(), 2) + self.V = dolfinx.fem.FunctionSpace(domain.fluid.msh, P2) + + # Stress (Tensor) + P3 = ufl.TensorElement("Lagrange", domain.fluid.msh.ufl_cell(), 2) + self.T = dolfinx.fem.FunctionSpace(domain.fluid.msh, P3) + + P4 = ufl.FiniteElement("DG", domain.fluid.msh.ufl_cell(), 0) + + self.DG = dolfinx.fem.FunctionSpace(domain.fluid.msh, P4) + + self.first_call_to_solver = True + self.first_call_to_surface_pressure = True + + # Store the dimension of the problem for convenience + self.ndim = domain.fluid.msh.topology.dim + self.facet_dim = self.ndim - 1 + + # find hmin in mesh + num_cells = domain.fluid.msh.topology.index_map(self.ndim).size_local + h = dolfinx.cpp.mesh.h(domain.fluid.msh, self.ndim, range(num_cells)) + + # This value of hmin is local to the mesh portion owned by the process + hmin_local = np.amin(h) + + # collect the minimum hmin from all processes + self.hmin = np.zeros(1) + self.comm.Allreduce(hmin_local, self.hmin, op=MPI.MIN) + self.hmin = self.hmin[0] + + self.num_Q_dofs = ( + self.Q.dofmap.index_map_bs * self.Q.dofmap.index_map.size_global + ) + self.num_V_dofs = ( + self.V.dofmap.index_map_bs * self.V.dofmap.index_map.size_global + ) + + if self.rank == 0: + print(f"hmin = {self.hmin}") + print(f"Total num dofs = {self.num_Q_dofs + self.num_V_dofs}") + + def build_boundary_conditions(self, domain, params): + """Build the boundary conditions + + A method to manage the building of boundary conditions, including the + steps of identifying entities on the boundary, marking those degrees + of freedom either by the identified facets or a gmsh marker function, + and finally assembling a list of Boundary objects that enforce the + correct value. + + Args: + domain (:obj:`pvade.geometry.MeshManager.Domain`): A Domain object + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + + """ + + self.bcu, self.inflow_profile = build_velocity_boundary_conditions( + domain, params, self.V + ) + + self.bcp = build_pressure_boundary_conditions(domain, params, self.Q) + + def build_forms(self, domain, params): + """Builds all variational statements + + This method creates all the functions, expressions, and variational + forms that will be needed for the numerical solution of Navier Stokes + using a fractional step method. This includes the calculation of a + tentative velocity, the calculation of the change in pressure + required to correct the tentative velocity to enforce continuity, and + the update to the velocity field to reflect this change in + pressure. + + Args: + domain (:obj:`pvade.geometry.MeshManager.Domain`): A Domain object + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + + """ + # Define fluid properties + self.dpdx = dolfinx.fem.Constant(domain.fluid.msh, (0.0, 0.0, 0.0)) + self.dt_c = dolfinx.fem.Constant(domain.fluid.msh, (params.solver.dt)) + nu = dolfinx.fem.Constant(domain.fluid.msh, (params.fluid.nu)) + + # Define trial and test functions for velocity + self.u = ufl.TrialFunction(self.V) + self.v = ufl.TestFunction(self.V) + + # Define trial and test functions for pressure + self.p = ufl.TrialFunction(self.Q) + self.q = ufl.TestFunction(self.Q) + + # Define functions for velocity solutions + self.u_k = dolfinx.fem.Function(self.V, name="velocity") + self.u_k1 = dolfinx.fem.Function(self.V, name="velocity") + self.u_k2 = dolfinx.fem.Function(self.V) + + # Define functions for pressure solutions + self.p_k = dolfinx.fem.Function(self.Q, name="pressure") + self.p_k1 = dolfinx.fem.Function(self.Q, name="pressure") + + initialize_flow = True + + if initialize_flow: + self.u_k1.x.array[:] = self.inflow_profile.x.array + self.u_k2.x.array[:] = self.inflow_profile.x.array + self.u_k.x.array[:] = self.inflow_profile.x.array + + # Define expressions used in variational forms + # Crank-Nicolson velocity + U_CN = 0.5 * (self.u + self.u_k1) + + # Adams-Bashforth velocity + U_AB = 1.5 * self.u_k1 - 0.5 * self.u_k2 + + use_eddy_viscosity = params.fluid.use_eddy_viscosity + + if use_eddy_viscosity: + # By default, don't use any eddy viscosity + filter_scale = ufl.CellVolume(domain.fluid.msh) ** ( + 1.0 / domain.fluid.msh.topology.dim + ) + + # Strain rate tensor, 0.5*(du_i/dx_j + du_j/dx_i) + Sij = ufl.sym(ufl.nabla_grad(U_AB)) + + # ufl.sqrt(Sij*Sij) + strainMag = (2.0 * ufl.inner(Sij, Sij)) ** 0.5 + + # Smagorinsky dolfinx.fem.constant, typically close to 0.17 + Cs = 0.17 + + # Eddy viscosity + self.nu_T = Cs**2 * filter_scale**2 * strainMag + + else: + self.nu_T = dolfinx.fem.Constant(domain.fluid.msh, 0.0) + + # ================================================================# + # DEFINE VARIATIONAL FORMS + # ================================================================# + U = 0.5 * (self.u_k1 + self.u) + + def epsilon(u): + """Convenience expression for ufl.sym(ufl.nabla_grad(u)) + + Args: + u (dolfinx.fem.Function): A dolfinx function + + Returns: + ufl.dolfinx.fem.form: ufl.sym(ufl.nabla_grad(u)) + """ + return ufl.sym(ufl.nabla_grad(u)) + + # Define stress tensor + def sigma(u, p, nu): + """Convenience expression for fluid stress, sigma + + Args: + u (dolfinx.fem.Function): Velocity + p (dolfinx.fem.Function): Pressure + nu (float, dolfinx.fem.Function): Viscosity + + Returns: + ufl.dolfinx.fem.form: Stress in fluid, $2\nu \epsilon (u)$ + """ + return 2 * nu * epsilon(u) - p * ufl.Identity(len(u)) + + fractional_step_scheme = "IPCS" + U = 0.5 * (self.u_k1 + self.u) + n = ufl.FacetNormal(domain.fluid.msh) + f = dolfinx.fem.Constant( + domain.fluid.msh, + (PETSc.ScalarType(0), PETSc.ScalarType(0), PETSc.ScalarType(0)), + ) + # Define variational problem for step 1: tentative velocity + self.F1 = ( + (1.0 / self.dt_c) * ufl.inner(self.u - self.u_k1, self.v) * ufl.dx + + ufl.inner(ufl.dot(U_AB, ufl.nabla_grad(U_CN)), self.v) * ufl.dx + + (nu + self.nu_T) * ufl.inner(ufl.grad(U_CN), ufl.grad(self.v)) * ufl.dx + + ufl.inner(ufl.grad(self.p_k1), self.v) * ufl.dx + - ufl.inner(self.dpdx, self.v) * ufl.dx + ) + + self.a1 = dolfinx.fem.form(ufl.lhs(self.F1)) + self.L1 = dolfinx.fem.form(ufl.rhs(self.F1)) + + # Define variational problem for step 2: pressure correction + self.a2 = dolfinx.fem.form( + ufl.dot(ufl.nabla_grad(self.p), ufl.nabla_grad(self.q)) * ufl.dx + ) + self.L2 = dolfinx.fem.form( + ufl.dot(ufl.nabla_grad(self.p_k1), ufl.nabla_grad(self.q)) * ufl.dx + - (1.0 / self.dt_c) * ufl.div(self.u_k) * self.q * ufl.dx + ) + + # Define variational problem for step 3: velocity update + self.a3 = dolfinx.fem.form(ufl.dot(self.u, self.v) * ufl.dx) + self.L3 = dolfinx.fem.form( + ufl.dot(self.u_k, self.v) * ufl.dx + - self.dt_c * ufl.dot(ufl.nabla_grad(self.p_k - self.p_k1), self.v) * ufl.dx + ) + + # Define a function and the dolfinx.fem.form of the stress + self.panel_stress = dolfinx.fem.Function(self.T) + self.stress = sigma(self.u_k, self.p_k, nu + self.nu_T) + + # Define mesh normals + # self.n = FacetNormal(domain.mesh) + + # Compute traction vector + # self.traction = ufl.dot(self.stress, self.n) + + # Create a dolfinx.fem.form for projecting stress onto a tensor function space + # e.g., panel_stress.assign(project(stress, T)) + # self.a4 = ufl.inner(ufl.TrialFunction(self.T), ufl.TestFunction(self.T))*ufl.dx + # self.L4 = ufl.inner(self.stress, ufl.TestFunction(self.T))*ufl.dx + + # Create a dolfinx.fem.form for projecting CFL calculation onto DG function space + cell_diam = ufl.CellDiameter(domain.fluid.msh) + cfl_form = ufl.sqrt(ufl.inner(self.u_k, self.u_k)) * self.dt_c / cell_diam + + self.cfl_vec = dolfinx.fem.Function(self.DG) + self.a5 = dolfinx.fem.form( + ufl.inner(ufl.TrialFunction(self.DG), ufl.TestFunction(self.DG)) * ufl.dx + ) + self.L5 = dolfinx.fem.form( + ufl.inner(cfl_form, ufl.TestFunction(self.DG)) * ufl.dx + ) + + # Set up the functions to ensure a dolfinx.fem.constant flux through the outlet + outlet_cells = dolfinx.mesh.locate_entities( + domain.fluid.msh, self.ndim, lambda x: x[0] > params.domain.x_max - 1 + ) + self.flux_plane = dolfinx.fem.Function(self.V) + self.flux_plane.interpolate( + lambda x: (np.ones(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])), + outlet_cells, + ) + + # self.flux_plane = Expression('x[0] < cutoff ? 0.0 : 1.0', cutoff=domain.x_range[1]-1.0, degree=1) + self.flux_dx = ufl.Measure("dx", domain=domain.fluid.msh) # not needed ? + + # self.vol = dolfinx.fem.petsc.assemble_matrix(self.flux_plane*self.flux_dx) + # form1 = dolfinx.fem.form(self.flux_plane*self.flux_dx) + # self.vol = dolfinx.fem.petsc.assemble_vector(form1) + # self.J_initial = float(fem.petsc.assemble_matrix(self.flux_plane*self.u_k[0]*self.flux_dx)/self.vol) + + # if mpi_info['rank'] == 0: + # print('J_initial = %f' % (self.J_initial)) + + # self.J_history = [self.J_initial] + self.dpdx_history = [0.0] + + def _assemble_system(self, params): + """Pre-assemble all LHS matrices and RHS vectors + + Here we pre-assemble all the forms corresponding to the left-hand side + matrices and right-hand side vectors once outside the time loop. This + will enable us to re-use certain features like the sparsity pattern + during the timestepping without any modification of the function + calls. + + Args: + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + + self.A1 = dolfinx.fem.petsc.assemble_matrix(self.a1, bcs=self.bcu) + self.A2 = dolfinx.fem.petsc.assemble_matrix(self.a2, bcs=self.bcp) + self.A3 = dolfinx.fem.petsc.assemble_matrix(self.a3) + self.A5 = dolfinx.fem.petsc.assemble_matrix(self.a5) + + self.A1.assemble() + self.A2.assemble() + self.A3.assemble() + self.A5.assemble() + + self.b1 = dolfinx.fem.petsc.assemble_vector(self.L1) + self.b2 = dolfinx.fem.petsc.assemble_vector(self.L2) + self.b3 = dolfinx.fem.petsc.assemble_vector(self.L3) + self.b5 = dolfinx.fem.petsc.assemble_vector(self.L5) + + self.solver_1 = PETSc.KSP().create(self.comm) + self.solver_1.setOperators(self.A1) + self.solver_1.setType(params.solver.solver1_ksp) + self.solver_1.getPC().setType(params.solver.solver1_pc) + self.solver_1.setFromOptions() + + self.solver_2 = PETSc.KSP().create(self.comm) + self.solver_2.setOperators(self.A2) + self.solver_2.setType(params.solver.solver2_ksp) + self.solver_2.getPC().setType(params.solver.solver2_pc) + self.solver_2.setFromOptions() + + self.solver_3 = PETSc.KSP().create(self.comm) + self.solver_3.setOperators(self.A3) + self.solver_3.setType(params.solver.solver3_ksp) + self.solver_3.getPC().setType(params.solver.solver2_pc) + self.solver_3.setFromOptions() + + self.solver_5 = PETSc.KSP().create(self.comm) + self.solver_5.setOperators(self.A5) + self.solver_5.setType("cg") + self.solver_5.getPC().setType("jacobi") + self.solver_5.setFromOptions() + + def solve(self, params): + """Solve for a single timestep advancement + + Here we perform the three-step solution process (tentative velocity, + pressure correction, velocity update) to advance the fluid simulation + a single timestep. Additionally, we calculate the new CFL number + associated with the latest velocity solution. + + Args: + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + if self.first_call_to_solver: + if self.rank == 0: + print("Starting Fluid Solution") + + self._assemble_system(params) + + # Calculate the tentative velocity + self._solver_step_1(params) + + # Calculate the change in pressure to enforce incompressibility + self._solver_step_2(params) + + # Update the velocity according to the pressure field + self._solver_step_3(params) + + # Compute the CFL number + self.compute_cfl() + + # Update new -> old variables + self.u_k2.x.array[:] = self.u_k1.x.array + self.u_k1.x.array[:] = self.u_k.x.array + self.p_k1.x.array[:] = self.p_k.x.array + + if self.first_call_to_solver: + self.first_call_to_solver = False + + def _solver_step_1(self, params): + """Solve step 1: tentative velocity + + Here we calculate the tentative velocity which, not guaranteed to be + divergence free. + + Args: + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + # Step 0: Re-assemble A1 since using an implicit convective term + self.A1.zeroEntries() + self.A1 = dolfinx.fem.petsc.assemble_matrix(self.A1, self.a1, bcs=self.bcu) + self.A1.assemble() + self.solver_1.setOperators(self.A1) + + # Step 1: Tentative velocity step + with self.b1.localForm() as loc: + loc.set(0) + + self.b1 = dolfinx.fem.petsc.assemble_vector(self.b1, self.L1) + + dolfinx.fem.petsc.apply_lifting(self.b1, [self.a1], [self.bcu]) + + self.b1.ghostUpdate( + addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE + ) + + dolfinx.fem.petsc.set_bc(self.b1, self.bcu) + + self.solver_1.solve(self.b1, self.u_k.vector) + self.u_k.x.scatter_forward() + + def _solver_step_2(self, params): + """Solve step 2: pressure correction + + Here we calculate the pressure field that would be required to correct + the tentative velocity such that it is divergence free. + + Args: + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + # Step 0: Re-assemble A1 since using an implicit convective term + # Step 2: Pressure correction step + with self.b2.localForm() as loc: + loc.set(0) + + self.b2 = dolfinx.fem.petsc.assemble_vector(self.b2, self.L2) + + dolfinx.fem.petsc.apply_lifting(self.b2, [self.a2], [self.bcp]) + + self.b2.ghostUpdate( + addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE + ) + + dolfinx.fem.petsc.set_bc(self.b2, self.bcp) + + self.solver_2.solve(self.b2, self.p_k.vector) + self.p_k.x.scatter_forward() + + def _solver_step_3(self, params): + """Solve step 3: velocity update + + Here we update the tentative velocity with the effect of the modified, + continuity-enforcing pressure field. + + Args: + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + # Step 3: Velocity correction step + with self.b3.localForm() as loc: + loc.set(0) + + self.b3 = dolfinx.fem.petsc.assemble_vector(self.b3, self.L3) + + self.b3.ghostUpdate( + addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE + ) + + self.solver_3.solve(self.b3, self.u_k.vector) + self.u_k.x.scatter_forward() + + def compute_cfl(self): + """Solve for the CFL number + + Using the velocity, timestep size, and cell sizes, we calculate a CFL + number at every mesh cell. From that, we select the single highest + value of CFL number and record it for the purposes of monitoring + simulation stability. + + Args: + None + """ + with self.b5.localForm() as loc: + loc.set(0) + + self.b5 = dolfinx.fem.petsc.assemble_vector(self.b5, self.L5) + + self.b5.ghostUpdate( + addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE + ) + + self.solver_5.solve(self.b5, self.cfl_vec.vector) + self.cfl_vec.x.scatter_forward() + + cfl_max_local = np.amax(self.cfl_vec.vector.array) + + # collect the minimum hmin from all processes + self.cfl_max = np.zeros(1) + self.comm.Allreduce(cfl_max_local, self.cfl_max, op=MPI.MAX) + self.cfl_max = self.cfl_max[0] + + def adjust_dpdx_for_constant_flux(self): + """Adjust the forcing term, ``dpdx``, to maintain flowrate + + Here we calculate what the value of the dolfinx.fem.constant driving + pressure/force, ``dpdx``, should be adjusted to to maintain a + dolfinx.fem.constant flux through the domain. This is a useful option if + performing a periodic simulation in which the flux will slowly decay + due to viscous dissipation. The amount of change in ``dpdx`` is + calculated by comparing the current flux to the flux measured in the + initial condition and then employing a PID controller to adjust the + driving force to seek the target defined by the initial condition. + + Args: + None + """ + + def pid_controller(J_init, J_history, dt): + """Summary + + Args: + J_init (TYPE): Description + J_history (TYPE): Description + dt (TYPE): Description + + Returns: + TYPE: Description + """ + assert type(J_history) is list + + # K_p = 0.1 + # K_i = 0.4 + # K_d = 0.001 + # K_p = 0.4 + # K_i = 0.4 + # K_d = 0.01 + K_p = 0.6 + K_i = 0.4 + K_d = 0.05 + + err = J_init - np.array(J_history) + + c_p = K_p * err[-1] + c_i = K_i * dt * np.sum(err) + + try: + c_d = K_d * (err[-1] - err[-2]) / dt + except: + c_d = 0.0 + + output = c_p + c_i + c_d + + return output + + # Compute the flow through the outflow flux plane + # self.A1 = dolfinx.fem.petsc.assemble_matrix(self.a1 + J = float( + dolfinx.fem.petsc.assemble_matrix( + self.flux_plane * self.u_k[0] * self.flux_dx + ) + / self.vol + ) + self.J_history.append(J) + + dpdx_val = pid_controller(self.J_initial, self.J_history, float(self.dt_c)) + self.dpdx_history.append(dpdx_val) + + if dpdx_val < 0.0: + print("WARNING: dpdx_val = %f" % (dpdx_val)) + + self.dpdx.assign(dolfinx.fem.Constant((dpdx_val, 0.0, 0.0))) diff --git a/pvopt/geometry/__init__.py b/pvade/fluid/__init__.py similarity index 100% rename from pvopt/geometry/__init__.py rename to pvade/fluid/__init__.py diff --git a/pvade/fluid/boundary_conditions.py b/pvade/fluid/boundary_conditions.py new file mode 100644 index 00000000..4b9b7b6a --- /dev/null +++ b/pvade/fluid/boundary_conditions.py @@ -0,0 +1,267 @@ +import dolfinx +from petsc4py import PETSc + +import numpy as np + +import warnings + + +def get_facet_dofs_by_gmsh_tag(domain, functionspace, location): + """Associate degrees of freedom with marker functions + + This function uses the marker information in the gmsh specification to + find the corresponding degrees of freedom for use in the actual + boundary condition specification. Note that this method does not + require access to the facet information computed with + ``_locate_boundary_entities``. + + Args: + domain (:obj:`pvade.geometry.MeshManager.Domain`): A Domain object + functionspace (dolfinx.fem.FunctionSpace): The dolfinx function space on which the boundary condition will be applied + location (str, int): A string identifier of the boundary condition location, e.g., "x_min" or the equivalent integer tag, e.g., 1 + + Returns: + list: A list of the degrees of freedom associated with this tag and functionspace + + """ + + facet_dim = domain.ndim - 1 + + if isinstance(location, str): + found_entities = domain.fluid.facet_tags.find( + domain.domain_markers[location]["idx"] + ) + elif isinstance(location, int): + found_entities = domain.fluid.facet_tags.find(location) + + # if len(found_entities) == 0: + # warnings.warn(f"Found no facets using location = {location}.") + + dofs = dolfinx.fem.locate_dofs_topological(functionspace, facet_dim, found_entities) + + return dofs + + +def build_vel_bc_by_type(bc_type, domain, functionspace, bc_location): + """Build a dirichlet bc + + Args: + bc_type (str): The type of simple boundary condition to apply, "slip" or "noslip" + domain (:obj:`pvade.geometry.MeshManager.Domain`): A Domain object + functionspace (dolfinx.fem.FunctionSpace): The dolfinx function space on which the boundary condition will be applied + bc_location (str, int): A string identifier of the boundary condition location, e.g., "x_min" or the equivalent integer tag, e.g., 1 + + Returns: + dolfinx.fem.dirichletbc: A dolfinx dirichlet boundary condition + """ + + if domain.rank == 0: + print(f"Setting '{bc_type}' BC on {bc_location}") + + if bc_type == "noslip": + zero_vec = dolfinx.fem.Constant( + domain.fluid.msh, PETSc.ScalarType((0.0, 0.0, 0.0)) + ) + + dofs = get_facet_dofs_by_gmsh_tag(domain, functionspace, bc_location) + + bc = dolfinx.fem.dirichletbc(zero_vec, dofs, functionspace) + + elif bc_type == "slip": + zero_scalar = dolfinx.fem.Constant(domain.fluid.msh, PETSc.ScalarType(0.0)) + + if bc_location in ["x_min", "x_max"]: + sub_id = 0 + elif bc_location in ["y_min", "y_max"]: + sub_id = 1 + elif bc_location in ["z_min", "z_max"]: + sub_id = 2 + + dofs = get_facet_dofs_by_gmsh_tag( + domain, functionspace.sub(sub_id), bc_location + ) + + bc = dolfinx.fem.dirichletbc(zero_scalar, dofs, functionspace.sub(sub_id)) + + else: + if domain.rank == 0: + raise ValueError(f"{bc_type} BC not recognized") + + return bc + + +class InflowVelocity: + def __init__(self, geom_dim, params): + """Inflow velocity object + + Args: + geom_dim (int): The geometric dimension (as opposed to the topological dimension, usually this is 3) + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + self.geom_dim = geom_dim + self.params = params + + def __call__(self, x): + """Define an inflow expression for use as boundary condition + + Args: + x (np.ndarray): Array of coordinates + + Returns: + np.ndarray: Value of velocity at each coordinate in input array + """ + + z0 = 0.05 + d0 = 0.5 + + inflow_values = np.zeros((3, x.shape[1]), dtype=PETSc.ScalarType) + + H = 0.41 + inflow_dy = H - x[1] + inflow_dz = H - x[2] + + u_hub = self.params.fluid.u_ref + z_hub = self.params.pv_array.elevation + + if self.params.general.example == "cylinder3d": + inflow_values[0] = ( + 16.0 + * self.params.fluid.u_ref + * x[1] + * x[2] + * inflow_dy + * inflow_dz + / H**4 + ) + elif self.params.general.example == "cylinder2d": + inflow_values[0] = ( + 4 + * (self.params.fluid.u_ref) + * np.sin(np.pi / 8) + * x[1] + * (0.41 - x[1]) + / (0.41**2) + ) + elif self.params.general.example == "panels3d": + inflow_values[0] = ( + (self.params.fluid.u_ref) + * np.log(((x[2]) - d0) / z0) + / (np.log((z_hub - d0) / z0)) + ) + elif self.params.general.example == "panels2d": + inflow_values[0] = ( + (self.params.fluid.u_ref) + * np.log(((x[1]) - d0) / z0) + / (np.log((z_hub - d0) / z0)) + ) + + return inflow_values + + +def get_inflow_profile_function(domain, params, functionspace): + ndim = domain.ndim + + # IMPORTANT: this is distinct from ndim because a mesh can + # be 2D but have 3 componenets of position, e.g., a mesh + # with triangular elements can still have a z dimension = 0. + # This only gets used in making sure the inflow value + # numpy array has the correct size, e.g.: + # + # np.zeros((geom_dim, x.shape[1]). dtype=...) + geom_dim = domain.fluid.msh.geometry.dim + + inflow_function = dolfinx.fem.Function(functionspace) + + inflow_velocity = InflowVelocity(geom_dim, params) + + if params.general.example in ["cylinder3d", "cylinder2d"]: + inflow_function.interpolate(inflow_velocity) + + else: + z0 = 0.05 + d0 = 0.5 + if ndim == 3: + upper_cells = dolfinx.mesh.locate_entities( + domain.fluid.msh, ndim, lambda x: x[2] > d0 + z0 + ) + lower_cells = dolfinx.mesh.locate_entities( + domain.fluid.msh, ndim, lambda x: x[2] <= d0 + z0 + ) + elif ndim == 2: + upper_cells = dolfinx.mesh.locate_entities( + domain.fluid.msh, ndim, lambda x: x[1] > d0 + z0 + ) + lower_cells = dolfinx.mesh.locate_entities( + domain.fluid.msh, ndim, lambda x: x[1] <= d0 + z0 + ) + + inflow_function.interpolate( + lambda x: np.zeros((geom_dim, x.shape[1]), dtype=PETSc.ScalarType) + ) + + inflow_function.interpolate(inflow_velocity, upper_cells) + + return inflow_function + + +def build_velocity_boundary_conditions(domain, params, functionspace): + """Build all boundary conditions on velocity + + This method builds all the boundary conditions associated with velocity and stores in a list, ``bcu``. + + Args: + domain (:obj:`pvade.geometry.MeshManager.Domain`): A Domain object + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + + """ + # Define velocity boundary conditions + + ndim = domain.ndim + + bcu = [] + + # Generate list of locations to loop over + if ndim == 2: + bc_location_list = ["y_min", "y_max"] + + elif ndim == 3: + bc_location_list = ["y_min", "y_max", "z_min", "z_max"] + + # Iterate over all user-prescribed boundaries + for bc_location in bc_location_list: + bc_type = getattr(params.fluid, f"bc_{bc_location}") + + bc = build_vel_bc_by_type(bc_type, domain, functionspace, bc_location) + + bcu.append(bc) + + # Set all interior surfaces to no slip + bc = build_vel_bc_by_type("noslip", domain, functionspace, "internal_surface") + bcu.append(bc) + + # Set the inflow boundary condition + inflow_function = get_inflow_profile_function(domain, params, functionspace) + dofs = get_facet_dofs_by_gmsh_tag(domain, functionspace, "x_min") + bcu.append(dolfinx.fem.dirichletbc(inflow_function, dofs)) + + return bcu, inflow_function + + +def build_pressure_boundary_conditions(domain, params, functionspace): + """Build all boundary conditions on pressure + + This method builds all the boundary conditions associated with pressure and stores in a list, ``bcp``. + + Args: + domain (:obj:`pvade.geometry.MeshManager.Domain`): A Domain object + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + # Define pressure boundary conditions + bcp = [] + + zero_scalar = dolfinx.fem.Constant(domain.fluid.msh, PETSc.ScalarType(0.0)) + + dofs = get_facet_dofs_by_gmsh_tag(domain, functionspace, "x_max") + bcp.append(dolfinx.fem.dirichletbc(zero_scalar, dofs, functionspace)) + + return bcp diff --git a/pvade/geometry/MeshManager.py b/pvade/geometry/MeshManager.py new file mode 100644 index 00000000..59057689 --- /dev/null +++ b/pvade/geometry/MeshManager.py @@ -0,0 +1,388 @@ +import gmsh +import numpy as np +import os +import time +import ufl +import dolfinx +import meshio + +from importlib import import_module + + +# from pvade.geometry.panels.DomainCreation import * +class FSIDomain: + """ + This class creates the computational domain for 3D examples(3D panels, 3D cylinder) + """ + + def __init__(self, params): + """The class is initialised here + Depending on the example we are solving, we import the corresponding DomainCrearion file + We define markers which will be used for boundary tag asignment + + Args: + params (input parameters): input parameters available in the input file + """ + + # Get MPI communicators + self.comm = params.comm + self.rank = params.rank + self.num_procs = params.num_procs + + self._get_domain_markers() + + def _get_domain_markers(self): + self.domain_markers = {} + + # Facet Markers + self.domain_markers["x_min"] = {"idx": 1, "entity": "facet", "gmsh_tags": []} + self.domain_markers["x_max"] = {"idx": 2, "entity": "facet", "gmsh_tags": []} + self.domain_markers["y_min"] = {"idx": 3, "entity": "facet", "gmsh_tags": []} + self.domain_markers["y_max"] = {"idx": 4, "entity": "facet", "gmsh_tags": []} + self.domain_markers["z_min"] = {"idx": 5, "entity": "facet", "gmsh_tags": []} + self.domain_markers["z_max"] = {"idx": 6, "entity": "facet", "gmsh_tags": []} + + self.domain_markers["internal_surface"] = { + "idx": 7, + "entity": "facet", + "gmsh_tags": [], + } + + # Cell Markers + self.domain_markers["fluid"] = {"idx": 8, "entity": "cell", "gmsh_tags": []} + self.domain_markers["structure"] = {"idx": 9, "entity": "cell", "gmsh_tags": []} + + def build(self, params): + """This function call builds the geometry, marks the boundaries and creates a mesh using Gmsh.""" + + domain_creation_module = ( + f"pvade.geometry.{params.general.example}.DomainCreation" + ) + try: + # This is equivalent to "import pvade.geometry.{params.general.example}.DomainCreation as dcm" + dcm = import_module(domain_creation_module) + except: + raise ValueError(f"Could not import {domain_creation_module}") + + self.geometry = dcm.DomainCreation(params) + + # Only rank 0 builds the geometry and meshes the domain + if self.rank == 0: + self.geometry.build(params) + self.domain_markers = self.geometry.mark_surfaces( + params, self.domain_markers + ) + self.geometry.set_length_scales(params, self.domain_markers) + + if params.fluid.periodic: + self._enforce_periodicity() + + self._generate_mesh() + + # All ranks receive their portion of the mesh from rank 0 (like an MPI scatter) + self.msh, self.cell_tags, self.facet_tags = dolfinx.io.gmshio.model_to_mesh( + self.geometry.gmsh_model, self.comm, 0 + ) + + self.ndim = self.msh.topology.dim + + # Specify names for the mesh elements + self.msh.name = "mesh_total" + self.cell_tags.name = "cell_tags" + self.facet_tags.name = "facet_tags" + + self._create_submeshes_from_parent() + + self._transfer_mesh_tags_to_submeshes(params) + + def _create_submeshes_from_parent(self): + """Create submeshes from a parent mesh by cell tags. + + This function uses the cell tags to identify meshes that + are part of the structure or fluid domain and saves them + as separate mesh entities inside a new object such that meshes + are accessed like domain.fluid.msh or domain.structure.msh. + + """ + + for sub_domain_name in ["fluid", "structure"]: + if self.rank == 0: + print(f"Creating {sub_domain_name} submesh") + + # Get the idx associated with either "fluid" or "structure" + marker_id = self.domain_markers[sub_domain_name]["idx"] + + # Find all cells where cell tag = marker_id + submesh_cells = self.cell_tags.find(marker_id) + + # Use those found cells to construct a new mesh + submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh( + self.msh, self.ndim, submesh_cells + ) + + class FSISubDomain: + pass + + sub_domain = FSISubDomain() + + sub_domain.msh = submesh + sub_domain.entity_map = entity_map + sub_domain.vertex_map = vertex_map + sub_domain.geom_map = geom_map + + setattr(self, sub_domain_name, sub_domain) + + def _transfer_mesh_tags_to_submeshes(self, params): + facet_dim = self.ndim - 1 + + f_map = self.msh.topology.index_map(facet_dim) + + # Get the total number of facets in the parent mesh + num_facets = f_map.size_local + f_map.num_ghosts + all_values = np.zeros(num_facets, dtype=np.int32) + + # Assign non-zero facet tags using the facet tag indices + all_values[self.facet_tags.indices] = self.facet_tags.values + + cell_to_facet = self.msh.topology.connectivity(self.ndim, facet_dim) + + for sub_domain_name in ["fluid", "structure"]: + # Get the sub-domain object + sub_domain = getattr(self, sub_domain_name) + + # Initialize facets and create connectivity between cells and facets + sub_domain.msh.topology.create_entities(facet_dim) + sub_f_map = sub_domain.msh.topology.index_map(facet_dim) + sub_domain.msh.topology.create_connectivity(self.ndim, facet_dim) + + sub_cell_to_facet = sub_domain.msh.topology.connectivity( + self.ndim, facet_dim + ) + + sub_num_facets = sub_f_map.size_local + sub_f_map.num_ghosts + sub_values = np.empty(sub_num_facets, dtype=np.int32) + + for k, entity in enumerate(sub_domain.entity_map): + child_facets = sub_cell_to_facet.links(k) + parent_facets = cell_to_facet.links(entity) + + for child, parent in zip(child_facets, parent_facets): + sub_values[child] = all_values[parent] + + sub_domain.facet_tags = dolfinx.mesh.meshtags( + sub_domain.msh, + sub_domain.msh.topology.dim - 1, + np.arange(sub_num_facets, dtype=np.int32), + sub_values, + ) + + sub_domain.cell_tags = dolfinx.mesh.meshtags( + sub_domain.msh, + sub_domain.msh.topology.dim, + np.arange(sub_num_facets, dtype=np.int32), + np.ones(sub_num_facets), + ) + + def read(self, read_path, params): + """Read the mesh from an external file. + The User can load an existing mesh file (mesh.xdmf) + and use it to solve the CFD/CSD problem + """ + if self.rank == 0: + print("Reading the mesh from file ...") + + path_to_mesh = os.path.join(read_path, "mesh.xdmf") + + with dolfinx.io.XDMFFile(self.comm, path_to_mesh, "r") as xdmf: + self.msh = xdmf.read_mesh(name="mesh_total") + self.cell_tags = xdmf.read_meshtags(self.msh, name="cell_tags") + + self.ndim = self.msh.topology.dim + self.msh.topology.create_connectivity(self.ndim - 1, self.ndim) + + self.facet_tags = xdmf.read_meshtags(self.msh, "facet_tags") + + self._create_submeshes_from_parent() + + self._transfer_mesh_tags_to_submeshes(params) + + if self.rank == 0: + print("Done.") + + def _enforce_periodicity(self): + # TODO: Make this a generic mapping depending on which walls are marked for peridic BCs + # TODO: Copy code to enforce periodicity from old generate_and_convert_3d_meshes.py + + # Mark the front/back walls for periodic mapping + front_back_translation = [ + 1, + 0, + 0, + 0, + 0, + 1, + 0, + self.y_span, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 1, + ] + + self.geometry.gmsh_model.mesh.setPeriodic( + 2, self.dom_tags["back"], self.dom_tags["front"], front_back_translation + ) + + # Mark the left/right walls for periodic mapping + left_right_translation = [ + 1, + 0, + 0, + self.x_span, + 0, + 1, + 0, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 0, + 1, + ] + + self.geometry.gmsh_model.mesh.setPeriodic( + 2, self.dom_tags["right"], self.dom_tags["left"], left_right_translation + ) + + def _generate_mesh(self): + """This function call generates the mesh.""" + print("Starting mesh generation... ", end="") + + # Generate the mesh + tic = time.time() + + if self.geometry.ndim == 3: + self._generate_mesh_3d() + elif self.geometry.ndim == 2: + self._generate_mesh_2d() + + toc = time.time() + + if self.rank == 0: + print("Finished.") + print(f"Total meshing time = {toc-tic:.1f} s") + + def _generate_mesh_3d(self): + # Mesh.Algorithm 2D mesh algorithm + # (1: MeshAdapt, 2: Automatic, 3: Initial mesh only, 5: Delaunay, 6: Frontal-Delaunay, 7: BAMG, 8: Frontal-Delaunay for Quads, 9: Packing of Parallelograms) + # Default value: 6 + gmsh.option.setNumber("Mesh.Algorithm", 6) + + # 3D mesh algorithm + # (1: Delaunay, 3: Initial mesh only, 4: Frontal, 7: MMG3D, 9: R-tree, 10: HXT) + # Default value: 1 + gmsh.option.setNumber("Mesh.Algorithm3D", 1) + + # Mesh recombination algorithm + # (0: simple, 1: blossom, 2: simple full-quad, 3: blos- som full-quad) + # Default value: 1 + gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 0) + + # Apply recombination algorithm to all surfaces, ignoring per-surface spec Default value: 0 + # gmsh.option.setNumber("Mesh.RecombineAll", 1) + + self.geometry.gmsh_model.mesh.generate(3) + self.geometry.gmsh_model.mesh.setOrder(2) + + self.geometry.gmsh_model.mesh.optimize("Relocate3D") + self.geometry.gmsh_model.mesh.generate(3) + + def _generate_mesh_2d(self): + # Mesh.Algorithm 2D mesh algorithm + # (1: MeshAdapt, 2: Automatic, 3: Initial mesh only, 5: Delaunay, 6: Frontal-Delaunay, 7: BAMG, 8: Frontal-Delaunay for Quads, 9: Packing of Parallelograms) + # Default value: 6 + gmsh.option.setNumber("Mesh.Algorithm", 2) + + gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) + # gmsh.option.setNumber("Mesh.RecombineAll", 1) + gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) + + self.geometry.gmsh_model.mesh.generate(3) + self.geometry.gmsh_model.mesh.setOrder(2) + + self.geometry.gmsh_model.mesh.optimize("Netgen") + self.geometry.gmsh_model.mesh.generate(3) + + def write_mesh_file(self, params): + """TODO: when saving a mesh file using only dolfinx functions + it's possible certain elements of the data aren't preserved + and that the mesh won't be able to be properly read in later + on. MAKE SURE YOU CAN SAVE A MESH USING ONLY DOLFINX FUNCTIONS + AND THEN READ IN THAT SAME MESH WITHOUT A LOSS OF CAPABILITY. + """ + + if self.rank == 0: + # Save the *.msh file and *.vtk file (the latter is for visualization only) + print(f"Writing Mesh to {params.general.output_dir_mesh}...") + + if os.path.exists(params.general.output_dir_mesh) == False: + os.makedirs(params.general.output_dir_mesh) + + mesh_filename = os.path.join(params.general.output_dir_mesh, "mesh.xdmf") + + with dolfinx.io.XDMFFile(self.comm, mesh_filename, "w") as fp: + fp.write_mesh(self.msh) + fp.write_meshtags(self.cell_tags) + fp.write_meshtags(self.facet_tags) + + # TODO: This will fail if no structure mesh exists + # for sub_domain_name in ["fluid", "structure"]: + # mesh_filename = os.path.join( + # params.general.output_dir_mesh, f"mesh_{sub_domain_name}.xdmf" + # ) + + # sub_domain = getattr(self, sub_domain_name) + + # with dolfinx.io.XDMFFile(self.comm, mesh_filename, "w") as fp: + # fp.write_mesh(sub_domain.msh) + # sub_domain.msh.topology.create_connectivity( + # self.ndim - 1, self.ndim + # ) + # fp.write_meshtags(sub_domain.facet_tags) + + if self.rank == 0: + print("Done.") + + def test_mesh_functionspace(self): + P2 = ufl.VectorElement("Lagrange", self.msh.ufl_cell(), 2) + P1 = ufl.FiniteElement("Lagrange", self.msh.ufl_cell(), 1) + V = dolfinx.fem.FunctionSpace(self.msh, P2) + Q = dolfinx.fem.FunctionSpace(self.msh, P1) + + local_rangeV = V.dofmap.index_map.local_range + dofsV = np.arange(*local_rangeV) + + local_rangeQ = Q.dofmap.index_map.local_range + dofsQ = np.arange(*local_rangeQ) + + nodes_dim = 0 + self.msh.topology.create_connectivity(nodes_dim, 0) + num_nodes_owned_by_proc = self.msh.topology.index_map(nodes_dim).size_local + geometry_entitites = dolfinx.cpp.mesh.entities_to_geometry( + self.msh, + nodes_dim, + np.arange(num_nodes_owned_by_proc, dtype=np.int32), + False, + ) + points = self.msh.geometry.x + + coords = points[:] + + print(f"Rank {self.rank} owns {num_nodes_owned_by_proc} nodes\n{coords}") diff --git a/pvopt/geometry/cylinder2d/__init__.py b/pvade/geometry/__init__.py similarity index 100% rename from pvopt/geometry/cylinder2d/__init__.py rename to pvade/geometry/__init__.py diff --git a/pvade/geometry/cylinder2d/DomainCreation.py b/pvade/geometry/cylinder2d/DomainCreation.py new file mode 100644 index 00000000..e982d0d3 --- /dev/null +++ b/pvade/geometry/cylinder2d/DomainCreation.py @@ -0,0 +1,47 @@ +import gmsh + +import numpy as np + +from pvade.geometry.template.TemplateDomainCreation import TemplateDomainCreation + + +class DomainCreation(TemplateDomainCreation): + """_summary_ test + + Args: + TemplateDomainCreation (_type_): _description_ + """ + + def __init__(self, params): + """Initialize the DomainCreation object + This initializes an object that creates the computational domain. + + Args: + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + super().__init__(params) + + def build(self, params): + """This function creates the computational domain for a flow around a 2D cylinder. + + Returns: + The function returns gmsh.model which contains the geometric description of the computational domain + """ + + # All ranks create a Gmsh model object + c_x = c_y = 0.2 + r = params.domain.cyld_radius + + rectangle = self.gmsh_model.occ.addRectangle( + params.domain.x_min, + params.domain.y_min, + 0, + params.domain.x_max - params.domain.x_min, + params.domain.y_max - params.domain.y_min, + ) + + obstacle = self.gmsh_model.occ.addDisk(c_x, c_y, 0, r, r) + + self.gmsh_model.occ.cut([(2, rectangle)], [(2, obstacle)]) + + self.gmsh_model.occ.synchronize() diff --git a/pvopt/geometry/cylinder3d/__init__.py b/pvade/geometry/cylinder2d/__init__.py similarity index 100% rename from pvopt/geometry/cylinder3d/__init__.py rename to pvade/geometry/cylinder2d/__init__.py diff --git a/pvade/geometry/cylinder3d/DomainCreation.py b/pvade/geometry/cylinder3d/DomainCreation.py new file mode 100644 index 00000000..af3e5367 --- /dev/null +++ b/pvade/geometry/cylinder3d/DomainCreation.py @@ -0,0 +1,111 @@ +import gmsh +import numpy as np + +from pvade.geometry.template.TemplateDomainCreation import TemplateDomainCreation + + +class DomainCreation(TemplateDomainCreation): + """_summary_ test + + Args: + TemplateDomainCreation (_type_): _description_ + """ + + def __init__(self, params): + """Initialize the DomainCreation object + This initializes an object that creates the computational domain. + + Args: + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + super().__init__(params) + + def build(self, params): + """This function creates the computational domain for a flow around a 3D cylinder. + + Returns: + The function returns gmsh.model which contains the geometric description of the computational domain + """ + # Compute and store some useful geometric quantities + self.x_span = params.domain.x_max - params.domain.x_min + self.y_span = params.domain.y_max - params.domain.y_min + self.z_span = params.domain.z_max - params.domain.z_min + + # Create the fluid domain extent + domain = self.gmsh_model.occ.addBox( + params.domain.x_min, + params.domain.y_min, + params.domain.z_min, + self.x_span, + self.y_span, + self.z_span, + ) + + self.cyld_radius = 0.05 + height = 100.0 + + cylinder = self.gmsh_model.occ.addCylinder( + 0.5, -height / 2.0, 0.2, 0.0, height, 0.0, self.cyld_radius + ) + + # Remove each panel from the overall domain + self.gmsh_model.occ.cut([(3, domain)], [(3, cylinder)]) + + self.gmsh_model.occ.synchronize() + + def set_length_scales(self, params, domain_markers): + res_min = params.domain.l_char + + # Define a distance field from the immersed panels + distance = self.gmsh_model.mesh.field.add("Distance", 1) + self.gmsh_model.mesh.field.setNumbers( + distance, "FacesList", domain_markers["internal_surface"]["gmsh_tags"] + ) + + threshold = self.gmsh_model.mesh.field.add("Threshold") + self.gmsh_model.mesh.field.setNumber(threshold, "IField", distance) + + factor = params.domain.l_char + + self.cyld_radius = params.domain.cyld_radius + resolution = factor * self.cyld_radius / 10 + self.gmsh_model.mesh.field.setNumber(threshold, "LcMin", resolution) + self.gmsh_model.mesh.field.setNumber(threshold, "LcMax", 20 * resolution) + self.gmsh_model.mesh.field.setNumber( + threshold, "DistMin", 0.5 * self.cyld_radius + ) + self.gmsh_model.mesh.field.setNumber(threshold, "DistMax", self.cyld_radius) + + # Define a distance field from the immersed panels + zmin_dist = self.gmsh_model.mesh.field.add("Distance") + self.gmsh_model.mesh.field.setNumbers( + zmin_dist, "FacesList", domain_markers["z_min"]["gmsh_tags"] + ) + + zmin_thre = self.gmsh_model.mesh.field.add("Threshold") + self.gmsh_model.mesh.field.setNumber(zmin_thre, "IField", zmin_dist) + self.gmsh_model.mesh.field.setNumber(zmin_thre, "LcMin", 2 * resolution) + self.gmsh_model.mesh.field.setNumber(zmin_thre, "LcMax", 5 * resolution) + self.gmsh_model.mesh.field.setNumber(zmin_thre, "DistMin", 0.1) + self.gmsh_model.mesh.field.setNumber(zmin_thre, "DistMax", 0.5) + + xy_dist = self.gmsh_model.mesh.field.add("Distance") + self.gmsh_model.mesh.field.setNumbers( + xy_dist, "FacesList", domain_markers["x_min"]["gmsh_tags"] + ) + self.gmsh_model.mesh.field.setNumbers( + xy_dist, "FacesList", domain_markers["x_max"]["gmsh_tags"] + ) + + xy_thre = self.gmsh_model.mesh.field.add("Threshold") + self.gmsh_model.mesh.field.setNumber(xy_thre, "IField", xy_dist) + self.gmsh_model.mesh.field.setNumber(xy_thre, "LcMin", 2 * resolution) + self.gmsh_model.mesh.field.setNumber(xy_thre, "LcMax", 5 * resolution) + self.gmsh_model.mesh.field.setNumber(xy_thre, "DistMin", 0.1) + self.gmsh_model.mesh.field.setNumber(xy_thre, "DistMax", 0.5) + + minimum = self.gmsh_model.mesh.field.add("Min") + self.gmsh_model.mesh.field.setNumbers( + minimum, "FieldsList", [threshold, xy_thre, zmin_thre] + ) + self.gmsh_model.mesh.field.setAsBackgroundMesh(minimum) diff --git a/pvopt/geometry/panels/__init__.py b/pvade/geometry/cylinder3d/__init__.py similarity index 100% rename from pvopt/geometry/panels/__init__.py rename to pvade/geometry/cylinder3d/__init__.py diff --git a/pvade/geometry/panels2d/DomainCreation.py b/pvade/geometry/panels2d/DomainCreation.py new file mode 100644 index 00000000..563c4ae2 --- /dev/null +++ b/pvade/geometry/panels2d/DomainCreation.py @@ -0,0 +1,72 @@ +import gmsh +import numpy as np + +from pvade.geometry.template.TemplateDomainCreation import TemplateDomainCreation + + +class DomainCreation(TemplateDomainCreation): + """_summary_ test + + Args: + TemplateDomainCreation (_type_): _description_ + """ + + def __init__(self, params): + """Initialize the DomainCreation object + This initializes an object that creates the computational domain. + + Args: + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + super().__init__(params) + + def build(self, params): + """This function creates the computational domain for a 2d simulation involving N panels. + The panels are set at a distance apart, rotated at an angle theta and are elevated with a distance H from the ground. + + Returns: + The function returns gmsh.model which contains the geometric description of the computational domain + """ + + # Compute and store some useful geometric quantities + self.x_span = params.domain.x_max - params.domain.x_min + self.y_span = params.domain.y_max - params.domain.y_min + # self.z_span = params.domain.z_max - params.domain.z_min + tracker_angle_rad = np.radians(params.pv_array.tracker_angle) + + # Create the fluid domain extent + domain = self.gmsh_model.occ.addRectangle( + params.domain.x_min, + params.domain.y_min, + 0, # params.domain.z_min, + self.x_span, + self.y_span + # self.z_span, + ) + + for panel_id in range(params.pv_array.num_rows): + panel_box = self.gmsh_model.occ.addRectangle( + -0.5 * params.pv_array.panel_length, + -0.5 * params.pv_array.panel_thickness, + 0, + params.pv_array.panel_length, + params.pv_array.panel_thickness, + ) + + # Rotate the panel currently centered at (0, 0, 0) + self.gmsh_model.occ.rotate( + [(2, panel_box)], 0, 0, 0, 0, 0, 1, tracker_angle_rad + ) + + # Translate the panel [panel_loc, 0, elev] + self.gmsh_model.occ.translate( + [(2, panel_box)], + panel_id * params.pv_array.spacing[0], + params.pv_array.elevation, + 0, + ) + + # Remove each panel from the overall domain + self.gmsh_model.occ.cut([(2, domain)], [(2, panel_box)]) + + self.gmsh_model.occ.synchronize() diff --git a/pvopt/geometry/panels2d/__init__.py b/pvade/geometry/panels2d/__init__.py similarity index 100% rename from pvopt/geometry/panels2d/__init__.py rename to pvade/geometry/panels2d/__init__.py diff --git a/pvade/geometry/panels3d/DomainCreation.py b/pvade/geometry/panels3d/DomainCreation.py new file mode 100644 index 00000000..d349268e --- /dev/null +++ b/pvade/geometry/panels3d/DomainCreation.py @@ -0,0 +1,144 @@ +import gmsh +import numpy as np + +from pvade.geometry.template.TemplateDomainCreation import TemplateDomainCreation + + +class DomainCreation(TemplateDomainCreation): + """_summary_ test + + Args: + TemplateDomainCreation (_type_): _description_ + """ + + def __init__(self, params): + """Initialize the DomainCreation object + This initializes an object that creates the computational domain. + + Args: + params (:obj:`pvade.Parameters.SimParams`): A SimParams object + """ + super().__init__(params) + + def build(self, params): + """This function creates the computational domain for a 3d simulation involving N panels. + The panels are set at a distance apart, rotated at an angle theta and are elevated with a distance H from the ground. + + Returns: + The function returns gmsh.model which contains the geometric description of the computational domain + """ + # Compute and store some useful geometric quantities + self.x_span = params.domain.x_max - params.domain.x_min + self.y_span = params.domain.y_max - params.domain.y_min + self.z_span = params.domain.z_max - params.domain.z_min + tracker_angle_rad = np.radians(params.pv_array.tracker_angle) + + ndim = 3 + + # Create the fluid domain extent + domain_id = self.gmsh_model.occ.addBox( + params.domain.x_min, + params.domain.y_min, + params.domain.z_min, + self.x_span, + self.y_span, + self.z_span, + 0, + ) + + domain_tag = (ndim, domain_id) + + domain_tag_list = [] + domain_tag_list.append(domain_tag) + + panel_tag_list = [] + + for k in range(params.pv_array.num_rows): + panel_id = self.gmsh_model.occ.addBox( + -0.5 * params.pv_array.panel_length, + 0.0, + -0.5 * params.pv_array.panel_thickness, + params.pv_array.panel_length, + params.pv_array.panel_width, + params.pv_array.panel_thickness, + ) + + panel_tag = (ndim, panel_id) + panel_tag_list.append(panel_tag) + + # Rotate the panel currently centered at (0, 0, 0) + self.gmsh_model.occ.rotate([panel_tag], 0, 0, 0, 0, 1, 0, tracker_angle_rad) + + # Translate the panel [panel_loc, 0, elev] + self.gmsh_model.occ.translate( + [panel_tag], + k * params.pv_array.spacing[0], + 0, + params.pv_array.elevation, + ) + + # Fragment all panels from the overall domain + self.gmsh_model.occ.fragment(domain_tag_list, panel_tag_list) + + self.gmsh_model.occ.synchronize() + + def set_length_scales(self, params, domain_markers): + res_min = params.domain.l_char + + # Define a distance field from the immersed panels + distance = self.gmsh_model.mesh.field.add("Distance", 1) + self.gmsh_model.mesh.field.setNumbers( + distance, "FacesList", domain_markers["internal_surface"]["gmsh_tags"] + ) + + threshold = self.gmsh_model.mesh.field.add("Threshold") + self.gmsh_model.mesh.field.setNumber(threshold, "IField", distance) + + factor = params.domain.l_char + + resolution = factor * 10 * params.pv_array.panel_thickness / 2 + half_panel = params.pv_array.panel_length * np.cos( + params.pv_array.tracker_angle + ) + self.gmsh_model.mesh.field.setNumber(threshold, "LcMin", resolution * 0.5) + self.gmsh_model.mesh.field.setNumber(threshold, "LcMax", 5 * resolution) + self.gmsh_model.mesh.field.setNumber( + threshold, "DistMin", params.pv_array.spacing[0] + ) + self.gmsh_model.mesh.field.setNumber( + threshold, "DistMax", params.pv_array.spacing + half_panel + ) + + # Define a distance field from the immersed panels + zmin_dist = self.gmsh_model.mesh.field.add("Distance") + self.gmsh_model.mesh.field.setNumbers( + zmin_dist, "FacesList", domain_markers["z_min"]["gmsh_tags"] + ) + + zmin_thre = self.gmsh_model.mesh.field.add("Threshold") + self.gmsh_model.mesh.field.setNumber(zmin_thre, "IField", zmin_dist) + self.gmsh_model.mesh.field.setNumber(zmin_thre, "LcMin", 2 * resolution) + self.gmsh_model.mesh.field.setNumber(zmin_thre, "LcMax", 5 * resolution) + self.gmsh_model.mesh.field.setNumber(zmin_thre, "DistMin", 0.1) + self.gmsh_model.mesh.field.setNumber(zmin_thre, "DistMax", 0.5) + + xy_dist = self.gmsh_model.mesh.field.add("Distance") + self.gmsh_model.mesh.field.setNumbers( + xy_dist, "FacesList", domain_markers["x_min"]["gmsh_tags"] + ) + self.gmsh_model.mesh.field.setNumbers( + xy_dist, "FacesList", domain_markers["x_max"]["gmsh_tags"] + ) + + xy_thre = self.gmsh_model.mesh.field.add("Threshold") + self.gmsh_model.mesh.field.setNumber(xy_thre, "IField", xy_dist) + self.gmsh_model.mesh.field.setNumber(xy_thre, "LcMin", 2 * resolution) + self.gmsh_model.mesh.field.setNumber(xy_thre, "LcMax", 5 * resolution) + self.gmsh_model.mesh.field.setNumber(xy_thre, "DistMin", 0.1) + self.gmsh_model.mesh.field.setNumber(xy_thre, "DistMax", 0.5) + + minimum = self.gmsh_model.mesh.field.add("Min") + self.gmsh_model.mesh.field.setNumbers( + minimum, "FieldsList", [threshold, xy_thre, zmin_thre] + ) + self.gmsh_model.mesh.field.setAsBackgroundMesh(minimum) diff --git a/pvopt/geometry/template/__init__.py b/pvade/geometry/panels3d/__init__.py similarity index 100% rename from pvopt/geometry/template/__init__.py rename to pvade/geometry/panels3d/__init__.py diff --git a/pvade/geometry/template/TemplateDomainCreation.py b/pvade/geometry/template/TemplateDomainCreation.py new file mode 100644 index 00000000..daf4ebf4 --- /dev/null +++ b/pvade/geometry/template/TemplateDomainCreation.py @@ -0,0 +1,129 @@ +import gmsh +import numpy as np + + +class TemplateDomainCreation: + """This class creates the geometry used for a given example. + Gmsh is used to create the computational domain + + """ + + def __init__(self, params): + """The class is initialised here + + Args: + params (_type_): _description_ + """ + + # Get MPI communicators + self.comm = params.comm + self.rank = params.rank + self.num_procs = params.num_procs + + # Initialize Gmsh options + gmsh.initialize() + gmsh.option.setNumber("General.Terminal", 0) + gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) + + # All ranks create a Gmsh model object + self.gmsh_model = gmsh.model() + self.gmsh_model.add("domain") + self.gmsh_model.setCurrent("domain") + + def build(self, params): + """ + panels: This function creates the computational domain for a 3d simulation involving N panels. + The panels are set at a distance apart, rotated at an angle theta and are elevated with a distance H from the ground. + panels2d: This function creates the computational domain for a 2d simulation involving N panels. + The panels are set at a distance apart, rotated at an angle theta and are elevated with a distance H from the ground. + cylinder3d: This function creates the computational domain for a flow around a 3D cylinder. + cylinder2d: This function creates the computational domain for a flow around a 2D cylinder. + Returns: + The function returns gmsh.model which contains the geometric description of the computational domain + """ + pass + + def set_length_scales(self, params, domain_markers): + """This function call defines the characteristic length for the mesh in locations of interst + LcMin,LcMax,DistMin and DistMax are used to create a refined mesh in specific locations + which results in a high fidelity mesh without using a uniform element size in the whole mesh. + """ + if self.rank == 0: + all_pts = self.gmsh_model.occ.getEntities(0) + self.gmsh_model.mesh.setSize(all_pts, params.domain.l_char) + + def mark_surfaces(self, params, domain_markers): + """This function call iterates over all boundaries and assigns tags for each boundary. + The Tags are being used when appying boundary condition. + """ + + self.ndim = self.gmsh_model.get_dimension() + + # Surfaces are the entities with dimension 1 less than the mesh dimension + # i.e., surfaces have dim=2 (facets) on a 3d mesh + # and dim=1 (lines) on a 2d mesh + surf_tag_list = self.gmsh_model.occ.getEntities(self.ndim - 1) + + for surf_tag in surf_tag_list: + surf_id = surf_tag[1] + + com = self.gmsh_model.occ.getCenterOfMass(self.ndim - 1, surf_id) + + if np.isclose(com[0], params.domain.x_min): + domain_markers["x_min"]["gmsh_tags"].append(surf_id) + + elif np.allclose(com[0], params.domain.x_max): + domain_markers["x_max"]["gmsh_tags"].append(surf_id) + + elif np.allclose(com[1], params.domain.y_min): + domain_markers["y_min"]["gmsh_tags"].append(surf_id) + + elif np.allclose(com[1], params.domain.y_max): + domain_markers["y_max"]["gmsh_tags"].append(surf_id) + + elif self.ndim == 3 and np.allclose(com[2], params.domain.z_min): + domain_markers["z_min"]["gmsh_tags"].append(surf_id) + + elif self.ndim == 3 and np.allclose(com[2], params.domain.z_max): + domain_markers["z_max"]["gmsh_tags"].append(surf_id) + + else: + domain_markers["internal_surface"]["gmsh_tags"].append(surf_id) + + # Volumes are the entities with dimension equal to the mesh dimension + vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim) + + if len(vol_tag_list) > 1: + for vol_tag in vol_tag_list: + vol_id = vol_tag[1] + + if vol_id <= params.pv_array.num_rows: + # This is a panel volume, vol_id = [1, 2, ..., num_panels] + domain_markers["structure"]["gmsh_tags"].append(vol_id) + + else: + # This is the fluid around the panels, vol_id = num_panels+1 + domain_markers["fluid"]["gmsh_tags"].append(vol_id) + + else: + vol_tag = vol_tag_list[0] + vol_id = vol_tag[1] + domain_markers["fluid"]["gmsh_tags"].append(vol_id) + + for key, data in domain_markers.items(): + if len(data["gmsh_tags"]) > 0: + # Cells (i.e., entities of dim = msh.topology.dim) + if data["entity"] == "cell": + self.gmsh_model.addPhysicalGroup( + self.ndim, data["gmsh_tags"], data["idx"] + ) + self.gmsh_model.setPhysicalName(self.ndim, data["idx"], key) + + # Facets (i.e., entities of dim = msh.topology.dim - 1) + if data["entity"] == "facet": + self.gmsh_model.addPhysicalGroup( + self.ndim - 1, data["gmsh_tags"], data["idx"] + ) + self.gmsh_model.setPhysicalName(self.ndim - 1, data["idx"], key) + + return domain_markers diff --git a/pvade/geometry/template/__init__.py b/pvade/geometry/template/__init__.py new file mode 100644 index 00000000..e8edaa2f --- /dev/null +++ b/pvade/geometry/template/__init__.py @@ -0,0 +1 @@ +# Init file for geometry diff --git a/pvopt/input_schema.yaml b/pvade/input_schema.yaml similarity index 100% rename from pvopt/input_schema.yaml rename to pvade/input_schema.yaml diff --git a/pvade/tests/inputs_test/2d_cyld.yaml b/pvade/tests/inputs_test/2d_cyld.yaml new file mode 100644 index 00000000..854dda9a --- /dev/null +++ b/pvade/tests/inputs_test/2d_cyld.yaml @@ -0,0 +1,55 @@ +general: + example: cylinder2d + output_dir_mesh: output/cylinder2d/mesh + output_dir_sol: output/cylinder2d/solution + create_mesh: True + read_mesh: false + mesh_only: False +domain: + x_min: 0.0 + x_max: 2.5 + y_min: 0.0 + y_max: 0.41 + l_char: .05 + cyld_radius: 0.05 +pv_array: + num_rows: 1 # not used + elevation: 1.5 # not used + spacing: [7.0] # not used + panel_length: 2.0 # not used + panel_width: 7.0 # not used + panel_thickness: 0.1 # not used + tracker_angle: 30.0 # not used +solver: + dt: 0.000625 + t_final: 3 + save_text_every: 0.1 + save_vtk_every: 0.1 + solver1_ksp: gmres + solver2_ksp: cg + solver3_ksp: cg + solver1_pc: jacobi + solver2_pc: jacobi + solver3_pc: jacobi + save_text_interval: 0.1 + save_xdmf_interval: 0.1 +fluid: + use_eddy_viscosity: False + u_ref: 2.5 + nu: 0.001 # Establish re = 20 with diam = 0.1 and u_bar = u_ref * (4/9) + dpdx: 0.0 + turbulence_model: none + periodic: false + # c_s: 0.325 + bc_y_max: noslip # slip noslip free + bc_y_min: noslip # slip noslip free +structure: + youngs: 1.0e+09 +# Boundary_conditions: +# top: noslip +# bottom: noslip +# left: inflow +# right: free +# front: noslip +# back: noslip + diff --git a/pvade/tests/inputs_test/3d_cyld.yaml b/pvade/tests/inputs_test/3d_cyld.yaml new file mode 100644 index 00000000..d4d5c436 --- /dev/null +++ b/pvade/tests/inputs_test/3d_cyld.yaml @@ -0,0 +1,59 @@ +general: + example: cylinder3d + output_dir_mesh: output/cylinder3d/mesh + output_dir_sol: output/cylinder3d/solution + create_mesh: True + read_mesh: False + mesh_only: True +domain: + x_min: 0.0 + x_max: 2.5 + y_min: 0.0 + y_max: 0.41 + z_min: 0.0 + z_max: 0.41 + l_char: 20 #.02 + cyld_radius: 0.05 +pv_array: + num_rows: 1 # not used + elevation: 1.5 # not used + spacing: [7.0] # not used + panel_length: 2.0 # not used + panel_width: 7.0 # not used + panel_thickness: 0.1 # not used + tracker_angle: 30.0 # not used +solver: + dt: 0.025 + t_final: 0.1 + save_text_every: 0.1 + save_vtk_every: 0.1 + solver1_ksp: cg + solver2_ksp: cg + solver3_ksp: cg + solver1_pc: jacobi + solver2_pc: jacobi + solver3_pc: jacobi + save_text_interval: 0.1 + save_xdmf_interval: 0.1 +fluid: + use_eddy_viscosity: True + u_ref: 0.45 + nu: 0.001 # Establish re = 20 with diam = 0.1 and u_bar = u_ref * (4/9) + dpdx: 0.0 + turbulence_model: none + periodic: false + # c_s: 0.325 + bc_y_max: noslip # slip noslip free + bc_y_min: noslip # slip noslip free + bc_z_max: noslip # slip noslip free + bc_z_min: noslip # slip noslip free +structure: + youngs: 1.0e+09 +# Boundary_conditions: +# top: noslip +# bottom: noslip +# left: inflow +# right: free +# front: noslip +# back: noslip + diff --git a/pvade/tests/inputs_test/sim_params_alt.yaml b/pvade/tests/inputs_test/sim_params_alt.yaml new file mode 100644 index 00000000..938738fb --- /dev/null +++ b/pvade/tests/inputs_test/sim_params_alt.yaml @@ -0,0 +1,50 @@ +general: + example: panels3d + output_dir_mesh: output/panels/mesh + output_dir_sol: output/panels/solution + create_mesh: True + read_mesh: False + mesh_only: True +domain: + x_min: -10 + x_max: 50 + y_min: -20 + y_max: 27 + z_min: 0 + z_max: 20 + l_char: 20 +pv_array: + num_rows: 7 + elevation: 1.5 + spacing: [7.0] + panel_length: 2.0 + panel_width: 7.0 + panel_thickness: 0.1 + tracker_angle: 30.0 +solver: + dt: 0.005 + t_final: 10.0 + save_text_every: 0.01 + save_vtk_every: 0.01 + solver1_ksp: cg + solver2_ksp: cg + solver3_ksp: cg + solver1_pc: jacobi + solver2_pc: jacobi + solver3_pc: jacobi + save_text_interval: 0.1 + save_xdmf_interval: 0.1 +fluid: + use_eddy_viscosity: True + u_ref: 8.0 + nu: 0.01 + dpdx_val: 0.0 + turbulence_model: none + periodic: false + bc_y_max: slip # slip noslip free + bc_y_min: slip # slip noslip free + bc_z_max: slip # slip noslip free + bc_z_min: noslip # slip noslip free + # c_s: 0.325 +structure: + youngs: 1.0e+09 diff --git a/pvade/tests/inputs_test/sim_params_alt_2D.yaml b/pvade/tests/inputs_test/sim_params_alt_2D.yaml new file mode 100644 index 00000000..1a251a74 --- /dev/null +++ b/pvade/tests/inputs_test/sim_params_alt_2D.yaml @@ -0,0 +1,46 @@ +general: + example: panels2d + output_dir_mesh: output/panels2d/mesh + output_dir_sol: output/panels2d/solution + create_mesh: True + read_mesh: False + mesh_only: False +domain: + x_min: -10 + x_max: 50 + y_min: 0 + y_max: 20 + l_char: 2.0 +pv_array: + num_rows: 3 + elevation: 1.5 + spacing: [7.0] + panel_length: 2.0 + panel_width: 7.0 + panel_thickness: 0.1 + tracker_angle: 30.0 +solver: + dt: 0.005 + t_final: 0.1 + save_text_every: 0.01 + save_vtk_every: 0.01 + solver1_ksp: gmres + solver2_ksp: cg + solver3_ksp: cg + solver1_pc: jacobi + solver2_pc: jacobi + solver3_pc: jacobi + save_text_interval: 0.1 + save_xdmf_interval: 0.1 +fluid: + use_eddy_viscosity: False + u_ref: 0.8 + nu: 0.001 + dpdx_val: 0.0 + turbulence_model: none + periodic: false + bc_y_max: slip # slip noslip free + bc_y_min: noslip # slip noslip free + # c_s: 0.325 +structure: + youngs: 1.0e+09 diff --git a/pvade/tests/test_Parameters.py b/pvade/tests/test_Parameters.py new file mode 100644 index 00000000..a372c0f0 --- /dev/null +++ b/pvade/tests/test_Parameters.py @@ -0,0 +1,62 @@ +import pytest +import os + +from pvade.Parameters import SimParams + +filename = "pvade/tests/test_inputs.yaml" + + +@pytest.fixture() +def create_params(): + # Create the file + with open(filename, "w") as fp: + fp.write("general:\n") + fp.write(" example: test_name\n") + fp.write(" output_dir_sol: path/to/solution/dir\n") + fp.write("domain:\n") + fp.write(" x_min: 0.0\n") + fp.write(" x_max: 2.5\n") + fp.write("solver:\n") + fp.write(" t_final: 1.0\n") + fp.write(" dt: 0.1\n") + fp.write(" save_xdmf_interval: 0.5\n") + fp.write(" save_text_interval: 0.5\n") + + # Run the actual test + yield + + # Teardown + os.remove(filename) + + +@pytest.fixture() +def create_bad_params(): + # Create the file + with open(filename, "w") as fp: + fp.write("generalTYPO:\n") + fp.write(" example: test_name\n") + fp.write(" output_dir_sol: path/to/solution/dir\n") + + # Run the actual test + yield + + # Teardown + os.remove(filename) + + +@pytest.mark.unit +def test_parameters(create_params): + params = SimParams(filename) + + assert params.general.output_dir_sol == "path/to/solution/dir" + + assert params.domain.x_min == pytest.approx(0.0) + + assert params.solver.t_steps == 10 + assert params.solver.save_xdmf_interval_n == 5 + + +@pytest.mark.unit +def test_bad_parameters(create_bad_params): + with pytest.raises(Exception): + params = SimParams(filename) diff --git a/pvade/tests/test_fsi_mesh.py b/pvade/tests/test_fsi_mesh.py new file mode 100644 index 00000000..5e349d4a --- /dev/null +++ b/pvade/tests/test_fsi_mesh.py @@ -0,0 +1,74 @@ +import pytest +import os +import sys + +from pvade.geometry.MeshManager import FSIDomain +from pvade.DataStream import DataStream +from pvade.Parameters import SimParams +from pvade.Utilities import get_input_file, write_metrics + + +@pytest.mark.unit +def test_meshing_cylinder2d(): + input_path = "pvade/tests/inputs_test/" + + # Get the path to the input file from the command line + input_file = input_path + "2d_cyld.yaml" # get_input_file() + + # Load the parameters object specified by the input file + params = SimParams(input_file) + + # Initialize the domain and construct the initial mesh + domain = FSIDomain(params) + domain.build(params) + domain.write_mesh_file(params) + + +@pytest.mark.unit +def test_meshing_2dpanels(): + input_path = "pvade/tests/inputs_test/" + + # Get the path to the input file from the command line + input_file = input_path + "sim_params_alt_2D.yaml" # get_input_file() + + # Load the parameters object specified by the input file + params = SimParams(input_file) + + # Initialize the domain and construct the initial mesh + domain = FSIDomain(params) + domain.build(params) + domain.write_mesh_file(params) + domain.test_mesh_functionspace() + + +@pytest.mark.unit() +def test_meshing_cylinder3d(): + input_path = "pvade/tests/inputs_test/" + + # get the path to the input file from the command line + input_file = input_path + "3d_cyld.yaml" # get_input_file() + + # load the parameters object specified by the input file + params = SimParams(input_file) + + # initialize the domain and construct the initial mesh + domain = FSIDomain(params) + domain.build(params) + domain.write_mesh_file(params) + + +@pytest.mark.unit +def test_meshing_3dpanels(): + input_path = "pvade/tests/inputs_test/" + + # Get the path to the input file from the command line + input_file = input_path + "sim_params_alt.yaml" # get_input_file() + + # Load the parameters object specified by the input file + params = SimParams(input_file) + + # Initialize the domain and construct the initial mesh + domain = FSIDomain(params) + domain.build(params) + domain.write_mesh_file(params) + domain.test_mesh_functionspace() diff --git a/pvade/tests/test_mesh/cylinder2d/mesh.h5 b/pvade/tests/test_mesh/cylinder2d/mesh.h5 new file mode 100644 index 00000000..83aabe9d Binary files /dev/null and b/pvade/tests/test_mesh/cylinder2d/mesh.h5 differ diff --git a/pvade/tests/test_mesh/cylinder2d/mesh.xdmf b/pvade/tests/test_mesh/cylinder2d/mesh.xdmf new file mode 100644 index 00000000..b1fe1e95 --- /dev/null +++ b/pvade/tests/test_mesh/cylinder2d/mesh.xdmf @@ -0,0 +1,32 @@ + + + + + + + mesh.h5:/Mesh/mesh_total/topology + + + mesh.h5:/Mesh/mesh_total/geometry + + + + + + mesh.h5:/MeshTags/cell_tags/topology + + + mesh.h5:/MeshTags/cell_tags/Values + + + + + + mesh.h5:/MeshTags/facet_tags/topology + + + mesh.h5:/MeshTags/facet_tags/Values + + + + diff --git a/pvade/tests/test_mesh/cylinder3d/mesh.h5 b/pvade/tests/test_mesh/cylinder3d/mesh.h5 new file mode 100644 index 00000000..82891f6f Binary files /dev/null and b/pvade/tests/test_mesh/cylinder3d/mesh.h5 differ diff --git a/pvade/tests/test_mesh/cylinder3d/mesh.xdmf b/pvade/tests/test_mesh/cylinder3d/mesh.xdmf new file mode 100644 index 00000000..3c3b6403 --- /dev/null +++ b/pvade/tests/test_mesh/cylinder3d/mesh.xdmf @@ -0,0 +1,32 @@ + + + + + + + mesh.h5:/Mesh/mesh_total/topology + + + mesh.h5:/Mesh/mesh_total/geometry + + + + + + mesh.h5:/MeshTags/cell_tags/topology + + + mesh.h5:/MeshTags/cell_tags/Values + + + + + + mesh.h5:/MeshTags/facet_tags/topology + + + mesh.h5:/MeshTags/facet_tags/Values + + + + diff --git a/pvade/tests/test_mesh/panels2d/mesh.h5 b/pvade/tests/test_mesh/panels2d/mesh.h5 new file mode 100644 index 00000000..dfd6d0c6 Binary files /dev/null and b/pvade/tests/test_mesh/panels2d/mesh.h5 differ diff --git a/pvade/tests/test_mesh/panels2d/mesh.xdmf b/pvade/tests/test_mesh/panels2d/mesh.xdmf new file mode 100644 index 00000000..358dd44a --- /dev/null +++ b/pvade/tests/test_mesh/panels2d/mesh.xdmf @@ -0,0 +1,32 @@ + + + + + + + mesh.h5:/Mesh/mesh_total/topology + + + mesh.h5:/Mesh/mesh_total/geometry + + + + + + mesh.h5:/MeshTags/cell_tags/topology + + + mesh.h5:/MeshTags/cell_tags/Values + + + + + + mesh.h5:/MeshTags/facet_tags/topology + + + mesh.h5:/MeshTags/facet_tags/Values + + + + diff --git a/pvade/tests/test_mesh/panels3d/mesh.h5 b/pvade/tests/test_mesh/panels3d/mesh.h5 new file mode 100644 index 00000000..832639e7 Binary files /dev/null and b/pvade/tests/test_mesh/panels3d/mesh.h5 differ diff --git a/pvade/tests/test_mesh/panels3d/mesh.xdmf b/pvade/tests/test_mesh/panels3d/mesh.xdmf new file mode 100644 index 00000000..9a8fce8c --- /dev/null +++ b/pvade/tests/test_mesh/panels3d/mesh.xdmf @@ -0,0 +1,32 @@ + + + + + + + mesh.h5:/Mesh/mesh_total/topology + + + mesh.h5:/Mesh/mesh_total/geometry + + + + + + mesh.h5:/MeshTags/cell_tags/topology + + + mesh.h5:/MeshTags/cell_tags/Values + + + + + + mesh.h5:/MeshTags/facet_tags/topology + + + mesh.h5:/MeshTags/facet_tags/Values + + + + diff --git a/pvade/tests/test_solve.py b/pvade/tests/test_solve.py new file mode 100644 index 00000000..e0e46b8a --- /dev/null +++ b/pvade/tests/test_solve.py @@ -0,0 +1,157 @@ +import pytest +import os + +import numpy as np + +from pvade.fluid.FlowManager import Flow +from pvade.DataStream import DataStream +from pvade.Parameters import SimParams +from pvade.Utilities import get_input_file, write_metrics +from pvade.geometry.MeshManager import FSIDomain + +from dolfinx.common import TimingType, list_timings +import cProfile +import sys +import tqdm.autonotebook + +input_path = "pvade/tests/inputs_test/" + +solve_iter = 10 + +rtol = 3.0e-5 + + +@pytest.mark.unit +def test_flow_3dpanels(): + # Get the path to the input file from the command line + input_file = input_path + "sim_params_alt.yaml" # get_input_file() + + # Load the parameters object specified by the input file + params = SimParams(input_file) + + # Initialize the domain and construct the initial mesh + domain = FSIDomain(params) + domain.read("pvade/tests/test_mesh/panels3d", params) + # Initialize the function spaces for the flow + flow = Flow(domain) + # # # Specify the boundary conditions + flow.build_boundary_conditions(domain, params) + # # # Build the fluid forms + flow.build_forms(domain, params) + dataIO = DataStream(domain, flow, params) + + for t_step in range(solve_iter): + flow.solve(params) + + max_velocity = np.amax(flow.u_k.x.array) + max_pressure = np.amax(flow.p_k.x.array) + + print("max_velocity = ", max_velocity) + print("max_pressure = ", max_pressure) + + max_velocity_truth = 16.78961717598599 + max_pressure_truth = 161.00087027677822 + assert np.isclose(max_velocity, max_velocity_truth, rtol=rtol) + assert np.isclose(max_pressure, max_pressure_truth, rtol=rtol) + + +@pytest.mark.unit +def test_flow_2dpanels(): + # Get the path to the input file from the command line + input_file = input_path + "sim_params_alt_2D.yaml" # get_input_file() + + # Load the parameters object specified by the input file + params = SimParams(input_file) + + # Initialize the domain and construct the initial mesh + domain = FSIDomain(params) + domain.read("pvade/tests/test_mesh/panels2d", params) + # Initialize the function spaces for the flow + flow = Flow(domain) + # # # Specify the boundary conditions + flow.build_boundary_conditions(domain, params) + # # # Build the fluid forms + flow.build_forms(domain, params) + dataIO = DataStream(domain, flow, params) + + for t_step in range(solve_iter): + flow.solve(params) + + max_velocity = np.amax(flow.u_k.x.array) + max_pressure = np.amax(flow.p_k.x.array) + + print("max_velocity = ", max_velocity) + print("max_pressure = ", max_pressure) + + max_velocity_truth = 3.4734894184978726 + max_pressure_truth = 1.698213865642233 + assert np.isclose(max_velocity, max_velocity_truth, rtol=rtol) + assert np.isclose(max_pressure, max_pressure_truth, rtol=rtol) + + +@pytest.mark.unit +def test_flow_2dcylinder(): + # Get the path to the input file from the command line + input_file = input_path + "2d_cyld.yaml" # get_input_file() + + # Load the parameters object specified by the input file + params = SimParams(input_file) + + # Initialize the domain and construct the initial mesh + domain = FSIDomain(params) + domain.read("pvade/tests/test_mesh/cylinder2d", params) + # Initialize the function spaces for the flow + flow = Flow(domain) + # # # Specify the boundary conditions + flow.build_boundary_conditions(domain, params) + # # # Build the fluid forms + flow.build_forms(domain, params) + dataIO = DataStream(domain, flow, params) + + for t_step in range(solve_iter): + flow.solve(params) + + max_velocity = np.amax(flow.u_k.x.array) + max_pressure = np.amax(flow.p_k.x.array) + + print("max_velocity = ", max_velocity) + print("max_pressure = ", max_pressure) + + max_velocity_truth = 1.8113852701695827 + max_pressure_truth = 1.3044593668958533 + assert np.isclose(max_velocity, max_velocity_truth, rtol=rtol) + assert np.isclose(max_pressure, max_pressure_truth, rtol=rtol) + + +@pytest.mark.unit +def test_flow_3dcylinder(): + # Get the path to the input file from the command line + input_file = input_path + "3d_cyld.yaml" # get_input_file() + + # Load the parameters object specified by the input file + params = SimParams(input_file) + + # Initialize the domain and construct the initial mesh + domain = FSIDomain(params) + domain.read("pvade/tests/test_mesh/cylinder3d", params) + # Initialize the function spaces for the flow + flow = Flow(domain) + # # # Specify the boundary conditions + flow.build_boundary_conditions(domain, params) + # # # Build the fluid forms + flow.build_forms(domain, params) + dataIO = DataStream(domain, flow, params) + + for t_step in range(solve_iter): + flow.solve(params) + + max_velocity = np.amax(flow.u_k.x.array) + max_pressure = np.amax(flow.p_k.x.array) + + print("max_velocity = ", max_velocity) + print("max_pressure = ", max_pressure) + + max_velocity_truth = 0.6242970092279582 + max_pressure_truth = 0.30929163498498147 + assert np.isclose(max_velocity, max_velocity_truth, rtol=rtol) + assert np.isclose(max_pressure, max_pressure_truth, rtol=rtol) diff --git a/pvade/tests/test_transfer_facet_tags.py b/pvade/tests/test_transfer_facet_tags.py new file mode 100644 index 00000000..b459d869 --- /dev/null +++ b/pvade/tests/test_transfer_facet_tags.py @@ -0,0 +1,112 @@ +import pytest +import dolfinx +import numpy as np + +from pvade.Parameters import SimParams +from pvade.geometry.MeshManager import FSIDomain + + +@pytest.mark.unit +@pytest.mark.parametrize( + "sub_domain_name, marker_name", + [ + ("fluid", "x_min"), + ("structure", "x_min"), + ("fluid", "x_max"), + ("structure", "x_max"), + ("fluid", "y_min"), + ("structure", "y_min"), + ("fluid", "y_max"), + ("structure", "y_max"), + ("fluid", "z_min"), + ("structure", "z_min"), + ("fluid", "z_max"), + ("structure", "z_max"), + ("fluid", "internal_surface"), + ("structure", "internal_surface"), + ], +) +def test_transfer_facet_tags(sub_domain_name, marker_name): + input_path = "pvade/tests/inputs_test/" + + # Get the path to the input file from the command line + input_file = input_path + "sim_params_alt.yaml" # get_input_file() + + # Load the parameters object specified by the input file + params = SimParams(input_file) + + # Initialize the domain and construct the initial mesh + domain = FSIDomain(params) + domain.build(params) + + # Get the sub-domain object + sub_domain = getattr(domain, sub_domain_name) + + idx = domain.domain_markers[marker_name]["idx"] + + facets_from_tag = sub_domain.facet_tags.find(idx) + + def x_min_wall(x): + return np.isclose(x[0], params.domain.x_min) + + def x_max_wall(x): + return np.isclose(x[0], params.domain.x_max) + + def y_min_wall(x): + return np.isclose(x[1], params.domain.y_min) + + def y_max_wall(x): + return np.isclose(x[1], params.domain.y_max) + + def z_min_wall(x): + return np.isclose(x[2], params.domain.z_min) + + def z_max_wall(x): + return np.isclose(x[2], params.domain.z_max) + + def internal_surface(x): + tol = 1e-3 + + x_mid = np.logical_and( + params.domain.x_min + tol < x[0], x[0] < params.domain.x_max - tol + ) + y_mid = np.logical_and( + params.domain.y_min + tol < x[1], x[1] < params.domain.y_max - tol + ) + z_mid = np.logical_and( + params.domain.z_min + tol < x[2], x[2] < params.domain.z_max - tol + ) + + return np.logical_and(x_mid, np.logical_and(y_mid, z_mid)) + + if marker_name == "x_min": + facets_from_dolfinx = dolfinx.mesh.locate_entities_boundary( + sub_domain.msh, domain.ndim - 1, x_min_wall + ) + if marker_name == "x_max": + facets_from_dolfinx = dolfinx.mesh.locate_entities_boundary( + sub_domain.msh, domain.ndim - 1, x_max_wall + ) + if marker_name == "y_min": + facets_from_dolfinx = dolfinx.mesh.locate_entities_boundary( + sub_domain.msh, domain.ndim - 1, y_min_wall + ) + if marker_name == "y_max": + facets_from_dolfinx = dolfinx.mesh.locate_entities_boundary( + sub_domain.msh, domain.ndim - 1, y_max_wall + ) + if marker_name == "z_min": + facets_from_dolfinx = dolfinx.mesh.locate_entities_boundary( + sub_domain.msh, domain.ndim - 1, z_min_wall + ) + if marker_name == "z_max": + facets_from_dolfinx = dolfinx.mesh.locate_entities_boundary( + sub_domain.msh, domain.ndim - 1, z_max_wall + ) + if marker_name == "internal_surface": + facets_from_dolfinx = dolfinx.mesh.locate_entities_boundary( + sub_domain.msh, domain.ndim - 1, internal_surface + ) + + assert len(facets_from_tag) == len(facets_from_dolfinx) + assert np.allclose(facets_from_tag, facets_from_dolfinx) diff --git a/pvopt/.DS_Store b/pvopt/.DS_Store deleted file mode 100644 index 0bdc1e32..00000000 Binary files a/pvopt/.DS_Store and /dev/null differ diff --git a/pvopt/DataStream.py b/pvopt/DataStream.py deleted file mode 100644 index b7492511..00000000 --- a/pvopt/DataStream.py +++ /dev/null @@ -1,62 +0,0 @@ -from dolfinx import * -import numpy as np -import time -import os -import shutil -from dolfinx.io import XDMFFile, VTKFile -from mpi4py import MPI -from pathlib import Path -import pytest - - -class DataStream: - def __init__(self, domain, flow, params): - - self.comm = params.comm - self.rank = params.rank - self.num_procs = params.num_procs - - self.ndim = domain.msh.topology.dim - - self.results_filename = f"{params.general.output_dir_sol}/solution.xdmf" - - # if self.num_procs > 1: - # encoding = [XDMFFile.Encoding.HDF5] - # else: - # encoding = [XDMFFile.Encoding.HDF5, XDMFFile.Encoding.ASCII] - - # # @pytest.mark.parametrize("encoding", encodings) - - with XDMFFile(self.comm, self.results_filename, "w") as xdmf_file: - tt = 0.0 - xdmf_file.write_mesh(domain.msh) - xdmf_file.write_function(flow.u_k, 0.0) - xdmf_file.write_function(flow.p_k, 0.0) - - - def save_XDMF_files(self, flow, tt): - with XDMFFile(self.comm, self.results_filename, "a") as xdmf_file: - xdmf_file.write_function(flow.u_k, tt) - xdmf_file.write_function(flow.p_k, tt) - - # def save_XDMF_files(self, flow, params, step, xdmf_file): - # current_time = step * params["dt"] - - # xdmf_file.write(flow.u_k, current_time) - # xdmf_file.write(flow.p_k, current_time) - - # try: - # nu_T_fn = project(flow.nu_T, flow.Q, solver_type="cg") - # except: - # print("Could not project nu_T form") - # else: - # nu_T_fn.rename("nu_T", "nu_T") - # xdmf_file.write(nu_T_fn, current_time) - - # def finalize(self, mpi_info): - # if mpi_info["rank"] == 0: - # print("Solver Finished.") - # self.toc = time.time() - # self.fp_log = open(self.fp_log_name, "a") - # self.fp_log.write("\nTotal Solver Time = %f (s)\n" % (self.toc - self.tic)) - # self.fp_log.close() diff --git a/pvopt/FlowManager.py b/pvopt/FlowManager.py deleted file mode 100644 index 6de4890b..00000000 --- a/pvopt/FlowManager.py +++ /dev/null @@ -1,878 +0,0 @@ -import numpy as np -from math import log -import ufl -from dolfinx import cpp as _cpp -from dolfinx import fem -import dolfinx -from dolfinx.fem import ( - Constant, - Function, - FunctionSpace, - dirichletbc, - extract_function_spaces, - form, - locate_dofs_geometrical, - locate_dofs_topological, -) -from dolfinx.io import XDMFFile -from dolfinx.mesh import ( - CellType, - GhostMode, - create_rectangle, - create_unit_cube, - locate_entities, - locate_entities_boundary, -) -from ufl import ( - div, - dx, - grad, - inner, - dot, - nabla_grad, - lhs, - rhs, - transpose, - sym, - Identity, - ds, - CellDiameter, - sqrt, - CellVolume, -) - -from mpi4py import MPI -from petsc4py import PETSc -from petsc4py.PETSc import ScalarType - -try: - import scipy.interpolate as interp - - found_scipy = True -except: - found_scipy = False -from pathlib import Path - - -class Flow: - """This class solves the CFD problem""" - - def __init__(self, domain): - """The class is initialised here - - Args: - domain (_type_): all variables related to the computational domain - """ - - # Pressure (Scalar) - P1 = ufl.FiniteElement("Lagrange", domain.msh.ufl_cell(), 1) - self.Q = fem.FunctionSpace(domain.msh, P1) - - # Velocity (Vector) - P2 = ufl.VectorElement("Lagrange", domain.msh.ufl_cell(), 2) - self.V = fem.FunctionSpace(domain.msh, P2) - - # Stress (Tensor) - P3 = ufl.TensorElement("Lagrange", domain.msh.ufl_cell(), 2) - self.T = fem.FunctionSpace(domain.msh, P3) - - P4 = ufl.FiniteElement("DG", domain.msh.ufl_cell(), 0) - - self.DG = fem.FunctionSpace(domain.msh, P4) - - self.first_call_to_solver = True - self.first_call_to_surface_pressure = True - - # Store the dimension of the problem for convenience - self.ndim = domain.msh.topology.dim - - # Store the comm and mpi info for convenience - self.comm = domain.comm - self.rank = domain.rank - self.num_procs = domain.num_procs - - # find hmin in mesh - num_cells = domain.msh.topology.index_map(self.ndim).size_local - h = _cpp.mesh.h(domain.msh, self.ndim, range(num_cells)) - - # This value of hmin is local to the mesh portion owned by the process - hmin_local = np.amin(h) - - # collect the minimum hmin from all processes - self.hmin = np.zeros(1) - self.comm.Allreduce(hmin_local, self.hmin, op=MPI.MIN) - self.hmin = self.hmin[0] - - self.num_Q_dofs = ( - self.Q.dofmap.index_map_bs * self.Q.dofmap.index_map.size_global - ) - self.num_V_dofs = ( - self.V.dofmap.index_map_bs * self.V.dofmap.index_map.size_global - ) - - if self.rank == 0: - print(f"hmin = {self.hmin}") - print(f"Total num dofs = {self.num_Q_dofs + self.num_V_dofs}") - - def build_boundary_conditions(self, domain, params): - """define and apply boundary conditions - - Args: - domain (_type_): computational domain - params (_type_): input parameters - """ - self.facet_dim = self.ndim - 1 - - self.zero_scalar = Constant(domain.msh, PETSc.ScalarType(0.0)) - self.zero_vec = Constant(domain.msh, PETSc.ScalarType((0.0, 0.0, 0.0))) - - self._locate_boundary_entities(domain, params) - - # self._locate_boundary_dofs() - self._locate_boundary_dofs_tags(domain) - - self._build_velocity_boundary_conditions(domain, params) - - self._build_pressure_boundary_conditions(domain, params) - - def _locate_boundary_entities(self, domain, params): - """functions that defines facets using coordinates - - Args: - domain (_type_): computational domain - params (_type_): input parameters - """ - - def x_min_wall(x): - return np.isclose(x[0], params.domain.x_min) - - def x_max_wall(x): - return np.isclose(x[0], params.domain.x_max) - - def y_min_wall(x): - return np.isclose(x[1], params.domain.y_min) - - def y_max_wall(x): - return np.isclose(x[1], params.domain.y_max) - - if self.ndim == 3: - def z_min_wall(x): - return np.isclose(x[2], params.domain.z_min) - - def z_max_wall(x): - return np.isclose(x[2], params.domain.z_max) - def internal_surface(x): - x_mid = np.logical_and( - params.domain.x_min < x[0], x[0] < params.domain.x_max - ) - y_mid = np.logical_and( - params.domain.y_min < x[1], x[1] < params.domain.y_max - ) - if self.ndim == 3: - z_mid = np.logical_and( - params.domain.z_min < x[2], x[2] < params.domain.z_max - ) - return np.logical_and(x_mid, np.logical_and(y_mid, z_mid)) - elif self.ndim == 2: - return np.logical_and(x_mid, y_mid) - - self.x_min_facets = locate_entities_boundary( - domain.msh, self.facet_dim, x_min_wall - ) - self.x_max_facets = locate_entities_boundary( - domain.msh, self.facet_dim, x_max_wall - ) - self.y_min_facets = locate_entities_boundary( - domain.msh, self.facet_dim, y_min_wall - ) - self.y_max_facets = locate_entities_boundary( - domain.msh, self.facet_dim, y_max_wall - ) - if self.ndim == 3: - self.z_min_facets = locate_entities_boundary( - domain.msh, self.facet_dim, z_min_wall - ) - self.z_max_facets = locate_entities_boundary( - domain.msh, self.facet_dim, z_max_wall - ) - - self.internal_surface_facets = locate_entities_boundary( - domain.msh, self.facet_dim, internal_surface - ) - - def _locate_boundary_dofs(self): - """Locate dof on boundary for a given function space""" - self.x_min_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, self.x_min_facets - ) - self.x_max_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, self.x_max_facets - ) - self.y_min_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, self.y_min_facets - ) - self.y_max_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, self.y_max_facets - ) - self.z_min_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, self.z_min_facets - ) - self.z_max_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, self.z_max_facets - ) - self.internal_surface_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, self.internal_surface_facets - ) - - self.x_min_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, self.x_min_facets - ) - self.x_max_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, self.x_max_facets - ) - self.y_min_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, self.y_min_facets - ) - self.y_max_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, self.y_max_facets - ) - self.z_min_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, self.z_min_facets - ) - self.z_max_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, self.z_max_facets - ) - self.internal_surface_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, self.internal_surface_facets - ) - - def _locate_boundary_dofs_tags(self, domain): - """Locate DOFS using boundary tags - - Args: - domain (_type_): computational domain - params (_type_): input parameters - """ - self.x_min_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, domain.ft.find(domain.x_min_marker) - ) - - self.x_max_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, domain.ft.find(domain.x_max_marker) - ) - - self.y_min_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, domain.ft.find(domain.y_min_marker) - ) - - self.y_max_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, domain.ft.find(domain.y_max_marker) - ) - if self.ndim == 3: - self.z_min_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, domain.ft.find(domain.z_min_marker)) - - self.z_max_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, domain.ft.find(domain.z_max_marker)) - - self.internal_surface_V_dofs = locate_dofs_topological( - self.V, self.facet_dim, domain.ft.find(domain.internal_surface_marker) - ) - - self.x_min_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, domain.ft.find(domain.x_min_marker) - ) - - self.x_max_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, domain.ft.find(domain.x_max_marker) - ) - - self.y_min_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, domain.ft.find(domain.y_min_marker) - ) - - self.y_max_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, domain.ft.find(domain.y_max_marker) - ) - - if self.ndim == 3: - self.z_min_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, domain.ft.find(domain.z_min_marker)) - - self.z_max_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, domain.ft.find(domain.z_max_marker)) - - self.internal_surface_Q_dofs = locate_dofs_topological( - self.Q, self.facet_dim, domain.ft.find(domain.internal_surface_marker) - ) - - def _applybc(self, value, domain, V, marker): - """Function that applies a type of boundary - condition for given inputs - - Args: - value : scalar or function set on the dof - domain : computational domain - V : function space - marker : boundary tag - - Returns: - dirichlet boundary conditions - """ - dofs = locate_dofs_topological(V, self.facet_dim, domain.ft.find(marker)) - - # print(value, dofs, V) - - bc = dirichletbc(value, dofs, V) - - return bc - - def _build_velocity_boundary_conditions(self, domain, params): - """Sets boundary condition on the velocity - - Args: - domain (_type_): computational domain - params (_type_): input parameters - """ - # Define velocity boundary conditions - self.bcu = [] - - # x_min and x_max wall BCs - # self.bcu.append(dirichletbc(self.zero_vec, self.x_min_V_dofs, self.V)) - # self.bcu.append(dirichletbc(self.zero_vec, self.x_max_V_dofs, self.V)) - - # y_min and y_max wall BCs - # self.bcu.append(dirichletbc(self.zero_vec, self.y_min_V_dofs, self.V)) - # self.bcu.append(dirichletbc(self.zero_vec, self.y_max_V_dofs, self.V)) - # self.bcu.append(dirichletbc(self.zero_vec, self.z_min_V_dofs, self.V)) - # self.bcu.append(dirichletbc(self.zero_vec, self.z_max_V_dofs, self.V)) - - def _apply_type_bc(bc_name, functionspace, marker, bcname): - if bc_name == "free": - if domain.rank == 0: - print("free bc_name on ", bcname) - elif bc_name == "noslip": - bc = self._applybc(self.zero_vec,domain,functionspace,marker) - if domain.rank == 0: - print("noslip bc_name on ", bcname) - elif bc_name == "slip": - if "y" in bcname: - temp_func = functionspace.sub(1) - elif "z" in bcname: - temp_func = functionspace.sub(2) - bc = self._applybc(self.zero_scalar, domain, temp_func, marker) - if domain.rank == 0: - print("Slip bc_name on ", bcname, "boundary normal velocity in ", bcname, "are set to 0.") - else: - if domain.rank == 0: - print("bc_name not recognized") - - exit() - - return bc - - if self.ndim == 3: - # iterate over all noudaries - for bclocation in ( - "bc_ywall_max", - "bc_ywall_min", - "bc_zwall_max", - "bc_zwall_min", - ): - if bclocation == "bc_ywall_max": - temp_marker = domain.y_max_marker - temp_bcname = "y max" - bc_value = params.fluid.bc_ywall_max - elif bclocation == "bc_ywall_min": - temp_marker = domain.y_min_marker - temp_bcname = "y min" - bc_value = params.fluid.bc_ywall_min - elif bclocation == "bc_zwall_max": - temp_marker = domain.z_max_marker - temp_bcname = "z max" - bc_value = params.fluid.bc_zwall_max - elif bclocation == "bc_zwall_min": - temp_marker = domain.z_min_marker - temp_bcname = "z min" - bc_value = params.fluid.bc_zwall_min - - self.bcu.append( - _apply_type_bc(bc_value, self.V, temp_marker, temp_bcname) - ) - - elif self.ndim == 2: - # iterate over all noudaries - for bclocation in "bc_ywall_min", "bc_ywall_max": - if bclocation == "bc_ywall_max": - temp_marker = domain.y_max_marker - temp_bcname = "y max" - bc_value = params.fluid.bc_ywall_max - elif bclocation == "bc_ywall_min": - temp_marker = domain.y_min_marker - temp_bcname = "y min" - bc_value = params.fluid.bc_ywall_min - - self.bcu.append( - _apply_type_bc(bc_value, self.V, temp_marker, temp_bcname) - ) - - # self.bcu.append(self._noslip(self.zero_scalar,domain,self.V.sub(1),domain.geometry.y_min_marker)) - # self.bcu.append(self._noslip(self.zero_scalar,domain,self.V.sub(1),domain.geometry.y_max_marker)) - - # # z_min and z_max wall BCs - # self.bcu.append(dirichletbc(self.zero_vec,domain,self.V,domain.geometry.z_min_marker)) - # self.bcu.append(dirichletbc(self.zero_vec,domain,self.V,domain.geometry.z_max_marker)) - - # interior surface (obstacle) - self.bcu.append( - dirichletbc(self.zero_vec, self.internal_surface_V_dofs, self.V) - ) - - def inflow_profile_expression(x): - inflow_values = np.zeros( - (domain.msh.geometry.dim, x.shape[1]), dtype=PETSc.ScalarType - ) - - H = 0.41 - inflow_dy = H - x[1] - inflow_dz = H - x[2] - - u_hub = params.fluid.u_ref - z_hub = params.pv_array.elevation - - if params.general.example == 'cylinder3d': - inflow_values[0] = ( - 16.0 - * params.fluid.u_ref - * x[1] - * x[2] - * inflow_dy - * inflow_dz - / H**4 - ) - elif params.general.example == 'cylinder2d': - inflow_values[0] = ( - 4 * (params.fluid.u_ref) * np.sin( np.pi/8) * x[1] * (0.41 - x[1])/(0.41**2) - # 16.0 * params.fluid.u_ref * x[1] * inflow_dy / H**4 - ) - elif params.general.example == 'panels': - inflow_values[0] = ( - (params.fluid.u_ref) * np.log(((x[2])-d0)/z0) / (np.log((z_hub-d0)/z0)) - # 4 * params.fluid.u_ref * np.sin(x[2]* np.pi/params.domain.z_max) * x[2] * (params.domain.z_max - x[2])/(params.domain.z_max**2) - ) - elif params.general.example == 'panels2d': - inflow_values[0] = ( - (params.fluid.u_ref) * np.log(((x[1])-d0)/z0) / (np.log((z_hub-d0)/z0)) - # 4 * params.fluid.u_ref * np.sin(x[2]* np.pi/params.domain.z_max) * x[2] * (params.domain.z_max - x[2])/(params.domain.z_max**2) - ) - - return inflow_values - - self.inflow_profile = fem.Function(self.V) - - if params.general.example in ['cylinder3d', 'cylinder2d']: - self.inflow_profile.interpolate(inflow_profile_expression) - - else: - z0 = 0.05 - d0 = 0.5 - if self.ndim == 3: - upper_cells = locate_entities( - domain.msh, self.ndim, lambda x: x[2] > d0+z0 - ) - lower_cells = locate_entities( - domain.msh, self.ndim, lambda x: x[2] <= d0+z0 - ) - elif self.ndim == 2: - upper_cells = locate_entities( - domain.msh, self.ndim, lambda x: x[1] > d0+z0 - ) - lower_cells = locate_entities( - domain.msh, self.ndim, lambda x: x[1] <= d0+z0 - ) - - - self.inflow_profile.interpolate(lambda x: np.zeros( - (domain.msh.geometry.dim, x.shape[1]), dtype=PETSc.ScalarType - )) - self.inflow_profile.interpolate(inflow_profile_expression, upper_cells) - - self.bcu.append(dirichletbc(self.inflow_profile, self.x_min_V_dofs)) - - def _build_pressure_boundary_conditions(self, domain, params): - """Sets boundary condition on the velocity - - Args: - domain (_type_): computational domain - params (_type_): input parameters - """ - # Define pressure boundary conditions - self.bcp = [] - - self.bcp.append(dirichletbc(self.zero_scalar, self.x_max_Q_dofs, self.Q)) - - def build_forms(self, domain, params): - """Builds variational statements - - Args: - domain (_type_): _description_ - params (_type_): _description_ - - Returns: - _type_: _description_ - """ - # Define fluid properties - self.dpdx = fem.Constant(domain.msh, (0.0, 0.0, 0.0)) - self.dt_c = fem.Constant(domain.msh, (params.solver.dt)) - nu = fem.Constant(domain.msh, (params.fluid.nu)) - - # Define trial and test functions for velocity - self.u = ufl.TrialFunction(self.V) - self.v = ufl.TestFunction(self.V) - - # Define trial and test functions for pressure - self.p = ufl.TrialFunction(self.Q) - self.q = ufl.TestFunction(self.Q) - - # Define functions for velocity solutions - self.u_k = fem.Function(self.V, name="velocity") - self.u_k1 = fem.Function(self.V, name="velocity") - self.u_k2 = fem.Function(self.V) - - # Define functions for pressure solutions - self.p_k = fem.Function(self.Q, name="pressure") - self.p_k1 = fem.Function(self.Q, name="pressure") - - initialize_flow = True - - if initialize_flow: - self.u_k1.x.array[:] = self.inflow_profile.x.array - self.u_k2.x.array[:] = self.inflow_profile.x.array - self.u_k.x.array[:] = self.inflow_profile.x.array - - # Define expressions used in variational forms - # Crank-Nicolson velocity - U_CN = 0.5 * (self.u + self.u_k1) - - # Adams-Bashforth velocity - U_AB = 1.5 * self.u_k1 - 0.5 * self.u_k2 - - use_eddy_viscosity = True - - if use_eddy_viscosity: - # By default, don't use any eddy viscosity - filter_scale = CellVolume(domain.msh) ** (1.0 / domain.msh.topology.dim) - - # Strain rate tensor, 0.5*(du_i/dx_j + du_j/dx_i) - Sij = sym(nabla_grad(U_AB)) - - # sqrt(Sij*Sij) - strainMag = (2.0 * inner(Sij, Sij)) ** 0.5 - - # Smagorinsky constant, typically close to 0.17 - Cs = 0.17 - - # Eddy viscosity - self.nu_T = Cs**2 * filter_scale**2 * strainMag - - else: - self.nu_T = fem.Constant(domain.msh, 0.0) - - - # ================================================================# - # DEFINE VARIATIONAL FORMS - # ================================================================# - U = 0.5 * (self.u_k1 + self.u) - - def epsilon(u): - return sym(nabla_grad(u)) - - # Define stress tensor - def sigma(u, p, nu): - return 2 * nu * epsilon(u) - p * Identity(len(u)) - - fractional_step_scheme = "IPCS" - U = 0.5 * (self.u_k1 + self.u) - n = ufl.FacetNormal(domain.msh) - f = Constant( - domain.msh, (PETSc.ScalarType(0), PETSc.ScalarType(0), PETSc.ScalarType(0)) - ) - # Define variational problem for step 1: tentative velocity - self.F1 = ( - (1.0 / self.dt_c) * inner(self.u - self.u_k1, self.v) * dx - + inner(dot(U_AB, nabla_grad(U_CN)), self.v) * dx - + (nu + self.nu_T) * inner(grad(U_CN), grad(self.v)) * dx - + inner(grad(self.p_k1), self.v) * dx - - inner(self.dpdx, self.v) * dx - ) - # self.F1 = ( - # dot((self.u - self.u_k1) / self.dt_c, self.v) * dx - # + dot(dot(self.u_k1, nabla_grad(self.u_k1)), self.v) * dx - # + inner(sigma(U, self.p_k1, nu + self.nu_T), epsilon(self.v)) * dx - # + dot(self.p_k1 * n, self.v) * ds - # - dot(nu * nabla_grad(U) * n, self.v) * ds - # - dot(f, self.v) * dx - # ) - - self.a1 = form(lhs(self.F1)) - self.L1 = form(rhs(self.F1)) - - # Define variational problem for step 2: pressure correction - self.a2 = form(dot(nabla_grad(self.p), nabla_grad(self.q)) * dx) - self.L2 = form( - dot(nabla_grad(self.p_k1), nabla_grad(self.q)) * dx - - (1.0 / self.dt_c) * div(self.u_k) * self.q * dx - ) - - # Define variational problem for step 3: velocity update - self.a3 = form(dot(self.u, self.v) * dx) - self.L3 = form( - dot(self.u_k, self.v) * dx - - self.dt_c * dot(nabla_grad(self.p_k - self.p_k1), self.v) * dx - ) - - # Define a function and the form of the stress - self.panel_stress = Function(self.T) - self.stress = sigma(self.u_k, self.p_k, nu + self.nu_T) - - # Define mesh normals - # self.n = FacetNormal(domain.mesh) - - # Compute traction vector - # self.traction = dot(self.stress, self.n) - - # Create a form for projecting stress onto a tensor function space - # e.g., panel_stress.assign(project(stress, T)) - # self.a4 = inner(ufl.TrialFunction(self.T), ufl.TestFunction(self.T))*dx - # self.L4 = inner(self.stress, ufl.TestFunction(self.T))*dx - - # Create a form for projecting CFL calculation onto DG function space - cell_diam = CellDiameter(domain.msh) - cfl_form = sqrt(inner(self.u_k, self.u_k)) * self.dt_c / cell_diam - - self.cfl_vec = Function(self.DG) - self.a5 = form( - inner(ufl.TrialFunction(self.DG), ufl.TestFunction(self.DG)) * dx - ) - self.L5 = form(inner(cfl_form, ufl.TestFunction(self.DG)) * dx) - - # Set up the functions to ensure a constant flux through the outlet - outlet_cells = locate_entities( - domain.msh, self.ndim, lambda x: x[0] > params.domain.x_max - 1 - ) - self.flux_plane = fem.Function(self.V) - self.flux_plane.interpolate( - lambda x: (np.ones(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])), - outlet_cells, - ) - - # self.flux_plane = Expression('x[0] < cutoff ? 0.0 : 1.0', cutoff=domain.x_range[1]-1.0, degree=1) - self.flux_dx = ufl.Measure("dx", domain=domain.msh) # not needed ? - - # self.vol = fem.petsc.assemble_matrix(self.flux_plane*self.flux_dx) - # form1 = form(self.flux_plane*self.flux_dx) - # self.vol = fem.petsc.assemble_vector(form1) - # self.J_initial = float(fem.petsc.assemble_matrix(self.flux_plane*self.u_k[0]*self.flux_dx)/self.vol) - - # if mpi_info['rank'] == 0: - # print('J_initial = %f' % (self.J_initial)) - - # self.J_history = [self.J_initial] - self.dpdx_history = [0.0] - - def _assemble_system(self,params): - """Assemble left-hand side matrices outside the time loop and set solver options""" - - self.A1 = fem.petsc.assemble_matrix(self.a1, bcs=self.bcu) - self.A2 = fem.petsc.assemble_matrix(self.a2, bcs=self.bcp) - self.A3 = fem.petsc.assemble_matrix(self.a3) - self.A5 = fem.petsc.assemble_matrix(self.a5) - - self.A1.assemble() - self.A2.assemble() - self.A3.assemble() - self.A5.assemble() - - self.b1 = fem.petsc.assemble_vector(self.L1) - self.b2 = fem.petsc.assemble_vector(self.L2) - self.b3 = fem.petsc.assemble_vector(self.L3) - self.b5 = fem.petsc.assemble_vector(self.L5) - - self.solver_1 = PETSc.KSP().create(self.comm) - self.solver_1.setOperators(self.A1) - self.solver_1.setType(params.solver.solver1_ksp) - self.solver_1.getPC().setType(params.solver.solver1_pc) - self.solver_1.setFromOptions() - - self.solver_2 = PETSc.KSP().create(self.comm) - self.solver_2.setOperators(self.A2) - self.solver_2.setType(params.solver.solver2_ksp) - self.solver_2.getPC().setType(params.solver.solver2_pc) - self.solver_2.setFromOptions() - - self.solver_3 = PETSc.KSP().create(self.comm) - self.solver_3.setOperators(self.A3) - self.solver_3.setType(params.solver.solver3_ksp) - self.solver_3.getPC().setType(params.solver.solver2_pc) - self.solver_3.setFromOptions() - - self.solver_5 = PETSc.KSP().create(self.comm) - self.solver_5.setOperators(self.A5) - self.solver_5.setType("cg") - self.solver_5.getPC().setType("jacobi") - self.solver_5.setFromOptions() - - def solve(self, params): - """solve the variational forms - - Args: - params (_type_): _description_ - """ - if self.first_call_to_solver: - if self.rank == 0: - print("Starting Fluid Solution") - - self._assemble_system(params) - - # Calculate the tentative velocity - self._solver_step_1(params) - - # Calculate the change in pressure to enforce incompressibility - self._solver_step_2(params) - - # Update the velocity according to the pressure field - self._solver_step_3(params) - - # Compute the CFL number - self.compute_cfl() - - # Update new -> old variables - self.u_k2.x.array[:] = self.u_k1.x.array - self.u_k1.x.array[:] = self.u_k.x.array - self.p_k1.x.array[:] = self.p_k.x.array - - if self.first_call_to_solver: - self.first_call_to_solver = False - - def _solver_step_1(self, params): - - # Step 0: Re-assemble A1 since using an implicit convective term - self.A1.zeroEntries() - self.A1 = fem.petsc.assemble_matrix(self.A1, self.a1, bcs=self.bcu) - self.A1.assemble() - self.solver_1.setOperators(self.A1) - - # Step 1: Tentative velocity step - with self.b1.localForm() as loc: - loc.set(0) - - self.b1 = fem.petsc.assemble_vector(self.b1, self.L1) - - dolfinx.fem.petsc.apply_lifting(self.b1, [self.a1], [self.bcu]) - - self.b1.ghostUpdate( - addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE - ) - - fem.petsc.set_bc(self.b1, self.bcu) - - self.solver_1.solve(self.b1, self.u_k.vector) - self.u_k.x.scatter_forward() - - def _solver_step_2(self, params): - # Step 2: Pressure correction step - with self.b2.localForm() as loc: - loc.set(0) - - self.b2 = fem.petsc.assemble_vector(self.b2, self.L2) - - dolfinx.fem.petsc.apply_lifting(self.b2, [self.a2], [self.bcp]) - - self.b2.ghostUpdate( - addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE - ) - - fem.petsc.set_bc(self.b2, self.bcp) - - self.solver_2.solve(self.b2, self.p_k.vector) - self.p_k.x.scatter_forward() - - def _solver_step_3(self, params): - # Step 3: Velocity correction step - with self.b3.localForm() as loc: - loc.set(0) - - self.b3 = fem.petsc.assemble_vector(self.b3, self.L3) - - self.b3.ghostUpdate( - addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE - ) - - self.solver_3.solve(self.b3, self.u_k.vector) - self.u_k.x.scatter_forward() - - def compute_cfl(self): - # Compute the CFL number - # TODO: only do this on save/print steps? - with self.b5.localForm() as loc: - loc.set(0) - - self.b5 = fem.petsc.assemble_vector(self.b5, self.L5) - - self.b5.ghostUpdate( - addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE - ) - - self.solver_5.solve(self.b5, self.cfl_vec.vector) - self.cfl_vec.x.scatter_forward() - - cfl_max_local = np.amax(self.cfl_vec.vector.array) - - # collect the minimum hmin from all processes - self.cfl_max = np.zeros(1) - self.comm.Allreduce(cfl_max_local, self.cfl_max, op=MPI.MAX) - self.cfl_max = self.cfl_max[0] - - def adjust_dpdx_for_constant_flux(self, mpi_info): - def pid_controller(J_init, J_history, dt): - assert type(J_history) is list - - # K_p = 0.1 - # K_i = 0.4 - # K_d = 0.001 - # K_p = 0.4 - # K_i = 0.4 - # K_d = 0.01 - K_p = 0.6 - K_i = 0.4 - K_d = 0.05 - - err = J_init - np.array(J_history) - - c_p = K_p * err[-1] - c_i = K_i * dt * np.sum(err) - - try: - c_d = K_d * (err[-1] - err[-2]) / dt - except: - c_d = 0.0 - - output = c_p + c_i + c_d - - return output - - # Compute the flow through the outflow flux plane - # self.A1 = fem.petsc.assemble_matrix(self.a1 - J = float( - fem.petsc.assemble_matrix(self.flux_plane * self.u_k[0] * self.flux_dx) - / self.vol - ) - self.J_history.append(J) - - dpdx_val = pid_controller(self.J_initial, self.J_history, float(self.dt_c)) - self.dpdx_history.append(dpdx_val) - - if dpdx_val < 0.0 and mpi_info["rank"] == 0: - print("WARNING: dpdx_val = %f" % (dpdx_val)) - - self.dpdx.assign(Constant((dpdx_val, 0.0, 0.0))) diff --git a/pvopt/Parameters.py b/pvopt/Parameters.py deleted file mode 100644 index 58d33a80..00000000 --- a/pvopt/Parameters.py +++ /dev/null @@ -1,225 +0,0 @@ -import os -import yaml -import argparse - -from mpi4py import MPI -from pandas import json_normalize -from jsonschema import validate - - -class SimParams: - def __init__(self, input_file_path=None): - # Get MPI communicators - self.comm = MPI.COMM_WORLD - self.rank = self.comm.Get_rank() - self.num_procs = self.comm.Get_size() - - # Open the schema file for reading and load its contents into a dictionary - pvopt_dir = os.path.dirname(os.path.abspath(__file__)) - with open(os.path.join(pvopt_dir, "input_schema.yaml"), "r") as fp: - self.schema_dict = yaml.safe_load(fp) - - # Flatten the schema dictionary to make navigation and CLI parsing easier - self._flatten_schema_dict() - - # Initialize the single dict that will hold and combine ALL inputs - self.input_dict = {} - - if input_file_path is not None: - # Assert that the file exists - assert os.path.isfile( - input_file_path - ), f"Could not find file '{input_file_path}', check it exists" - - # Open the input file for reading and load its contents into a dictionary - with open(input_file_path, "r") as fp: - self.input_file_dict = yaml.safe_load(fp) - - # Set values as specified by the input yaml file - self._set_user_params_from_file() - - # Notes from WALID: - # - Print a warning if the value set from the file isn't in the schema - # - Print the values that don't fit - # - everything you expect to store in parameters should be in the yaml schema - - else: - # Assign all default parameters using the flattened schema - - # We can still support this, but maybe not the official getting started strategy - # (provide a yaml file for demos) - self._initialize_params_to_default() - - # Override any previous value with a command line input - self._set_user_params_from_cli() - - # Check that the complete parameter spec conforms to the schema - self._validate_inputs() - - # Store the nested dictionary as attributes on this object for easy access - # e.g., params.domain.x_max instead of params['domain']['x_max'] - self._store_dict_as_attrs() - - self._add_derived_quantities() - - def _flatten_schema_dict(self): - - flat_schema_raw = json_normalize(self.schema_dict, sep=".").to_dict() - - self.flat_schema_dict = {} - - for key in flat_schema_raw.keys(): - if "default" in key: - root = key.split(".default")[0] - - short_key = root.split(".") - short_key = [sk for sk in short_key if sk != "properties"] - short_key = ".".join(short_key) - - self.flat_schema_dict[short_key] = {} - - for subkey in ["default", "type", "units", "description"]: - try: - val = flat_schema_raw[f"{root}.{subkey}"][0] - except: - val = None - - self.flat_schema_dict[short_key][subkey] = val - - def _set_user_params_from_file(self): - - flat_input_file_dict = json_normalize(self.input_file_dict, sep=".").to_dict() - - for key, val in flat_input_file_dict.items(): - path_to_input = key.split(".") - self._set_nested_dict_value( - self.input_dict, path_to_input, val[0], error_on_missing_key=False - ) - - def _initialize_params_to_default(self): - - for key, val in self.flat_schema_dict.items(): - path_to_input = key.split(".") - self._set_nested_dict_value( - self.input_dict, - path_to_input, - val["default"], - error_on_missing_key=False, - ) - - def _set_user_params_from_cli(self): - - parser = argparse.ArgumentParser() - - for key, value in self.flat_schema_dict.items(): - - help_message = f"{value['description']} (data type = {value['type']}, Units = {value['units']})" - - if value["type"] == "string": - cli_type = str - elif value["type"] == "number": - cli_type = float - elif value["type"] == "integer": - cli_type = int - else: - cli_type = None - - parser.add_argument( - f"--{key}", metavar="", type=cli_type, help=help_message - ) - - command_line_inputs = parser.parse_args() - - # Find any command line arguments that were used and replace those entries in params - for key, value in vars(command_line_inputs).items(): - if value is not None: - path_to_input = key.split(".") - self._set_nested_dict_value( - self.input_dict, path_to_input, value, error_on_missing_key=False - ) - - def _validate_inputs(self): - # This compares the input dictionary against the yaml schema - # to ensure all values are set correctly - validate(self.input_dict, self.schema_dict) - - # self._pprint_dict(self.input_dict) - - def _store_dict_as_attrs(self): - flat_input_dict = json_normalize(self.input_dict, sep=".").to_dict() - - for key, val in flat_input_dict.items(): - path_to_input = key.split(".") - self._rec_settattr(self, path_to_input, val[0], error_on_missing_key=False) - - def _set_nested_dict_value( - self, parent_obj, path_to_input, value, error_on_missing_key=True - ): - - key = path_to_input[-1] - path_to_input = path_to_input[:-1] - - for p in path_to_input: - if error_on_missing_key: - assert ( - p in parent_obj - ), f"Could not find option '{p}' in set of valid inputs." - parent_obj = parent_obj.setdefault(p, {}) - - if error_on_missing_key: - assert ( - key in parent_obj - ), f"Could not find option '{key}' in set of valid inputs." - parent_obj[key] = value - - def _rec_settattr( - self, parent_obj, path_to_input, value, error_on_missing_key=True - ): - class ParamGroup: - pass - - if len(path_to_input) == 1: - if error_on_missing_key: - assert hasattr( - parent_obj, path_to_input[0] - ), f"Attribute '{path_to_input[0]}' is invalid." - - setattr(parent_obj, path_to_input[0], value) - - else: - if error_on_missing_key: - child_obj = getattr(parent_obj, path_to_input[0]) - else: - child_obj = getattr(parent_obj, path_to_input[0], ParamGroup()) - - setattr(parent_obj, path_to_input[0], child_obj) - parent_obj = child_obj - self._rec_settattr( - parent_obj, - path_to_input[1:], - value, - error_on_missing_key=error_on_missing_key, - ) - - def _pprint_dict(self, d, indent=0): - for key, val in d.items(): - for k in range(indent): - print("| ", end="") - - if isinstance(val, dict): - print(f"{key}:") - indent += 1 - self._pprint_dict(val, indent) - indent -= 1 - - else: - print(f"{key}: {val} ({type(val)})") - - def _add_derived_quantities(self): - self.solver.t_steps = int(self.solver.t_final / self.solver.dt) - self.solver.save_xdmf_interval_n = int( - self.solver.save_xdmf_interval / self.solver.dt - ) - self.solver.save_text_interval_n = int( - self.solver.save_text_interval / self.solver.dt - ) diff --git a/pvopt/__pycache__/DataStream.cpython-310.pyc b/pvopt/__pycache__/DataStream.cpython-310.pyc deleted file mode 100644 index 16d8803a..00000000 Binary files a/pvopt/__pycache__/DataStream.cpython-310.pyc and /dev/null differ diff --git a/pvopt/__pycache__/FlowManager.cpython-310.pyc b/pvopt/__pycache__/FlowManager.cpython-310.pyc deleted file mode 100644 index 1101dc7d..00000000 Binary files a/pvopt/__pycache__/FlowManager.cpython-310.pyc and /dev/null differ diff --git a/pvopt/__pycache__/FlowManager_copy.cpython-310.pyc b/pvopt/__pycache__/FlowManager_copy.cpython-310.pyc deleted file mode 100644 index 658b7402..00000000 Binary files a/pvopt/__pycache__/FlowManager_copy.cpython-310.pyc and /dev/null differ diff --git a/pvopt/__pycache__/IOManager.cpython-310.pyc b/pvopt/__pycache__/IOManager.cpython-310.pyc deleted file mode 100644 index 0292b8cf..00000000 Binary files a/pvopt/__pycache__/IOManager.cpython-310.pyc and /dev/null differ diff --git a/pvopt/__pycache__/MeshManager.cpython-310.pyc b/pvopt/__pycache__/MeshManager.cpython-310.pyc deleted file mode 100644 index 35c0bf09..00000000 Binary files a/pvopt/__pycache__/MeshManager.cpython-310.pyc and /dev/null differ diff --git a/pvopt/__pycache__/MeshManager.cpython-39.pyc b/pvopt/__pycache__/MeshManager.cpython-39.pyc deleted file mode 100644 index 73a9460a..00000000 Binary files a/pvopt/__pycache__/MeshManager.cpython-39.pyc and /dev/null differ diff --git a/pvopt/__pycache__/Parameters.cpython-310.pyc b/pvopt/__pycache__/Parameters.cpython-310.pyc deleted file mode 100644 index 37b4fa3e..00000000 Binary files a/pvopt/__pycache__/Parameters.cpython-310.pyc and /dev/null differ diff --git a/pvopt/__pycache__/__init__.cpython-310.pyc b/pvopt/__pycache__/__init__.cpython-310.pyc deleted file mode 100644 index ab189671..00000000 Binary files a/pvopt/__pycache__/__init__.cpython-310.pyc and /dev/null differ diff --git a/pvopt/__pycache__/__init__.cpython-39.pyc b/pvopt/__pycache__/__init__.cpython-39.pyc deleted file mode 100644 index 2b0f3514..00000000 Binary files a/pvopt/__pycache__/__init__.cpython-39.pyc and /dev/null differ diff --git a/pvopt/enforce_periodicity.py b/pvopt/enforce_periodicity.py deleted file mode 100644 index a8c84657..00000000 --- a/pvopt/enforce_periodicity.py +++ /dev/null @@ -1,142 +0,0 @@ -from dolfin import * - -parameters["mesh_partitioner"] = "ParMETIS" -parameters["partitioning_approach"] = "REPARTITION" -parameters["ParMETIS_repartitioning_weight"] = 100000000.0 -import numpy as np - -comm = MPI.comm_world -rank = comm.Get_rank() -num_procs = comm.Get_size() - - -x_range = [-10.5, 10.5] -y_range = [-3.5, 3.5] -z_range = [0.0, 15.0] - -method = 3 - -if method == 1: - - mesh_name = "periodic_3_panel/theta_30_coarse/mesh.xdmf" - - mesh = Mesh() - with XDMFFile(MPI.comm_world, mesh_name) as infile: - infile.read(mesh) - - # if rank == 0: - # print('Successfully read XDMF file') - - # mvc = MeshValueCollection('size_t', mesh, 1) - # with XDMFFile(MPI.comm_world, mesh_name.split('.xdmf')[0]+'_mf.xdmf') as infile: - # infile.read(mvc, 'name_to_read') - - # mf = MeshFunction('size_t', mesh, mvc) - - # hdfile = HDF5File(MPI.comm_world, 'converted_mesh.h5', "w") - # hdfile.write(mesh,'mesh') - # hdfile.write(mf,'boundary_markers') - # hdfile.close() - -elif method == 2: - mesh = Mesh(MPI.comm_world) - h5file = HDF5File( - MPI.comm_world, "periodic_3_panel/theta_30_coarse/h5_mesh.h5", "r" - ) - h5file.read(mesh, "/mesh", False) - # mf = MeshFunction("size_t", mesh, mesh.geometry().dim()-1) - # h5file.read(mf, '/boundary_markers') - -else: - nn = 25 - pt1 = Point(x_range[0], y_range[0], z_range[0]) - pt2 = Point(x_range[1], y_range[1], z_range[1]) - mesh = BoxMesh(pt1, pt2, nn * 3, nn, nn * 3) - -print("Finished reading mesh.") - - -class PeriodicBoundary(SubDomain): - # Left boundary is "target domain" G - def inside(self, x, on_boundary): - # return on_boundary and near(x[0], x_range[0]) or on_boundary and near(x[1], y_range[0]) - return ( - near(x[0], x_range[0]) - and x[1] < y_range[1] - DOLFIN_EPS - and on_boundary - or near(x[1], y_range[0]) - and x[0] < x_range[1] - DOLFIN_EPS - and on_boundary - ) - - # Map right boundary to left boundary and back boundary to front boundary - def map(self, x, y): - if near(x[0], x_range[1]) and near(x[1], y_range[1]): - y[0] = x[0] - (x_range[1] - x_range[0]) - y[1] = x[1] - (y_range[1] - y_range[0]) - y[2] = x[2] - - elif near(x[0], x_range[1]): - y[0] = x[0] - (x_range[1] - x_range[0]) - y[1] = x[1] - y[2] = x[2] - - elif near(x[1], y_range[1]): - y[0] = x[0] - y[1] = x[1] - (y_range[1] - y_range[0]) - y[2] = x[2] - - else: - y[0] = x[0] - 1000 - y[1] = x[1] - 1000 - y[2] = x[2] - 1000 - - -# This gives wiggle room when mapping one periodic surface to another -periodic_map_tolerance = 1e-10 - -# Velocity (Vector) -V = VectorFunctionSpace( - mesh, "P", 2, constrained_domain=PeriodicBoundary(map_tol=periodic_map_tolerance) -) -Q = FunctionSpace( - mesh, "P", 1, constrained_domain=PeriodicBoundary(map_tol=periodic_map_tolerance) -) - -interp_exp = Expression(("x[0]", "x[1]", "x[2]"), degree=2) -fn = project(interp_exp, V, solver_type="cg", preconditioner_type="jacobi") - -fn_rank = Function(Q) -fn_rank.vector()[:] = rank -fn_rank.rename("rank", "rank") - -fp_rank = File("enforce_periodicity/rank.pvd") -fp_rank << fn_rank -proj_method = 2 - -if proj_method == 1: - interp_exp = Expression(("x[0]", "x[1]", "x[2]"), degree=2) - fn = project(interp_exp, V, solver_type="cg") - -elif proj_method == 2: - - fn = Function(V) - coords = V.tabulate_dof_coordinates()[::3] - vec = np.zeros(np.size(coords)) - eps = 1e-6 - - for k, coords in enumerate(coords): - if coords[0] < x_range[0] + eps: - vec[3 * k] = -1.0 - # pass - elif coords[0] > x_range[1] - eps: - vec[3 * k] = 1.0 - # pass - - fn.vector()[:] = vec - -fn.rename("function", "function") -fp = File("enforce_periodicity/fn.pvd") -fp << fn - -print("Finished") diff --git a/pvopt/generate_and_convert_3d_meshes.py b/pvopt/generate_and_convert_3d_meshes.py deleted file mode 100644 index adfe607d..00000000 --- a/pvopt/generate_and_convert_3d_meshes.py +++ /dev/null @@ -1,386 +0,0 @@ -# purpose : generate mesh using gmsh api and convert to format readable by fenics -# boundaries are defined here using markers - -# mesh generated use adaptive somehow, might need more refinement outflow inflow - -# input: variable from file (bypassed) -# hard coded: every variable in yaml -# output: xdmf file for mesh and cdmf for boundary - -# objective: use yaml input, use adaptivity inflow out-flow, periodic vs non-periodic -# longterm terrain - -# from dolfin import * -# import meshio -# import gmsh - -# import numpy as np -# import os -# import time -# import sys -# import yaml - - -# import sys -import os -import numpy as np -import time - -# try: -# import gmsh -# except ImportError: -# print("This demo requires gmsh to be installed") -# sys.exit(0) -import yaml -import gmsh -from dolfinx.io import XDMFFile, gmshio - -from dolfinx import * -from mpi4py import MPI - - -# ================================================================ - - -def get_domain_tags(model, x_range, y_range, z_range): - - # Loop through all surfaces to find periodic tags - surf_ids = model.occ.getEntities(2) - - xm = 0.5 * (x_range[0] + x_range[1]) - ym = 0.5 * (y_range[0] + y_range[1]) - zm = 0.5 * (z_range[0] + z_range[1]) - - dom_tags = {} - - for surf in surf_ids: - tag = surf[1] - - com = model.occ.getCenterOfMass(2, tag) - - if np.isclose(com[0], x_range[0]): - dom_tags["left"] = [tag] - - elif np.allclose(com[0], x_range[1]): - dom_tags["right"] = [tag] - - elif np.allclose(com[1], y_range[0]): - dom_tags["front"] = [tag] - - elif np.allclose(com[1], y_range[1]): - dom_tags["back"] = [tag] - - elif np.allclose(com[2], z_range[0]): - dom_tags["bottom"] = [tag] - - elif np.allclose(com[2], z_range[1]): - dom_tags["top"] = [tag] - - else: - if "panel_surface" in dom_tags: - dom_tags["panel_surface"].append(tag) - else: - dom_tags["panel_surface"] = [tag] - print(dom_tags) - return dom_tags - - -# ================================================================ - - -def write_xdmf_mesh_file(msh, mt, ft, proj_dir, save_boundary_mesh=True, verbose=False): - output_mesh_name = "%s/mesh.xdmf" % (proj_dir) - with XDMFFile(msh.comm, output_mesh_name, "w") as file: - file.write_mesh(msh) - msh.topology.create_connectivity(2, 3) - file.write_meshtags( - mt, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry" - ) - file.write_meshtags( - ft, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry" - ) - - -# ================================================================ -def generate_periodic_panel_gmsh( - panel, dom, num_panel_rows, theta_deg, proj_dir, comm, lc_fac=1.0 -): - num_procs = comm.Get_size() - rank = comm.Get_rank() - model_rank = 0 - # print('rank',rank) - - if rank == model_rank: - # model = generate_periodic_panel_gmsh(panel, dom, num_panel_rows, theta_deg, proj_dir, lc_fac=3.0) - - # lc_fac=3.0 - lc_fac = 3.0 - - print("characteristic length set to", lc_fac) - - gmsh.initialize() - gmsh.option.setNumber("General.Terminal", 0) - # Set the file format, version 2.2 - gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) - - # Name the *.geo file - gmsh.model.add("multi_panel_mesh") - - # Convert to radians - theta_rad = np.radians(theta_deg) - - # Calculate some useful dimensions - x_range = [-0.5 * dom["length"], 0.5 * dom["length"]] - y_range = [-0.5 * dom["span"], 0.5 * dom["span"]] - z_range = [0, dom["height"]] - panel_loc = panel["spacing"] * np.arange(num_panel_rows) - panel_loc = panel_loc - np.mean(panel_loc) - - # Create the computational domain - print("Create the computational domain") - domain = gmsh.model.occ.addBox( - x_range[0], - y_range[0], - z_range[0], - dom["length"], - dom["span"], - dom["height"], - ) - - for z in range(num_panel_rows): - panel_box = gmsh.model.occ.addBox( - -0.5 * panel["length"], - y_range[0], - -0.5 * panel["thickness"], - panel["length"], - dom["span"], - panel["thickness"], - ) - - # Rotate the panel currently centered at (0, 0, 0) - gmsh.model.occ.rotate([(3, panel_box)], 0, 0, 0, 0, -1, 0, theta_rad) - - # Translate the panel [panel_loc, 0, elev] - gmsh.model.occ.translate( - [(3, panel_box)], panel_loc[z], 0, panel["elevation"] - ) - - # Remove each panel from the overall domain - gmsh.model.occ.cut([(3, domain)], [(3, panel_box)]) - - gmsh.model.occ.synchronize() - - # Get the surface tags for the domain and panel walls - print("Get the surface tags for the domain and panel walls") - dom_tags = get_domain_tags(gmsh.model, x_range, y_range, z_range) - - gmsh.model.addPhysicalGroup(3, [1], 11) - gmsh.model.setPhysicalName(3, 11, "fluid") - - gmsh.model.addPhysicalGroup(2, dom_tags["left"], 1) - gmsh.model.setPhysicalName(2, 1, "test") - gmsh.model.addPhysicalGroup(2, dom_tags["right"], 2) - gmsh.model.setPhysicalName(2, 2, "right") - gmsh.model.addPhysicalGroup(2, dom_tags["front"], 3) - gmsh.model.setPhysicalName(2, 3, "front") - gmsh.model.addPhysicalGroup(2, dom_tags["back"], 4) - gmsh.model.setPhysicalName(2, 4, "back") - gmsh.model.addPhysicalGroup(2, dom_tags["bottom"], 5) - gmsh.model.setPhysicalName(2, 5, "bottom") - gmsh.model.addPhysicalGroup(2, dom_tags["top"], 6) - gmsh.model.setPhysicalName(2, 6, "top") - gmsh.model.addPhysicalGroup(2, dom_tags["panel_surface"], 7) - gmsh.model.setPhysicalName(2, 7, "panel_surface") - - # Mark the front/back walls for periodic mapping - front_back_translation = [ - 1, - 0, - 0, - 0, - 0, - 1, - 0, - dom["span"], - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - gmsh.model.mesh.setPeriodic( - 2, dom_tags["back"], dom_tags["front"], front_back_translation - ) - - # Mark the left/right walls for periodic mapping - left_right_translation = [ - 1, - 0, - 0, - dom["length"], - 0, - 1, - 0, - 0, - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - gmsh.model.mesh.setPeriodic( - 2, dom_tags["right"], dom_tags["left"], left_right_translation - ) - - # Set size scales for the mesh - lc = lc_fac * panel["thickness"] - eps = 0.1 - - # Set the mesh size at each point on the panel - # panel_bbox = gmsh.model.getEntitiesInBoundingBox(x_range[0]+eps, y_range[0]-eps, z_range[0]+eps, - # x_range[1]-eps, y_range[1]+eps, z_range[1]-eps) - # gmsh.model.mesh.setSize(panel_bbox, lc) - - print("creation of the mesh") - # Define a distance field from the bottom of the domain - gmsh.model.mesh.field.add("Distance", 1) - gmsh.model.mesh.field.setNumbers(1, "FacesList", dom_tags["bottom"]) - - gmsh.model.mesh.field.add("Threshold", 2) - gmsh.model.mesh.field.setNumber(2, "IField", 1) - gmsh.model.mesh.field.setNumber(2, "LcMin", 2.0 * lc) - # gmsh.model.mesh.field.setNumber(2, 'LcMin', lc) - gmsh.model.mesh.field.setNumber(2, "LcMax", 6.0 * lc) - gmsh.model.mesh.field.setNumber(2, "DistMin", 4.5) - gmsh.model.mesh.field.setNumber(2, "DistMax", 1.0 * z_range[1]) - - # Define a distance field from the immersed panels - gmsh.model.mesh.field.add("Distance", 3) - gmsh.model.mesh.field.setNumbers(3, "FacesList", dom_tags["panel_surface"]) - - gmsh.model.mesh.field.add("Threshold", 4) - gmsh.model.mesh.field.setNumber(4, "IField", 3) - gmsh.model.mesh.field.setNumber(4, "LcMin", lc) - gmsh.model.mesh.field.setNumber(4, "LcMax", 6.0 * lc) - gmsh.model.mesh.field.setNumber(4, "DistMin", 0.5) - gmsh.model.mesh.field.setNumber(4, "DistMax", 0.6 * z_range[1]) - - gmsh.model.mesh.field.add("Min", 5) - gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2, 4]) - gmsh.model.mesh.field.setAsBackgroundMesh(5) - - gmsh.model.mesh.field.setAsBackgroundMesh(5) - - print("Meshing... ", end="") - - # Generate the mesh - tic = time.time() - gmsh.model.mesh.generate(3) - toc = time.time() - - print("Finished.") - print("Total Meshing Time = %.1f s" % (toc - tic)) - - return gmsh - - -# ================================================================ - - -def main(theta_deg=30.0, mod=None): - comm = MPI.COMM_WORLD - num_procs = comm.Get_size() - rank = comm.Get_rank() - - # Set the dimensions of the panel (all in meters) - panel = {} - with open("sim_params.yaml") as fp: - output_dic = yaml.load(fp, Loader=yaml.FullLoader) - - # Specify the number of panels to mesh - num_panel_rows = output_dic["num_panel_rows"] - - # panel['spacing'] = 7.0 - panel["spacing"] = output_dic["spacing"] - panel["elevation"] = output_dic["elevation"] - panel["length"] = output_dic["panel_length"] - panel["thickness"] = output_dic["thickness"] - - # Set the dimensions of the domain (all in meters) - dom = {} - dom["length"] = num_panel_rows * panel["spacing"] - dom["span"] = 1.0 * panel["spacing"] - dom["height"] = 10.0 * panel["elevation"] - - # Create a project directory for this run - proj_dir = "periodic_3_panel_alt/theta_%04.1f" % (theta_deg) - os.makedirs("%s" % (proj_dir), exist_ok=True) - - # earlier code generated the mesh when np = 1 and computed dofs when np ~= 1 - # with dolfinx we can use any number of np, dolfinx should generate the mesh - # on rank 1 and pass the information to all other ranks without having to write and read the mesh - - gmsh = generate_periodic_panel_gmsh( - panel, dom, num_panel_rows, theta_deg, proj_dir, comm, lc_fac=3.0 - ) - - msh, mt, ft = gmshio.model_to_mesh(gmsh.model, comm, 0) - - # print(rank,msh) - - msh.name = "cmpt_domain" - mt.name = f"{msh.name}_cells" - ft.name = f"{msh.name}_facets" - - write_xdmf_mesh_file(msh, mt, ft, proj_dir) - - mesh = msh - - # V = VectorFunctionSpace(mesh, 'P', 2) - V = fem.FunctionSpace(msh, ("Lagrange", 2)) - Q = fem.FunctionSpace(msh, ("Lagrange", 1)) - - local_rangeV = V.dofmap.index_map.local_range - dofsV = np.arange(*local_rangeV) - - local_rangeQ = Q.dofmap.index_map.local_range - dofsQ = np.arange(*local_rangeQ) - - size = np.size(dofsQ) - - print("Finished.") - print("Total Dofs = ", (size)) - print("Total Dofs = ", (dofsV)) - print("Total Dofs = ", (dofsQ)) - - -# mod would be the perturbation around .1 -def run_all_angles(mod=None): # run multiple angle - for k in range(10): - print("Found mod", mod) - - theta_deg_list = np.array( - [-70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60, 70], - dtype=np.float64, - ) - if mod is not None: - theta_deg_list += mod - - for theta_deg in theta_deg_list: - # Specify the rotation of the panels - print(theta_deg) - - main(theta_deg=theta_deg, mod=mod) - - -# run_all_angles() -# run_all_angles(mod=float(sys.argv[1])) -main() diff --git a/pvopt/geometry/MeshManager.py b/pvopt/geometry/MeshManager.py deleted file mode 100644 index 030a25ff..00000000 --- a/pvopt/geometry/MeshManager.py +++ /dev/null @@ -1,232 +0,0 @@ -from dolfinx.io import XDMFFile, gmshio -from dolfinx.fem import VectorFunctionSpace, FunctionSpace -from dolfinx.cpp import mesh as cppmesh - -from mpi4py import MPI -import gmsh -import numpy as np -import os -import time -import ufl -import dolfinx -import meshio - -# from pvopt.geometry.panels.DomainCreation import * -class Domain: - """This class creates the computational domain""" - - def __init__(self, params): - """The class is initialised here - - Args: - params (input parameters): input parameters available in the input file - """ - - # Get MPI communicators - self.comm = params.comm - self.rank = params.rank - self.num_procs = params.num_procs - - # Store a full copy of params on this object - self.params = params - - problem = self.params.general.example - - if problem == "panels": - from pvopt.geometry.panels.DomainCreation import DomainCreation - - elif problem == "cylinder3d": - from pvopt.geometry.cylinder3d.DomainCreation import DomainCreation - - elif problem == "cylinder2d": - from pvopt.geometry.cylinder2d.DomainCreation import DomainCreation - - self.geometry = DomainCreation(self.params) - - def build(self): - - # Only rank 0 builds the geometry and meshes the domain - if self.rank == 0: - self.geometry.build() - self.geometry.mark_surfaces() - self.geometry.set_length_scales() - - if self.params.fluid.periodic: - self._enforce_periodicity(self.geometry.gmsh_model) - - self._generate_mesh(self.geometry.gmsh_model) - - # All ranks receive their portion of the mesh from rank 0 (like an MPI scatter) - self.mesh, self.mt, self.ft = gmshio.model_to_mesh( - self.geometry.gmsh_model, self.comm, 0 - ) - - self.ndim = self.mesh.topology.dim - - # Specify names for the mesh elements - self.mesh.name = self.params.general.example - self.mt.name = f"{self.mesh.name}_cells" - self.ft.name = f"{self.mesh.name}_facets" - - def read(self): - """Read the mesh from external file located in output/mesh""" - if self.rank == 0: - print("Reading the mesh from file ...") - with dolfinx.io.XDMFFile( - MPI.COMM_WORLD, self.params.general.output_dir_mesh + "/mesh.xdmf", "r" - ) as xdmf: - self.mesh = xdmf.read_mesh(name="Grid") - - self.mesh.topology.create_connectivity( - self.mesh.topology.dim - 1, self.mesh.topology.dim - ) - with XDMFFile( - MPI.COMM_WORLD, self.params.general.output_dir_mesh + "/mesh_mf.xdmf", "r" - ) as infile: - self.ft = infile.read_meshtags(self.mesh, "Grid") - if self.rank == 0: - print("Done.") - - def _enforce_periodicity(self, gmsh_model): - - # TODO: Make this a generic mapping depending on which walls are marked for peridic BCs - # TODO: Copy code to enforce periodicity from old generate_and_convert_3d_meshes.py - - # Mark the front/back walls for periodic mapping - front_back_translation = [ - 1, - 0, - 0, - 0, - 0, - 1, - 0, - self.y_span, - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - gmsh_model.mesh.setPeriodic( - 2, self.dom_tags["back"], self.dom_tags["front"], front_back_translation - ) - - # Mark the left/right walls for periodic mapping - left_right_translation = [ - 1, - 0, - 0, - self.x_span, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - gmsh_model.mesh.setPeriodic( - 2, self.dom_tags["right"], self.dom_tags["left"], left_right_translation - ) - - def _generate_mesh(self, gmsh_model): - print("Starting mesh generation... ", end="") - - # Generate the mesh - tic = time.time() - # gmsh.option.setNumber("Mesh.Algorithm", 8) - # gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) - # # gmsh.option.setNumber("Mesh.RecombineAll", 1) - # gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) - # gmsh_model.mesh.generate(3) - # gmsh_model.mesh.setOrder(2) - # gmsh_model.mesh.optimize("Netgen") - gmsh_model.mesh.generate(3) - toc = time.time() - print("Finished.") - print(f"Total meshing time = {toc-tic:.1f} s") - - def write_mesh_file(self): - """ - TODO: when saving a mesh file using only dolfinx functions - it's possible certain elements of the data aren't preserved - and that the mesh won't be able to be properly read in later - on. MAKE SURE YOU CAN SAVE A MESH USING ONLY DOLFINX FUNCTIONS - AND THEN READ IN THAT SAME MESH WITHOUT A LOSS OF CAPABILITY. - """ - - if self.rank == 0: - # Save the *.msh file and *.vtk file (the latter is for visualization only) - print( - "Writing Mesh to %s... " % (self.params.general.output_dir_mesh), end="" - ) - - if os.path.exists(self.params.general.output_dir_mesh) == False: - os.mkdir(self.params.general.output_dir_mesh) - gmsh.write("%s/mesh.msh" % (self.params.general.output_dir_mesh)) - gmsh.write("%s/mesh.vtk" % (self.params.general.output_dir_mesh)) - - def create_mesh(mesh, clean_points, cell_type): - cells = mesh.get_cells_type(cell_type) - cell_data = mesh.get_cell_data("gmsh:physical", cell_type) - - out_mesh = meshio.Mesh( - points=clean_points, - cells={cell_type: cells}, - cell_data={"name_to_read": [cell_data]}, - ) - return out_mesh - - mesh_from_file = meshio.read( - f"{self.params.general.output_dir_mesh}/mesh.msh" - ) - pts = mesh_from_file.points - tetra_mesh = create_mesh(mesh_from_file, pts, "tetra") - tri_mesh = create_mesh(mesh_from_file, pts, "triangle") - - meshio.write(f"{self.params.general.output_dir_mesh}/mesh.xdmf", tetra_mesh) - meshio.write( - f"{self.params.general.output_dir_mesh}/mesh_mf.xdmf", tri_mesh - ) - print("Done.") - - def test_mesh_functionspace(self): - - P2 = ufl.VectorElement("Lagrange", self.mesh.ufl_cell(), 2) - P1 = ufl.FiniteElement("Lagrange", self.mesh.ufl_cell(), 1) - V = FunctionSpace(self.mesh, P2) - Q = FunctionSpace(self.mesh, P1) - - local_rangeV = V.dofmap.index_map.local_range - dofsV = np.arange(*local_rangeV) - - local_rangeQ = Q.dofmap.index_map.local_range - dofsQ = np.arange(*local_rangeQ) - - # coords = self.mesh.coordinates() - - nodes_dim = 0 - self.mesh.topology.create_connectivity(nodes_dim, 0) - num_nodes_owned_by_proc = self.mesh.topology.index_map(nodes_dim).size_local - geometry_entitites = cppmesh.entities_to_geometry( - self.mesh, - nodes_dim, - np.arange(num_nodes_owned_by_proc, dtype=np.int32), - False, - ) - points = self.mesh.geometry.x - - coords = points[:] - - print(f"Rank {self.rank} owns {num_nodes_owned_by_proc} nodes\n{coords}") diff --git a/pvopt/geometry/MeshManager2d.py b/pvopt/geometry/MeshManager2d.py deleted file mode 100644 index e2b54232..00000000 --- a/pvopt/geometry/MeshManager2d.py +++ /dev/null @@ -1,376 +0,0 @@ -from dolfinx.io import XDMFFile, gmshio -from dolfinx.fem import VectorFunctionSpace, FunctionSpace -from dolfinx.cpp import mesh as cppmesh - -from mpi4py import MPI -import gmsh -import numpy as np -import os -import time -import ufl -import dolfinx -import meshio - -# from pvopt.geometry.panels.DomainCreation import * -class FSIDomain: - """This class creates the computational domain - """ - def __init__(self, params): - """The class is initialised here - - Args: - params (input parameters): input parameters available in the input file - """ - - # Get MPI communicators - self.comm = MPI.COMM_WORLD - self.rank = self.comm.Get_rank() - self.num_procs = self.comm.Get_size() - - # Store a full copy of params on this object - self.params = params - - problem = self.params.general.example - - if problem == "panels": - from pvopt.geometry.panels.DomainCreation import DomainCreation - dim = 3 - elif problem == "panels2d": - from pvopt.geometry.panels2d.DomainCreation import DomainCreation - dim = 2 - elif problem == "cylinder3d": - from pvopt.geometry.cylinder3d.DomainCreation import DomainCreation - dim = 3 - elif problem == "cylinder2d": - from pvopt.geometry.cylinder2d.DomainCreation import DomainCreation - dim = 2 - - # define markers for boundaries - if dim == 3: - self.x_min_marker = 1 - self.y_min_marker = 2 - self.z_min_marker = 3 - self.x_max_marker = 4 - self.y_max_marker = 5 - self.z_max_marker = 6 - self.internal_surface_marker = 7 - self.fluid_marker = 8 - elif dim ==2: - self.x_min_marker = 1 - self.y_min_marker = 2 - self.x_max_marker = 4 - self.y_max_marker = 5 - self.internal_surface_marker = 7 - self.fluid_marker = 8 - - def build(self): - """This function call builds the geometry using Gmsh - """ - self.mesh_comm = MPI.COMM_WORLD - self.model_rank = 0 - self.gdim = 3 - - - problem = self.params.general.example - - if problem == "panels": - from pvopt.geometry.panels.DomainCreation import DomainCreation - elif problem == "panels2d": - from pvopt.geometry.panels2d.DomainCreation import DomainCreation - elif problem == "cylinder3d": - from pvopt.geometry.cylinder3d.DomainCreation import DomainCreation - elif problem == "cylinder2d": - from pvopt.geometry.cylinder2d.DomainCreation import DomainCreation - - - geometry = DomainCreation(self.params) - - # Only rank 0 builds the geometry and meshes the domain - # if self.rank == 0: - self.pv_model = geometry.build() - - - self._mark_surfaces() - self._set_length_scales_mod() # TODO: this should probably be a method of the domain creation class - - if self.params.fluid.periodic: - self._enforce_periodicity() - - self._generate_mesh() - - # All ranks receive their portion of the mesh from rank 0 (like an MPI scatter) - self.msh, self.mt, self.ft = gmshio.model_to_mesh(self.pv_model, self.comm, 0) - - self.ndim = self.msh.topology.dim - - # Specify names for the mesh elements - self.msh.name = self.params.general.example - self.mt.name = f"{self.msh.name}_cells" - self.ft.name = f"{self.msh.name}_facets" - - def read(self): - """Read the mesh from external file located in output/mesh - """ - if self.rank == 0: - print("Reading the mesh from file ...") - with dolfinx.io.XDMFFile(MPI.COMM_WORLD, self.params.general.output_dir_mesh+"/mesh.xdmf", "r") as xdmf: - self.msh = xdmf.read_mesh(name="Grid") - - self.msh.topology.create_connectivity(self.msh.topology.dim-1, self.msh.topology.dim) - with XDMFFile(MPI.COMM_WORLD,self.params.general.output_dir_mesh+"/mesh_mf.xdmf",'r') as infile: - self.ft = infile.read_meshtags(self.msh, "Grid") - if self.rank == 0: - print("Done.") - - - def _mark_surfaces(self): - """Creates boundary tags using gmsh - """ - # Loop through all surfaces to find periodic tags - surf_ids = self.pv_model.occ.getEntities(1) - - self.dom_tags = {} - - for surf in surf_ids: - tag = surf[1] - - com = self.pv_model.occ.getCenterOfMass(1, tag) - - if np.isclose(com[0], self.params.domain.x_min): - self.dom_tags["left"] = [tag] - - elif np.allclose(com[0], self.params.domain.x_max): - self.dom_tags["right"] = [tag] - - elif np.allclose(com[1], self.params.domain.y_min): - self.dom_tags["bottom"] = [tag] - - elif np.allclose(com[1], self.params.domain.y_max): - self.dom_tags["top"] = [tag] - - else: - if "panel_surface" in self.dom_tags: - self.dom_tags["panel_surface"].append(tag) - else: - self.dom_tags["panel_surface"] = [tag] - - - self.pv_model.addPhysicalGroup(2, [1], 11) - self.pv_model.setPhysicalName(2, 11, "fluid") - - self.pv_model.addPhysicalGroup(1, self.dom_tags["left"], self.x_min_marker) - self.pv_model.setPhysicalName(1, self.x_min_marker, "left") - self.pv_model.addPhysicalGroup(1, self.dom_tags["right"], self.x_max_marker) - self.pv_model.setPhysicalName(1, self.x_max_marker, "right") - self.pv_model.addPhysicalGroup(1, self.dom_tags["bottom"], self.y_min_marker) - self.pv_model.setPhysicalName(1, self.y_min_marker, "bottom") - self.pv_model.addPhysicalGroup(1, self.dom_tags["top"], self.y_max_marker) - self.pv_model.setPhysicalName(1, self.y_max_marker, "top") - - self.pv_model.addPhysicalGroup(1, self.dom_tags["panel_surface"], self.internal_surface_marker) - self.pv_model.setPhysicalName(1, self.internal_surface_marker, "panel_surface") - - def _set_length_scales_mod(self): - res_min = self.params.domain.l_char - if self.mesh_comm.rank == self.model_rank: - # Define a distance field from the immersed panels - self.pv_model.mesh.field.add("Distance", 1) - self.pv_model.mesh.field.setNumbers( - 1, "FacesList", self.dom_tags["panel_surface"] - ) - self.pv_model.mesh.field.add("Threshold", 2) - self.pv_model.mesh.field.setNumber(2, "IField", 1) - if 'cylinder2d' in self.params.general.example: - self.cyld_radius = self.params.domain.cyld_radius - self.pv_model.mesh.field.setNumber(2, "LcMin", self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(2, "LcMax", 3.0 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(2, "DistMin", 2.0 * self.cyld_radius) - self.pv_model.mesh.field.setNumber(2, "DistMax", 10.0 * self.cyld_radius) - - else: - self.pv_model.mesh.field.setNumber(2, "LcMin", 0.2*self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(2, "LcMax", 3.0 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(2, "DistMin", .5 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(2, "DistMax", 1. * self.params.domain.l_char) - - # self.pv_model.mesh.field.setNumber(4, "LcMin", self.params.pv_array.panel_thickness) - # self.pv_model.mesh.field.setNumber(4, "LcMax", 3.0 * self.params.pv_array.panel_thickness) - # self.pv_model.mesh.field.setNumber(4, "DistMin", 2* self.params.pv_array.panel_length) - # self.pv_model.mesh.field.setNumber(4, "DistMax", 10.* self.params.pv_array.panel_length) - - self.pv_model.mesh.field.add("Distance", 3) - self.pv_model.mesh.field.setNumbers(3, "FacesList", self.dom_tags["top"]) - self.pv_model.mesh.field.setNumbers(3, "FacesList", self.dom_tags["bottom"]) - - self.pv_model.mesh.field.add("Threshold", 4) - self.pv_model.mesh.field.setNumber(4, "IField", 3) - self.pv_model.mesh.field.setNumber(4, "LcMin", 5.0 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(4, "LcMax", 10.0 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(4, "DistMin", 0.1) - self.pv_model.mesh.field.setNumber(4, "DistMax", 0.5) - - self.pv_model.mesh.field.add("Distance", 5) - self.pv_model.mesh.field.setNumbers(5, "FacesList", self.dom_tags["left"]) - self.pv_model.mesh.field.setNumbers(5, "FacesList", self.dom_tags["right"]) - - self.pv_model.mesh.field.add("Threshold", 6) - self.pv_model.mesh.field.setNumber(6, "IField", 5) - self.pv_model.mesh.field.setNumber(6, "LcMin", 2.0 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(6, "LcMax", 5.0 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(6, "DistMin", 0.1) - self.pv_model.mesh.field.setNumber(6, "DistMax", 0.5) - - - self.pv_model.mesh.field.add("Min", 7) - self.pv_model.mesh.field.setNumbers(7, "FieldsList", [2,4,6]) - self.pv_model.mesh.field.setAsBackgroundMesh(7) - - - - if self.mesh_comm.rank == self.model_rank: - gmsh.option.setNumber("Mesh.Algorithm", 8) - gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) - gmsh.option.setNumber("Mesh.RecombineAll", 1) - gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) - - gmsh.model.mesh.generate(2) - gmsh.model.mesh.setOrder(2) - gmsh.model.mesh.optimize("Netgen") - - - def _enforce_periodicity(self): - - # TODO: Make this a generic mapping depending on which walls are marked for peridic BCs - # TODO: Copy code to enforce periodicity from old generate_and_convert_3d_meshes.py - - # Mark the front/back walls for periodic mapping - front_back_translation = [ - 1, - 0, - 0, - 0, - 0, - 1, - 0, - self.y_span, - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - self.pv_model.mesh.setPeriodic( - 2, self.dom_tags["back"], self.dom_tags["front"], front_back_translation - ) - - # Mark the left/right walls for periodic mapping - left_right_translation = [ - 1, - 0, - 0, - self.x_span, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - self.pv_model.mesh.setPeriodic( - 2, self.dom_tags["right"], self.dom_tags["left"], left_right_translation - ) - - def _generate_mesh(self): - if self.rank == 0: - print("Starting mesh generation... ", end="") - - # Generate the mesh - tic = time.time() - gmsh.option.setNumber("Mesh.Algorithm", 8) - gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) - # gmsh.option.setNumber("Mesh.RecombineAll", 1) - gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) - self.pv_model.mesh.generate(3) - self.pv_model.mesh.setOrder(2) - self.pv_model.mesh.optimize("Netgen") - self.pv_model.mesh.generate(3) - toc = time.time() - if self.rank == 0: - print("Finished.") - print(f"Total meshing time = {toc-tic:.1f} s") - - - def write_mesh_file(self): - ''' - TODO: when saving a mesh file using only dolfinx functions - it's possible certain elements of the data aren't preserved - and that the mesh won't be able to be properly read in later - on. MAKE SURE YOU CAN SAVE A MESH USING ONLY DOLFINX FUNCTIONS - AND THEN READ IN THAT SAME MESH WITHOUT A LOSS OF CAPABILITY. - ''' - - if self.rank == 0: - # Save the *.msh file and *.vtk file (the latter is for visualization only) - print('Writing Mesh to %s... ' % (self.params.general.output_dir_mesh), end='') - - if os.path.exists(self.params.general.output_dir_mesh) == False: - os.mkdir(self.params.general.output_dir_mesh) - gmsh.write('%s/mesh.msh' % (self.params.general.output_dir_mesh)) - gmsh.write('%s/mesh.vtk' % (self.params.general.output_dir_mesh)) - def create_mesh(mesh, clean_points, cell_type): - cells = mesh.get_cells_type(cell_type) - cell_data = mesh.get_cell_data("gmsh:physical", cell_type) - - out_mesh = meshio.Mesh(points=clean_points, cells={ - cell_type: cells}, cell_data={"name_to_read": [cell_data]}) - return out_mesh - - mesh_from_file = meshio.read(f'{self.params.general.output_dir_mesh}/mesh.msh') - pts = mesh_from_file.points - tetra_mesh = create_mesh(mesh_from_file, pts, "quad") - tri_mesh = create_mesh(mesh_from_file, pts, "line") - - meshio.write(f'{self.params.general.output_dir_mesh}/mesh.xdmf', tetra_mesh) - meshio.write(f'{self.params.general.output_dir_mesh}/mesh_mf.xdmf', tri_mesh) - print("Done.") - - def test_mesh_functionspace(self): - - P2 = ufl.VectorElement("Lagrange", self.msh.ufl_cell(), 2) - P1 = ufl.FiniteElement("Lagrange", self.msh.ufl_cell(), 1) - V = FunctionSpace(self.msh, P2) - Q = FunctionSpace(self.msh, P1) - - local_rangeV = V.dofmap.index_map.local_range - dofsV = np.arange(*local_rangeV) - - local_rangeQ = Q.dofmap.index_map.local_range - dofsQ = np.arange(*local_rangeQ) - - # coords = self.mesh.coordinates() - - nodes_dim = 0 - self.msh.topology.create_connectivity(nodes_dim, 0) - num_nodes_owned_by_proc = self.msh.topology.index_map(nodes_dim).size_local - geometry_entitites = cppmesh.entities_to_geometry( - self.msh, - nodes_dim, - np.arange(num_nodes_owned_by_proc, dtype=np.int32), - False, - ) - points = self.msh.geometry.x - - coords = points[:] - - print(f"Rank {self.rank} owns {num_nodes_owned_by_proc} nodes\n{coords}") diff --git a/pvopt/geometry/MeshManager3d.py b/pvopt/geometry/MeshManager3d.py deleted file mode 100644 index bcd7b31f..00000000 --- a/pvopt/geometry/MeshManager3d.py +++ /dev/null @@ -1,434 +0,0 @@ -from dolfinx.io import XDMFFile, gmshio -from dolfinx.fem import VectorFunctionSpace, FunctionSpace -from dolfinx.cpp import mesh as cppmesh - -from mpi4py import MPI -import gmsh -import numpy as np -import os -import time -import ufl -import dolfinx -import meshio - -# from pvopt.geometry.panels.DomainCreation import * -class FSIDomain: - """This class creates the computational domain - """ - def __init__(self, params): - """The class is initialised here - - Args: - params (input parameters): input parameters available in the input file - """ - - # Get MPI communicators - self.comm = MPI.COMM_WORLD - self.rank = self.comm.Get_rank() - self.num_procs = self.comm.Get_size() - - # Store a full copy of params on this object - self.params = params - - problem = self.params.general.example - - if problem == "panels": - from pvopt.geometry.panels.DomainCreation import DomainCreation - dim = 3 - elif problem == "panels2d": - from pvopt.geometry.panels2d.DomainCreation import DomainCreation - dim = 2 - elif problem == "cylinder3d": - from pvopt.geometry.cylinder3d.DomainCreation import DomainCreation - dim = 3 - elif problem == "cylinder2d": - from pvopt.geometry.cylinder2d.DomainCreation import DomainCreation - dim = 2 - - # define markers for boundaries - if dim == 3: - self.x_min_marker = 1 - self.y_min_marker = 2 - self.z_min_marker = 3 - self.x_max_marker = 4 - self.y_max_marker = 5 - self.z_max_marker = 6 - self.internal_surface_marker = 7 - self.fluid_marker = 8 - elif dim ==2: - self.x_min_marker = 1 - self.y_min_marker = 2 - self.x_max_marker = 4 - self.y_max_marker = 5 - self.internal_surface_marker = 7 - self.fluid_marker = 8 - - def build(self): - """This function call builds the geometry using Gmsh - """ - self.mesh_comm = MPI.COMM_WORLD - self.model_rank = 0 - self.gdim = 3 - - - problem = self.params.general.example - - if problem == "panels": - from pvopt.geometry.panels.DomainCreation import DomainCreation - elif problem == "panels2d": - from pvopt.geometry.panels2d.DomainCreation import DomainCreation - elif problem == "cylinder3d": - from pvopt.geometry.cylinder3d.DomainCreation import DomainCreation - elif problem == "cylinder2d": - from pvopt.geometry.cylinder2d.DomainCreation import DomainCreation - - - geometry = DomainCreation(self.params) - - # Only rank 0 builds the geometry and meshes the domain - # if self.rank == 0: - self.pv_model = geometry.build() - - - self._mark_surfaces() - self._set_length_scales_mod() # TODO: this should probably be a method of the domain creation class - - if self.params.fluid.periodic: - self._enforce_periodicity() - - self._generate_mesh() - - # All ranks receive their portion of the mesh from rank 0 (like an MPI scatter) - self.msh, self.mt, self.ft = gmshio.model_to_mesh(self.pv_model, self.comm, 0) - - self.ndim = self.msh.topology.dim - - # Specify names for the mesh elements - self.msh.name = self.params.general.example - self.mt.name = f"{self.msh.name}_cells" - self.ft.name = f"{self.msh.name}_facets" - - def read(self): - """Read the mesh from external file located in output/mesh - """ - if self.rank == 0: - print("Reading the mesh from file ...") - with dolfinx.io.XDMFFile(MPI.COMM_WORLD, self.params.general.output_dir_mesh+"/mesh.xdmf", "r") as xdmf: - self.msh = xdmf.read_mesh(name="Grid") - - self.msh.topology.create_connectivity(self.msh.topology.dim-1, self.msh.topology.dim) - with XDMFFile(MPI.COMM_WORLD,self.params.general.output_dir_mesh+"/mesh_mf.xdmf",'r') as infile: - self.ft = infile.read_meshtags(self.msh, "Grid") - if self.rank == 0: - print("Done.") - - - def _mark_surfaces(self): - """Creates boundary tags using gmsh - """ - # Loop through all surfaces to find periodic tags - surf_ids = self.pv_model.occ.getEntities(2) - - self.dom_tags = {} - - for surf in surf_ids: - tag = surf[1] - - com = self.pv_model.occ.getCenterOfMass(2, tag) - - if np.isclose(com[0], self.params.domain.x_min): - self.dom_tags["left"] = [tag] - - elif np.allclose(com[0], self.params.domain.x_max): - self.dom_tags["right"] = [tag] - - elif np.allclose(com[1], self.params.domain.y_min): - self.dom_tags["front"] = [tag] - - elif np.allclose(com[1], self.params.domain.y_max): - self.dom_tags["back"] = [tag] - - elif np.allclose(com[2], self.params.domain.z_min): - self.dom_tags["bottom"] = [tag] - - elif np.allclose(com[2], self.params.domain.z_max): - self.dom_tags["top"] = [tag] - - else: - if "panel_surface" in self.dom_tags: - self.dom_tags["panel_surface"].append(tag) - else: - self.dom_tags["panel_surface"] = [tag] - - - self.pv_model.addPhysicalGroup(3, [1], 11) - self.pv_model.setPhysicalName(3, 11, "fluid") - - self.pv_model.addPhysicalGroup(2, self.dom_tags["left"], self.x_min_marker) - self.pv_model.setPhysicalName(2, self.x_min_marker, "left") - self.pv_model.addPhysicalGroup(2, self.dom_tags["right"], self.x_max_marker) - self.pv_model.setPhysicalName(2, self.x_max_marker, "right") - self.pv_model.addPhysicalGroup(2, self.dom_tags["front"], self.y_min_marker) - self.pv_model.setPhysicalName(2, self.y_min_marker, "front") - self.pv_model.addPhysicalGroup(2, self.dom_tags["back"], self.y_max_marker) - self.pv_model.setPhysicalName(2, self.y_max_marker, "back") - self.pv_model.addPhysicalGroup(2, self.dom_tags["bottom"], self.z_min_marker) - self.pv_model.setPhysicalName(2, self.z_min_marker, "bottom") - self.pv_model.addPhysicalGroup(2, self.dom_tags["top"], self.z_max_marker) - self.pv_model.setPhysicalName(2, self.z_max_marker, "top") - self.pv_model.addPhysicalGroup(2, self.dom_tags["panel_surface"], self.internal_surface_marker) - self.pv_model.setPhysicalName(2, self.internal_surface_marker, "panel_surface") - - def _set_length_scales_mod(self): - res_min = self.params.domain.l_char - if self.mesh_comm.rank == self.model_rank: - # Define a distance field from the immersed panels - distance = self.pv_model.mesh.field.add("Distance", 1) - self.pv_model.mesh.field.setNumbers(distance, "FacesList", self.dom_tags["panel_surface"]) - - threshold = self.pv_model.mesh.field.add("Threshold") - self.pv_model.mesh.field.setNumber(threshold, "IField", distance) - - - factor = self.params.domain.l_char - if 'cylinder3d' in self.params.general.example: - self.cyld_radius = self.params.domain.cyld_radius - resolution = factor * self.cyld_radius / 10 - self.pv_model.mesh.field.setNumber(threshold, "LcMin", resolution) - self.pv_model.mesh.field.setNumber(threshold, "LcMax", 20 * resolution) - self.pv_model.mesh.field.setNumber(threshold, "DistMin", .5 * self.cyld_radius) - self.pv_model.mesh.field.setNumber(threshold, "DistMax", self.cyld_radius) - - else: - resolution = factor * 10*self.params.pv_array.panel_thickness/2 - half_panel = self.params.pv_array.panel_length * np.cos(self.params.pv_array.tracker_angle) - self.pv_model.mesh.field.setNumber(threshold, "LcMin", resolution*0.5) - self.pv_model.mesh.field.setNumber(threshold, "LcMax", 5*resolution) - self.pv_model.mesh.field.setNumber(threshold, "DistMin", self.params.pv_array.spacing[0]) - self.pv_model.mesh.field.setNumber(threshold, "DistMax", self.params.pv_array.spacing+half_panel) - - - # # Define a distance field from the immersed panels - # zmin_dist = self.pv_model.mesh.field.add("Distance") - # self.pv_model.mesh.field.setNumbers(zmin_dist, "FacesList", self.dom_tags["bottom"]) - - # zmin_thre = self.pv_model.mesh.field.add("Threshold") - # self.pv_model.mesh.field.setNumber(zmin_thre, "IField", zmin_dist) - # self.pv_model.mesh.field.setNumber(zmin_thre, "LcMin", .5*resolution) - # self.pv_model.mesh.field.setNumber(zmin_thre, "LcMax", 5*resolution) - # self.pv_model.mesh.field.setNumber(zmin_thre, "DistMin", 1.5) - # self.pv_model.mesh.field.setNumber(zmin_thre, "DistMax", 1.5+half_panel) - - - - - # other_dist = self.pv_model.mesh.field.add("Distance") - # self.pv_model.mesh.field.setNumbers(other_dist, "FacesList", self.dom_tags["front"]) - # self.pv_model.mesh.field.setNumbers(other_dist, "FacesList", self.dom_tags["back"]) - # self.pv_model.mesh.field.setNumbers(other_dist, "FacesList", self.dom_tags["top"]) - - # other_thre = self.pv_model.mesh.field.add("Threshold") - # self.pv_model.mesh.field.setNumber(other_thre, "IField", other_dist) - # self.pv_model.mesh.field.setNumber(other_thre, "LcMin", 5.0 * self.params.domain.l_char) - # self.pv_model.mesh.field.setNumber(other_thre, "LcMax", 10.0 * self.params.domain.l_char) - # self.pv_model.mesh.field.setNumber(other_thre, "DistMin", 0.1) - # self.pv_model.mesh.field.setNumber(other_thre, "DistMax", 0.5) - - - xy_dist = self.pv_model.mesh.field.add("Distance") - self.pv_model.mesh.field.setNumbers(xy_dist, "FacesList", self.dom_tags["left"]) - self.pv_model.mesh.field.setNumbers(xy_dist, "FacesList", self.dom_tags["right"]) - - xy_thre = self.pv_model.mesh.field.add("Threshold") - self.pv_model.mesh.field.setNumber(xy_thre, "IField", xy_dist) - self.pv_model.mesh.field.setNumber(xy_thre, "LcMin", 2 * resolution) - self.pv_model.mesh.field.setNumber(xy_thre, "LcMax", 5* resolution) - self.pv_model.mesh.field.setNumber(xy_thre, "DistMin", 0.1) - self.pv_model.mesh.field.setNumber(xy_thre, "DistMax", 0.5) - - - minimum = self.pv_model.mesh.field.add("Min") - self.pv_model.mesh.field.setNumbers(minimum, "FieldsList", [threshold, xy_thre ]) - self.pv_model.mesh.field.setAsBackgroundMesh(minimum) - - # self.pv_model.mesh.field.setAsBackgroundMesh(7) - - # gmsh.model.mesh.field.add("Distance", 1) - # gmsh.model.mesh.field.setNumbers(1, "FacesList", self.dom_tags["panel_surface"]) - # r = res_min - # resolution = r - # gmsh.model.mesh.field.add("Threshold", 2) - # gmsh.model.mesh.field.setNumber(2, "IField", 1) - # gmsh.model.mesh.field.setNumber(2, "LcMin", resolution) - # gmsh.model.mesh.field.setNumber(2, "LcMax", 2*resolution) - # gmsh.model.mesh.field.setNumber(2, "DistMin", 0.5*r) - # gmsh.model.mesh.field.setNumber(2, "DistMax", r) - - # # We add a similar threshold at the inlet to ensure fully resolved inlet flow - - # gmsh.model.mesh.field.add("Distance", 3) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['back'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['front'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['top'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['bottom'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['left'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['right'] ) - # gmsh.model.mesh.field.add("Threshold", 4) - # gmsh.model.mesh.field.setNumber(4, "IField", 3) - # gmsh.model.mesh.field.setNumber(4, "LcMin", 5*resolution) - # gmsh.model.mesh.field.setNumber(4, "LcMax", 10*resolution) - # gmsh.model.mesh.field.setNumber(4, "DistMin", 0.1) - # gmsh.model.mesh.field.setNumber(4, "DistMax", 0.5) - - # # We combine these fields by using the minimum field - # gmsh.model.mesh.field.add("Min", 5) - # gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2, 4]) - # gmsh.model.mesh.field.setAsBackgroundMesh(5) - - - # if self.mesh_comm.rank == self.model_rank: - # gmsh.option.setNumber("Mesh.Algorithm", 5) - # gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) - # gmsh.option.setNumber("Mesh.RecombineAll", 1) - # gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) - # gmsh.model.mesh.generate(3) - # gmsh.model.mesh.setOrder(2) - # gmsh.model.mesh.optimize("Netgen") - - - def _enforce_periodicity(self): - - # TODO: Make this a generic mapping depending on which walls are marked for peridic BCs - # TODO: Copy code to enforce periodicity from old generate_and_convert_3d_meshes.py - - # Mark the front/back walls for periodic mapping - front_back_translation = [ - 1, - 0, - 0, - 0, - 0, - 1, - 0, - self.y_span, - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - self.pv_model.mesh.setPeriodic( - 2, self.dom_tags["back"], self.dom_tags["front"], front_back_translation - ) - - # Mark the left/right walls for periodic mapping - left_right_translation = [ - 1, - 0, - 0, - self.x_span, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - self.pv_model.mesh.setPeriodic( - 2, self.dom_tags["right"], self.dom_tags["left"], left_right_translation - ) - - def _generate_mesh(self): - if self.rank == 0: - print("Starting mesh generation... ", end="") - - # Generate the mesh - tic = time.time() - gmsh.option.setNumber("Mesh.Algorithm", 8) - gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) - # gmsh.option.setNumber("Mesh.RecombineAll", 1) - gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) - self.pv_model.mesh.generate(3) - self.pv_model.mesh.setOrder(2) - self.pv_model.mesh.optimize("Netgen") - self.pv_model.mesh.generate(3) - toc = time.time() - if self.rank == 0: - print("Finished.") - print(f"Total meshing time = {toc-tic:.1f} s") - - - def write_mesh_file(self): - ''' - TODO: when saving a mesh file using only dolfinx functions - it's possible certain elements of the data aren't preserved - and that the mesh won't be able to be properly read in later - on. MAKE SURE YOU CAN SAVE A MESH USING ONLY DOLFINX FUNCTIONS - AND THEN READ IN THAT SAME MESH WITHOUT A LOSS OF CAPABILITY. - ''' - - if self.rank == 0: - # Save the *.msh file and *.vtk file (the latter is for visualization only) - print('Writing Mesh to %s... ' % (self.params.general.output_dir_mesh), end='') - - if os.path.exists(self.params.general.output_dir_mesh) == False: - os.mkdir(self.params.general.output_dir_mesh) - gmsh.write('%s/mesh.msh' % (self.params.general.output_dir_mesh)) - gmsh.write('%s/mesh.vtk' % (self.params.general.output_dir_mesh)) - def create_mesh(mesh, clean_points, cell_type): - cells = mesh.get_cells_type(cell_type) - cell_data = mesh.get_cell_data("gmsh:physical", cell_type) - - out_mesh = meshio.Mesh(points=clean_points, cells={ - cell_type: cells}, cell_data={"name_to_read": [cell_data]}) - return out_mesh - - mesh_from_file = meshio.read(f'{self.params.general.output_dir_mesh}/mesh.msh') - pts = mesh_from_file.points - tetra_mesh = create_mesh(mesh_from_file, pts, "tetra") - tri_mesh = create_mesh(mesh_from_file, pts, "triangle") - - meshio.write(f'{self.params.general.output_dir_mesh}/mesh.xdmf', tetra_mesh) - meshio.write(f'{self.params.general.output_dir_mesh}/mesh_mf.xdmf', tri_mesh) - print("Done.") - - def test_mesh_functionspace(self): - - P2 = ufl.VectorElement("Lagrange", self.msh.ufl_cell(), 2) - P1 = ufl.FiniteElement("Lagrange", self.msh.ufl_cell(), 1) - V = FunctionSpace(self.msh, P2) - Q = FunctionSpace(self.msh, P1) - - local_rangeV = V.dofmap.index_map.local_range - dofsV = np.arange(*local_rangeV) - - local_rangeQ = Q.dofmap.index_map.local_range - dofsQ = np.arange(*local_rangeQ) - - # coords = self.mesh.coordinates() - - nodes_dim = 0 - self.msh.topology.create_connectivity(nodes_dim, 0) - num_nodes_owned_by_proc = self.msh.topology.index_map(nodes_dim).size_local - geometry_entitites = cppmesh.entities_to_geometry( - self.msh, - nodes_dim, - np.arange(num_nodes_owned_by_proc, dtype=np.int32), - False, - ) - points = self.msh.geometry.x - - coords = points[:] - - print(f"Rank {self.rank} owns {num_nodes_owned_by_proc} nodes\n{coords}") diff --git a/pvopt/geometry/__pycache__/MeshManager.cpython-310.pyc b/pvopt/geometry/__pycache__/MeshManager.cpython-310.pyc deleted file mode 100644 index e98e516e..00000000 Binary files a/pvopt/geometry/__pycache__/MeshManager.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/__pycache__/MeshManager.cpython-39.pyc b/pvopt/geometry/__pycache__/MeshManager.cpython-39.pyc deleted file mode 100644 index 531765d1..00000000 Binary files a/pvopt/geometry/__pycache__/MeshManager.cpython-39.pyc and /dev/null differ diff --git a/pvopt/geometry/__pycache__/MeshManager2d.cpython-310.pyc b/pvopt/geometry/__pycache__/MeshManager2d.cpython-310.pyc deleted file mode 100644 index deeb1272..00000000 Binary files a/pvopt/geometry/__pycache__/MeshManager2d.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/__pycache__/MeshManager3d.cpython-310.pyc b/pvopt/geometry/__pycache__/MeshManager3d.cpython-310.pyc deleted file mode 100644 index d9ae344c..00000000 Binary files a/pvopt/geometry/__pycache__/MeshManager3d.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/__pycache__/__init__.cpython-310.pyc b/pvopt/geometry/__pycache__/__init__.cpython-310.pyc deleted file mode 100644 index 21352ede..00000000 Binary files a/pvopt/geometry/__pycache__/__init__.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/__pycache__/__init__.cpython-39.pyc b/pvopt/geometry/__pycache__/__init__.cpython-39.pyc deleted file mode 100644 index 6a2cb97a..00000000 Binary files a/pvopt/geometry/__pycache__/__init__.cpython-39.pyc and /dev/null differ diff --git a/pvopt/geometry/cylinder2d/DomainCreation.py b/pvopt/geometry/cylinder2d/DomainCreation.py deleted file mode 100644 index ee8f7cb1..00000000 --- a/pvopt/geometry/cylinder2d/DomainCreation.py +++ /dev/null @@ -1,51 +0,0 @@ -from dolfinx.io import XDMFFile, gmshio -from dolfinx.fem import VectorFunctionSpace, FunctionSpace -from dolfinx.cpp import mesh as cppmesh -from mpi4py import MPI -import gmsh -import numpy as np -import os -import time -import ufl -import dolfinx -import meshio -# from pvopt.geometry.panels.domain_creation import * - - -class FSIDomain: - def __init__(self, params): - - # Get MPI communicators - self.comm = MPI.COMM_WORLD - self.rank = self.comm.Get_rank() - self.num_procs = self.comm.Get_size() - - # Store a full copy of params on this object - self.params = params - - def build(self): - gmsh.initialize() - gmsh.option.setNumber("General.Terminal", 0) - gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) - - # All ranks create a Gmsh model object - self.pv_model = gmsh.model() - c_x = c_y = 0.2 - r = self.params.domain.cyld_radius - gdim = 2 - mesh_comm = MPI.COMM_WORLD - model_rank = 0 - if mesh_comm.rank == model_rank: - rectangle = self.pv_model.occ.addRectangle(self.params.domain.x_min , - self.params.domain.y_min , - 0, - self.params.domain.x_max , - self.params.domain.y_max , tag=1) - obstacle = self.pv_model.occ.addDisk(c_x, c_y, 0, r, r) - - if mesh_comm.rank == model_rank: - self.pv_model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)]) - self.pv_model.occ.synchronize() - return self.pv_model - - diff --git a/pvopt/geometry/cylinder2d/DomainCreation_temp.py b/pvopt/geometry/cylinder2d/DomainCreation_temp.py deleted file mode 100644 index a93cd099..00000000 --- a/pvopt/geometry/cylinder2d/DomainCreation_temp.py +++ /dev/null @@ -1,679 +0,0 @@ -from dolfinx.io import XDMFFile, gmshio -from dolfinx.fem import VectorFunctionSpace, FunctionSpace -from dolfinx.cpp import mesh as cppmesh - -from mpi4py import MPI -import gmsh -import numpy as np -import os -import time -import ufl -import dolfinx -import meshio -# from pvopt.geometry.panels.domain_creation import * -class FSIDomain: - def __init__(self, params): - - # Get MPI communicators - self.comm = MPI.COMM_WORLD - self.rank = self.comm.Get_rank() - self.num_procs = self.comm.Get_size() - - # Store a full copy of params on this object - self.params = params - - self.x_min_marker = 1 - self.y_min_marker = 2 - self.z_min_marker = 3 - self.x_max_marker = 4 - self.y_max_marker = 5 - self.z_max_marker = 6 - self.internal_surface_marker = 7 - self.fluid_marker = 8 - - def build(self): - test = 1 - if test == 0: - gmsh.initialize() - - L = 2.2 - H = 0.41 - c_x = c_y =0.2 - r = 0.05 - gdim = 2 - mesh_comm = MPI.COMM_WORLD - model_rank = 0 - if mesh_comm.rank == model_rank: - rectangle = gmsh.model.occ.addRectangle(0,0,0, L, H, tag=1) - obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r) - - if mesh_comm.rank == model_rank: - fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)]) - gmsh.model.occ.synchronize() - fluid_marker = 1 - if mesh_comm.rank == model_rank: - volumes = gmsh.model.getEntities(dim=gdim) - assert(len(volumes) == 1) - gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker) - gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid") - - self.inlet_marker, self.outlet_marker, self.wall_marker, self.obstacle_marker = 2, 3, 4, 5 - inflow, outflow, walls, obstacle = [], [], [], [] - if mesh_comm.rank == model_rank: - boundaries = gmsh.model.getBoundary(volumes, oriented=False) - for boundary in boundaries: - center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1]) - if np.allclose(center_of_mass, [0, H/2, 0]): - inflow.append(boundary[1]) - elif np.allclose(center_of_mass, [L, H/2, 0]): - outflow.append(boundary[1]) - elif np.allclose(center_of_mass, [L/2, H, 0]) or np.allclose(center_of_mass, [L/2, 0, 0]): - walls.append(boundary[1]) - else: - obstacle.append(boundary[1]) - gmsh.model.addPhysicalGroup(1, walls, self.wall_marker) - gmsh.model.setPhysicalName(1, self.wall_marker, "Walls") - gmsh.model.addPhysicalGroup(1, inflow, self.inlet_marker) - gmsh.model.setPhysicalName(1, self.inlet_marker, "Inlet") - gmsh.model.addPhysicalGroup(1, outflow, self.outlet_marker) - gmsh.model.setPhysicalName(1, self.outlet_marker, "Outlet") - gmsh.model.addPhysicalGroup(1, obstacle, self.obstacle_marker) - gmsh.model.setPhysicalName(1, self.obstacle_marker, "Obstacle") - - # Create distance field from obstacle. - # Add threshold of mesh sizes based on the distance field - # LcMax - /-------- - # / - # LcMin -o---------/ - # | | | - # Point DistMin DistMax - res_min = r / 3 - if mesh_comm.rank == model_rank: - distance_field = gmsh.model.mesh.field.add("Distance") - gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle) - threshold_field = gmsh.model.mesh.field.add("Threshold") - gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field) - gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min) - gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * H) - gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r) - gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * H) - min_field = gmsh.model.mesh.field.add("Min") - gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field]) - gmsh.model.mesh.field.setAsBackgroundMesh(min_field) - - if mesh_comm.rank == model_rank: - gmsh.option.setNumber("Mesh.Algorithm", 8) - gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) - gmsh.option.setNumber("Mesh.RecombineAll", 1) - gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) - gmsh.model.mesh.generate(gdim) - gmsh.model.mesh.setOrder(2) - gmsh.model.mesh.optimize("Netgen") - - self.msh, self._, self.ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim) - # Specify names for the mesh elements - self.msh.name = "pv_domain" - self._.name = f"{self.msh.name}_cells" - self.ft.name = f"{self.msh.name}_facets" - with XDMFFile(self.msh.comm, "output_mesh_name.pvd", "w") as fp: - fp.write_mesh(self.msh) - - - - else: - self.mesh_comm = MPI.COMM_WORLD - self.model_rank = 0 - self.gdim = 3 - - # Initialize Gmsh options - gmsh.initialize() - gmsh.option.setNumber("General.Terminal", 0) - gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) - - # All ranks create a Gmsh model object - self.pv_model = gmsh.model() - - # Only rank 0 builds the geometry and meshes the domain - # if self.rank == 0: - self._construct_geometry() - self._mark_surfaces() - self._set_length_scales_mod() - - if self.params.fluid.periodic: - self._enforce_periodicity() - - self._generate_mesh() - - # All ranks receive their portion of the mesh from rank 0 (like an MPI scatter) - self.msh, self.mt, self.ft = gmshio.model_to_mesh(self.pv_model, self.comm, 0) - - # Specify names for the mesh elements - self.msh.name = "pv_domain" - self.mt.name = f"{self.msh.name}_cells" - self.ft.name = f"{self.msh.name}_facets" - - def read(self): - if self.rank == 0: - print("Reading the mesh from file ...") - with dolfinx.io.XDMFFile(MPI.COMM_WORLD, self.params.general.output_dir_mesh+"/mesh.xdmf", "r") as xdmf: - self.msh = xdmf.read_mesh(name="Grid") - - self.msh.topology.create_connectivity(self.msh.topology.dim-1, self.msh.topology.dim) - with XDMFFile(MPI.COMM_WORLD,self.params.general.output_dir_mesh+"/mesh_mf.xdmf",'r') as infile: - self.ft = infile.read_meshtags(self.msh, "Grid") - if self.rank == 0: - print("Done.") - - def _construct_geometry(self): - - self.pv_model.add("pv_domain") - self.pv_model.setCurrent("pv_domain") - - # Compute and store some useful geometric quantities - self.x_span = self.params.domain.x_max - self.params.domain.x_min - self.y_span = self.params.domain.y_max - self.params.domain.y_min - self.z_span = self.params.domain.z_max - self.params.domain.z_min - tracker_angle_rad = np.radians(self.params.pv_array.tracker_angle) - - # Create the fluid domain extent - domain = self.pv_model.occ.addBox( - self.params.domain.x_min, - self.params.domain.y_min, - self.params.domain.z_min, - self.x_span, - self.y_span, - self.z_span, - ) - - for panel_id in range(self.params.pv_array.num_rows): - panel_box = self.pv_model.occ.addBox( - -0.5 * self.params.pv_array.panel_length, - 0.0, - -0.5 * self.params.pv_array.panel_thickness, - self.params.pv_array.panel_length, - self.params.pv_array.panel_width, - self.params.pv_array.panel_thickness, - ) - - # Rotate the panel currently centered at (0, 0, 0) - self.pv_model.occ.rotate( - [(3, panel_box)], 0, 0, 0, 0, -1, 0, tracker_angle_rad - ) - - # Translate the panel [panel_loc, 0, elev] - self.pv_model.occ.translate( - [(3, panel_box)], - panel_id * self.params.pv_array.spacing[0], - 0, - self.params.pv_array.elevation, - ) - - # Remove each panel from the overall domain - self.pv_model.occ.cut([(3, domain)], [(3, panel_box)]) - - self.pv_model.occ.synchronize() - - def _construct_geometry_mod(self): - - if self.mesh_comm.rank == self.model_rank: - # Initialize Gmsh options - gmsh.initialize() - gmsh.option.setNumber("General.Terminal", 0) - gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) - gmsh.model.add("pv_domain") - gmsh.model.setCurrent("pv_domain") - - # Compute and store some useful geometric quantities - self.x_span = self.params.domain.x_max - self.params.domain.x_min - self.y_span = self.params.domain.y_max - self.params.domain.y_min - self.z_span = self.params.domain.z_max - self.params.domain.z_min - tracker_angle_rad = np.radians(self.params.pv_array.tracker_angle) - - # Create the fluid domain extent - domain = gmsh.model.occ.addBox( - self.params.domain.x_min, - self.params.domain.y_min, - self.params.domain.z_min, - self.x_span, - self.y_span, - self.z_span, - ) - - for panel_id in range(self.params.pv_array.num_rows): - panel_box = gmsh.model.occ.addBox( - -0.5 * self.params.pv_array.panel_length, - 0.0, - -0.5 * self.params.pv_array.panel_thickness, - self.params.pv_array.panel_length, - self.params.pv_array.panel_width, - self.params.pv_array.panel_thickness, - ) - - # Rotate the panel currently centered at (0, 0, 0) - gmsh.model.occ.rotate( - [(3, panel_box)], 0, 0, 0, 0, -1, 0, tracker_angle_rad - ) - - # Translate the panel [panel_loc, 0, elev] - gmsh.model.occ.translate( - [(3, panel_box)], - panel_id * self.params.pv_array.spacing[0], - 0, - self.params.pv_array.elevation, - ) - - # Remove each panel from the overall domain - gmsh.model.occ.cut([(3, domain)], [(3, panel_box)]) - - gmsh.model.occ.synchronize() - - def _mark_surfaces(self): - - # Loop through all surfaces to find periodic tags - surf_ids = self.pv_model.occ.getEntities(2) - - self.dom_tags = {} - - for surf in surf_ids: - tag = surf[1] - - com = self.pv_model.occ.getCenterOfMass(2, tag) - - if np.isclose(com[0], self.params.domain.x_min): - self.dom_tags["left"] = [tag] - - elif np.allclose(com[0], self.params.domain.x_max): - self.dom_tags["right"] = [tag] - - elif np.allclose(com[1], self.params.domain.y_min): - self.dom_tags["front"] = [tag] - - elif np.allclose(com[1], self.params.domain.y_max): - self.dom_tags["back"] = [tag] - - elif np.allclose(com[2], self.params.domain.z_min): - self.dom_tags["bottom"] = [tag] - - elif np.allclose(com[2], self.params.domain.z_max): - self.dom_tags["top"] = [tag] - - else: - if "panel_surface" in self.dom_tags: - self.dom_tags["panel_surface"].append(tag) - else: - self.dom_tags["panel_surface"] = [tag] - - - self.pv_model.addPhysicalGroup(3, [1], 11) - self.pv_model.setPhysicalName(3, 11, "fluid") - - self.pv_model.addPhysicalGroup(2, self.dom_tags["left"], self.x_min_marker) - self.pv_model.setPhysicalName(2, self.x_min_marker, "left") - self.pv_model.addPhysicalGroup(2, self.dom_tags["right"], self.x_max_marker) - self.pv_model.setPhysicalName(2, self.x_max_marker, "right") - self.pv_model.addPhysicalGroup(2, self.dom_tags["front"], self.y_min_marker) - self.pv_model.setPhysicalName(2, self.y_min_marker, "front") - self.pv_model.addPhysicalGroup(2, self.dom_tags["back"], self.y_max_marker) - self.pv_model.setPhysicalName(2, self.y_max_marker, "back") - self.pv_model.addPhysicalGroup(2, self.dom_tags["bottom"], self.z_min_marker) - self.pv_model.setPhysicalName(2, self.z_min_marker, "bottom") - self.pv_model.addPhysicalGroup(2, self.dom_tags["top"], self.z_max_marker) - self.pv_model.setPhysicalName(2, self.z_max_marker, "top") - self.pv_model.addPhysicalGroup(2, self.dom_tags["panel_surface"], self.internal_surface_marker) - self.pv_model.setPhysicalName(2, self.internal_surface_marker, "panel_surface") - - def _mark_surfaces_mod(self): - - fluid_marker = 1 - if self.mesh_comm.rank == self.model_rank: - volumes = gmsh.model.getEntities(dim=3) - assert(len(volumes) == 1) - gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker) - gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid") - - self.inlet_marker, self.outlet_marker, self.wall_z_marker, self.obstacle_marker, self.wall_y_marker = 2, 3, 4, 5, 6 - self.inflow, self.outflow, self.walls_z, self.walls_y, self.obstacle = [], [], [], [],[] - inflow_b_found,outflow_b_found,wallz_b_found, wally_b_found, obstacle_b_found = False,False,False,False,False - if self.mesh_comm.rank == self.model_rank: - boundaries = gmsh.model.getBoundary(volumes, oriented=False) - for boundary in boundaries: - center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1]) - center_x = self.params.domain.x_min + abs(self.params.domain.x_max-self.params.domain.x_min)/2 - center_y = self.params.domain.y_min + abs(self.params.domain.y_max-self.params.domain.y_min)/2 - center_z = self.params.domain.z_min + abs(self.params.domain.z_max-self.params.domain.z_min)/2 - # if np.allclose(center_of_mass, [self.params.domain.x_min, center_y, center_z]): - # self.inflow.append(boundary[1]) - # inflow_b_found = True - if np.allclose(center_of_mass[0], [self.params.domain.x_min]): - self.inflow.append(boundary[1]) - inflow_b_found = True - elif np.allclose(center_of_mass[0], [self.params.domain.x_max]): - self.outflow.append(boundary[1]) - outflow_b_found = True - elif np.allclose(center_of_mass[1], [self.params.domain.y_max]) or np.allclose(center_of_mass[1], [self.params.domain.y_min]): - self.walls_y.append(boundary[1]) - wally_b_found = True - elif np.allclose(center_of_mass[2], [self.params.domain.z_max]) or np.allclose(center_of_mass[2], [self.params.domain.z_min]): - self.walls_z.append(boundary[1]) - wallz_b_found = True - else: - self.obstacle.append(boundary[1]) - obstacle_b_found = True - - if (inflow_b_found and outflow_b_found and wallz_b_found and wally_b_found and obstacle_b_found) == True: - print("all boundaries are found ") - else: - print("boundary missing") - quit - gmsh.model.addPhysicalGroup(1, self.walls_z, self.wall_z_marker) - gmsh.model.addPhysicalGroup(1, self.walls_y, self.wall_y_marker) - gmsh.model.setPhysicalName(1, self.wall_z_marker, "Walls_z") - gmsh.model.setPhysicalName(1, self.wall_y_marker, "Walls_y") - gmsh.model.addPhysicalGroup(1, self.inflow, self.inlet_marker) - gmsh.model.setPhysicalName(1, self.inlet_marker, "Inlet") - gmsh.model.addPhysicalGroup(1, self.outflow, self.outlet_marker) - gmsh.model.setPhysicalName(1, self.outlet_marker, "Outlet") - gmsh.model.addPhysicalGroup(1, self.obstacle, self.obstacle_marker) - gmsh.model.setPhysicalName(1, self.obstacle_marker, "Obstacle") - - def _set_length_scales(self): - - # Set size scales for the mesh - # eps = 0.1 - - # Set the mesh size at each point on the panel - # panel_bbox = self.pv_model.getEntitiesInBoundingBox(self.params.domain.x_min+eps, self.params.domain.y_min-eps, self.params.domain.z_min+eps, - # self.params.domain.x_max-eps, self.params.domain.y_max+eps, self.params.domain.z_max-eps) - # self.pv_model.mesh.setSize(panel_bbox, self.params.domain.l_char) - - # Define a distance field from the bottom of the domain - self.pv_model.mesh.field.add("Distance", 1) - self.pv_model.mesh.field.setNumbers(1, "FacesList", self.obstacle) - # self.pv_model.mesh.field.setNumbers(1, "FacesList", self.dom_tags["bottom"]) - - self.pv_model.mesh.field.add("Threshold", 2) - self.pv_model.mesh.field.setNumber(2, "IField", 1) - self.pv_model.mesh.field.setNumber(2, "LcMin", 2.0 * self.params.domain.l_char) - # self.pv_model.mesh.field.setNumber(2, 'LcMin', self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(2, "LcMax", 6.0 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(2, "DistMin", 4.5) - self.pv_model.mesh.field.setNumber(2, "DistMax", 1.0 * self.params.domain.z_max) - - # Define a distance field from the immersed panels - self.pv_model.mesh.field.add("Distance", 3) - self.pv_model.mesh.field.setNumbers( - 3, "FacesList", self.obstacle - ) - - self.pv_model.mesh.field.add("Threshold", 4) - self.pv_model.mesh.field.setNumber(4, "IField", 3) - self.pv_model.mesh.field.setNumber(4, "LcMin", 0.2*self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(4, "LcMax", 6.0 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(4, "DistMin", 0.5) - self.pv_model.mesh.field.setNumber(4, "DistMax", 0.6 * self.params.domain.z_max) - - self.pv_model.mesh.field.add("Min", 5) - self.pv_model.mesh.field.setNumbers(5, "FieldsList", [2, 4]) - - threshold_field = gmsh.model.mesh.field.add("Threshold") - distance_field = gmsh.model.mesh.field.add("Distance") - self.pv_model.mesh.field.setNumbers(distance_field, "EdgesList", self.obstacle) - self.pv_model.mesh.field.setNumber(threshold_field, "IField", distance_field) - self.pv_model.mesh.field.setNumber(threshold_field, "LcMin", self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * self.params.domain.z_max) - self.pv_model.mesh.field.setNumber(threshold_field, "DistMin", self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(threshold_field, "DistMax", 2 * self.params.domain.z_max) - min_field = gmsh.model.mesh.field.add("Min") - self.pv_model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field]) - self.pv_model.mesh.field.setAsBackgroundMesh(min_field) - - # self.pv_model.mesh.field.setAsBackgroundMesh(5) - - # self.pv_model.mesh.field.setAsBackgroundMesh(5) - - def _set_length_scales_mod(self): - res_min = self.params.domain.l_char - if self.mesh_comm.rank == self.model_rank: - # Define a distance field from the bottom of the domain - self.pv_model.mesh.field.add("Distance", 1) - self.pv_model.mesh.field.setNumbers(1, "FacesList", self.dom_tags["bottom"]) - - self.pv_model.mesh.field.add("Threshold", 2) - self.pv_model.mesh.field.setNumber(2, "IField", 1) - self.pv_model.mesh.field.setNumber(2, "LcMin", 2.0 * self.params.domain.l_char) - # self.pv_model.mesh.field.setNumber(2, 'LcMin', self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(2, "LcMax", 6.0 * self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(2, "DistMin", 4.5) - self.pv_model.mesh.field.setNumber(2, "DistMax", 1.0 * self.params.domain.z_max) - - # Define a distance field from the immersed panels - self.pv_model.mesh.field.add("Distance", 3) - self.pv_model.mesh.field.setNumbers( - 3, "FacesList", self.dom_tags["panel_surface"] - ) - - self.pv_model.mesh.field.add("Threshold", 4) - self.pv_model.mesh.field.setNumber(4, "IField", 3) - - self.pv_model.mesh.field.setNumber(4, "LcMin", self.params.domain.l_char*0.5) - self.pv_model.mesh.field.setNumber(4, "LcMax", 5*self.params.domain.l_char) - self.pv_model.mesh.field.setNumber(4, "DistMin", 0.5) - self.pv_model.mesh.field.setNumber(4, "DistMax", 0.6 * self.params.domain.z_max) - - self.pv_model.mesh.field.add("Min", 5) - self.pv_model.mesh.field.setNumbers(5, "FieldsList", [2, 4]) - self.pv_model.mesh.field.setAsBackgroundMesh(5) - - - # gmsh.model.mesh.field.add("Distance", 1) - # gmsh.model.mesh.field.setNumbers(1, "FacesList", self.dom_tags["panel_surface"]) - # r = res_min - # resolution = r - # gmsh.model.mesh.field.add("Threshold", 2) - # gmsh.model.mesh.field.setNumber(2, "IField", 1) - # gmsh.model.mesh.field.setNumber(2, "LcMin", resolution) - # gmsh.model.mesh.field.setNumber(2, "LcMax", 2*resolution) - # gmsh.model.mesh.field.setNumber(2, "DistMin", 0.5*r) - # gmsh.model.mesh.field.setNumber(2, "DistMax", r) - - # # We add a similar threshold at the inlet to ensure fully resolved inlet flow - - # gmsh.model.mesh.field.add("Distance", 3) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['back'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['front'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['top'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['bottom'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['left'] ) - # gmsh.model.mesh.field.setNumbers(3, "FacesList", self.dom_tags['right'] ) - # gmsh.model.mesh.field.add("Threshold", 4) - # gmsh.model.mesh.field.setNumber(4, "IField", 3) - # gmsh.model.mesh.field.setNumber(4, "LcMin", 5*resolution) - # gmsh.model.mesh.field.setNumber(4, "LcMax", 10*resolution) - # gmsh.model.mesh.field.setNumber(4, "DistMin", 0.1) - # gmsh.model.mesh.field.setNumber(4, "DistMax", 0.5) - - # # We combine these fields by using the minimum field - # gmsh.model.mesh.field.add("Min", 5) - # gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2, 4]) - # gmsh.model.mesh.field.setAsBackgroundMesh(5) - - - if self.mesh_comm.rank == self.model_rank: - gmsh.option.setNumber("Mesh.Algorithm", 5) - gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) - # gmsh.option.setNumber("Mesh.RecombineAll", 1) - gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) - gmsh.model.mesh.generate(3) - gmsh.model.mesh.setOrder(2) - gmsh.model.mesh.optimize("Netgen") - - - - - # self.msh, self._, self.ft = gmshio.model_to_mesh(gmsh.model, self.mesh_comm, self.model_rank, gdim=2) - # Set size scales for the mesh - # eps = 0.1 - - # Set the mesh size at each point on the panel - # panel_bbox = self.pv_model.getEntitiesInBoundingBox(self.params.domain.x_min+eps, self.params.domain.y_min-eps, self.params.domain.z_min+eps, - # self.params.domain.x_max-eps, self.params.domain.y_max+eps, self.params.domain.z_max-eps) - # self.pv_model.mesh.setSize(panel_bbox, self.params.domain.l_char) - - # Define a distance field from the bottom of the domain - # self.pv_model.mesh.field.add("Distance", 1) - # self.pv_model.mesh.field.setNumbers(1, "FacesList", self.dom_tags["bottom"]) - - # self.pv_model.mesh.field.add("Threshold", 2) - # self.pv_model.mesh.field.setNumber(2, "IField", 1) - # self.pv_model.mesh.field.setNumber(2, "LcMin", 2.0 * self.params.domain.l_char) - # # self.pv_model.mesh.field.setNumber(2, 'LcMin', self.params.domain.l_char) - # self.pv_model.mesh.field.setNumber(2, "LcMax", 6.0 * self.params.domain.l_char) - # self.pv_model.mesh.field.setNumber(2, "DistMin", 4.5) - # self.pv_model.mesh.field.setNumber(2, "DistMax", 1.0 * self.params.domain.z_max) - - # # Define a distance field from the immersed panels - # self.pv_model.mesh.field.add("Distance", 3) - # self.pv_model.mesh.field.setNumbers( - # 3, "FacesList", self.dom_tags["panel_surface"] - # ) - - # self.pv_model.mesh.field.add("Threshold", 4) - # self.pv_model.mesh.field.setNumber(4, "IField", 3) - # self.pv_model.mesh.field.setNumber(4, "LcMin", self.params.domain.l_char) - # self.pv_model.mesh.field.setNumber(4, "LcMax", 6.0 * self.params.domain.l_char) - # self.pv_model.mesh.field.setNumber(4, "DistMin", 0.5) - # self.pv_model.mesh.field.setNumber(4, "DistMax", 0.6 * self.params.domain.z_max) - - # self.pv_model.mesh.field.add("Min", 5) - # self.pv_model.mesh.field.setNumbers(5, "FieldsList", [2, 4]) - # self.pv_model.mesh.field.setAsBackgroundMesh(5) - - # self.pv_model.mesh.field.setAsBackgroundMesh(5) - - def _enforce_periodicity(self): - - # TODO: Make this a generic mapping depending on which walls are marked for peridic BCs - # TODO: Copy code to enforce periodicity from old generate_and_convert_3d_meshes.py - - # Mark the front/back walls for periodic mapping - front_back_translation = [ - 1, - 0, - 0, - 0, - 0, - 1, - 0, - self.y_span, - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - self.pv_model.mesh.setPeriodic( - 2, self.dom_tags["back"], self.dom_tags["front"], front_back_translation - ) - - # Mark the left/right walls for periodic mapping - left_right_translation = [ - 1, - 0, - 0, - self.x_span, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - 0, - 0, - 0, - 0, - 1, - ] - - self.pv_model.mesh.setPeriodic( - 2, self.dom_tags["right"], self.dom_tags["left"], left_right_translation - ) - - def _generate_mesh(self): - if self.rank == 0: - print("Starting mesh generation... ", end="") - - # Generate the mesh - tic = time.time() - gmsh.option.setNumber("Mesh.Algorithm", 8) - gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) - # gmsh.option.setNumber("Mesh.RecombineAll", 1) - gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) - self.pv_model.mesh.generate(3) - self.pv_model.mesh.setOrder(2) - self.pv_model.mesh.optimize("Netgen") - self.pv_model.mesh.generate(3) - toc = time.time() - if self.rank == 0: - print("Finished.") - print(f"Total meshing time = {toc-tic:.1f} s") - - def write_mesh_file(self): - if self.rank == 0: - # Save the *.msh file and *.vtk file (the latter is for visualization only) - print('Writing Mesh to %s... ' % (self.params.general.output_dir_mesh), end='') - gmsh.write('%s/mesh.msh' % (self.params.general.output_dir_mesh)) - gmsh.write('%s/mesh.vtk' % (self.params.general.output_dir_mesh)) - def create_mesh(mesh, clean_points, cell_type): - cells = mesh.get_cells_type(cell_type) - cell_data = mesh.get_cell_data("gmsh:physical", cell_type) - - out_mesh = meshio.Mesh(points=clean_points, cells={ - cell_type: cells}, cell_data={"name_to_read": [cell_data]}) - return out_mesh - - mesh_from_file = meshio.read(f'{self.params.general.output_dir_mesh}/mesh.msh') - pts = mesh_from_file.points - tetra_mesh = create_mesh(mesh_from_file, pts, "tetra") - tri_mesh = create_mesh(mesh_from_file, pts, "triangle") - - meshio.write(f'{self.params.general.output_dir_mesh}/mesh.xdmf', tetra_mesh) - meshio.write(f'{self.params.general.output_dir_mesh}/mesh_mf.xdmf', tri_mesh) - print("Done.") - - def test_mesh_functionspace(self): - - P2 = ufl.VectorElement("Lagrange", self.msh.ufl_cell(), 2) - P1 = ufl.FiniteElement("Lagrange", self.msh.ufl_cell(), 1) - V = FunctionSpace(self.msh, P2) - Q = FunctionSpace(self.msh, P1) - - local_rangeV = V.dofmap.index_map.local_range - dofsV = np.arange(*local_rangeV) - - local_rangeQ = Q.dofmap.index_map.local_range - dofsQ = np.arange(*local_rangeQ) - - self.ndim = self.msh.topology.dim - - # coords = self.mesh.coordinates() - - nodes_dim = 0 - self.msh.topology.create_connectivity(nodes_dim, 0) - num_nodes_owned_by_proc = self.msh.topology.index_map(nodes_dim).size_local - geometry_entitites = cppmesh.entities_to_geometry( - self.msh, - nodes_dim, - np.arange(num_nodes_owned_by_proc, dtype=np.int32), - False, - ) - points = self.msh.geometry.x - - coords = points[:] - - print(f"Rank {self.rank} owns {num_nodes_owned_by_proc} nodes\n{coords}") - diff --git a/pvopt/geometry/cylinder2d/__pycache__/DomainCreation.cpython-310.pyc b/pvopt/geometry/cylinder2d/__pycache__/DomainCreation.cpython-310.pyc deleted file mode 100644 index 80c1f5bd..00000000 Binary files a/pvopt/geometry/cylinder2d/__pycache__/DomainCreation.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/cylinder2d/__pycache__/__init__.cpython-310.pyc b/pvopt/geometry/cylinder2d/__pycache__/__init__.cpython-310.pyc deleted file mode 100644 index c98d444e..00000000 Binary files a/pvopt/geometry/cylinder2d/__pycache__/__init__.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/cylinder3d/.DS_Store b/pvopt/geometry/cylinder3d/.DS_Store deleted file mode 100644 index 3606724a..00000000 Binary files a/pvopt/geometry/cylinder3d/.DS_Store and /dev/null differ diff --git a/pvopt/geometry/cylinder3d/DomainCreation.py b/pvopt/geometry/cylinder3d/DomainCreation.py deleted file mode 100644 index 8a212414..00000000 --- a/pvopt/geometry/cylinder3d/DomainCreation.py +++ /dev/null @@ -1,60 +0,0 @@ -from dolfinx.io import XDMFFile, gmshio -from dolfinx.fem import VectorFunctionSpace, FunctionSpace -from dolfinx.cpp import mesh as cppmesh -from mpi4py import MPI -from pvopt.geometry.template.TemplateDomainCreation import TemplateDomainCreation - -import gmsh -import numpy as np -import os -import time - - -class DomainCreation(TemplateDomainCreation): - def __init__(self, params): - super().__init__(params) - - def build(self): - - # Compute and store some useful geometric quantities - self.x_span = self.params.domain.x_max - self.params.domain.x_min - self.y_span = self.params.domain.y_max - self.params.domain.y_min - self.z_span = self.params.domain.z_max - self.params.domain.z_min - - # Create the fluid domain extent - domain = self.gmsh_model.occ.addBox( - self.params.domain.x_min, - self.params.domain.y_min, - self.params.domain.z_min, - self.x_span, - self.y_span, - self.z_span, - ) - - self.cyld_radius = 0.05 - height = 100.0 - - cylinder = self.gmsh_model.occ.addCylinder( - 0.5, -height / 2.0, 0.2, 0.0, height, 0.0, self.cyld_radius - ) - - # Remove each panel from the overall domain - self.gmsh_model.occ.cut([(3, domain)], [(3, cylinder)]) - - self.gmsh_model.occ.synchronize() - - def set_length_scales(self): - - # Define a distance field from the cylinder - self.gmsh_model.mesh.field.add("Distance", 1) - self.gmsh_model.mesh.field.setNumbers(1, "FacesList", self.dom_tags["internal_surface"]) - - self.gmsh_model.mesh.field.add("Threshold", 2) - self.gmsh_model.mesh.field.setNumber(2, "IField", 1) - self.gmsh_model.mesh.field.setNumber(2, "LcMin", self.params.domain.l_char) - self.gmsh_model.mesh.field.setNumber(2, "LcMax", 6.0 * self.params.domain.l_char) - self.gmsh_model.mesh.field.setNumber(2, "DistMin", 2.0 * self.cyld_radius) - self.gmsh_model.mesh.field.setNumber(2, "DistMax", 4.0*self.cyld_radius) - - self.gmsh_model.mesh.field.setAsBackgroundMesh(2) - diff --git a/pvopt/geometry/cylinder3d/__pycache__/DomainCreation.cpython-310.pyc b/pvopt/geometry/cylinder3d/__pycache__/DomainCreation.cpython-310.pyc deleted file mode 100644 index 1307a8d7..00000000 Binary files a/pvopt/geometry/cylinder3d/__pycache__/DomainCreation.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/cylinder3d/__pycache__/__init__.cpython-310.pyc b/pvopt/geometry/cylinder3d/__pycache__/__init__.cpython-310.pyc deleted file mode 100644 index 50c46856..00000000 Binary files a/pvopt/geometry/cylinder3d/__pycache__/__init__.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/panels/.DS_Store b/pvopt/geometry/panels/.DS_Store deleted file mode 100644 index fd74a931..00000000 Binary files a/pvopt/geometry/panels/.DS_Store and /dev/null differ diff --git a/pvopt/geometry/panels/DomainCreation.py b/pvopt/geometry/panels/DomainCreation.py deleted file mode 100644 index d4124793..00000000 --- a/pvopt/geometry/panels/DomainCreation.py +++ /dev/null @@ -1,67 +0,0 @@ -from dolfinx.io import XDMFFile, gmshio -from dolfinx.fem import VectorFunctionSpace, FunctionSpace -from dolfinx.cpp import mesh as cppmesh -from mpi4py import MPI -from pvopt.geometry.template.TemplateDomainCreation import TemplateDomainCreation - -import gmsh -import numpy as np -import os -import time -import ufl -import dolfinx - -# import meshio - - -class DomainCreation(TemplateDomainCreation): - def __init__(self, params): - super().__init__(params) - - def build(self): - - # Compute and store some useful geometric quantities - self.x_span = self.params.domain.x_max - self.params.domain.x_min - self.y_span = self.params.domain.y_max - self.params.domain.y_min - self.z_span = self.params.domain.z_max - self.params.domain.z_min - tracker_angle_rad = np.radians(self.params.pv_array.tracker_angle) - - # Create the fluid domain extent - domain = self.gmsh_model.occ.addBox( - self.params.domain.x_min, - self.params.domain.y_min, - self.params.domain.z_min, - self.x_span, - self.y_span, - self.z_span, - ) - - for panel_id in range(self.params.pv_array.num_rows): - panel_box = self.gmsh_model.occ.addBox( - -0.5 * self.params.pv_array.panel_length, - 0.0, - -0.5 * self.params.pv_array.panel_thickness, - self.params.pv_array.panel_length, - self.params.pv_array.panel_width, - self.params.pv_array.panel_thickness, - ) - - # Rotate the panel currently centered at (0, 0, 0) - self.gmsh_model.occ.rotate( - [(3, panel_box)], 0, 0, 0, 0, -1, 0, tracker_angle_rad - ) - - # Translate the panel [panel_loc, 0, elev] - self.gmsh_model.occ.translate( - [(3, panel_box)], - panel_id * self.params.pv_array.spacing[0], - 0, - self.params.pv_array.elevation, - ) - - # Remove each panel from the overall domain - self.gmsh_model.occ.cut([(3, domain)], [(3, panel_box)]) - - self.gmsh_model.occ.synchronize() - return self.gmsh_model - diff --git a/pvopt/geometry/panels/__pycache__/DomainCreation.cpython-310.pyc b/pvopt/geometry/panels/__pycache__/DomainCreation.cpython-310.pyc deleted file mode 100644 index fb7d0495..00000000 Binary files a/pvopt/geometry/panels/__pycache__/DomainCreation.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/panels/__pycache__/__init__.cpython-310.pyc b/pvopt/geometry/panels/__pycache__/__init__.cpython-310.pyc deleted file mode 100644 index 5a8b96b7..00000000 Binary files a/pvopt/geometry/panels/__pycache__/__init__.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/panels/__pycache__/domain_creation.cpython-310.pyc b/pvopt/geometry/panels/__pycache__/domain_creation.cpython-310.pyc deleted file mode 100644 index 184edc33..00000000 Binary files a/pvopt/geometry/panels/__pycache__/domain_creation.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/panels2d/DomainCreation.py b/pvopt/geometry/panels2d/DomainCreation.py deleted file mode 100644 index 67d8b315..00000000 --- a/pvopt/geometry/panels2d/DomainCreation.py +++ /dev/null @@ -1,99 +0,0 @@ -from dolfinx.io import XDMFFile, gmshio -from dolfinx.fem import VectorFunctionSpace, FunctionSpace -from dolfinx.cpp import mesh as cppmesh -from mpi4py import MPI -import gmsh -import numpy as np -import os -import time -import ufl -import dolfinx -import meshio - -class DomainCreation: - def __init__(self, params): - - # Get MPI communicators - self.comm = MPI.COMM_WORLD - self.rank = self.comm.Get_rank() - self.num_procs = self.comm.Get_size() - - # Store a full copy of params on this object - self.params = params - - self.x_min_marker = 1 - self.y_min_marker = 2 - self.z_min_marker = 3 - self.x_max_marker = 4 - self.y_max_marker = 5 - self.z_max_marker = 6 - self.internal_surface_marker = 7 - self.fluid_marker = 8 - - def build(self): - - self.mesh_comm = MPI.COMM_WORLD - self.model_rank = 0 - self.gdim = 3 - - # Initialize Gmsh options - gmsh.initialize() - gmsh.option.setNumber("General.Terminal", 0) - gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) - - # All ranks create a Gmsh model object - self.pv_model = gmsh.model() - - # Only rank 0 builds the geometry and meshes the domain - # if self.rank == 0: - - self.pv_model.add("pv_domain") - self.pv_model.setCurrent("pv_domain") - - # Compute and store some useful geometric quantities - self.x_span = self.params.domain.x_max - self.params.domain.x_min - self.y_span = self.params.domain.y_max - self.params.domain.y_min - # self.z_span = self.params.domain.z_max - self.params.domain.z_min - tracker_angle_rad = np.radians(self.params.pv_array.tracker_angle) - - # Create the fluid domain extent - domain = self.pv_model.occ.addRectangle( - self.params.domain.x_min, - self.params.domain.y_min, - 0,#self.params.domain.z_min, - self.x_span, - self.y_span - #self.z_span, - ) - - for panel_id in range(self.params.pv_array.num_rows): - panel_box = self.pv_model.occ.addRectangle( - -0.5 * self.params.pv_array.panel_length, - -0.5 * self.params.pv_array.panel_thickness, - 0, - self.params.pv_array.panel_length, - self.params.pv_array.panel_thickness - ) - - # Rotate the panel currently centered at (0, 0, 0) - self.pv_model.occ.rotate( - [(2, panel_box)], 0, 0, 0, 0, 0, 1, tracker_angle_rad - ) - - # Translate the panel [panel_loc, 0, elev] - self.pv_model.occ.translate( - [(2, panel_box)], - panel_id * self.params.pv_array.spacing[0], - self.params.pv_array.elevation, - 0 - ) - - # Remove each panel from the overall domain - self.pv_model.occ.cut([(2, domain)], [(2, panel_box)]) - - - - - self.pv_model.occ.synchronize() - return self.pv_model - \ No newline at end of file diff --git a/pvopt/geometry/panels2d/__pycache__/DomainCreation.cpython-310.pyc b/pvopt/geometry/panels2d/__pycache__/DomainCreation.cpython-310.pyc deleted file mode 100644 index 800cd538..00000000 Binary files a/pvopt/geometry/panels2d/__pycache__/DomainCreation.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/panels2d/__pycache__/__init__.cpython-310.pyc b/pvopt/geometry/panels2d/__pycache__/__init__.cpython-310.pyc deleted file mode 100644 index 9e1de0b7..00000000 Binary files a/pvopt/geometry/panels2d/__pycache__/__init__.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/template/TemplateDomainCreation.py b/pvopt/geometry/template/TemplateDomainCreation.py deleted file mode 100644 index a13e71d5..00000000 --- a/pvopt/geometry/template/TemplateDomainCreation.py +++ /dev/null @@ -1,108 +0,0 @@ -from dolfinx.io import XDMFFile, gmshio -from dolfinx.fem import VectorFunctionSpace, FunctionSpace -from dolfinx.cpp import mesh as cppmesh -from mpi4py import MPI -import gmsh -import numpy as np -import os -import time - - -class TemplateDomainCreation: - def __init__(self, params): - - # Store a full copy of params on this object - self.params = params - - # Get MPI communicators - self.comm = self.params.comm - self.rank = self.params.rank - self.num_procs = self.params.num_procs - - self.x_min_marker = 1 - self.x_max_marker = 2 - self.y_min_marker = 3 - self.y_max_marker = 4 - self.z_min_marker = 5 - self.z_max_marker = 6 - self.internal_surface_marker = 7 - self.fluid_marker = 8 - self.structure_marker = 8 - - # Initialize Gmsh options - gmsh.initialize() - gmsh.option.setNumber("General.Terminal", 0) - gmsh.option.setNumber("Mesh.MshFileVersion", 2.2) - - # All ranks create a Gmsh model object - self.gmsh_model = gmsh.model() - self.gmsh_model.add("domain") - self.gmsh_model.setCurrent("domain") - - def build(self): - pass - - def set_length_scales(self): - - if self.rank == 0: - all_pts = self.gmsh_model.occ.getEntities(0) - self.gmsh_model.mesh.setSize(all_pts, self.params.domain.l_char) - - def mark_surfaces(self): - """Creates boundary tags using gmsh""" - # Loop through all surfaces to find periodic tags - surf_ids = self.gmsh_model.occ.getEntities(2) - - self.dom_tags = {} - - for surf in surf_ids: - tag = surf[1] - - com = self.gmsh_model.occ.getCenterOfMass(2, tag) - - if np.isclose(com[0], self.params.domain.x_min): - self.dom_tags["x_min"] = [tag] - - elif np.allclose(com[0], self.params.domain.x_max): - self.dom_tags["x_max"] = [tag] - - elif np.allclose(com[1], self.params.domain.y_min): - self.dom_tags["y_min"] = [tag] - - elif np.allclose(com[1], self.params.domain.y_max): - self.dom_tags["y_max"] = [tag] - - elif np.allclose(com[2], self.params.domain.z_min): - self.dom_tags["z_min"] = [tag] - - elif np.allclose(com[2], self.params.domain.z_max): - self.dom_tags["z_max"] = [tag] - - else: - if "internal_surface" in self.dom_tags: - self.dom_tags["internal_surface"].append(tag) - else: - self.dom_tags["internal_surface"] = [tag] - print(self.dom_tags) - - self.gmsh_model.addPhysicalGroup(3, [1], self.fluid_marker) - self.gmsh_model.setPhysicalName(3, self.fluid_marker, "fluid") - - self.gmsh_model.addPhysicalGroup(2, self.dom_tags["x_min"], self.x_min_marker) - self.gmsh_model.setPhysicalName(2, self.x_min_marker, "x_min") - self.gmsh_model.addPhysicalGroup(2, self.dom_tags["x_max"], self.x_max_marker) - self.gmsh_model.setPhysicalName(2, self.x_max_marker, "x_max") - self.gmsh_model.addPhysicalGroup(2, self.dom_tags["y_min"], self.y_min_marker) - self.gmsh_model.setPhysicalName(2, self.y_min_marker, "y_min") - self.gmsh_model.addPhysicalGroup(2, self.dom_tags["y_max"], self.y_max_marker) - self.gmsh_model.setPhysicalName(2, self.y_max_marker, "y_max") - self.gmsh_model.addPhysicalGroup(2, self.dom_tags["z_min"], self.z_min_marker) - self.gmsh_model.setPhysicalName(2, self.z_min_marker, "z_min") - self.gmsh_model.addPhysicalGroup(2, self.dom_tags["z_max"], self.z_max_marker) - self.gmsh_model.setPhysicalName(2, self.z_max_marker, "z_max") - self.gmsh_model.addPhysicalGroup( - 2, self.dom_tags["internal_surface"], self.internal_surface_marker - ) - self.gmsh_model.setPhysicalName( - 2, self.internal_surface_marker, "internal_surface" - ) diff --git a/pvopt/geometry/template/__pycache__/TemplateDomainCreation.cpython-310.pyc b/pvopt/geometry/template/__pycache__/TemplateDomainCreation.cpython-310.pyc deleted file mode 100644 index ea44cbc1..00000000 Binary files a/pvopt/geometry/template/__pycache__/TemplateDomainCreation.cpython-310.pyc and /dev/null differ diff --git a/pvopt/geometry/template/__pycache__/__init__.cpython-310.pyc b/pvopt/geometry/template/__pycache__/__init__.cpython-310.pyc deleted file mode 100644 index 463d4226..00000000 Binary files a/pvopt/geometry/template/__pycache__/__init__.cpython-310.pyc and /dev/null differ diff --git a/pvopt/save_all_net_traction.py b/pvopt/save_all_net_traction.py deleted file mode 100644 index 70d8e8c2..00000000 --- a/pvopt/save_all_net_traction.py +++ /dev/null @@ -1,65 +0,0 @@ -import numpy as np -import glob -import sys -import time -from mpi4py import MPI - - -def get_net_traction(path_to_files): - - try: - fp = open("%s/sim_log.txt" % (path_to_files), "r") - except: - raise ValueError("Cannot find %s" % (path_to_files)) - else: - fp.close() - - csv_list = sorted(glob.glob("%s/time_series/*.csv" % (path_to_files))) - - data_out = [] - - for csv in csv_list: - raw_data = np.genfromtxt(csv, delimiter=",", skip_header=1) - p_top = raw_data[:, 3] - p_bot = raw_data[:, 7] - - data_out.append(p_top - p_bot) - - data_out = np.array(data_out) - - output_filename = "%s/net_traction" % (path_to_files) - np.savetxt(output_filename + ".csv", data_out, delimiter=",") - np.save(output_filename + ".npy", data_out) - - -comm = MPI.COMM_WORLD -rank = comm.Get_rank() -num_procs = comm.Get_size() - -try: - prefix = sys.argv[1] -except: - prefix = "output/param_study" - -run_list_global = sorted(glob.glob("%s/theta_*/uref_*" % (prefix))) -run_list_local = np.array_split(run_list_global, num_procs)[rank] - -if rank == 0: - for k in run_list_global: - print(k) - print("================================================================") - print("Found %d runs to process" % (len(run_list_global))) - print( - "Dividing work between %d cores (~%d runs per core)" - % (num_procs, len(run_list_local)) - ) - print("================================================================") - -for path_to_files in run_list_local: - tic = time.time() - get_net_traction(path_to_files) - toc = time.time() - print( - "Rank %d finished processing %s in %f seconds." - % (rank, path_to_files, toc - tic) - ) diff --git a/pytest.ini b/pytest.ini new file mode 100644 index 00000000..93ad3d1a --- /dev/null +++ b/pytest.ini @@ -0,0 +1,4 @@ +[pytest] +markers = + unit: A short test that relies mostly on quick pass/fail outcomes. + regression: A longer test that requires solving a small problem and comparing floating point results. diff --git a/setup.py b/setup.py new file mode 100644 index 00000000..697eadac --- /dev/null +++ b/setup.py @@ -0,0 +1,8 @@ +from setuptools import setup + +setup( + name="PV Aerodynamic Design Engineering, PVade", + version="0.0", + packages=["pvade"], + zip_safe=False, +)