First version for a python package

.eggs/README.txt

+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.
+

.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/LICENSE

+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.

.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/PKG-INFO

+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.

.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/RECORD

+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,,

.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/WHEEL

+Wheel-Version: 1.0
+Generator: bdist_wheel (0.42.0)
+Root-Is-Purelib: true
+Tag: py3-none-any
+

.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/entry_points.txt

+[distutils.commands]
+ptr = ptr:PyTest
+pytest = ptr:PyTest

.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/requires.txt

+
+[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

.eggs/pytest_runner-6.0.1-py3.12.egg/EGG-INFO/top_level.txt

+ptr

.eggs/pytest_runner-6.0.1-py3.12.egg/ptr/__init__.py

+"""
+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)

BscElectrolyteModels.egg-info/PKG-INFO

+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)
+
+[![Pipeline Status](https://git.rwth-aachen.de/Jan.Habscheid/bsc-electrolytemodels/badges/main/pipeline.svg)](https://git.rwth-aachen.de/Jan.Habscheid/bsc-electrolytemodels/pipelines)
+[![Documentation](https://img.shields.io/badge/docs-latest-blue)](https://janhab.pages.rwth-aachen.de/bsc-electrolytemodels/)
+[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.13645296.svg)](https://doi.org/10.5281/zenodo.13645296)
+[![GitLab Version](https://img.shields.io/badge/version-1.0-blue.svg)](https://git.rwth-aachen.de/jan.habscheid/bsc-electrolytemodels/-/tags)
+[![License](https://img.shields.io/badge/license-GPLv3-blue)](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

BscElectrolyteModels.egg-info/SOURCES.txt

+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

BscElectrolyteModels.egg-info/dependency_links.txt

+

BscElectrolyteModels.egg-info/entry_points.txt

+[console_scripts]
+my-example-utility = example.example_module:main

BscElectrolyteModels.egg-info/requires.txt

+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

BscElectrolyteModels.egg-info/top_level.txt

+Bsc-ElectrolyteModels

BscElectrolyteModels.egg-info/zip-safe

+

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

build/lib/Bsc-ElectrolyteModels/ElectrolyticDiode.py

+'''
+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

build/lib/Bsc-ElectrolyteModels/Eq02.py

+'''
+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()

build/lib/Bsc-ElectrolyteModels/Eq04.py

+'''
+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()