Merge pull request #213 from jotelha/joss-review-issue-185

CLI and FEniCS integration test cases, addressing Joss review issue 185

matscipy/cli/electrochemistry/poisson_nernst_planck_solver_cli.py

@@ -39,7 +39,6 @@ def main():
 
     Specify quantities in SI units at this command line interface."""
 
-
     # in order to have both:
     # * preformatted help text and ...
     # * automatic display of defaults
@@ -222,6 +221,12 @@ class ArgumentDefaultsAndRawDescriptionHelpFormatter(
         raise ValueError("Boundary conditions '{}' not implemented!".format(
                          args.boundary_conditions))
 
+    if not args.outfile:
+        # fenics writes some log messages to stdout. If we pipe the output,
+        # we don't want that, hence here we have to suppress fenics logging
+        # if no output file has been specified
+        dolfin.cpp.log.set_log_active(False)
+
     pnp.solve()
 
     extra_kwargs = {}

tests/test_c2d.py → tests/test_electrochemistry_cli.py

@@ -31,8 +31,9 @@
 
 from looseversion import LooseVersion
 
-class c2dCliTest(matscipytest.MatSciPyTestCase):
-    """Tests c2d and pnp command line interfaces"""
+
+class ElectrochemistryCliTest(matscipytest.MatSciPyTestCase):
+    """Tests continuous2discrete and poisson-nernst-planck command line interfaces"""
 
     def setUp(self):
         """Reads reference data files"""
@@ -60,9 +61,9 @@ def setUp(self):
         # command line interface scripts are expected to reside within
         # ../matscipy/cli/electrochemistry relative to this test directory
         self.pnp_cli = os.path.join(
-            self.test_path,os.path.pardir,'matscipy','cli','electrochemistry','pnp.py')
+            self.test_path,os.path.pardir,'matscipy','cli','electrochemistry','poisson_nernst_planck_solver_cli.py')
         self.c2d_cli = os.path.join(
-            self.test_path,os.path.pardir,'matscipy','cli','electrochemistry','c2d.py')
+            self.test_path,os.path.pardir,'matscipy','cli','electrochemistry','continuous2discrete_cli.py')
 
         self.assertTrue(os.path.exists(self.pnp_cli))
         self.assertTrue(os.path.exists(self.c2d_cli))
@@ -86,7 +87,7 @@ def tearDown(self):
         self.bin_path.cleanup()
 
     def test_c2d_input_format_npz_output_format_xyz(self):
-        """c2d NaCl.npz NaCl.xyz"""
+        """continuous2discrete NaCl.npz NaCl.xyz"""
         print("  RUN test_c2d_input_format_npz_output_format_xyz")
         with tempfile.TemporaryDirectory() as tmpdir:
             subprocess.run(
@@ -101,7 +102,7 @@ def test_c2d_input_format_npz_output_format_xyz(self):
             self.assertTrue( ( xyz.cell == self.ref_xyz.cell ).all() )
 
     def test_c2d_input_format_txt_output_format_xyz(self):
-        """c2d NaCl.txt NaCl.xyz"""
+        """continuous2discrete NaCl.txt NaCl.xyz"""
         print("  RUN test_c2d_input_format_txt_output_format_xyz")
         with tempfile.TemporaryDirectory() as tmpdir:
             subprocess.run(
@@ -119,7 +120,7 @@ def test_c2d_input_format_txt_output_format_xyz(self):
         """ LAMMPS data file won't work for ASE version up until 3.18.1,
             LAMMPS data file input broken in ASE 3.19.0, skipped""")
     def test_c2d_input_format_npz_output_format_lammps(self):
-        """c2d NaCl.npz NaCl.lammps"""
+        """continuous2discrete NaCl.npz NaCl.lammps"""
         print("  RUN test_c2d_input_format_npz_output_format_lammps")
         with tempfile.TemporaryDirectory() as tmpdir:
             subprocess.run(
@@ -147,13 +148,14 @@ def test_pnp_output_format_npz(self):
 
             test_npz = np.load(os.path.join(tmpdir,'out.npz'))
 
-            # depending on which solver is used for the tests, absolute values might deviate more than the default tolerance
+            # depending on which solver is used for the tests,
+            # absolute values might deviate more than the default tolerance
             self.assertArrayAlmostEqual(test_npz['x'], self.ref_npz['x'], tol=1e-7)
             self.assertArrayAlmostEqual(test_npz['u'], self.ref_npz['u'], tol=1e-5)
             self.assertArrayAlmostEqual(test_npz['c'], self.ref_npz['c'], tol=1e-3)
 
     def test_pnp_output_format_txt(self):
-        """pnp -c 0.1 0.1 -u 0.05 -l 1.0e-7 -bc cell NaCl.txt"""
+        """poisson-nernst-planck -c 0.1 0.1 -u 0.05 -l 1.0e-7 -bc cell NaCl.txt"""
         print("  RUN test_pnp_output_format_txt")
         with tempfile.TemporaryDirectory() as tmpdir:
             subprocess.run(
@@ -164,13 +166,14 @@ def test_pnp_output_format_txt(self):
 
             test_txt = np.loadtxt(os.path.join(tmpdir,'out.txt'), unpack=True)
 
-            # depending on which solver is used for the tests, absolute values might deviate more than the default tolerance
+            # depending on which solver is used for the tests,
+            # absolute values might deviate more than the default tolerance
             self.assertArrayAlmostEqual(test_txt[0,:], self.ref_txt[0,:], tol=1e-5) # x
             self.assertArrayAlmostEqual(test_txt[1,:],self.ref_txt[1,:], tol=1e-5) # u
             self.assertArrayAlmostEqual(test_txt[2:,:],self.ref_txt[2:,:], tol=1e-3) # c
 
     def test_pnp_c2d_pipeline_mode(self):
-        """pnp -c 0.1 0.1 -u 0.05 -l 1.0e-7 -bc cell | c2d > NaCl.xyz"""
+        """poisson-nernst-planck -c 0.1 0.1 -u 0.05 -l 1.0e-7 -bc cell | continuous2discrete > NaCl.xyz"""
         print("  RUN test_pnp_c2d_pipeline_mode")
         with tempfile.TemporaryDirectory() as tmpdir:
             pnp = subprocess.Popen(

tests/test_poisson_nernst_planck_solver.py

@@ -24,6 +24,7 @@
 
 from matscipy.electrochemistry import PoissonNernstPlanckSystem
 
+
 class PoissonNernstPlanckSolverTest(matscipytest.MatSciPyTestCase):
 
     def setUp(self):
@@ -45,5 +46,6 @@ def test_poisson_nernst_planck_solver_std_interface_bc(self):
         self.assertArrayAlmostEqual(pnp.potential, self.ref_data ['u'])
         self.assertArrayAlmostEqual(pnp.concentration, self.ref_data ['c'], 1e-6)
 
+
 if __name__ == '__main__':
     unittest.main()

tests/test_poisson_nernst_planck_solver_fenics.py

@@ -0,0 +1,63 @@
+#
+# Copyright 2019-2020 Johannes Hoermann (U. Freiburg)
+#
+# matscipy - Materials science with Python at the atomic-scale
+# https://github.com/libAtoms/matscipy
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+import matscipytest
+import numpy as np
+import os.path
+import sys
+import unittest
+
+
+try:
+    import fenics
+except ImportError:
+    print("fenics not found: skipping fenics-dependent tests")
+
+from matscipy.electrochemistry.poisson_nernst_planck_solver_fenics \
+    import PoissonNernstPlanckSystemFEniCS as PoissonNernstPlanckSystem
+
+
+class PoissonNernstPlanckSolverTest(matscipytest.MatSciPyTestCase):
+
+    def setUp(self):
+        """Provides 0.1 mM NaCl solution at 0.05 V across 100 nm open half space reference data from binary npz file"""
+        self.test_path = os.path.dirname(os.path.abspath(__file__))
+        self.ref_data = np.load(
+            os.path.join(self.test_path, 'electrochemistry_data',
+            'NaCl_c_0.1_mM_0.1_mM_z_+1_-1_L_1e-7_u_0.05_V_seg_200_interface.npz') )
+
+    @unittest.skipIf("fenics" not in sys.modules,
+                     "fenics required")
+    def test_poisson_nernst_planck_solver_fenics_std_interface_bc(self):
+        """Tests PNP solver against simple interfacial BC"""
+        pnp = PoissonNernstPlanckSystem(
+            c=[0.1,0.1], z=[1,-1], L=1e-7, delta_u=0.05,
+            N=200, e=1e-12, maxit=20)
+        pnp.useStandardInterfaceBC()
+        pnp.solve()
+
+        # Reference data has been generated with controlled-volume solver and is slightly off the FEM results,
+        # hence the generous tolerances below
+        self.assertArrayAlmostEqual(pnp.grid, self.ref_data ['x'])
+        self.assertArrayAlmostEqual(pnp.potential, self.ref_data ['u'], 1e-6)
+        self.assertArrayAlmostEqual(pnp.concentration, self.ref_data ['c'], 1e-5)
+
+
+if __name__ == '__main__':
+    unittest.main()