diff --git a/.eggs/README.txt b/.eggs/README.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5d01668824f45c3a6683e12d1b9048bb1d273041
--- /dev/null
+++ b/.eggs/README.txt
@@ -0,0 +1,6 @@
+This directory contains eggs that were downloaded by setuptools to build, test, and run plug-ins.
+
+This directory caches those eggs to prevent repeated downloads.
+
+However, it is safe to delete this directory.
+
diff --git a/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/LICENSE b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/LICENSE
new file mode 100644
index 0000000000000000000000000000000000000000..353924be0e59b9ad7e6c22848c2189398481821d
--- /dev/null
+++ b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/LICENSE
@@ -0,0 +1,19 @@
+Copyright Jason R. Coombs
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to
+deal in the Software without restriction, including without limitation the
+rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+sell copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
+IN THE SOFTWARE.
diff --git a/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/PKG-INFO b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/PKG-INFO
new file mode 100644
index 0000000000000000000000000000000000000000..75c36d65e0bc03a2c5e8d7995f786956ebc4823e
--- /dev/null
+++ b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/PKG-INFO
@@ -0,0 +1,185 @@
+Metadata-Version: 2.1
+Name: pytest-runner
+Version: 6.0.1
+Summary: Invoke py.test as distutils command with dependency resolution
+Home-page: https://github.com/pytest-dev/pytest-runner/
+Author: Jason R. Coombs
+Author-email: jaraco@jaraco.com
+Classifier: Development Status :: 7 - Inactive
+Classifier: Intended Audience :: Developers
+Classifier: License :: OSI Approved :: MIT License
+Classifier: Programming Language :: Python :: 3
+Classifier: Programming Language :: Python :: 3 :: Only
+Classifier: Framework :: Pytest
+Requires-Python: >=3.7
+License-File: LICENSE
+Provides-Extra: docs
+Requires-Dist: sphinx ; extra == 'docs'
+Requires-Dist: jaraco.packaging >=9 ; extra == 'docs'
+Requires-Dist: rst.linker >=1.9 ; extra == 'docs'
+Requires-Dist: jaraco.tidelift >=1.4 ; extra == 'docs'
+Provides-Extra: testing
+Requires-Dist: pytest >=6 ; extra == 'testing'
+Requires-Dist: pytest-checkdocs >=2.4 ; extra == 'testing'
+Requires-Dist: pytest-flake8 ; extra == 'testing'
+Requires-Dist: pytest-cov ; extra == 'testing'
+Requires-Dist: pytest-enabler >=1.0.1 ; extra == 'testing'
+Requires-Dist: pytest-virtualenv ; extra == 'testing'
+Requires-Dist: types-setuptools ; extra == 'testing'
+Requires-Dist: pytest-black >=0.3.7 ; (platform_python_implementation != "PyPy") and extra == 'testing'
+Requires-Dist: pytest-mypy >=0.9.1 ; (platform_python_implementation != "PyPy") and extra == 'testing'
+
+.. image:: https://img.shields.io/pypi/v/pytest-runner.svg
+ :target: `PyPI link`_
+
+.. image:: https://img.shields.io/pypi/pyversions/pytest-runner.svg
+ :target: `PyPI link`_
+
+.. _PyPI link: https://pypi.org/project/pytest-runner
+
+.. image:: https://github.com/pytest-dev/pytest-runner/workflows/tests/badge.svg
+ :target: https://github.com/pytest-dev/pytest-runner/actions?query=workflow%3A%22tests%22
+ :alt: tests
+
+.. image:: https://img.shields.io/badge/code%20style-black-000000.svg
+ :target: https://github.com/psf/black
+ :alt: Code style: Black
+
+.. .. image:: https://readthedocs.org/projects/skeleton/badge/?version=latest
+.. :target: https://skeleton.readthedocs.io/en/latest/?badge=latest
+
+.. image:: https://img.shields.io/badge/skeleton-2022-informational
+ :target: https://blog.jaraco.com/skeleton
+
+.. image:: https://tidelift.com/badges/package/pypi/pytest-runner
+ :target: https://tidelift.com/subscription/pkg/pypi-pytest-runner?utm_source=pypi-pytest-runner&utm_medium=readme
+
+Setup scripts can use pytest-runner to add setup.py test support for pytest
+runner.
+
+Deprecation Notice
+==================
+
+pytest-runner depends on deprecated features of setuptools and relies on features that break security
+mechanisms in pip. For example 'setup_requires' and 'tests_require' bypass ``pip --require-hashes``.
+See also `pypa/setuptools#1684 <https://github.com/pypa/setuptools/issues/1684>`_.
+
+It is recommended that you:
+
+- Remove ``'pytest-runner'`` from your ``setup_requires``, preferably removing the ``setup_requires`` option.
+- Remove ``'pytest'`` and any other testing requirements from ``tests_require``, preferably removing the ``tests_requires`` option.
+- Select a tool to bootstrap and then run tests such as tox.
+
+Usage
+=====
+
+- Add 'pytest-runner' to your 'setup_requires'. Pin to '>=2.0,<3dev' (or
+ similar) to avoid pulling in incompatible versions.
+- Include 'pytest' and any other testing requirements to 'tests_require'.
+- Invoke tests with ``setup.py pytest``.
+- Pass ``--index-url`` to have test requirements downloaded from an alternate
+ index URL (unnecessary if specified for easy_install in setup.cfg).
+- Pass additional py.test command-line options using ``--addopts``.
+- Set permanent options for the ``python setup.py pytest`` command (like ``index-url``)
+ in the ``[pytest]`` section of ``setup.cfg``.
+- Set permanent options for the ``py.test`` run (like ``addopts`` or ``pep8ignore``) in the ``[pytest]``
+ section of ``pytest.ini`` or ``tox.ini`` or put them in the ``[tool:pytest]``
+ section of ``setup.cfg``. See `pytest issue 567
+ <https://github.com/pytest-dev/pytest/issues/567>`_.
+- Optionally, set ``test=pytest`` in the ``[aliases]`` section of ``setup.cfg``
+ to cause ``python setup.py test`` to invoke pytest.
+
+Example
+=======
+
+The most simple usage looks like this in setup.py::
+
+ setup(
+ setup_requires=[
+ 'pytest-runner',
+ ],
+ tests_require=[
+ 'pytest',
+ ],
+ )
+
+Additional dependencies require to run the tests (e.g. mock or pytest
+plugins) may be added to tests_require and will be downloaded and
+required by the session before invoking pytest.
+
+Follow `this search on github
+<https://github.com/search?utf8=%E2%9C%93&q=filename%3Asetup.py+pytest-runner&type=Code&ref=searchresults>`_
+for examples of real-world usage.
+
+Standalone Example
+==================
+
+This technique is deprecated - if you have standalone scripts
+you wish to invoke with dependencies, `use pip-run
+<https://pypi.org/project/pip-run>`_.
+
+Although ``pytest-runner`` is typically used to add pytest test
+runner support to maintained packages, ``pytest-runner`` may
+also be used to create standalone tests. Consider `this example
+failure <https://gist.github.com/jaraco/d979a558bc0bf2194c23>`_,
+reported in `jsonpickle #117
+<https://github.com/jsonpickle/jsonpickle/issues/117>`_
+or `this MongoDB test
+<https://gist.github.com/jaraco/0b9e482f5c0a1300dc9a>`_
+demonstrating a technique that works even when dependencies
+are required in the test.
+
+Either example file may be cloned or downloaded and simply run on
+any system with Python and Setuptools. It will download the
+specified dependencies and run the tests. Afterward, the the
+cloned directory can be removed and with it all trace of
+invoking the test. No other dependencies are needed and no
+system configuration is altered.
+
+Then, anyone trying to replicate the failure can do so easily
+and with all the power of pytest (rewritten assertions,
+rich comparisons, interactive debugging, extensibility through
+plugins, etc).
+
+As a result, the communication barrier for describing and
+replicating failures is made almost trivially low.
+
+Considerations
+==============
+
+Conditional Requirement
+-----------------------
+
+Because it uses Setuptools setup_requires, pytest-runner will install itself
+on every invocation of setup.py. In some cases, this causes delays for
+invocations of setup.py that will never invoke pytest-runner. To help avoid
+this contingency, consider requiring pytest-runner only when pytest
+is invoked::
+
+ needs_pytest = {'pytest', 'test', 'ptr'}.intersection(sys.argv)
+ pytest_runner = ['pytest-runner'] if needs_pytest else []
+
+ # ...
+
+ setup(
+ #...
+ setup_requires=[
+ #... (other setup requirements)
+ ] + pytest_runner,
+ )
+
+For Enterprise
+==============
+
+Available as part of the Tidelift Subscription.
+
+This project and the maintainers of thousands of other packages are working with Tidelift to deliver one enterprise subscription that covers all of the open source you use.
+
+`Learn more <https://tidelift.com/subscription/pkg/pypi-PROJECT?utm_source=pypi-PROJECT&utm_medium=referral&utm_campaign=github>`_.
+
+Security Contact
+================
+
+To report a security vulnerability, please use the
+`Tidelift security contact <https://tidelift.com/security>`_.
+Tidelift will coordinate the fix and disclosure.
diff --git a/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/RECORD b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/RECORD
new file mode 100644
index 0000000000000000000000000000000000000000..d7ff24dbeef0d3d0fe2e158334c12c4a101718ee
--- /dev/null
+++ b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/RECORD
@@ -0,0 +1,7 @@
+ptr/__init__.py,sha256=0UfzhCooVgCNTBwVEOPOVGEPck4pnl_6PTfsC-QzNGM,6730
+pytest_runner-6.0.1.dist-info/LICENSE,sha256=2z8CRrH5J48VhFuZ_sR4uLUG63ZIeZNyL4xuJUKF-vg,1050
+pytest_runner-6.0.1.dist-info/METADATA,sha256=Ho3FvAFjFHeY5OQ64WFzkLigFaIpuNr4G3uSmOk3nho,7319
+pytest_runner-6.0.1.dist-info/WHEEL,sha256=oiQVh_5PnQM0E3gPdiz09WCNmwiHDMaGer_elqB3coM,92
+pytest_runner-6.0.1.dist-info/entry_points.txt,sha256=BqezBqeO63XyzSYmHYE58gKEFIjJUd-XdsRQkXHy2ig,58
+pytest_runner-6.0.1.dist-info/top_level.txt,sha256=DPzHbWlKG8yq8EOD5UgEvVNDWeJRPyimrwfShwV6Iuw,4
+pytest_runner-6.0.1.dist-info/RECORD,,
diff --git a/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/WHEEL b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/WHEEL
new file mode 100644
index 0000000000000000000000000000000000000000..98c0d20b7a64f4f998d7913e1d38a05dba20916c
--- /dev/null
+++ b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/WHEEL
@@ -0,0 +1,5 @@
+Wheel-Version: 1.0
+Generator: bdist_wheel (0.42.0)
+Root-Is-Purelib: true
+Tag: py3-none-any
+
diff --git a/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/entry_points.txt b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/entry_points.txt
new file mode 100644
index 0000000000000000000000000000000000000000..08606705d825d4ad98f0a7acfdc654d1e6c1b8fa
--- /dev/null
+++ b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/entry_points.txt
@@ -0,0 +1,3 @@
+[distutils.commands]
+ptr = ptr:PyTest
+pytest = ptr:PyTest
diff --git a/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/requires.txt b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/requires.txt
new file mode 100644
index 0000000000000000000000000000000000000000..153518886b4e238ea0fecf91de9d43e7f4ff31f9
--- /dev/null
+++ b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/requires.txt
@@ -0,0 +1,17 @@
+
+[docs]
+sphinx
+jaraco.packaging>=9
+rst.linker>=1.9
+jaraco.tidelift>=1.4
+
+[testing]
+pytest>=6
+pytest-checkdocs>=2.4
+pytest-flake8
+pytest-cov
+pytest-enabler>=1.0.1
+pytest-virtualenv
+types-setuptools
+pytest-black>=0.3.7
+pytest-mypy>=0.9.1
diff --git a/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/top_level.txt b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/top_level.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e9148ae387cc078aec4d97be87acde54320781de
--- /dev/null
+++ b/.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/top_level.txt
@@ -0,0 +1 @@
+ptr
diff --git a/.eggs/pytest_runner-6.0.1-py3.12.egg/ptr/__init__.py b/.eggs/pytest_runner-6.0.1-py3.12.egg/ptr/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..41192fa83ea9f3194d5a31f429e7fd5515c7529d
--- /dev/null
+++ b/.eggs/pytest_runner-6.0.1-py3.12.egg/ptr/__init__.py
@@ -0,0 +1,216 @@
+"""
+Implementation
+"""
+
+import os as _os
+import shlex as _shlex
+import contextlib as _contextlib
+import sys as _sys
+import operator as _operator
+import itertools as _itertools
+import warnings as _warnings
+
+import pkg_resources
+import setuptools.command.test as orig
+from setuptools import Distribution
+
+
+@_contextlib.contextmanager
+def _save_argv(repl=None):
+ saved = _sys.argv[:]
+ if repl is not None:
+ _sys.argv[:] = repl
+ try:
+ yield saved
+ finally:
+ _sys.argv[:] = saved
+
+
+class CustomizedDist(Distribution):
+
+ allow_hosts = None
+ index_url = None
+
+ def fetch_build_egg(self, req):
+ """Specialized version of Distribution.fetch_build_egg
+ that respects respects allow_hosts and index_url."""
+ from setuptools.command.easy_install import easy_install
+
+ dist = Distribution({'script_args': ['easy_install']})
+ dist.parse_config_files()
+ opts = dist.get_option_dict('easy_install')
+ keep = (
+ 'find_links',
+ 'site_dirs',
+ 'index_url',
+ 'optimize',
+ 'site_dirs',
+ 'allow_hosts',
+ )
+ for key in list(opts):
+ if key not in keep:
+ del opts[key] # don't use any other settings
+ if self.dependency_links:
+ links = self.dependency_links[:]
+ if 'find_links' in opts:
+ links = opts['find_links'][1].split() + links
+ opts['find_links'] = ('setup', links)
+ if self.allow_hosts:
+ opts['allow_hosts'] = ('test', self.allow_hosts)
+ if self.index_url:
+ opts['index_url'] = ('test', self.index_url)
+ install_dir_func = getattr(self, 'get_egg_cache_dir', _os.getcwd)
+ install_dir = install_dir_func()
+ cmd = easy_install(
+ dist,
+ args=["x"],
+ install_dir=install_dir,
+ exclude_scripts=True,
+ always_copy=False,
+ build_directory=None,
+ editable=False,
+ upgrade=False,
+ multi_version=True,
+ no_report=True,
+ user=False,
+ )
+ cmd.ensure_finalized()
+ return cmd.easy_install(req)
+
+
+class PyTest(orig.test):
+ """
+ >>> import setuptools
+ >>> dist = setuptools.Distribution()
+ >>> cmd = PyTest(dist)
+ """
+
+ user_options = [
+ ('extras', None, "Install (all) setuptools extras when running tests"),
+ (
+ 'index-url=',
+ None,
+ "Specify an index url from which to retrieve dependencies",
+ ),
+ (
+ 'allow-hosts=',
+ None,
+ "Whitelist of comma-separated hosts to allow "
+ "when retrieving dependencies",
+ ),
+ (
+ 'addopts=',
+ None,
+ "Additional options to be passed verbatim to the pytest runner",
+ ),
+ ]
+
+ def initialize_options(self):
+ self.extras = False
+ self.index_url = None
+ self.allow_hosts = None
+ self.addopts = []
+ self.ensure_setuptools_version()
+
+ @staticmethod
+ def ensure_setuptools_version():
+ """
+ Due to the fact that pytest-runner is often required (via
+ setup-requires directive) by toolchains that never invoke
+ it (i.e. they're only installing the package, not testing it),
+ instead of declaring the dependency in the package
+ metadata, assert the requirement at run time.
+ """
+ pkg_resources.require('setuptools>=27.3')
+
+ def finalize_options(self):
+ if self.addopts:
+ self.addopts = _shlex.split(self.addopts)
+
+ @staticmethod
+ def marker_passes(marker):
+ """
+ Given an environment marker, return True if the marker is valid
+ and matches this environment.
+ """
+ return (
+ not marker
+ or not pkg_resources.invalid_marker(marker)
+ and pkg_resources.evaluate_marker(marker)
+ )
+
+ def install_dists(self, dist):
+ """
+ Extend install_dists to include extras support
+ """
+ return _itertools.chain(
+ orig.test.install_dists(dist), self.install_extra_dists(dist)
+ )
+
+ def install_extra_dists(self, dist):
+ """
+ Install extras that are indicated by markers or
+ install all extras if '--extras' is indicated.
+ """
+ extras_require = dist.extras_require or {}
+
+ spec_extras = (
+ (spec.partition(':'), reqs) for spec, reqs in extras_require.items()
+ )
+ matching_extras = (
+ reqs
+ for (name, sep, marker), reqs in spec_extras
+ # include unnamed extras or all if self.extras indicated
+ if (not name or self.extras)
+ # never include extras that fail to pass marker eval
+ and self.marker_passes(marker)
+ )
+ results = list(map(dist.fetch_build_eggs, matching_extras))
+ return _itertools.chain.from_iterable(results)
+
+ @staticmethod
+ def _warn_old_setuptools():
+ msg = (
+ "pytest-runner will stop working on this version of setuptools; "
+ "please upgrade to setuptools 30.4 or later or pin to "
+ "pytest-runner < 5."
+ )
+ ver_str = pkg_resources.get_distribution('setuptools').version
+ ver = pkg_resources.parse_version(ver_str)
+ if ver < pkg_resources.parse_version('30.4'):
+ _warnings.warn(msg)
+
+ def run(self):
+ """
+ Override run to ensure requirements are available in this session (but
+ don't install them anywhere).
+ """
+ self._warn_old_setuptools()
+ dist = CustomizedDist()
+ for attr in 'allow_hosts index_url'.split():
+ setattr(dist, attr, getattr(self, attr))
+ for attr in (
+ 'dependency_links install_requires tests_require extras_require '
+ ).split():
+ setattr(dist, attr, getattr(self.distribution, attr))
+ installed_dists = self.install_dists(dist)
+ if self.dry_run:
+ self.announce('skipping tests (dry run)')
+ return
+ paths = map(_operator.attrgetter('location'), installed_dists)
+ with self.paths_on_pythonpath(paths):
+ with self.project_on_sys_path():
+ return self.run_tests()
+
+ @property
+ def _argv(self):
+ return ['pytest'] + self.addopts
+
+ def run_tests(self):
+ """
+ Invoke pytest, replacing argv. Return result code.
+ """
+ with _save_argv(_sys.argv[:1] + self.addopts):
+ result_code = __import__('pytest').main()
+ if result_code:
+ raise SystemExit(result_code)
diff --git a/BscElectrolyteModels.egg-info/PKG-INFO b/BscElectrolyteModels.egg-info/PKG-INFO
new file mode 100644
index 0000000000000000000000000000000000000000..4600c3ca15b76a50af6f556a017f94b85d207b6d
--- /dev/null
+++ b/BscElectrolyteModels.egg-info/PKG-INFO
@@ -0,0 +1,117 @@
+Metadata-Version: 2.1
+Name: BscElectrolyteModels
+Version: 1.0
+Summary: A description is yet to follow
+Home-page: https://git.rwth-aachen.de/JanHab/Bsc-ElectrolyteModels
+Author: Jan Habscheid, Lambert Theisen, Manuel Torrilhon
+Author-email: Jan.Habscheid@rwth-aachen.de, lambert.theisen@rwth-aachen.de, mt@mathcces.rwth-aachen.de
+License: GNU General Public License v3.0 or later
+Keywords: fenics-project,Battery Simulation,finite element method
+Classifier: License :: OSI Approved :: GNU General Public License v3 (GPLv3)
+Classifier: Programming Language :: Python :: 3
+Requires-Python: ==3.12.3
+Description-Content-Type: text/markdown
+License-File: LICENSE
+Requires-Dist: pyvista==0.43.10
+Requires-Dist: numpy==1.26.4
+Requires-Dist: scipy==1.14.0
+Provides-Extra: notebook
+Requires-Dist: jupyter>=1.0.0; extra == "notebook"
+Provides-Extra: tests
+Requires-Dist: pytest==8.3.3; extra == "tests"
+Provides-Extra: docs
+Requires-Dist: sphinx==7.3.7; extra == "docs"
+Requires-Dist: myst-parser==4.0.0; extra == "docs"
+Requires-Dist: sphinx-copybotton==0.5.2; extra == "docs"
+Requires-Dist: sphinx_rtd_theme==3.0.1; extra == "docs"
+Provides-Extra: build
+Requires-Dist: build>=0.7.0; extra == "build"
+Provides-Extra: dev
+Requires-Dist: BScElectrolytemodels[tests]; extra == "dev"
+Requires-Dist: BScElectrolytemodels[docs]; extra == "dev"
+Requires-Dist: BScElectrolytemodels[build]; extra == "dev"
+
+# Reproducibility Repository for Numerical Treatment of a Thermodynamically Consistent Electrolyte Model (B.Sc. Thesis - Jan Habscheid)
+
+[](https://git.rwth-aachen.de/Jan.Habscheid/bsc-electrolytemodels/pipelines)
+[](https://janhab.pages.rwth-aachen.de/bsc-electrolytemodels/)
+[](https://doi.org/10.5281/zenodo.13645296)
+[](https://git.rwth-aachen.de/jan.habscheid/bsc-electrolytemodels/-/tags)
+[](https://git.rwth-aachen.de/Jan.Habscheid/bsc-electrolytemodels/-/blob/main/LICENSE?ref_type=heads)
+
+## Thesis
+
+This repository contains the code to reproduce the results presented in the bachelor thesis: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model
+Find the thesis at [https://doi.org/10.18154/RWTH-2024-09837](https://doi.org/10.18154/RWTH-2024-09837)
+
+### Abstract
+
+Batteries play a crucial role in the energy transition. The production of green energy depends on external factors. Storing energy in batteries is necessary to access green energy at any time.
+
+Better optimized batteries are essential for the future. Lifetime, loading time, and energy loss are just some aspects that must be improved to prepare for a greener future. Numerical simulations are crucial to understanding and optimizing batteries' behavior. Those simulations enable researchers to test many different materials without considerable additional expenses to, for example, find the best combination of anions and cations.
+
+The classical Nernst-Planck model for the ion transport in an electrolyte fails to predict the correct concentration in the boundaries of the electrolyte. This work will present and analyze a thermodynamically consistent electrolyte model with dimensionless units under isothermal conditions. A simplified version of the system for the one-dimensional equilibrium of an ideal mixture and the incompressible limit will be considered. The numerical implementation of the model with the open-source software FEniCSx will be discussed. Furthermore, the influence of different boundary conditions, material parameters, solvation, and compressibility on the electric potential, pressure, and ion concentration will be investigated, and the model will be compared with the classical Nernst-Planck model. Examples of the double layer capacity and electrolytic diode will be considered.
+
+## Installation
+
+As a numerical solver, mainly FEniCSx was used and installed via conda.
+All the calculations were performed on a Linux machine. According to the documentation, everything should work well on macOS, but this was not tested. FEniCSx offers some beta versions for Windows support, but it is recommended to use WSL2 instead.
+
+```
+conda create --name fenicsx-env python=3.12.3 -y
+conda activate fenicsx-env
+conda install -c conda-forge fenics-dolfinx=0.8.0 mpich=4.2.1 pyvista=0.43.10 matplotlib=3.8.4 numpy=1.26.4 scipy=1.14.0 pytest==8.3.3 -y
+```
+
+### Alternative installation
+
+Use the "environment.yml" file to install all necessary environments
+
+```
+conda env create -f environment.yml
+```
+
+### macOS installation using Docker
+
+```
+docker compose build
+docker compose run solver
+```
+
+### Testing
+
+Use pytest with
+```
+python -m pytest
+```
+
+to verify that everything was installed correctly.
+
+
+## Usage
+
+Find the visualizations from the thesis and some extra calculations in the "examples" folder.
+In the subfolder "ReproducableCode" is the code, to execute the calculations with some first visualizations.
+The subfolder "Data" stores the data for all the simulations in a *.npz file, which can be read with numpy `np.load(file.npz)`.
+"Visualizations" creates the necessary figures from the thesis and stores them in *.svg format in "Figures".
+
+In "src" there are the generic FEniCSx implementations, that were used to calculate the examples.
+
+
+## Contact
+
+**Author**
+- Jan Habscheid
+- Jan.Habscheid@rwth-aachen.de
+
+**Supervisor**
+- Dr. Lambert Theissen
+- ACoM - Applied and Computational Mathematics
+- RWTH Aachen University
+- theisen@acom.rwth-aachen.de
+
+**Supervisor**
+- Prof. Dr. Manuel Torrilhon
+- ACoM - Applied and Computational Mathematics
+- RWTH Aachen University
+- mt@acom.rwth-aachen.de
diff --git a/BscElectrolyteModels.egg-info/SOURCES.txt b/BscElectrolyteModels.egg-info/SOURCES.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ec989abfcf3359a0372724bcc4fc2f9cd9c274a0
--- /dev/null
+++ b/BscElectrolyteModels.egg-info/SOURCES.txt
@@ -0,0 +1,23 @@
+LICENSE
+README.md
+pyproject.toml
+setup.cfg
+setup.py
+BscElectrolyteModels.egg-info/PKG-INFO
+BscElectrolyteModels.egg-info/SOURCES.txt
+BscElectrolyteModels.egg-info/dependency_links.txt
+BscElectrolyteModels.egg-info/entry_points.txt
+BscElectrolyteModels.egg-info/requires.txt
+BscElectrolyteModels.egg-info/top_level.txt
+BscElectrolyteModels.egg-info/zip-safe
+src/ElectrolyticDiode.py
+src/Eq02.py
+src/Eq04.py
+src/EqN.py
+src/Helper_DoubleLayerCapacity.py
+src/__init__.py
+tests/test_ElectrolyticDiode.py
+tests/test_Eq02.py
+tests/test_Eq04.py
+tests/test_EqN.py
+tests/test_Helper_DoubleLayerCapacity.py
\ No newline at end of file
diff --git a/BscElectrolyteModels.egg-info/dependency_links.txt b/BscElectrolyteModels.egg-info/dependency_links.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8b137891791fe96927ad78e64b0aad7bded08bdc
--- /dev/null
+++ b/BscElectrolyteModels.egg-info/dependency_links.txt
@@ -0,0 +1 @@
+
diff --git a/BscElectrolyteModels.egg-info/entry_points.txt b/BscElectrolyteModels.egg-info/entry_points.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7bcc5d3ffd0275679d5e3d0d72b5a5c1dabcd4d6
--- /dev/null
+++ b/BscElectrolyteModels.egg-info/entry_points.txt
@@ -0,0 +1,2 @@
+[console_scripts]
+my-example-utility = example.example_module:main
diff --git a/BscElectrolyteModels.egg-info/requires.txt b/BscElectrolyteModels.egg-info/requires.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f08ddad2fbba14b117fd049223a2d5abadb233e2
--- /dev/null
+++ b/BscElectrolyteModels.egg-info/requires.txt
@@ -0,0 +1,23 @@
+pyvista==0.43.10
+numpy==1.26.4
+scipy==1.14.0
+
+[build]
+build>=0.7.0
+
+[dev]
+BScElectrolytemodels[tests]
+BScElectrolytemodels[docs]
+BScElectrolytemodels[build]
+
+[docs]
+sphinx==7.3.7
+myst-parser==4.0.0
+sphinx-copybotton==0.5.2
+sphinx_rtd_theme==3.0.1
+
+[notebook]
+jupyter>=1.0.0
+
+[tests]
+pytest==8.3.3
diff --git a/BscElectrolyteModels.egg-info/top_level.txt b/BscElectrolyteModels.egg-info/top_level.txt
new file mode 100644
index 0000000000000000000000000000000000000000..eac2d36c1ef1fda9a6f9c8786db84293051c1009
--- /dev/null
+++ b/BscElectrolyteModels.egg-info/top_level.txt
@@ -0,0 +1 @@
+Bsc-ElectrolyteModels
diff --git a/BscElectrolyteModels.egg-info/zip-safe b/BscElectrolyteModels.egg-info/zip-safe
new file mode 100644
index 0000000000000000000000000000000000000000..8b137891791fe96927ad78e64b0aad7bded08bdc
--- /dev/null
+++ b/BscElectrolyteModels.egg-info/zip-safe
@@ -0,0 +1 @@
+
diff --git a/README.md b/README.md
index 0067dace29cb3ac2028e1c91a5ced5e4d531a620..5c040e987b20b5020b42615183314375ae8a89f5 100644
--- a/README.md
+++ b/README.md
@@ -27,7 +27,7 @@ All the calculations were performed on a Linux machine. According to the documen
```
conda create --name fenicsx-env python=3.12.3 -y
conda activate fenicsx-env
-conda install -c conda-forge fenics-dolfinx=0.8.0 mpich=4.2.1 pyvista=0.43.10 matplotlib=3.8.4 numpy=1.26.4 scipy=1.14.0 pytest==8.3.3 -y
+conda install -c conda-forge fenics-dolfinx=0.8.0 mpich=4.2.1 pyvista=0.43.10 gcc=12.4.0 matplotlib=3.8.4 numpy=1.26.4 scipy=1.14.0 pytest==8.3.3 -y
```
### Alternative installation
diff --git a/build/lib/Bsc-ElectrolyteModels/ElectrolyticDiode.py b/build/lib/Bsc-ElectrolyteModels/ElectrolyticDiode.py
new file mode 100644
index 0000000000000000000000000000000000000000..836b48c19f548142f22d74f01577dde207448d40
--- /dev/null
+++ b/build/lib/Bsc-ElectrolyteModels/ElectrolyticDiode.py
@@ -0,0 +1,401 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+
+This module solves the dimensionless system of equations for the example of an electric diode
+
+# ! Disclaimer: The results conform with [A Numerical Strategy for Nernst–Planck Systems with Solvation Effect, J. Fuhrmann, DOI:10.1002/fuce.201500215]. As the size of the domain, the boundary conditions, and the parameters are chosen differently, it can not be expected to match the same results quantitatively but qualitatively. In [Fuhrmann], only the potential and the cations are visualized. The potential behaves the same, and the cations are pushed to the outside of the domain for the forward bias and to the center of the domain for the backward bias. Furthermore, a rectification effect in the concentrations of the ions cannot be seen, as the range of the concentrations is the same along the different biases. On the other hand, this rectification process can be observed in [Entropy and convergence analysis for two finite volume schemes for a Nernst-Planck-Poisson system with ion volume constraints, B. Gaudeaul, J. Fuhrmann, DOI:10.1007/s00211-022-01279-y]. The electric potential can be verified to match the behavior from [Gaudeaul & Fuhrmann]. Furthermore, the ion concentrations can also be validated qualitatively. However, in [Gaudeaul & Fuhrmann], a rectification effect can be observed, reducing the concentrations of anions and cations for the reverse bias and no bias compared to the forward bias.
+'''
+
+import numpy as np
+from mpi4py import MPI
+from dolfinx import mesh, fem, log, io
+from dolfinx.fem.petsc import NonlinearProblem
+from dolfinx.nls.petsc import NewtonSolver
+from ufl import TestFunctions, split, dot, grad, dx, inner, ln
+from basix.ufl import element, mixed_element
+import matplotlib.pyplot as plt
+from ufl import Measure
+from dolfinx.mesh import locate_entities, meshtags
+
+def ElectrolyticDiode(Bias_type:str, phi_bias:float, g_phi:float, z_A:float, z_C:float, y_A_bath:float, y_C_bath:float, K:float|str, Lambda2:float, a2:float, number_cells:list, solvation:float = 0, PoissonBoltzmann:bool=False, relax_param:float=None, Lx:float=2, Ly:float=10, rtol:float=1e-8, max_iter:float=500, return_type:str='Vector'):
+ '''
+ Solves the system of equations for the example of an electric diode
+
+ Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
+
+ ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
+
+ Parameters
+ ----------
+ Bias_type : str
+ ForwardBias, NoBias, BackwardBias
+ phi_bias : float
+ Bias in φ
+ g_phi : float
+ Neumann boundary condition for φ
+ z_A : float
+ Charge number of species A
+ z_C : float
+ Charge number of species C
+ y_A_bath : float
+ Bath concentration of anions
+ y_C_bath : float
+ Bath concentration of cations
+ K : float | str
+ Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
+ ! Only implemented for incompressible mixtures, yet
+ Lambda2 : float
+ Dimensionless parameter
+ a2 : float
+ Dimensionless parameter
+ number_cells : int
+ Number of cells in the mesh
+ solvation : float, optional
+ solvation number, by default 0
+ PoissonBoltzmann : bool, optional
+ Solve classical Nernst-Planck model with the use of the Poisson-Boltzmann formulation if True, else solve the presented model by Dreyer, Guhlke, Müller, by default False
+ relax_param : float, optional
+ Relaxation parameter for the Newton solver
+ xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
+ , by default None -> Determined automatically
+ Lx : float, optional
+ Length of domain in x-direction, by default 2
+ Ly : float, optional
+ Length of domain in y-direction, by default 10
+ rtol : float, optional
+ Relative tolerance for Newton solver, by default 1e-8
+ max_iter : float, optional
+ Maximum number of Newton iterations, by default 500
+ return_type : str, optional
+ Type of return value, by default 'Vector'
+ 'Vector' -> Returns the solution as a numpy array for triangle elements
+ 'Extended' -> Returns the solution as both, a numpy array for triangle elements and the fenicsx function u
+
+
+ Returns
+ -------
+ y_A_vector, y_C_vector, phi_vector, p_vector, x_vector
+ Returns atomic fractions for species A and C, electric potential, pressure, and the mesh as numpy arrays for triangle elements. x_vector is a list of both dimensions [x, y]
+ If return_type is 'Extended', also returns the fenicsx function u
+ '''
+ if K != 'incompressible':
+ raise NotImplementedError('Only incompressible electrolytes are implemented yet')
+ x0 = np.array([0, 0])
+ x1 = np.array([Lx, Ly])
+ match Bias_type:
+ case 'ForwardBias': phi_bias = phi_bias
+ case 'NoBias': phi_bias = 0
+ case 'BackwardBias': phi_bias = -phi_bias
+ case _: raise ValueError('Invalid Bias_type')
+
+ # Define boundaries
+ geom_tol = 1E-4 # ! Geometric tolerance, may need to be adjusted
+ def Left(x):
+ return np.isclose(x[0], x0[0], geom_tol, geom_tol)
+
+ def Right(x):
+ return np.isclose(x[0], x1[0], geom_tol, geom_tol)
+
+ def Top(x):
+ return np.isclose(x[1], x1[1], geom_tol, geom_tol)
+
+ def Bottom(x):
+ return np.isclose(x[1], x0[1], geom_tol, geom_tol)
+
+ def Right_Top(x):
+ NeumannPart = np.logical_and(1/2*Ly < x[1], x[1] < 3/4*Ly)
+ RightPart = np.isclose(x[0], x1[0], geom_tol, geom_tol)
+ return np.logical_and(RightPart, NeumannPart)
+
+ def Right_Bottom(x):
+ NeumannPart = np.logical_and(1/4*Ly < x[1], x[1] < 1/2*Ly)
+ RightPart = np.isclose(x[0], x1[0], geom_tol, geom_tol)
+ return np.logical_and(RightPart, NeumannPart)
+
+ def Left_Top(x):
+ NeumannPart = np.logical_and(1/2*Ly < x[1], x[1] < 3/4*Ly)
+ LeftPart = np.isclose(x[0], x0[0], geom_tol, geom_tol)
+ return np.logical_and(LeftPart, NeumannPart)
+
+ def Left_Bottom(x):
+ NeumannPart = np.logical_and(1/4*Ly < x[1], x[1] < 1/2*Ly)
+ LeftPart = np.isclose(x[0], x0[0], geom_tol, geom_tol)
+ return np.logical_and(LeftPart, NeumannPart)
+
+ # Create the mesh
+ msh = mesh.create_rectangle(MPI.COMM_WORLD, [x0, x1], number_cells)
+
+ # Define Finite Elements
+ CG1_elem = element('Lagrange', msh.basix_cell(), 1)
+ # from ufl import FiniteElement
+ # CG1_elem = FiniteElement("Lagrange", msh.ufl_cell(), 1)
+
+ # Define Mixed Function Space
+ W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
+ W = fem.functionspace(msh, W_elem)
+
+ # Define Trial- and Testfunctions
+ u = fem.Function(W)
+ y_A, y_C, phi, p = split(u)
+ (v_A, v_C, v_1, v_2) = TestFunctions(W)
+
+ # Define boundary conditions
+ W0, _ = W.sub(0).collapse()
+ W1, _ = W.sub(1).collapse()
+ W2, _ = W.sub(2).collapse()
+ W3, _ = W.sub(3).collapse()
+
+ def phi_bottom_(x):
+ return np.full_like(x[1], 0)
+ def phi_top_(x):
+ return np.full_like(x[1], phi_bias)
+ def y_A_bottom_(x):
+ return np.full_like(x[1], y_A_bath)
+ def y_A_top_(x):
+ return np.full_like(x[1], y_A_bath)
+ def y_C_bottom_(x):
+ return np.full_like(x[1], y_C_bath)
+ def y_C_top_(x):
+ return np.full_like(x[1], y_C_bath)
+
+ phi_bottom_bcs = fem.Function(W2)
+ phi_bottom_bcs.interpolate(phi_bottom_)
+ phi_top_bcs = fem.Function(W2)
+ phi_top_bcs.interpolate(phi_top_)
+ y_A_bottom_bcs = fem.Function(W0)
+ y_A_bottom_bcs.interpolate(y_A_bottom_)
+ y_A_top_bcs = fem.Function(W0)
+ y_A_top_bcs.interpolate(y_A_top_)
+ y_C_bottom_bcs = fem.Function(W1)
+ y_C_bottom_bcs.interpolate(y_C_bottom_)
+ y_C_top_bcs = fem.Function(W1)
+ y_C_top_bcs.interpolate(y_C_top_)
+
+ # Identify face tags and create dirichlet conditions
+ facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Bottom)
+ facet_top_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Top)
+ bc_bottom_phi = fem.dirichletbc(phi_bottom_bcs, facet_bottom_dofs, W.sub(2))
+ bc_top_phi = fem.dirichletbc(phi_top_bcs, facet_top_dofs, W.sub(2))
+
+ facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Bottom)
+ bc_bottom_y_A = fem.dirichletbc(y_A_bottom_bcs, facet_bottom_dofs, W.sub(0))
+ facet_top_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Top)
+ bc_top_y_A = fem.dirichletbc(y_A_top_bcs, facet_top_dofs, W.sub(0))
+
+ facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Bottom)
+ bc_bottom_y_C = fem.dirichletbc(y_C_bottom_bcs, facet_bottom_dofs, W.sub(1))
+ facet_top_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Top)
+ bc_top_y_C = fem.dirichletbc(y_C_top_bcs, facet_top_dofs, W.sub(1))
+
+ # Collect boundary conditions
+ bcs = [bc_bottom_phi, bc_top_phi, bc_bottom_y_A, bc_top_y_A, bc_bottom_y_C, bc_top_y_C]
+ # bcs = [bc_bottom_phi, bc_top_phi, bc_bottom_y_A, bc_bottom_y_C]
+
+
+ # def p_bottom_(x):
+ # return np.full_like(x[1], 0)
+ # def p_top_(x):
+ # return np.full_like(x[1], 0)
+ # p_bottom_bcs = fem.Function(W3)
+ # p_bottom_bcs.interpolate(p_bottom_)
+ # p_top_bcs = fem.Function(W3)
+ # p_top_bcs.interpolate(p_top_)
+ # facet_bottom_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Bottom)
+ # facet_top_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Top)
+ # bc_bottom_p = fem.dirichletbc(p_bottom_bcs, facet_bottom_dofs, W.sub(3))
+ # bc_top_p = fem.dirichletbc(p_top_bcs, facet_top_dofs, W.sub(3))
+ # bcs.append(bc_bottom_p)
+ # bcs.append(bc_top_p)
+
+ # Variational formulation
+ if K == 'incompressible':
+ # def nF(y_A, y_C):
+ # return (z_C * y_C + z_A * y_A) # n = 1
+
+ # A = (
+ # inner(grad(phi), grad(v_1)) * dx
+ # - 1 / Lambda2 * nF(y_A, y_C) * v_1 * dx
+ # ) + (
+ # inner(grad(p), grad(v_2)) * dx
+ # + 1 / a2 * nF(y_A, y_C) * dot(grad(phi), grad(v_2)) * dx
+ # ) + (
+ # inner(grad(ln(y_A) + a2 * solvation * p - ln(1-y_A-y_C) + z_A * phi), grad(v_A)) * dx
+ # + inner(grad(ln(y_C) + a2 * solvation * p- ln(1-y_A-y_C) + z_C * phi), grad(v_C)) * dx
+ # )
+ # def nF(y_A, y_C):
+ # # n = const = 1
+ # return (z_C * y_C + z_A * y_A)
+
+ # def J_A(y_A, y_C, phi, p):
+ # g_A = (solvation + 1) * a2 * (p - 1) # g_Aref, but constant and take gradient
+ # mu_A = g_A + ln(y_A)
+ # g_N = a2 * (p - 1) # solvation_N = 0, g_Sref, but constant and take gradient
+ # mu_N = g_N + ln(1 - y_A - y_C)
+ # return grad(mu_A - mu_N + z_A * phi)
+ # # return grad(ln(y_A) - ln(1 - y_A - y_C) + z_A * phi)
+
+ # def J_C(y_A, y_C, phi, p):
+ # g_C = (solvation + 1) * a2 * (p - 1) # g_Cref, but constant and take gradient
+ # mu_C = g_C + ln(y_C)
+ # g_N = a2 * (p - 1)
+ # mu_N = g_N + ln(1 - y_A - y_C)
+ # return grad(mu_C - mu_N + z_C * phi)
+ # # return grad(ln(y_C) - ln(1 - y_A - y_C) + z_C * phi)
+
+ def nF(y_A, y_C):
+ return (z_C * y_C + z_A * y_A)
+
+ # Diffusion fluxes for species A and C
+ def J_A(y_A, y_C, phi, p):
+ return grad(ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi)
+
+ def J_C(y_A, y_C, phi, p):
+ return grad(ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi)
+
+ # if PoissonBoltzmann:
+ # def J_A(y_A, y_C, phi, p):
+ # return grad(ln(y_A) + z_A * phi)
+ # def J_C(y_A, y_C, phi, p):
+ # return grad(ln(y_C) + z_C * phi)
+
+ A = (
+ inner(grad(phi), grad(v_1)) * dx
+ - 1 / Lambda2 * nF(y_A, y_C) * v_1 * dx
+ ) + (
+ inner(grad(p), grad(v_2)) * dx
+ + 1 / a2 * nF(y_A, y_C) * dot(grad(phi), grad(v_2)) * dx
+ ) + (
+ inner(J_A(y_A, y_C, phi, p), grad(v_A)) * dx
+ + inner(J_C(y_A, y_C, phi, p), grad(v_C)) * dx
+ )
+
+ # Define Neumann boundaries
+ boundaries = [(0, Right_Top), (1, Right_Bottom), (2, Left_Top), (3, Left_Bottom)]
+ facet_indices, facet_markers = [], []
+ fdim = msh.topology.dim - 1
+ for (marker, locator) in boundaries:
+ facets = locate_entities(msh, fdim, locator)
+ facet_indices.append(facets)
+ facet_markers.append(np.full_like(facets, marker))
+ facet_indices = np.hstack(facet_indices).astype(np.int32)
+ facet_markers = np.hstack(facet_markers).astype(np.int32)
+ sorted_facets = np.argsort(facet_indices)
+ facet_tag = meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
+ ds = Measure("ds", domain=msh, subdomain_data=facet_tag)
+ L = 0
+ L += (
+ inner(g_phi, v_1) * ds(0)
+ + inner(-g_phi, v_1) * ds(1)
+ # Uncomment, if you want to apply the Neumann bcs on both sides/ left side
+ # + inner(g_phi, v_1) * ds(2)
+ # + inner(-g_phi, v_1) * ds(3)
+ )
+
+ F = A + L
+
+ # Initialize initial guess for u
+ y_C_init = fem.Function(W1)
+ y_A_init = fem.Function(W0)
+ y_C_init.interpolate(lambda x: np.full_like(x[0], y_C_bath))
+ y_A_init.interpolate(lambda x: np.full_like(x[0], y_A_bath))
+
+ with u.vector.localForm() as u_loc:
+ u_loc.set(0)
+ u.sub(0).interpolate(y_A_init)
+ u.sub(1).interpolate(y_C_init)
+
+ # Define Nonlinear Problem
+ problem = NonlinearProblem(F, u, bcs=bcs)
+
+ solver = NewtonSolver(MPI.COMM_WORLD, problem)
+ solver.convergence_criterion = 'residual' #"incremental"
+ solver.rtol = rtol
+ solver.relaxation_parameter = relax_param
+ solver.max_it = max_iter
+ solver.report = True
+
+ # Solve the problem
+ log.set_log_level(log.LogLevel.INFO)
+ n, converged = solver.solve(u)
+ assert (converged)
+ print(f"Number of interations: {n:d}")
+
+ # Extract solution
+ y_A_vector = u.sub(0).collapse().x.array
+ y_C_vector = u.sub(1).collapse().x.array
+ phi_vector = u.sub(2).collapse().x.array
+ p_vector = u.sub(3).collapse().x.array
+
+ X = W.sub(0).collapse()[0].tabulate_dof_coordinates()
+ x_vector = [X[:, 0], X[:, 1]]
+
+ if return_type == 'Vector':
+ return y_A_vector, y_C_vector, phi_vector, p_vector, x_vector
+ elif return_type == 'Extended':
+ return y_A_vector, y_C_vector, phi_vector, p_vector, x_vector, u
+ else:
+ raise ValueError('Invalid return_type')
+
+if __name__ == '__main__':
+ phi_bias = 10#10
+ Bias_type = 'ForwardBias' # 'ForwardBias', 'NoBias', 'BackwardBias'
+ g_phi = 350#5
+ y_fixed = 0.01#0.01
+ z_A = -1.0
+ z_C = 1.0
+ K = 'incompressible'
+ Lambda2 = 8.553e-6 # ! Change back to 1e-6
+ # g_phi *= np.sqrt(Lambda2) # ! Unsure about this scaling
+ # g_phi *= Lambda2
+ # g_phi = g_phi / np.sqrt(Lambda2)
+ a2 = 7.5412e-4
+ number_cells = [20, 128]#[20, 100]
+ Lx = 0.02
+ Ly = 0.1
+ x0 = np.array([0, 0])
+ x1 = np.array([Lx, Ly])
+ refinement_style = 'uniform'
+ solvation = 5
+ PoissonBoltzmann = False
+ rtol = 1e-3 # ToDo: Change back to 1e-8, currently just for testing
+ relax_param = 0.15 # 0.1
+ max_iter = 15_000
+
+ # Solve the system
+ y_A, y_C, phi, p, X = ElectrolyticDiode(Bias_type, phi_bias, g_phi, z_A, z_C, y_fixed, y_fixed, K, Lambda2, a2, number_cells, solvation, PoissonBoltzmann, relax_param, Lx, Ly, rtol, max_iter, return_type='Vector')
+ x, y = X[0], X[1]
+ y_S = 1 - y_A - y_C
+
+ # Plot results
+ levelsf = 10
+ levels = 10
+ fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(8,10))
+
+ c = axs[0,0].tricontourf(x, y, y_A, levelsf)
+ fig.colorbar(c)
+ axs[0,0].tricontour(x, y, y_A, levels, colors='black')
+ axs[0,0].set_title('$y_A$')
+
+ c = axs[0,1].tricontourf(x, y, y_C, levelsf)
+ fig.colorbar(c)
+ axs[0,1].tricontour(x, y, y_C, levels, colors='black')
+ axs[0,1].set_title('$y_C$')
+
+ c = axs[0,2].tricontourf(x, y, y_S, levelsf)
+ fig.colorbar(c)
+ axs[0,2].tricontour(x, y, y_S, levels, colors='black')
+ axs[0,2].set_title('$y_S$')
+
+ c = axs[1,0].tricontourf(x, y, phi, levelsf)
+ fig.colorbar(c)
+ axs[1,0].tricontour(x, y, phi, levels, colors='black')
+ axs[1,0].set_title('$\\varphi$')
+
+ c = axs[1,1].tricontourf(x, y, p, levelsf)
+ fig.colorbar(c)
+ axs[1,1].tricontour(x, y, p, levels, colors='black')
+ axs[1,1].set_title('$p$')
+
+ fig.tight_layout()
+ fig.show()
\ No newline at end of file
diff --git a/build/lib/Bsc-ElectrolyteModels/Eq02.py b/build/lib/Bsc-ElectrolyteModels/Eq02.py
new file mode 100644
index 0000000000000000000000000000000000000000..0ec56e7491732d7519b83f8a79fdd3ffc8a34bcf
--- /dev/null
+++ b/build/lib/Bsc-ElectrolyteModels/Eq02.py
@@ -0,0 +1,309 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+
+This file implements the fenics solver for the reduced systme of equations to two equations for the electric potential and the pressure.
+'''
+
+import numpy as np
+from mpi4py import MPI
+from dolfinx import mesh, fem, log
+from dolfinx.fem.petsc import NonlinearProblem
+from dolfinx.nls.petsc import NewtonSolver
+from ufl import TestFunctions, split, dot, grad, dx, inner, Mesh, exp
+from basix.ufl import element, mixed_element
+import matplotlib.pyplot as plt
+
+# Define mesh
+def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
+ '''
+ Creates a one-dimensional mesh with a refined region at the left boundary
+
+ Parameters
+ ----------
+ refinement_style : str
+ How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
+ number_cells : int
+ Number of cells in the mesh
+
+ Returns
+ -------
+ Mesh
+ One-dimensional mesh, ready for use in FEniCSx
+ '''
+ if refinement_style == 'log':
+ coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
+ elif refinement_style == 'hard_log':
+ coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
+ coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
+ coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+ elif refinement_style == 'hard_hard_log':
+ coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
+ coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
+ coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+ num_vertices = len(coordinates_np)
+ num_cells = num_vertices - 1
+ cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
+ gdim = 1
+ shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
+ degree = 1
+ domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
+ coordinates_np_ = []
+ [coordinates_np_.append([coord]) for coord in coordinates_np]
+ msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
+ return msh
+
+def solve_System_2eq(phi_left:float, phi_right:float, p_right:float, z_A:float, z_C:float, y_A_R:float, y_C_R:float, K:float|str, Lambda2:float, a2:float, number_cells:int, solvation:float = 0, relax_param:float=None, x0:float=0, x1:float=1, refinement_style:str='uniform', return_type:str='Scalar', rtol:float=1e-8, max_iter:float=500):
+ '''
+ Solve the simplified dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
+
+ System of equations:
+ λ²Δ φ =−L²n^F
+
+ a²∇p=−n^F∇ φ
+
+ with the space charge
+ n^F = z_A y_A(φ, p) + z_C y_C(φ ,p)
+
+ if the mixture is compressible:
+ y_alpha = C_alpha * (K+p−1)^(−κ+1)a²K exp(−z_α φ)
+
+ if the mixture is incompressible:
+ y_alpha = D_alpha * exp(−(κ+1)a²p−z_α φ)
+
+ with φ the electric potential, p the pressure, n^F the total free charge density, J_α the diffusion fluxes of species α, λ² a dimensionless parameter, L²=1, a² a dimensionless parameter, N the number of species, and α the species index.
+
+ ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
+
+ Parameters
+ ----------
+ phi_left : float
+ Value of φ at the left boundary
+ phi_right : float
+ Value of φ at the right boundary
+ p_right : float
+ Value of p at the right boundary
+ z_A : float
+ Charge number of species A
+ z_C : float
+ Charge number of species C
+ y_A_R : float
+ Atomic fractions of species A at right boundary
+ y_C_R : float
+ Atomic fractions of species C at right boundary
+ K : float | str
+ Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
+ Lambda2 : float
+ Dimensionless parameter
+ a2 : float
+ Dimensionless parameter
+ number_cells : int
+ Number of cells in the mesh
+ solvation : float, optional
+ solvation number, by default 0
+ relax_param : float, optional
+ Relaxation parameter for the Newton solver
+ xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
+ , by default None -> Determined automatically
+ x0 : float, optional
+ Left boundary of the domain, by default 0
+ x1 : float, optional
+ Right boundary of the domain, by default 1
+ refinement_style : str, optional
+ Specify for refinement towards zero
+ Options are 'uniform', 'log', 'hard_log', 'hard_hard_log' by default 'uniform'
+ return_type : str, optional
+ 'Vector' or 'Scalar', 'Scalar' returns dolfinx.fem type and 'Vector' numpy arrays of the solution, by default 'Scalar'
+ rtol : float, optional
+ Relative tolerance for Newton solver, by default 1e-8
+ max_iter : float, optional
+ Maximum number of Newton iterations, by default 500
+
+ Returns
+ -------
+ y_A, y_C, phi, p, msh
+ Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
+ If return_type is 'Vector', the solution is returned as numpy arrays
+ '''
+ if return_type == 'Scalar':
+ raise NotImplementedError('Scalar return type is not implemented yet')
+ # Define boundaries of the domain
+ x0 = 0
+ x1 = 1
+
+ # Define boundaries for the boundary conditions
+ def Left(x):
+ return np.isclose(x[0], x0)
+
+ def Right(x):
+ return np.isclose(x[0], x1)
+
+ # Create mesh
+ if refinement_style == 'uniform':
+ msh = mesh.create_unit_interval(MPI.COMM_WORLD, number_cells, dtype=np.float64)
+ else:
+ msh = create_refined_mesh(refinement_style, number_cells)
+
+ # Define Finite Elements
+ CG1_elem = element('Lagrange', msh.basix_cell(), 1)
+
+ # Define Mixed Function Space
+ W_elem = mixed_element([CG1_elem, CG1_elem])#, CG1_elem, CG1_elem])
+ W = fem.functionspace(msh, W_elem)
+
+ # Define Trial- and Testfunctions
+ u = fem.Function(W)
+ phi, p = split(u)
+ (v_1, v_2) = TestFunctions(W)
+
+ # Collapse function space for bcs
+ W0, _ = W.sub(0).collapse()
+ W1, _ = W.sub(1).collapse()
+
+ # Define boundary conditions values
+ def phi_left_(x):
+ return np.full_like(x[0], phi_left)
+ def phi_right_(x):
+ return np.full_like(x[0], phi_right)
+ def p_right_(x):
+ return np.full_like(x[0], p_right)
+
+ # Interpolate bcs functions
+ phi_left_bcs = fem.Function(W0)
+ phi_left_bcs.interpolate(phi_left_)
+ phi_right_bcs = fem.Function(W0)
+ phi_right_bcs.interpolate(phi_right_)
+ p_right_bcs = fem.Function(W1)
+ p_right_bcs.interpolate(p_right_)
+
+ # Identify dofs for boundary conditions
+ # Define boundary conditions
+ facet_left_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Left)
+ facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
+ bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(0))
+ bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(0))
+
+ facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
+ bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(1))
+
+
+ # Combine boundary conditions into list
+ bcs = [bc_left_phi, bc_right_phi, bc_right_p]
+
+ def y_A(phi, p):
+ D_A = y_A_R / exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
+ return D_A * exp(-(solvation + 1) * a2 * p - z_A * phi)
+
+ def y_C(phi, p):
+ D_C = y_C_R / exp(-(solvation + 1) * a2 * p_right - z_C * phi_right)
+ return D_C * exp(-(solvation + 1) * a2 * p - z_C * phi)
+
+
+ # Define variational problem
+ if K == 'incompressible':
+ # total free charge density
+ def nF(y_A, y_C, p):
+ return (z_C * y_C + z_A * y_A)
+ else:
+ # total number density
+ def n(p):
+ return (p-1)/K + 1
+
+ # total free charge density
+ def nF(y_A, y_C, p):
+ return (z_C * y_C + z_A * y_A) * n(p)
+ # Variational Form
+ A = (
+ inner(grad(phi), grad(v_1)) * dx
+ - 1 / Lambda2 * nF(y_A(phi, p), y_C(phi, p), p) * v_1 * dx
+ ) + (
+ inner(grad(p), grad(v_2)) * dx
+ + 1 / a2 * nF(y_A(phi, p), y_C(phi, p), p) * dot(grad(phi), grad(v_2)) * dx
+ )
+ F = A
+
+ # Define Nonlinear Problem
+ problem = NonlinearProblem(F, u, bcs=bcs)
+
+ # Define Newton Solver and solver settings
+ solver = NewtonSolver(MPI.COMM_WORLD, problem)
+ solver.convergence_criterion = "incremental"
+ solver.rtol = rtol
+ if relax_param != None:
+ solver.relaxation_parameter = relax_param
+ else:
+ if phi_right == phi_left:
+ solver.relaxation_parameter = 1.0
+ else:
+ solver.relaxation_parameter = 1/(np.abs(phi_right-phi_left)**(5/4))
+ solver.max_it = max_iter
+ solver.report = True
+
+ # Solve the problem
+ log.set_log_level(log.LogLevel.INFO)
+ n, converged = solver.solve(u)
+ assert (converged)
+ print(f"Number of interations: {n:d}")
+
+ # Split the mixed function space into the individual components
+ phi, p = u.split()
+
+ # Return the solution
+ if return_type=='Vector':
+ x_vals = np.array(msh.geometry.x[:,0])
+ phi_vals = np.array(u.sub(0).collapse().x.array)
+ p_vals = np.array(u.sub(1).collapse().x.array)
+
+ # Calculate the atomic fractions
+ D_A = y_A_R / np.exp(-(solvation + 1) * a2 * p_right - z_A * phi_right)
+ y_A_vals = D_A * np.exp(-(solvation + 1) * a2 * p_vals - z_A * phi_vals)
+
+ D_C = y_C_R / np.exp(-(solvation + 1) * a2 * p_right - z_C * phi_right)
+ y_C_vals = D_C * np.exp(-(solvation + 1) * a2 * p_vals - z_C * phi_vals)
+
+ return y_A_vals, y_C_vals, phi_vals, p_vals, x_vals
+
+if __name__ == '__main__':
+ # Define the parameters
+ phi_left = 5.0
+ phi_right = 0.0
+ p_right = 0.0
+ y_A_R = 1/3
+ y_C_R = 1/3
+ z_A = -1.0
+ z_C = 1.0
+ K = 'incompressible'
+ Lambda2 = 8.553e-6
+ a2 = 7.5412e-4
+ number_cells = 1024
+ relax_param = .1
+ rtol = 1e-4
+ max_iter = 500
+
+ # Solve the system
+ y_A, y_C, phi, p, x = solve_System_2eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, relax_param=relax_param, x0=0, x1=1, refinement_style='uniform', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+ # Plot the solution
+ plt.plot(x, phi)
+ plt.xlim(0,0.05)
+ plt.grid()
+ plt.xlabel('x [-]')
+ plt.ylabel('$\\varphi$ [-]')
+ plt.show()
+
+ plt.plot(x, y_A, '--', color='tab:blue', label='$y_A$')
+ plt.plot(x, y_C, '-', color='tab:blue', label='$y_C$')
+ plt.plot(x, 1 - y_A - y_C, ':', color='tab:blue', label='$y_S$')
+ plt.xlim(0,0.05)
+ plt.legend()
+ plt.grid()
+ plt.xlabel('x [-]')
+ plt.ylabel('$y_\\alpha$ [-]')
+ plt.show()
+
+ plt.plot(x, p)
+ plt.xlim(0,0.05)
+ plt.grid()
+ plt.xlabel('x [-]')
+ plt.ylabel('$p$ [-]')
+ plt.show()
diff --git a/build/lib/Bsc-ElectrolyteModels/Eq04.py b/build/lib/Bsc-ElectrolyteModels/Eq04.py
new file mode 100644
index 0000000000000000000000000000000000000000..4856d848e3847f7229ac43754add17cc33cb8f61
--- /dev/null
+++ b/build/lib/Bsc-ElectrolyteModels/Eq04.py
@@ -0,0 +1,358 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+
+This script implements the fenics solver for a ternary electrolyte (N=3, A,C,S)
+'''
+
+import numpy as np
+from mpi4py import MPI
+from dolfinx import mesh, fem, log
+from dolfinx.fem.petsc import NonlinearProblem
+from dolfinx.nls.petsc import NewtonSolver
+from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Mesh
+from basix.ufl import element, mixed_element
+import matplotlib.pyplot as plt
+
+# Define mesh
+def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
+ '''
+ Creates a one-dimensional mesh with a refined region at the left boundary
+
+ Parameters
+ ----------
+ refinement_style : str
+ How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
+ number_cells : int
+ Number of cells in the mesh
+
+ Returns
+ -------
+ Mesh
+ One-dimensional mesh, ready for use in FEniCSx
+ '''
+ if refinement_style == 'log':
+ coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
+ elif refinement_style == 'hard_log':
+ coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
+ coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
+ coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+ elif refinement_style == 'hard_hard_log':
+ coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
+ coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
+ coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+ num_vertices = len(coordinates_np)
+ num_cells = num_vertices - 1
+ cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
+ gdim = 1
+ shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
+ degree = 1
+ domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
+ coordinates_np_ = []
+ [coordinates_np_.append([coord]) for coord in coordinates_np]
+ msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
+ return msh
+
+def solve_System_4eq(phi_left:float, phi_right:float, p_right:float, z_A:float, z_C:float, y_A_R:float, y_C_R:float, K:float|str, Lambda2:float, a2:float, number_cells:int, solvation:float = 0, PoissonBoltzmann:bool=False, relax_param:float=None, x0:float=0, x1:float=1, refinement_style:str='uniform', return_type:str='Scalar', rtol:float=1e-8, max_iter:float=500):
+ '''
+ Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
+
+ System of equations:
+ λ²Δ φ =−L²n^F
+
+ a²∇p=−n^F∇ φ
+
+ div(J_A)=0
+
+ div(J_C)=0
+
+ with φ the electric potential, p the pressure, n^F the total free charge density, J_α the diffusion fluxes of species α, λ² a dimensionless parameter, L²=1, a² a dimensionless parameter, N the number of species, and α the species index.
+
+ ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
+
+ Parameters
+ ----------
+ phi_left : float
+ Value of φ at the left boundary
+ phi_right : float
+ Value of φ at the right boundary
+ p_right : float
+ Value of p at the right boundary
+ z_A : float
+ Charge number of species A
+ z_C : float
+ Charge number of species C
+ y_A_R : float
+ Atomic fractions of species A at right boundary
+ y_C_R : float
+ Atomic fractions of species C at right boundary
+ K : float | str
+ Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
+ Lambda2 : float
+ Dimensionless parameter
+ a2 : float
+ Dimensionless parameter
+ number_cells : int
+ Number of cells in the mesh
+ solvation : float, optional
+ solvation number, by default 0
+ PoissonBoltzmann : bool, optional
+ Solve classical Nernst-Planck model with the use of the Poisson-Boltzmann formulation if True, else solve the presented model by Dreyer, Guhlke, Müller, by default False
+ relax_param : float, optional
+ Relaxation parameter for the Newton solver
+ xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
+ , by default None -> Determined automatically
+ x0 : float, optional
+ Left boundary of the domain, by default 0
+ x1 : float, optional
+ Right boundary of the domain, by default 1
+ refinement_style : str, optional
+ Specify for refinement towards zero
+ Options are 'uniform', 'log', 'hard_log', 'hard_hard_log' by default 'uniform'
+ return_type : str, optional
+ 'Vector' or 'Scalar', 'Scalar' returns dolfinx.fem type and 'Vector' numpy arrays of the solution, by default 'Scalar'
+ rtol : float, optional
+ Relative tolerance for Newton solver, by default 1e-8
+ max_iter : float, optional
+ Maximum number of Newton iterations, by default 500
+
+ Returns
+ -------
+ y_A, y_C, phi, p, msh
+ Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
+ If return_type is 'Vector', the solution is returned as numpy arrays
+ '''
+ # Define boundaries of the domain
+ x0 = 0
+ x1 = 1
+
+ # Define boundaries for the boundary conditions
+ def Left(x):
+ return np.isclose(x[0], x0)
+
+ def Right(x):
+ return np.isclose(x[0], x1)
+
+ # Create mesh
+ if refinement_style == 'uniform':
+ msh = mesh.create_unit_interval(MPI.COMM_WORLD, number_cells, dtype=np.float64)
+ else:
+ msh = create_refined_mesh(refinement_style, number_cells)
+
+ # Define Finite Elements
+ CG1_elem = element('Lagrange', msh.basix_cell(), 1)
+
+ # Define Mixed Function Space
+ W_elem = mixed_element([CG1_elem, CG1_elem, CG1_elem, CG1_elem])
+ W = fem.functionspace(msh, W_elem)
+
+ # Define Trial- and Testfunctions
+ u = fem.Function(W)
+ y_A, y_C, phi, p = split(u)
+ (v_A, v_C, v_1, v_2) = TestFunctions(W)
+
+ # Collapse function space for bcs
+ W0, _ = W.sub(0).collapse()
+ W1, _ = W.sub(1).collapse()
+ W2, _ = W.sub(2).collapse()
+ W3, _ = W.sub(3).collapse()
+
+ # Define boundary conditions values
+ def phi_left_(x):
+ return np.full_like(x[0], phi_left)
+ def phi_right_(x):
+ return np.full_like(x[0], phi_right)
+ def p_right_(x):
+ return np.full_like(x[0], p_right)
+ def y_A_right_(x):
+ return np.full_like(x[0], y_A_R)
+ def y_C_right_(x):
+ return np.full_like(x[0], y_C_R)
+
+ # Interpolate bcs functions
+ phi_left_bcs = fem.Function(W2)
+ phi_left_bcs.interpolate(phi_left_)
+ phi_right_bcs = fem.Function(W2)
+ phi_right_bcs.interpolate(phi_right_)
+ p_right_bcs = fem.Function(W3)
+ p_right_bcs.interpolate(p_right_)
+ y_A_right_bcs = fem.Function(W0)
+ y_A_right_bcs.interpolate(y_A_right_)
+ y_C_right_bcs = fem.Function(W1)
+ y_C_right_bcs.interpolate(y_C_right_)
+
+ # Identify dofs for boundary conditions
+ # Define boundary conditions
+ facet_left_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Left)
+ facet_right_dofs = fem.locate_dofs_geometrical((W.sub(2), W.sub(2).collapse()[0]), Right)
+ bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(2))
+ bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(2))
+
+ facet_right_dofs = fem.locate_dofs_geometrical((W.sub(3), W.sub(3).collapse()[0]), Right)
+ bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(3))
+
+ facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
+ bc_right_y_A = fem.dirichletbc(y_A_right_bcs, facet_right_dofs, W.sub(0))
+
+ facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
+ bc_right_y_C = fem.dirichletbc(y_C_right_bcs, facet_right_dofs, W.sub(1))
+
+ # Combine boundary conditions into list
+ bcs = [bc_left_phi, bc_right_phi, bc_right_p, bc_right_y_A, bc_right_y_C]
+
+
+
+ # Define variational problem
+ if K == 'incompressible':
+ # total free charge density
+ def nF(y_A, y_C):
+ return (z_C * y_C + z_A * y_A)
+
+ # Diffusion fluxes for species A and C
+ def J_A(y_A, y_C, phi, p):
+ return grad(ln(y_A) + a2 * (p - 1) * (solvation + 1) + z_A * phi)
+ # return grad(ln(y_A) + a2 * (p - 1) - solvation * ln(1-y_A-y_C) + z_A * phi)
+
+ def J_C(y_A, y_C, phi, p):
+ return grad(ln(y_C) + a2 * (p - 1) * (solvation + 1) + z_C * phi)
+ # return grad(ln(y_C) + a2 * (p - 1) - solvation * ln(1-y_A-y_C) + z_C * phi)
+
+ # Variational Form
+ A = (
+ inner(grad(phi), grad(v_1)) * dx
+ - 1 / Lambda2 * nF(y_A, y_C) * v_1 * dx
+ ) + (
+ inner(grad(p), grad(v_2)) * dx
+ + 1 / a2 * nF(y_A, y_C) * dot(grad(phi), grad(v_2)) * dx
+ ) + (
+ inner(J_A(y_A, y_C, phi, p), grad(v_A)) * dx
+ + inner(J_C(y_A, y_C, phi, p), grad(v_C)) * dx
+ )
+ if PoissonBoltzmann:
+ A += (
+ inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_A)) * dx
+ + inner(grad(- a2 * (p - 1) * (solvation + 1)), grad(v_C)) * dx
+ )
+ else:
+ # total number density
+ def n(p):
+ return (p-1)/K + 1
+
+ # total free charge density
+ def nF(y_A, y_C, p):
+ return (z_C * y_C + z_A * y_A) * n(p)
+
+ # Diffusion fluxes for species A and C
+ def J_A(y_A, y_C, phi, p):
+ return ln(y_A) + a2 * (solvation + 1) * K * ln(1 + 1/K * (p-1)) + z_A * phi
+
+ def J_C(y_A, y_C, phi, p):
+ return ln(y_C) + a2 * (solvation + 1)* K * ln(1 + 1/K * (p-1)) + z_C * phi
+
+ A = (
+ inner(grad(phi), grad(v_1)) * dx
+ - 1 / Lambda2 * nF(y_A, y_C, p) * v_1 * dx
+ ) + (
+ inner(grad(p), grad(v_2)) * dx
+ + 1 / a2 * nF(y_A, y_C, p) * dot(grad(phi), grad(v_2)) * dx
+ ) + (
+ inner(grad(J_A(y_A, y_C, phi, p)), grad(v_A)) * dx
+ + inner(grad(J_C(y_A, y_C, phi, p)), grad(v_C)) * dx
+ )
+ F = A
+
+ # Initialize initial guess for u
+ y_C_init = fem.Function(W1)
+ y_A_init = fem.Function(W0)
+ y_C_init.interpolate(lambda x: np.full_like(x[0], y_C_R))
+ y_A_init.interpolate(lambda x: np.full_like(x[0], y_A_R))
+
+ with u.vector.localForm() as u_loc:
+ u_loc.set(0)
+ u.sub(0).interpolate(y_A_init)
+ u.sub(1).interpolate(y_C_init)
+
+ # Define Nonlinear Problem
+ problem = NonlinearProblem(F, u, bcs=bcs)
+
+ # Define Newton Solver and solver settings
+ solver = NewtonSolver(MPI.COMM_WORLD, problem)
+ solver.convergence_criterion = "incremental"
+ solver.rtol = rtol
+ if relax_param != None:
+ solver.relaxation_parameter = relax_param
+ else:
+ if phi_right == phi_left:
+ solver.relaxation_parameter = 1.0
+ else:
+ solver.relaxation_parameter = 1/(np.abs(phi_right-phi_left)**(5/4))
+ solver.max_it = max_iter
+ solver.report = True
+
+ # Solve the problem
+ log.set_log_level(log.LogLevel.INFO)
+ n, converged = solver.solve(u)
+ assert (converged)
+ print(f"Number of interations: {n:d}")
+
+ # Split the mixed function space into the individual components
+ y_A, y_C, phi, p = u.split()
+
+ # Return the solution
+ if return_type=='Vector':
+ x_vals = np.array(msh.geometry.x[:,0])
+ y_A_vals = np.array(u.sub(0).collapse().x.array)
+ y_C_vals = np.array(u.sub(1).collapse().x.array)
+ phi_vals = np.array(u.sub(2).collapse().x.array)
+ p_vals = np.array(u.sub(3).collapse().x.array)
+
+ return y_A_vals, y_C_vals, phi_vals, p_vals, x_vals
+ elif return_type=='Scalar':
+ return y_A, y_C, phi, p, msh
+
+
+if __name__ == '__main__':
+ # Define the parameters
+ phi_left = 10.0
+ phi_right = 0.0
+ p_right = 0.0
+ y_A_R = 0.01#1/3
+ y_C_R = 0.01#1/3
+ z_A = -1.0
+ z_C = 1.0
+ K = 'incompressible'
+ Lambda2 = 8.553e-6
+ a2 = 7.5412e-4
+ solvation = 0
+ number_cells = 1024
+ relax_param = .1
+ rtol = 1e-4
+ max_iter = 500
+
+ # Solve the system
+ y_A, y_C, phi, p, x = solve_System_4eq(phi_left, phi_right, p_right, z_A, z_C, y_A_R, y_C_R, K, Lambda2, a2, number_cells, solvation=solvation, relax_param=relax_param, x0=0, x1=1, refinement_style='uniform', return_type='Vector', max_iter=max_iter, rtol=rtol)
+
+ # Plot the solution
+ plt.plot(x, phi)
+ plt.xlim(0,0.05)
+ plt.grid()
+ plt.xlabel('x [-]')
+ plt.ylabel('$\\varphi$ [-]')
+ plt.show()
+
+ plt.plot(x, y_A, '--', color='tab:blue', label='$y_A$')
+ plt.plot(x, y_C, '-', color='tab:blue', label='$y_C$')
+ plt.plot(x, 1 - y_A - y_C, ':', color='tab:blue', label='$y_S$')
+ plt.xlim(0,0.05)
+ plt.legend()
+ plt.grid()
+ plt.xlabel('x [-]')
+ plt.ylabel('$y_\\alpha$ [-]')
+ plt.show()
+
+ plt.plot(x, p)
+ plt.xlim(0,0.05)
+ plt.grid()
+ plt.xlabel('x [-]')
+ plt.ylabel('$p$ [-]')
+ plt.show()
diff --git a/build/lib/Bsc-ElectrolyteModels/EqN.py b/build/lib/Bsc-ElectrolyteModels/EqN.py
new file mode 100644
index 0000000000000000000000000000000000000000..9b48ed187e5cf869fa7338df76a185aa1c19b42a
--- /dev/null
+++ b/build/lib/Bsc-ElectrolyteModels/EqN.py
@@ -0,0 +1,324 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+
+This script implements the fenics solver for the generic system of equations for N species
+'''
+
+import numpy as np
+from mpi4py import MPI
+from dolfinx import mesh, fem, log
+from dolfinx.fem.petsc import NonlinearProblem
+from dolfinx.nls.petsc import NewtonSolver
+from ufl import TestFunctions, split, dot, grad, dx, inner, ln, Mesh
+from basix.ufl import element, mixed_element
+import matplotlib.pyplot as plt
+
+# Define mesh
+def create_refined_mesh(refinement_style:str, number_cells:int) -> Mesh:
+ '''
+ Creates a one-dimensional mesh with a refined region at the left boundary
+
+ Parameters
+ ----------
+ refinement_style : str
+ How the mesh should be refined. Options are 'log', 'hard_log', 'hard_hard_log'
+ number_cells : int
+ Number of cells in the mesh
+
+ Returns
+ -------
+ Mesh
+ One-dimensional mesh, ready for use in FEniCSx
+ '''
+ if refinement_style == 'log':
+ coordinates_np = (np.logspace(0, 1, number_cells+1) - 1) / 9
+ elif refinement_style == 'hard_log':
+ coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.1
+ coordinates_np2 = 0.1 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.9
+ coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+ elif refinement_style == 'hard_hard_log':
+ coordinates_np1 = (np.logspace(0,1,int(number_cells*0.9)+1,endpoint=False)-1)/9 * 0.004
+ coordinates_np2 = 0.004 + (np.logspace(0,1,int(number_cells*0.1)+1)-1)/9 * 0.996
+ coordinates_np = np.concatenate((coordinates_np1, coordinates_np2), axis=0)
+ num_vertices = len(coordinates_np)
+ num_cells = num_vertices - 1
+ cells_np = np.column_stack((np.arange(num_cells), np.arange(1, num_cells+1)))
+ gdim = 1
+ shape = 'interval' # 'interval', 'triangle', 'quadrilateral', 'tetrahedron', 'hexahedron'
+ degree = 1
+ domain = Mesh(element("Lagrange", shape, 1, shape=(1,)))
+ coordinates_np_ = []
+ [coordinates_np_.append([coord]) for coord in coordinates_np]
+ msh = mesh.create_mesh(MPI.COMM_WORLD, cells_np, coordinates_np_, domain)
+ return msh
+
+def solve_System_Neq(phi_left:float, phi_right:float, p_right:float, z_alpha:list, y_R:list, K:float|str, Lambda2:float, a2:float, number_cells:int, solvation:float = 0, PoissonBoltzmann:bool=False, relax_param:float=None, x0:float=0, x1:float=1, refinement_style:str='uniform', return_type:str='Vector', rtol:float=1e-8, max_iter:float=500):
+ '''
+ Solve the dimensionless system of equations presented in: Numerical Treatment of a Thermodynamically Consistent Electrolyte Model, B.Sc. Thesis Habscheid 2024
+
+ System of equations:
+ λ²Δ φ =−L²n^F
+
+ a²∇p=−n^F∇ φ
+
+ div(J_α)=0 α∈ {1,...,N−1}
+
+ with φ the electric potential, p the pressure, n^F the total free charge density, J_α the diffusion fluxes of species α, λ² a dimensionless parameter, L²=1, a² a dimensionless parameter, N the number of species, and α the species index.
+
+ ! If the Newton solver diverges, you may try to reduce the relaxation parameter.
+
+ Parameters
+ ----------
+ phi_left : float
+ Value of φ at the left boundary
+ phi_right : float
+ Value of φ at the right boundary
+ p_right : float
+ Value of p at the right boundary
+ z_alpha : list
+ Charge numbers for species α = 1,...,N-1
+ y_R : list
+ Atomic fractions at right boundary for species α = 1,...,N-1
+ K : float | str
+ Dimensioness bulk modulus of the electrolyte. If 'incompressible', the system is solved for an incompressible electrolyte
+ Lambda2 : float
+ Dimensionless parameter
+ a2 : float
+ Dimensionless parameter
+ number_cells : int
+ Number of cells in the mesh
+ solvation : float, optional
+ solvation number, not implemented yet, by default 0
+ PoissonBoltzmann : bool, optional
+ Solve classical Nernst-Planck model with the use of the Poisson-Boltzmann formulation if True, else solve the presented model by Dreyer, Guhlke, Müller, Not implemented yet, by default False
+ relax_param : float, optional
+ Relaxation parameter for the Newton solver
+ xₙ₊₁ = γ xₙ f(xₙ)/f'(xₙ) with γ the relaxation parameter
+ , by default None -> Determined automatically
+ x0 : float, optional
+ Left boundary of the domain, by default 0
+ x1 : float, optional
+ Right boundary of the domain, by default 1
+ refinement_style : str, optional
+ Specify for refinement towards zero
+ Options are 'uniform', 'log', 'hard_log', 'hard_hard_log' by default 'uniform'
+ return_type : str, optional
+ 'Vector' or 'Scalar' (not implemented yet, should be implemented in a later version), 'Scalar' returns dolfinx.fem type and 'Vector' numpy arrays of the solution, by default 'Vector'
+ rtol : float, optional
+ Relative tolerance for Newton solver, by default 1e-8
+ max_iter : float, optional
+ Maximum number of Newton iterations, by default 500
+
+ Returns
+ -------
+ y_A, y_C, phi, p, msh
+ Returns atomic fractions for species A and C, electric potential, pressure, and the mesh
+ If return_type is 'Vector', the solution is returned as numpy arrays
+ Only return_type 'Vector' is implemented yet
+ '''
+ if solvation != 0:
+ raise NotImplementedError('Solvation number not implemented yet')
+ if PoissonBoltzmann:
+ raise NotImplementedError('Poisson-Boltzmann not implemented yet')
+ # Define boundaries of the domain
+ x0 = 0
+ x1 = 1
+
+ # Define boundaries
+ def Left(x):
+ return np.isclose(x[0], x0)
+
+ def Right(x):
+ return np.isclose(x[0], x1)
+
+ # Create mesh
+ if refinement_style == 'uniform':
+ msh = mesh.create_unit_interval(MPI.COMM_WORLD, number_cells, dtype=np.float64)
+ else:
+ msh = create_refined_mesh(refinement_style, number_cells)
+
+ # Define Finite Elements
+ CG1_elem = element('Lagrange', msh.basix_cell(), 1)
+
+ # Define Mixed Function Space
+ Elem_list = [CG1_elem, CG1_elem]
+ [Elem_list.append(CG1_elem) for _ in range(len(z_alpha))]
+ W_elem = mixed_element(Elem_list)
+ W = fem.functionspace(msh, W_elem)
+
+ # Define Trial- and Testfunctions
+ u = fem.Function(W)
+ my_TrialFunctions = split(u)
+ my_TestFunctions = TestFunctions(W)
+
+ phi, p = my_TrialFunctions[0], my_TrialFunctions[1]
+ v_1, v_2 = my_TestFunctions[0], my_TestFunctions[1]
+ y_alpha = my_TrialFunctions[2:]
+ v_alpha = my_TestFunctions[2:]
+
+ # Collapse function space for bcs
+ W_ = []
+ [W_.append(W.sub(i).collapse()[0]) for i in range(len(z_alpha)+2)]
+
+ # Define boundary conditions values
+ def phi_left_(x):
+ return np.full_like(x[0], phi_left)
+ def phi_right_(x):
+ return np.full_like(x[0], phi_right)
+ def p_right_(x):
+ return np.full_like(x[0], p_right)
+
+ # Interpolate bcs functions
+ phi_left_bcs = fem.Function(W_[0])
+ phi_left_bcs.interpolate(phi_left_)
+ phi_right_bcs = fem.Function(W_[0])
+ phi_right_bcs.interpolate(phi_right_)
+ p_right_bcs = fem.Function(W_[1])
+ p_right_bcs.interpolate(p_right_)
+
+ # Identify dofs for boundary conditions
+ facet_left_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Left)
+ facet_right_dofs = fem.locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]), Right)
+ bc_left_phi = fem.dirichletbc(phi_left_bcs, facet_left_dofs, W.sub(0))
+ bc_right_phi = fem.dirichletbc(phi_right_bcs, facet_right_dofs, W.sub(0))
+
+ facet_right_dofs = fem.locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]), Right)
+ bc_right_p = fem.dirichletbc(p_right_bcs, facet_right_dofs, W.sub(1))
+
+ # Combine boundary conditions for electric potential and pressure into list
+ bcs = [bc_left_phi, bc_right_phi, bc_right_p]
+
+ # Repeat the same for the boundary conditoins for the atomic fractions
+ for i in range(len(z_alpha)):
+ y_right_bcs = fem.Function(W_[i+2])
+ def y_right_(x):
+ return np.full_like(x[0], y_R[i])
+ y_right_bcs.interpolate(y_right_)
+ facet_right_dofs = fem.locate_dofs_geometrical((W.sub(i+2), W.sub(i+2).collapse()[0]), Right)
+ bc_right_y = fem.dirichletbc(y_right_bcs, facet_right_dofs, W.sub(i+2))
+ bcs.append(bc_right_y)
+
+ # Define variational problem
+ if K == 'incompressible':
+ # total free charge density
+ def nF(y_alpha):
+ nF = 0
+ for i in range(len(z_alpha)):
+ nF += z_alpha[i] * y_alpha[i]
+ return nF
+
+ # Diffusion fluxes for species A and C
+ def J_alpha(y_alpha, alpha, phi, p):
+ mu_alpha = ln(y_alpha[alpha])
+ mu_S = ln(1 - sum(y_alpha))
+ return grad(mu_alpha - mu_S + z_alpha[alpha] * phi)
+
+ # Variational Form
+ A = (
+ inner(grad(phi), grad(v_1)) * dx
+ - 1 / Lambda2 * nF(y_alpha) * v_1 * dx
+ ) + (
+ inner(grad(p), grad(v_2)) * dx
+ + 1 / a2 * nF(y_alpha) * dot(grad(phi), grad(v_2)) * dx
+ )
+ for alpha in range(len(z_alpha)):
+ A += (
+ inner(J_alpha(y_alpha, alpha, phi, p), grad(v_alpha[alpha])) * dx
+ )
+ if PoissonBoltzmann:
+ raise ValueError('Poisson-Boltzmann not implemented for incompressible systems')
+ else:
+ raise ValueError('Only incompressible systems are implemented')
+ F = A
+
+ # Initialize initial guess for u
+ with u.vector.localForm() as u_loc:
+ u_loc.set(0)
+
+ # Initialize initial guess for u
+ for alpha in range(len(z_alpha)):
+ y_alpha_init = fem.Function(W_[alpha+2])
+ y_alpha_init.interpolate(lambda x: np.full_like(x[0], y_R[alpha]))
+ u.sub(alpha+2).interpolate(y_alpha_init)
+
+ # Define Nonlinear Problem
+ problem = NonlinearProblem(F, u, bcs=bcs)
+
+ # Define Newton Solver
+ solver = NewtonSolver(MPI.COMM_WORLD, problem)
+ solver.convergence_criterion = "incremental"
+ solver.rtol = rtol
+ if relax_param != None:
+ solver.relaxation_parameter = relax_param
+ else:
+ if phi_right == phi_left:
+ solver.relaxation_parameter = 1.0
+ else:
+ solver.relaxation_parameter = 1/(np.abs(phi_right-phi_left)**(5/4))
+ solver.max_it = max_iter
+ solver.report = True
+
+
+ log.set_log_level(log.LogLevel.INFO)
+ n, converged = solver.solve(u)
+ assert (converged)
+ print(f"Number of (interations: {n:d}")
+
+ # Return the solution
+ if return_type=='Vector':
+ x = np.array(msh.geometry.x[:,0])
+ phi = np.array(u.sub(0).collapse().x.array)
+ p = np.array(u.sub(1).collapse().x.array)
+ y = []
+ [y.append(u.sub(i+2).collapse().x.array) for i in range(len(z_alpha))]
+ y = np.array(y)
+ return y, phi, p, x
+
+
+if __name__ == '__main__':
+ # Define the parameters
+ phi_left = 8.0
+ phi_right = 0.0
+ p_right = 0.0
+ y_R = [3/6, 1/6, 1/6]
+ z_alpha = [-1.0, 1.0, 2.0]
+ K = 'incompressible'
+ Lambda2 = 8.553e-6
+ a2 = 7.5412e-4
+ number_cells = 1024
+ relax_param = .05
+ rtol = 1e-4
+ max_iter = 2_500
+ refinement_style = 'hard_log'
+ return_type = 'Vector'
+
+ # Solve the system
+ y, phi, p, x = solve_System_Neq(phi_left, phi_right, p_right, z_alpha, y_R, K, Lambda2, a2, number_cells, relax_param=relax_param, refinement_style=refinement_style, return_type=return_type, max_iter=max_iter, rtol=rtol)
+
+ # Plot the solution
+ plt.figure()
+ plt.plot(x, phi)
+ plt.xlim(0,0.05)
+ plt.grid()
+ plt.xlabel('x [-]')
+ plt.ylabel('$\\varphi$ [-]')
+ plt.show()
+
+ plt.figure()
+ plt.xlim(0,0.05)
+ plt.plot(x, p)
+ plt.grid()
+ plt.xlabel('x [-]')
+ plt.ylabel('$p$ [-]')
+ plt.show()
+
+ plt.figure()
+ for i in range(len(z_alpha)):
+ plt.plot(x, y[i], label=f'$y_{i}$')
+ plt.xlim(0,0.05)
+ plt.legend()
+ plt.grid()
+ plt.xlabel('x [-]')
+ plt.ylabel('$y_i$ [-]')
+ plt.show()
\ No newline at end of file
diff --git a/build/lib/Bsc-ElectrolyteModels/Helper_DoubleLayerCapacity.py b/build/lib/Bsc-ElectrolyteModels/Helper_DoubleLayerCapacity.py
new file mode 100644
index 0000000000000000000000000000000000000000..e47b37efa121e9c1b5ae89707ed25366f1a01fcc
--- /dev/null
+++ b/build/lib/Bsc-ElectrolyteModels/Helper_DoubleLayerCapacity.py
@@ -0,0 +1,411 @@
+'''
+Jan Habscheid
+Jan.Habscheid@rwth-aachen.de
+'''
+
+import numpy as np
+from scipy.optimize import fixed_point, fsolve
+
+
+# Helper functions
+def Phi_pot_center(Phi_pot:np.ndarray) -> np.ndarray:
+ '''
+ Returns vector with the center of the electric potential values.
+
+ Parameters
+ ----------
+ Phi_pot : np.ndarray
+ Input vector with electric potential values.
+
+ Returns
+ -------
+ np.ndarray
+ Vector with the center of the electric potential values.
+ '''
+ return (Phi_pot[1:] + Phi_pot[:-1]) / 2
+
+def dx(Phi_pot:np.ndarray) -> float:
+ '''
+ Returns the difference between the first two electric potential values.
+
+ Assumes that the electric potential values are equally spaced.
+
+ Parameters
+ ----------
+ Phi_pot : np.ndarray
+ Input vector with electric potential values.
+
+ Returns
+ -------
+ float
+ Difference between the first two electric potential values.
+ '''
+ return Phi_pot[1] - Phi_pot[0]
+
+def C_dl(Q_DL:np.ndarray, Phi_pot:np.ndarray) -> np.ndarray:
+ '''
+ Double Layer Capacity
+
+ Parameters
+ ----------
+ Q_DL : np.ndarray
+ Charge of the system
+ Phi_pot : np.ndarray
+ Electric potential values
+
+ Returns
+ -------
+ np.ndarray
+ Double Layer Capacity
+ '''
+ return (Q_DL[1:] - Q_DL[:-1]) / dx(Phi_pot)
+
+def n(p:np.ndarray, K:str|float) -> np.ndarray:
+ '''
+ Calculates the total number density
+
+ Parameters
+ ----------
+ p : np.ndarray
+ Pressure
+ K : str | float
+ Bulk modulus, use 'incompressible' for an incompressible mixture
+
+ Returns
+ -------
+ np.ndarray
+ Total number density
+ '''
+ if K == 'incompressible':
+ return np.ones_like(p)
+ n_Dimensionless = (p-1) / K + 1
+ return n_Dimensionless
+
+def Q_num_(y_A:np.ndarray, y_C:np.ndarray, n:np.ndarray, x:np.ndarray, z_A:float=-1.0, z_C:float=1.0) -> float:
+ '''
+ Calculates the charge of the system
+
+ Q = ∫_Ω n^F dΩ
+
+ Parameters
+ ----------
+ y_A : np.ndarray
+ Anion fraction
+ y_C : np.ndarray
+ Cation fraction
+ n : np.ndarray
+ Total number density
+ x : np.ndarray
+ Spatial discretization
+ z_A : float, optional
+ Charge number of anions, by default -1.0
+ z_C : float, optional
+ Charge number of cations, by default 1.0
+
+ Returns
+ -------
+ float
+ Charge of the system
+ '''
+ nF_dimensionless = (z_C * y_C + z_A * y_A) * n
+ nF_int = -np.trapz(nF_dimensionless, x)
+ return nF_int
+
+def Q_num_dim(y_A:np.ndarray, y_C:np.ndarray, n:np.ndarray, x:np.ndarray, z_A:float, z_C:float, nR_m:float, e0:float, LR:float) -> float:
+ '''
+ Calculates the charge of the system
+
+ Q = ∫_Ω n^F dΩ
+
+ Parameters
+ ----------
+ y_A : np.ndarray
+ Anion fraction
+ y_C : np.ndarray
+ Cation fraction
+ n : np.ndarray
+ Total number density
+ x : np.ndarray
+ Spatial discretization
+ z_A : float, optional
+ Charge number of anions, by default -1.0
+ z_C : float, optional
+ Charge number of cations, by default 1.0
+ nR_m : float
+ Reference number density in 1/m^3
+ e0 : float
+ Dielectric constant
+ LR : float
+ Reference length in m
+
+ Returns
+ -------
+ float
+ Charge of the system
+ '''
+ Q_DL = Q_num_(y_A, y_C, n, x, z_A, z_C)
+ Q_DL *= nR_m * e0 * LR
+ Q_DL *= 1e+6
+ Q_DL *= 1/(1e+4)
+ return Q_DL
+
+def Q_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, solvation:float) -> float:
+ '''
+ Calculates charge of the system using the analytical method in dimensionless units
+
+ Q = λ^2(∂ₓ φᴸ − ∂ₓ φᴿ)
+
+ Parameters
+ ----------
+ y_A_R : float
+ Value of anion fraction at the right boundary
+ y_C_R : float
+ Value of cation fraction at the right boundary
+ y_N_R : float
+ Value of neutral fraction at the right boundary
+ z_A : float
+ Charge number of anions
+ z_C : float
+ Charge number of cations
+ phi_L : float
+ Electric potential at the left boundary [V]
+ phi_R : float
+ Electric potential at the right boundary
+ p_R : float
+ Pressure at the right boundary
+ K : str | float
+ Bulk modulus, use 'incompressible' for an incompressible mixture
+ Lambda2 : float
+ Dimensionless parameter
+ a2 : float
+ Dimensionless parameter
+ solvation : float
+ Solvation number
+
+ Returns
+ -------
+ float
+ Charge of the system in dimensionless units
+ '''
+ if solvation != 0:
+ raise ValueError('Solvation number must be 0 for the analytical method')
+ z_N = 0
+ E_p = p_R
+ if K == 'incompressible':
+ # Assume E_p = p_R = 0
+ D_A = y_A_R / (np.exp(-a2*p_R-z_A*phi_R))
+ D_C = y_C_R / (np.exp(-a2*p_R-z_C*phi_R))
+ D_N = y_N_R / (np.exp(-a2*p_R))
+
+ # Analytical method
+ E_p = p_R
+ CLambda_L = np.log(D_A * np.exp(-z_A * phi_L - E_p) + D_C * np.exp(-z_C * phi_L - E_p) + D_N * np.exp(-z_N * phi_L - E_p))
+
+ Lambda = np.sqrt(Lambda2)
+ Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2) * (np.sqrt(CLambda_L))
+ return Q_DL
+ else:
+ C_A = y_A_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_A*phi_R))
+ C_C = y_C_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_C*phi_R))
+ C_N = y_N_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_N*phi_R))
+
+ CLambda_L = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N
+
+ Left = np.sqrt((1/CLambda_L)**(-1/(a2*K)) + 1 - K - E_p)
+
+ Lambda = np.sqrt(Lambda2)
+ a = np.sqrt(a2)
+
+ Q_DL = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * a * np.sqrt(2) * (Left)# - Right)
+ return Q_DL
+ raise ValueError('Invalid input for K')
+
+def Q_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, solvation:float) -> float:
+ '''
+ Calculates charge of the system using the analytical method in dimensionless units
+
+ Q = λ^2(∂ₓ φᴸ − ∂ₓ φᴿ)
+
+ Parameters
+ ----------
+ y_A_R : float
+ Value of anion fraction at the right boundary
+ y_C_R : float
+ Value of cation fraction at the right boundary
+ y_N_R : float
+ Value of neutral fraction at the right boundary
+ z_A : float
+ Charge number of anions
+ z_C : float
+ Charge number of cations
+ phi_L : float
+ Value of electric potential at the left boundary (dimensionless units)
+ phi_R : float
+ Value of electric potential at the right
+ p_R : float
+ Value of pressure at the right boundary
+ K : str | float
+ Bulk modulus, use 'incompressible' for an incompress
+ Lambda2 : float
+ Dimensionless parameter
+ a2 : float
+ Dimensionless parameter
+ nR_m : float
+ Reference number density in 1/m^3
+ e0 : float
+ Dielectric constant
+ LR : float
+ Reference length in m
+ solvation : float
+ Solvation number
+
+ Returns
+ -------
+ float
+ Charge of the system in µAs/cm³
+ '''
+ Q_DL = Q_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, phi_L, phi_R, p_R, K, Lambda2, a2, solvation)
+ Q_DL *= nR_m * e0 * LR
+ Q_DL *= 1e+6
+ Q_DL *= 1/(1e+4)
+ return Q_DL
+
+
+
+def C_DL_dimless_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, solvation:float) -> float:
+ '''
+ Calculates the double layer capacity of the system using the analytical method in dimensionless units
+
+ C_dl = ∂Q/∂ φᴸ
+
+ Parameters
+ ----------
+ y_A_R : float
+ Value of anion fraction at the right boundary
+ y_C_R : float
+ Value of cation fraction at the right boundary
+ y_N_R : float
+ Value of neutral fraction at the right boundary
+ z_A : float
+ Charge number of anions
+ z_C : float
+ Charge number of cations
+ phi_L : float
+ Electric potential at the left boundary [V]
+ phi_R : float
+ Electric potential at the right boundary
+ p_R : float
+ Pressure at the right boundary
+ K : str | float
+ Bulk modulus, use 'incompressible' for an incompressible mixture
+ Lambda2 : float
+ Dimensionless parameter
+ a2 : float
+ Dimensionless parameter
+ solvation : float
+ Solvation number
+
+ Returns
+ -------
+ float
+ Double layer capacity of the system in dimensionless units
+ '''
+ if solvation != 0:
+ raise ValueError('Solvation number must be 0 for the analytical method')
+ z_N = 0
+ E_p = p_R
+ if K == 'incompressible':
+ # Assume E_p = p_R = 0
+ D_A = y_A_R / (np.exp(-a2*p_R-z_A*phi_R))
+ D_C = y_C_R / (np.exp(-a2*p_R-z_C*phi_R))
+ D_N = y_N_R / (np.exp(-a2*p_R))
+
+ CLambda_L = np.log(D_A * np.exp(-z_A * phi_L - E_p * a2) + D_C * np.exp(-z_C * phi_L - E_p * a2) + D_N * np.exp(-z_N * phi_L - E_p * a2))
+
+ xi_CLambda_L = D_A * np.exp(-z_A * phi_L - E_p * a2) + D_C * np.exp(-z_C * phi_L - E_p * a2) + D_N * np.exp(-z_N * phi_L - E_p * a2)
+ dxi_CLambda_L = - D_A * z_A * np.exp(-z_A * phi_L - E_p * a2) - D_C * z_C * np.exp(-z_C * phi_L - E_p * a2) - D_N * z_N * np.exp(-z_N * phi_L - E_p * a2)
+
+ Lambda = np.sqrt(Lambda2)
+
+ PreTerm = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * np.sqrt(2)
+ DerivedTerm = dxi_CLambda_L / (2 * xi_CLambda_L * np.sqrt(np.log(xi_CLambda_L)))
+
+ return PreTerm * DerivedTerm
+ else:
+ C_A = y_A_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_A*phi_R))
+ C_C = y_C_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_C*phi_R))
+ C_N = y_N_R / ((K + p_R - 1)**(-1*a2*K)*np.exp(-z_N*phi_R))
+
+ z = 1 - K - E_p
+ y = 1/(a2*K)
+
+ CLambda = C_A * np.exp(-z_A*phi_L) + C_C * np.exp(-z_C*phi_L) + C_N * np.exp(-z_N*phi_L)
+ dClambda = -C_A * z_A * np.exp(-z_A*phi_L) - C_C * z_C * np.exp(-z_C*phi_L) - C_N * z_N * np.exp(-z_N*phi_L)
+
+ Lambda = np.sqrt(Lambda2)
+ a = np.sqrt(a2)
+
+ DerivedTerm = y * CLambda**(y-1) * dClambda / (2 * np.sqrt(CLambda**y + z))
+
+ PreTerm = (phi_L-phi_R) / np.abs(phi_L-phi_R) * Lambda * a * np.sqrt(2)
+
+ return PreTerm * DerivedTerm
+ raise ValueError('Invalid input for K')
+
+
+
+def C_DL_dim_ana(y_A_R:float, y_C_R:float, y_N_R:float, z_A:float, z_C:float, phi_L:float, phi_R:float, p_R:float, K:str|float, Lambda2:float, a2:float, nR_m:float, e0:float, LR:float, k:float, T:float, solvation:float) -> float:
+ '''
+ Calculates the double layer capacity of the system using the analytical method in dimensionless units
+
+ C_dl = ∂Q/∂ φᴸ
+
+ Parameters
+ ----------
+ y_A_R : float
+ Value of anion fraction at the right boundary
+ y_C_R : float
+ Value of cation fraction at the right boundary
+ y_N_R : float
+ Value of neutral fraction at the right boundary
+ z_A : float
+ Charge number of anions
+ z_C : float
+ Charge number of cations
+ phi_L : float
+ Value of electric potential at the left boundary (dimensionless units)
+ phi_R : float
+ Value of electric potential at the right
+ p_R : float
+ Value of pressure at the right boundary
+ K : str | float
+ Bulk modulus, use 'incompressible' for an incompress
+ Lambda2 : float
+ Dimensionless parameter
+ a2 : float
+ Dimensionless parameter
+ nR_m : float
+ Reference number density in 1/m^3
+ e0 : float
+ Dielectric constant
+ LR : float
+ Reference length in m
+ k : float
+ Boltzmann constant
+ T : float
+ Temperature
+ solvation : float
+ Solvation number
+
+ Returns
+ -------
+ float
+ Double layer capacity of the system in µAs/cm²
+ '''
+ C_DL = C_DL_dimless_ana(y_A_R, y_C_R, y_N_R, z_A, z_C, phi_L, phi_R, p_R, K, Lambda2, a2, solvation)
+ C_DL *= nR_m * e0 * LR
+ # ! ToDo
+ C_DL *= 1e+6
+ C_DL *= 1/(1e+4)
+ C_DL *= 1 /(k*T/e0)
+ return C_DL
\ No newline at end of file
diff --git a/build/lib/Bsc-ElectrolyteModels/__init__.py b/build/lib/Bsc-ElectrolyteModels/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/dist/BscElectrolyteModels-1.0-py3-none-any.whl b/dist/BscElectrolyteModels-1.0-py3-none-any.whl
new file mode 100644
index 0000000000000000000000000000000000000000..fd17b699cfa34864d347290b3d90d59d0c1e2287
Binary files /dev/null and b/dist/BscElectrolyteModels-1.0-py3-none-any.whl differ
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000000000000000000000000000000000000..7fd26b970b848120ad4f6ab398e1864ab0111e60
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,3 @@
+[build-system]
+requires = ["setuptools"]
+build-backend = "setuptools.build_meta"
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..d38a8f818a9483d56f32c832672a2ee394f41221
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,47 @@
+[metadata]
+name = BScElectrolytemodels
+version = attr: 1.0
+author = Jan Habscheid, Manuel Torrilhon, Lambert Theisen
+author_email = Jan.Habscheid@rwth-aachen.de, lambert.theisen@rwth-aachen.de, mt@mathcces.rwth-aachen.de
+description = A placeholder for the description
+long_description = file: README.md
+long_description_content_type = text/markdown
+url = https://git.rwth-aachen.de/JanHab/Bsc-Electrolytemodels
+keywords = fenics-project, Battery Simulation, finite element method
+license = GNU General Public License v3.0 or later
+classifiers =
+ License :: OSI Approved :: GNU General Public License v3 (GPLv3)
+ Programming Language :: Python :: 3
+
+[options]
+package_dir=
+ =src
+packages = find:
+zip_safe = True
+python_requires == 3.12.3
+include_package_data = True
+install_requires =
+ pyvista == 0.43.10
+ numpy == 1.26.4
+ scipy == 1.14.0
+
+[options.packages.find]
+where=src
+
+[options.entry_points]
+console_scripts =
+ my-example-utility = example.example_module:main
+
+[options.extras_require]
+notebook = jupyter >=1.0.0
+tests = pytest == 8.3.3
+docs =
+ sphinx == 7.3.7
+ myst-parser == 4.0.0
+ sphinx-copybotton == 0.5.2
+ sphinx_rtd_theme == 3.0.1
+build = build >= 0.7.0
+dev =
+ BScElectrolytemodels[tests]
+ BScElectrolytemodels[docs]
+ BScElectrolytemodels[build]
\ No newline at end of file
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000000000000000000000000000000000000..5d16151a1ae45c78235d6ce914fd6145a2b293e2
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,22 @@
+from setuptools import setup, find_packages
+
+setup(
+ name="BscElectrolyteModels",
+ version="1.0",
+ description="A description is yet to follow",
+ url="https://git.rwth-aachen.de/JanHab/Bsc-ElectrolyteModels",
+ author="Jan Habscheid, Lambert Theisen, Manuel Torrilhon",
+ author_email="Jan.Habscheid@rwth-aachen.de, lambert.theisen@rwth-aachen.de, mt@mathcces.rwth-aachen.de",
+ packages=["Bsc-ElectrolyteModels"],
+ package_dir={"Bsc-ElectrolyteModels":"src"},
+ package_data={"Bsc-ElectrolyteModels": ["tests/*"]},
+ python_requires="==3.12.3",
+ install_requires=[
+ "pyvista == 0.43.10",
+ "numpy == 1.26.4",
+ "scipy == 1.14.0"
+ ],
+ setup_requires=["pytest-runner"],
+ tests_require=["pytest>=7.4.0"],
+ test_suite="tests"
+)
diff --git a/tests/__pycache__/__init__.cpython-312.pyc b/tests/__pycache__/__init__.cpython-312.pyc
index 089662da0b56caba69f0501d02dac24686b4e462..651390ebc6f1c4a594785796f33eaf39a8f1de6a 100644
Binary files a/tests/__pycache__/__init__.cpython-312.pyc and b/tests/__pycache__/__init__.cpython-312.pyc differ
diff --git a/tests/__pycache__/test_ElectrolyticDiode.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_ElectrolyticDiode.cpython-312-pytest-8.3.3.pyc
index aebed3727349670ee4178aee4a98d1060dad63ca..75443412202ab6efa782e9dc013bf0d7a93ba147 100644
Binary files a/tests/__pycache__/test_ElectrolyticDiode.cpython-312-pytest-8.3.3.pyc and b/tests/__pycache__/test_ElectrolyticDiode.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/__pycache__/test_Eq02.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_Eq02.cpython-312-pytest-8.3.3.pyc
index df145f680fec694c6c644d8d8fdd81d57546a1f8..d6ac38f5546ba664255fcc939cd984fec4472e78 100644
Binary files a/tests/__pycache__/test_Eq02.cpython-312-pytest-8.3.3.pyc and b/tests/__pycache__/test_Eq02.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/__pycache__/test_Eq04.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_Eq04.cpython-312-pytest-8.3.3.pyc
index 389b1d06261f49dbfb4ffe9d76c7efeda9f2ad47..bfc8a082ed7595e3d98f6f3f974a4be8c0472035 100644
Binary files a/tests/__pycache__/test_Eq04.cpython-312-pytest-8.3.3.pyc and b/tests/__pycache__/test_Eq04.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/__pycache__/test_EqN.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_EqN.cpython-312-pytest-8.3.3.pyc
index 6316e86eb0c6a2b14ab1657d500b84b1b8dfdab3..48cd69033b4fac8207b44760ff99d7ef78e8235a 100644
Binary files a/tests/__pycache__/test_EqN.cpython-312-pytest-8.3.3.pyc and b/tests/__pycache__/test_EqN.cpython-312-pytest-8.3.3.pyc differ
diff --git a/tests/__pycache__/test_Helper_DoubleLayerCapacity.cpython-312-pytest-8.3.3.pyc b/tests/__pycache__/test_Helper_DoubleLayerCapacity.cpython-312-pytest-8.3.3.pyc
index 3f329ce080859a5eda5164959c5001f4000d54fa..109a7a873b2b538ea6e694f0006e8c569c1bab68 100644
Binary files a/tests/__pycache__/test_Helper_DoubleLayerCapacity.cpython-312-pytest-8.3.3.pyc and b/tests/__pycache__/test_Helper_DoubleLayerCapacity.cpython-312-pytest-8.3.3.pyc differ