Repository: pybamm-team/PyBaMM (python)

Issue:
Allow porosity to vary as function of x
**Summary**
 Specify porosity as a function of x in models

 **Motivation**
Model graded electrodes, etc. If step function, need to be careful with how we discretise (easiest thing is to make sure to have an edge at the step change)



Pull request description:
Issue 802 porosity
# Description

Allow porosity to vary in space

Fixes #802 

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/master/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [x] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [ ] Bug fix (non-breaking change which fixes an issue)


# Key checklist:

- [x] No style issues: `$ flake8`
- [x] All tests pass: `$ python run-tests.py --unit`
- [x] The documentation builds: `$ cd docs` and then `$ make clean; make html`

You can run all three at once, using `$ python run-tests.py --quick`.

## Further checks:

- [x] Code is commented, particularly in hard-to-understand areas
- [x] Tests added that prove fix is effective or that feature works


Patch:
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 8140fdeed..d6146c0fb 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -173,17 +173,6 @@ from .models.submodels import (
     tortuosity,
 )
 
-#
-# Parameters class and methods
-#
-from .parameters.parameter_values import ParameterValues
-from .parameters import geometric_parameters
-from .parameters import electrical_parameters
-from .parameters import thermal_parameters
-from .parameters import standard_parameters_lithium_ion, standard_parameters_lead_acid
-from .parameters.print_parameters import print_parameters, print_evaluated_parameters
-from .parameters import parameter_sets
-
 #
 # Geometry
 #
@@ -202,6 +191,18 @@ from .geometry.geometry import (
 from .expression_tree.independent_variable import KNOWN_SPATIAL_VARS, KNOWN_COORD_SYS
 from .geometry import standard_spatial_vars
 
+#
+# Parameters class and methods
+#
+from .parameters.parameter_values import ParameterValues
+from .parameters import geometric_parameters
+from .parameters import electrical_parameters
+from .parameters import thermal_parameters
+from .parameters import standard_parameters_lithium_ion, standard_parameters_lead_acid
+from .parameters.print_parameters import print_parameters, print_evaluated_parameters
+from .parameters import parameter_sets
+
+
 #
 # Mesh and Discretisation classes
 #
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index 6bac3d5d6..e1225bf03 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -38,6 +38,8 @@ class Concatenation(pybamm.Symbol):
         # combine domains from children
         domain = []
         for child in children:
+            if not isinstance(child, pybamm.Symbol):
+                raise TypeError("{} is not a pybamm symbol".format(child))
             child_domain = child.domain
             if set(domain).isdisjoint(child_domain):
                 domain += child_domain
diff --git a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
index f395877b0..4c3afaf42 100644
--- a/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/basic_dfn.py
@@ -100,9 +100,15 @@ class BasicDFN(BaseModel):
         # Porosity
         # Primary broadcasts are used to broadcast scalar quantities across a domain
         # into a vector of the right shape, for multiplying with other vectors
-        eps_n = pybamm.PrimaryBroadcast(param.epsilon_n, "negative electrode")
-        eps_s = pybamm.PrimaryBroadcast(param.epsilon_s, "separator")
-        eps_p = pybamm.PrimaryBroadcast(param.epsilon_p, "positive electrode")
+        eps_n = pybamm.PrimaryBroadcast(
+            pybamm.Parameter("Negative electrode porosity"), "negative electrode"
+        )
+        eps_s = pybamm.PrimaryBroadcast(
+            pybamm.Parameter("Separator porosity"), "separator"
+        )
+        eps_p = pybamm.PrimaryBroadcast(
+            pybamm.Parameter("Positive electrode porosity"), "positive electrode"
+        )
         eps = pybamm.Concatenation(eps_n, eps_s, eps_p)
 
         # Tortuosity
diff --git a/pybamm/models/submodels/porosity/constant_porosity.py b/pybamm/models/submodels/porosity/constant_porosity.py
index c487957f2..da07259ac 100644
--- a/pybamm/models/submodels/porosity/constant_porosity.py
+++ b/pybamm/models/submodels/porosity/constant_porosity.py
@@ -20,17 +20,10 @@ class Constant(BaseModel):
 
     def get_fundamental_variables(self):
 
-        eps_n_av = self.param.epsilon_n
-        eps_s_av = self.param.epsilon_s
-        eps_p_av = self.param.epsilon_p
+        eps_n = self.param.epsilon_n
+        eps_s = self.param.epsilon_s
+        eps_p = self.param.epsilon_p
 
-        eps_n = pybamm.FullBroadcast(
-            eps_n_av, "negative electrode", "current collector"
-        )
-        eps_s = pybamm.FullBroadcast(eps_s_av, "separator", "current collector")
-        eps_p = pybamm.FullBroadcast(
-            eps_p_av, "positive electrode", "current collector"
-        )
         eps = pybamm.Concatenation(eps_n, eps_s, eps_p)
 
         deps_n_dt = pybamm.FullBroadcast(0, "negative electrode", "current collector")
diff --git a/pybamm/parameters/standard_parameters_lithium_ion.py b/pybamm/parameters/standard_parameters_lithium_ion.py
index 7b0624287..88731b2b8 100644
--- a/pybamm/parameters/standard_parameters_lithium_ion.py
+++ b/pybamm/parameters/standard_parameters_lithium_ion.py
@@ -274,14 +274,11 @@ centre_y_tab_p = pybamm.geometric_parameters.centre_y_tab_p
 centre_z_tab_p = pybamm.geometric_parameters.centre_z_tab_p
 
 # Microscale geometry
-epsilon_n = pybamm.Parameter("Negative electrode porosity")
-epsilon_s = pybamm.Parameter("Separator porosity")
-epsilon_p = pybamm.Parameter("Positive electrode porosity")
-epsilon = pybamm.Concatenation(
-    pybamm.FullBroadcast(epsilon_n, ["negative electrode"], "current collector"),
-    pybamm.FullBroadcast(epsilon_s, ["separator"], "current collector"),
-    pybamm.FullBroadcast(epsilon_p, ["positive electrode"], "current collector"),
-)
+var = pybamm.standard_spatial_vars
+epsilon_n = pybamm.FunctionParameter("Negative electrode porosity", var.x_n)
+epsilon_s = pybamm.FunctionParameter("Separator porosity", var.x_s)
+epsilon_p = pybamm.FunctionParameter("Positive electrode porosity", var.x_p)
+epsilon = pybamm.Concatenation(epsilon_n, epsilon_s, epsilon_p)
 epsilon_s_n = pybamm.Parameter("Negative electrode active material volume fraction")
 epsilon_s_p = pybamm.Parameter("Positive electrode active material volume fraction")
 epsilon_inactive_n = 1 - epsilon_n - epsilon_s_n


Test patch:
diff --git a/tests/unit/test_expression_tree/test_concatenations.py b/tests/unit/test_expression_tree/test_concatenations.py
index 2a5a16114..eb18381dd 100644
--- a/tests/unit/test_expression_tree/test_concatenations.py
+++ b/tests/unit/test_expression_tree/test_concatenations.py
@@ -25,6 +25,10 @@ class TestConcatenations(unittest.TestCase):
         with self.assertRaises(TypeError):
             conc2.evaluate()
 
+        # trying to concatenate non-pybamm symbols
+        with self.assertRaises(TypeError):
+            pybamm.Concatenation(1, 2)
+
     def test_concatenation_domains(self):
         a = pybamm.Symbol("a", domain=["negative electrode"])
         b = pybamm.Symbol("b", domain=["separator", "positive electrode"])