Issue 765 compute grad fem

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added the gradient operation for the Finite Element Method ([#767](https://github.com/pybamm-team/PyBaMM/pull/767))
 -   Added `InputParameter` node for quickly changing parameter values ([#752](https://github.com/pybamm-team/PyBaMM/pull/752))
 -   Added submodels for operating modes other than current-controlled ([#751](https://github.com/pybamm-team/PyBaMM/pull/751))
 -   Changed finite volume discretisation to use exact values provided by Neumann boundary conditions when computing the gradient instead of adding ghost nodes([#748](https://github.com/pybamm-team/PyBaMM/pull/748))
@@ -35,6 +36,7 @@
 
 ## Bug fixes
 
+-   Fixed a bug which meant that the Ohmic heating in the current collectors was incorrect if using the Finite Element Method ([#767](https://github.com/pybamm-team/PyBaMM/pull/767))
 -   Improved automatic broadcasting ([#747](https://github.com/pybamm-team/PyBaMM/pull/747))
 -   Fixed bug with wrong temperature in initial conditions ([#737](https://github.com/pybamm-team/PyBaMM/pull/737))
 -   Improved flexibility of parameter values so that parameters (such as diffusivity or current) can be set as functions or scalars ([#723](https://github.com/pybamm-team/PyBaMM/pull/723))

examples/notebooks/change-input-current.ipynb

@@ -66,12 +66,12 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "118979747d4f418982a358d3769dc257",
+       "model_id": "000da8e81b0b4e91984553ba15ba179c",
        "version_major": 2,
        "version_minor": 0
       },
       "text/plain": [
-       "interactive(children=(FloatSlider(value=0.0, description='t', max=0.00787878787878788, step=0.005), Output()),…"
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=0.0036363636363636364, step=0.005), Output()…"
       ]
      },
      "metadata": {},
@@ -134,7 +134,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "397b778427504eb8ac88b53de8397683",
+       "model_id": "cc4fd2e7c9ea4129a7acf0e519fea268",
        "version_major": 2,
        "version_minor": 0
       },
@@ -177,12 +177,12 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "4acb95c01295406db377f5dd1593e014",
+       "model_id": "84f6cac9b14f43a1907322e650f3c9f8",
        "version_major": 2,
        "version_minor": 0
       },
       "text/plain": [
-       "interactive(children=(FloatSlider(value=0.0, description='t', max=0.057314594406781015, step=0.001), Output())…"
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=0.026550306076661374, step=0.001), Output())…"
       ]
      },
      "metadata": {},
@@ -301,12 +301,12 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "19720a31f3964e55b01905476f50393a",
+       "model_id": "09ced8a9e10b4e769a9b4cc8573b7fa2",
        "version_major": 2,
        "version_minor": 0
       },
       "text/plain": [
-       "interactive(children=(FloatSlider(value=0.0, description='t', max=0.0028657297203390506, step=0.00014328648601…"
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=0.0013275153038330688, step=6.63757651916534…"
       ]
      },
      "metadata": {},
@@ -324,7 +324,7 @@
     "# Example: simulate for 30 seconds\n",
     "simulation_time = 30  # end time in seconds\n",
     "tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge).evaluate(0)\n",
-    "npts = 50 * simulation_time * omega  # need enough timesteps to resolve output\n",
+    "npts = int(50 * simulation_time * omega)  # need enough timesteps to resolve output\n",
     "t_eval = np.linspace(0, simulation_time / tau, npts)\n",
     "solution = model.default_solver.solve(model, t_eval)\n",
     "label = [\"Frequency: {} Hz\".format(omega)]\n",
@@ -361,7 +361,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.4"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,

examples/scripts/compare_comsol/compare_comsol_DFN.py

@@ -63,7 +63,7 @@
 def get_interp_fun(variable_name, domain):
     """
     Create a :class:`pybamm.Function` object using the variable, to allow plotting with
-    :class:`'pybamm.QuickPlot'` (interpolate in space to match edges, and then create
+    :class:`pybamm.QuickPlot` (interpolate in space to match edges, and then create
     function to interpolate in time)
     """
     variable = comsol_variables[variable_name]
@@ -78,7 +78,17 @@ def get_interp_fun(variable_name, domain):
     variable = interp.interp1d(comsol_x, variable, axis=0)(pybamm_x)
 
     def myinterp(t):
-        return interp.interp1d(comsol_t, variable)(t)[:, np.newaxis]
+        try:
+            return interp.interp1d(
+                comsol_t, variable, fill_value="extrapolate", bounds_error=False,
+            )(t)[:, np.newaxis]
+        except ValueError as err:
+            raise ValueError(
+                """Failed to interpolate '{}' with time range [{}, {}] at time {}.
+                Original error: {}""".format(
+                    variable_name, comsol_t[0], comsol_t[-1], t, err
+                )
+            )
 
     # Make sure to use dimensional time
     fun = pybamm.Function(myinterp, pybamm.t * tau, name=variable_name + "_comsol")
@@ -92,7 +102,9 @@ def myinterp(t):
 comsol_phi_n = get_interp_fun("phi_n", ["negative electrode"])
 comsol_phi_e = get_interp_fun("phi_e", whole_cell)
 comsol_phi_p = get_interp_fun("phi_p", ["positive electrode"])
-comsol_voltage = interp.interp1d(comsol_t, comsol_variables["voltage"])
+comsol_voltage = interp.interp1d(
+    comsol_t, comsol_variables["voltage"], fill_value="extrapolate", bounds_error=False
+)
 
 # Create comsol model with dictionary of Matrix variables
 comsol_model = pybamm.BaseModel()

pybamm/expression_tree/unary_operators.py

@@ -327,7 +327,6 @@ class Gradient_Squared(SpatialOperator):
     operator with itself. In particular, this is useful in the finite element
     formualtion where we only require the (sclar valued) square of the gradient,
     and  not the gradient itself.
-
     **Extends:** :class:`SpatialOperator`
     """
 

pybamm/models/submodels/thermal/x_lumped/x_lumped_2D_current_collectors.py

@@ -49,8 +49,7 @@ def _current_collector_heating(self, variables):
         """Returns the heat source terms in the 2D current collector"""
         phi_s_cn = variables["Negative current collector potential"]
         phi_s_cp = variables["Positive current collector potential"]
-        # Note: grad not implemented in 2D weak form, but can compute grad squared
-        # directly
+
         Q_s_cn = self.param.sigma_cn_prime * pybamm.grad_squared(phi_s_cn)
         Q_s_cp = self.param.sigma_cp_prime * pybamm.grad_squared(phi_s_cp)
         return Q_s_cn, Q_s_cp

pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_2D_current_collector.py

@@ -25,8 +25,7 @@ def _current_collector_heating(self, variables):
         """Returns the heat source terms in the 2D current collector"""
         phi_s_cn = variables["Negative current collector potential"]
         phi_s_cp = variables["Positive current collector potential"]
-        # Note: grad not implemented in 2D weak form, but can compute grad squared
-        # directly
+
         Q_s_cn = self.param.sigma_cn_prime * pybamm.grad_squared(phi_s_cn)
         Q_s_cp = self.param.sigma_cp_prime * pybamm.grad_squared(phi_s_cp)
         return Q_s_cn, Q_s_cp

pybamm/spatial_methods/scikit_finite_element.py

@@ -3,7 +3,8 @@
 #
 import pybamm
 
-from scipy.sparse import csr_matrix
+from scipy.sparse import csr_matrix, csc_matrix
+from scipy.sparse.linalg import inv
 import numpy as np
 import skfem
 
@@ -67,10 +68,104 @@ def spatial_variable(self, symbol):
         return vector
 
     def gradient(self, symbol, discretised_symbol, boundary_conditions):
-        """Matrix-vector multiplication to implement the gradient operator.
-        See :meth:`pybamm.SpatialMethod.gradient`
+        """Matrix-vector multiplication to implement the gradient operator. The
+        gradient w of the function u is approximated by the finite element method
+        using the same function space as u, i.e. we solve w = grad(u), which
+        corresponds to the weak form w*v*dx = grad(u)*v*dx, where v is a suitable
+        test function.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol that we will take the laplacian of.
+        discretised_symbol: :class:`pybamm.Symbol`
+            The discretised symbol of the correct size
+        boundary_conditions : dict
+            The boundary conditions of the model
+            ({symbol.id: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
+
+        Returns
+        -------
+        :class: `pybamm.Concatenation`
+            A concatenation that contains the result of acting the discretised
+            gradient on the child discretised_symbol. The first column corresponds
+            to the y-component of the gradient and the second column corresponds
+            to the z component of the gradient.
         """
-        raise NotImplementedError
+        domain = symbol.domain[0]
+        mesh = self.mesh[domain][0]
+
+        # get gradient matrix
+        grad_y_matrix, grad_z_matrix = self.gradient_matrix(symbol, boundary_conditions)
+
+        # assemble mass matrix (there is no need to zero out entries here, since
+        # boundary conditions are already accounted for in the governing pde
+        # for the symbol we are taking the gradient of. we just want to get the
+        # correct weights)
+        @skfem.bilinear_form
+        def mass_form(u, du, v, dv, w):
+            return u * v
+
+        mass = skfem.asm(mass_form, mesh.basis)
+        # we need the inverse
+        mass_inv = pybamm.Matrix(inv(csc_matrix(mass)))
+
+        # compute gradient
+        grad_y = mass_inv @ (grad_y_matrix @ discretised_symbol)
+        grad_z = mass_inv @ (grad_z_matrix @ discretised_symbol)
+
+        # create concatenation
+        grad = pybamm.Concatenation(
+            grad_y, grad_z, check_domain=False, concat_fun=np.hstack
+        )
+        grad.domain = domain
+
+        return grad
+
+    def gradient_squared(self, symbol, discretised_symbol, boundary_conditions):
+        """Multiplication to implement the inner product of the gradient operator
+        with itself. See :meth:`pybamm.SpatialMethod.gradient_squared`
+        """
+        grad = self.gradient(symbol, discretised_symbol, boundary_conditions)
+        grad_y, grad_z = grad.orphans
+        return grad_y ** 2 + grad_z ** 2
+
+    def gradient_matrix(self, symbol, boundary_conditions):
+        """
+        Gradient matrix for finite elements in the appropriate domain.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol for which we want to calculate the gradient matrix
+        boundary_conditions : dict
+            The boundary conditions of the model
+            ({symbol.id: {"negative tab": neg. tab bc, "positive tab": pos. tab bc}})
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) finite element gradient matrix for the domain
+        """
+        # get primary domain mesh
+        domain = symbol.domain[0]
+        mesh = self.mesh[domain][0]
+
+        # make form for the gradient in the y direction
+        @skfem.bilinear_form
+        def gradient_dy(u, du, v, dv, w):
+            return du[0] * v[0]
+
+        # make form for the gradient in the z direction
+        @skfem.bilinear_form
+        def gradient_dz(u, du, v, dv, w):
+            return du[1] * v[1]
+
+        # assemble the matrices
+        grad_y = skfem.asm(gradient_dy, mesh.basis)
+        grad_z = skfem.asm(gradient_dz, mesh.basis)
+
+        return pybamm.Matrix(grad_y), pybamm.Matrix(grad_z)
 
     def divergence(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the divergence operator.
@@ -151,15 +246,6 @@ def unit_bc_load_form(v, dv, w):
 
         return -stiffness_matrix @ discretised_symbol + boundary_load
 
-    def gradient_squared(self, symbol, discretised_symbol, boundary_conditions):
-        """Matrix-vector multiplication to implement the inner product of the
-        gradient operator with itself.
-        See :meth:`pybamm.SpatialMethod.gradient_squared`
-        """
-        stiffness_matrix = self.stiffness_matrix(symbol, boundary_conditions)
-
-        return stiffness_matrix @ (discretised_symbol ** 2)
-
     def stiffness_matrix(self, symbol, boundary_conditions):
         """
         Laplacian (stiffness) matrix for finite elements in the appropriate domain.

pybamm/spatial_methods/spatial_method.py

@@ -185,7 +185,6 @@ def gradient_squared(self, symbol, discretised_symbol, boundary_conditions):
             The symbol that we will take the gradient of.
         discretised_symbol: :class:`pybamm.Symbol`
             The discretised symbol of the correct size
-
         boundary_conditions : dict
             The boundary conditions of the model
             ({symbol.id: {"left": left bc, "right": right bc}})

tests/unit/test_spatial_methods/test_scikit_finite_element.py

@@ -13,8 +13,6 @@ def test_not_implemented(self):
         spatial_method = pybamm.ScikitFiniteElement()
         spatial_method.build(mesh)
         self.assertEqual(spatial_method.mesh, mesh)
-        with self.assertRaises(NotImplementedError):
-            spatial_method.gradient(None, None, None)
         with self.assertRaises(NotImplementedError):
             spatial_method.divergence(None, None, None)
         with self.assertRaises(NotImplementedError):
@@ -54,7 +52,6 @@ def test_discretise_equations(self):
             pybamm.laplacian(var) - pybamm.source(unit_source, var, boundary=True),
             pybamm.laplacian(var)
             - pybamm.source(unit_source ** 2 + 1 / var, var, boundary=True),
-            pybamm.grad_squared(var),
         ]:
             # Check that equation can be evaluated in each case
             # Dirichlet
@@ -125,6 +122,40 @@ def test_discretise_equations(self):
         with self.assertRaises(pybamm.GeometryError):
             disc.process_symbol(x)
 
+    def test_gradient(self):
+        mesh = get_unit_2p1D_mesh_for_testing(ypts=32, zpts=32)
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume(),
+            "current collector": pybamm.ScikitFiniteElement(),
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # test gradient of 5*y + 6*z
+        var = pybamm.Variable("var", domain="current collector")
+        disc.set_variable_slices([var])
+
+        y = mesh["current collector"][0].coordinates[0, :]
+        z = mesh["current collector"][0].coordinates[1, :]
+
+        gradient = pybamm.grad(var)
+        grad_disc = disc.process_symbol(gradient)
+        grad_disc_y, grad_disc_z = grad_disc.children
+
+        np.testing.assert_array_almost_equal(
+            grad_disc_y.evaluate(None, 5 * y + 6 * z),
+            5 * np.ones_like(y)[:, np.newaxis],
+        )
+        np.testing.assert_array_almost_equal(
+            grad_disc_z.evaluate(None, 5 * y + 6 * z),
+            6 * np.ones_like(z)[:, np.newaxis],
+        )
+
+        # check grad_squared positive
+        eqn = pybamm.grad_squared(var)
+        eqn_disc = disc.process_symbol(eqn)
+        ans = eqn_disc.evaluate(None, 3 * y ** 2)
+        np.testing.assert_array_less(0, ans)
+
     def test_manufactured_solution(self):
         mesh = get_unit_2p1D_mesh_for_testing(ypts=32, zpts=32)
         spatial_methods = {