Skip to content

Automate source terms in FEM #1410

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Open
rtimms opened this issue Mar 2, 2021 · 6 comments · May be fixed by #4971
Open

Automate source terms in FEM #1410

rtimms opened this issue Mar 2, 2021 · 6 comments · May be fixed by #4971
Labels
difficulty: hard Will take several weeks priority: low No existing plans to resolve

Comments

@rtimms
Copy link
Contributor

rtimms commented Mar 2, 2021

When you solve an ODE (pointwise) using the FEM the variable on the lhs is multiplied by the mass matrix upon discretisation. If you wrap the rhs in pybamm.source then it also gets multiplied by the mass matrix so everything is fine. However, it would be easy to miss this. We should either NOT multiply the lhs by the mass matrix OR automate adding pybamm.source where it is needed.

See this example

import pybamm


class Model(pybamm.lithium_ion.BaseModel):
    def __init__(self, name="2+1D test"):
        super().__init__({"dimensionality": 2}, name)

        i = pybamm.PrimaryBroadcast(1, "current collector")
        a = pybamm.Variable("a", domain="current collector")
        b = pybamm.Variable(
            "b",
            domain="negative electrode",
            auxiliary_domains={"secondary": "current collector"},
        )
        c = pybamm.Variable("c", domain="current collector")

        self.rhs[a] = i
        self.rhs[b] = pybamm.PrimaryBroadcast(i, "negative electrode")
        self.rhs[c] = pybamm.source(i, c)

        self.initial_conditions[a] = pybamm.Scalar(0)
        self.initial_conditions[b] = pybamm.Scalar(0)
        self.initial_conditions[c] = pybamm.Scalar(0)

        self.variables = {
            "a": a,
            "b": b,
            "c": c,
            "b_av": pybamm.x_average(b),
            "true": i * pybamm.t,
        }

    def new_empty_copy(self):
        return pybamm.BaseModel.new_empty_copy(self)


model = Model()
var = pybamm.standard_spatial_vars
var_pts = {
    var.x_n: 3,
    var.x_s: 3,
    var.x_p: 3,
    var.r_n: 3,
    var.r_p: 3,
    var.y: 3,
    var.z: 3,
}
sim = pybamm.Simulation(model, var_pts=var_pts)
sim.solve([0, 3600])
sim.plot(["a", "b_av", "c", "true"])

a gives the wrong result, but b and c are correct (because b has domain negative electrode it is discretised by FVM and because c uses source)
@rtimms
Copy link
Contributor Author

rtimms commented Mar 2, 2021

This was the cause of the problem in #1399 where we had an ODE for the concentration at each point (y,z)

@valentinsulzer
Copy link
Member

Is the problem that the mass matrix always assumes everything is a laplacian unless you wrap it in source?

@rtimms
Copy link
Contributor Author

rtimms commented Mar 2, 2021 •
edited
Loading

It’s because you expect to end up with discretised equations of the form M*dc/dt = K*c + M*s, where the stiffness matrix K comes from the laplacian. We automatically put the M on the lhs but need to use source to get the M on the rhs.

It’s not a huge deal to do this manually at the moment, but it is confusing. If/when we add more spatial methods we’ll have to think a bit more carefully about issues like this.

@jaskiratsingh2000
Copy link

Hi! Seems interestihng. Should I take this up?

@valentinsulzer valentinsulzer added priority: low No existing plans to resolve difficulty: hard Will take several weeks labels May 8, 2023
@ArivoliR
Copy link

ArivoliR commented Apr 11, 2025 •
edited
Loading

Hii, I'd love to take my attempt at solving this issue.

I've defined a wrap_source() function in the class BaseModel, that runs a check if isinstance(rhs_term, pybamm.expression_tree.binary_operators.MatrixMultiplication). After making these changes I ran this

(v) ➜  ~ py test.py
<class 'pybamm.expression_tree.broadcasts.PrimaryBroadcast'> type of a
<class 'pybamm.expression_tree.broadcasts.PrimaryBroadcast'> type of b
<class 'pybamm.expression_tree.binary_operators.MatrixMultiplication'> type of c
2025-04-11 11:25:10.179 - [WARNING] base_model.default_solver(420): The default solver changed to IDAKLUSolver after the v25.4.0. release. You can swap back to the previous default by using `pybamm.CasadiSolver()` instead.
/home/blazevfx/PyBaMM/src/pybamm/models/base_model.py:162: ModelWarning: RHS term for variable 'a' was not wrapped in `pybamm.source()`. Auto-wrapping applied.
  wrapped_rhs[variable] = self.wrap_source(expression, variable)
Traceback (most recent call last):
  File "/home/blazevfx/test.py", line 52, in <module>
    sim.solve([0, 3600])
  File "/home/blazevfx/PyBaMM/src/pybamm/simulation.py", line 472, in solve
    self.build(initial_soc=initial_soc, inputs=inputs)
  File "/home/blazevfx/PyBaMM/src/pybamm/simulation.py", line 327, in build
    self._set_parameters()
  File "/home/blazevfx/PyBaMM/src/pybamm/simulation.py", line 255, in _set_parameters
    self._model_with_set_params = self._parameter_values.process_model(
                                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/blazevfx/PyBaMM/src/pybamm/parameters/parameter_values.py", line 486, in process_model
    model.rhs = new_rhs
    ^^^^^^^^^
  File "/home/blazevfx/PyBaMM/src/pybamm/models/base_model.py", line 162, in rhs
    wrapped_rhs[variable] = self.wrap_source(expression, variable)
                            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/blazevfx/PyBaMM/src/pybamm/models/base_model.py", line 1505, in wrap_source
    rhs_term = pybamm.source(rhs_term, variable)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/blazevfx/PyBaMM/src/pybamm/expression_tree/binary_operators.py", line 1549, in source
    raise pybamm.DomainError(
pybamm.expression_tree.exceptions.DomainError: 'source' only implemented in the 'current collector' domain, but symbols have domains ['negative electrode'] and ['negative electrode']

I ran the example you had sent, but I'm not able to test with it. Would love some guidance on how to proceed from here

/home/blazevfx/PyBaMM/src/pybamm/models/base_model.py:162: ModelWarning: RHS term for variable 'a' was not wrapped in `pybamm.source()`. Auto-wrapping applied.
  wrapped_rhs[variable] = self.wrap_source(expression, variable)

wrap function seems to work though :)

ps: still working on it, have made some progress

@ArivoliR ArivoliR linked a pull request Apr 11, 2025 that will close this issue
Draft
@ArivoliR
Copy link

ArivoliR commented Apr 12, 2025 •
edited
Loading

I have fixed the function. Am getting correct output now and all except one of the integration tests are passing

2025-04-12 13:18:36.795 - [WARNING] base_model.wrap_source(1504): [wrap_source] Wrapping RHS for 'a' in pybamm.source()

Can someone please review this?

the only integration test that failed was

__________________________________ TestLogger.test_logger __________________________________
[gw6] linux -- Python 3.12.2 /home/blazevfx/PyBaMM/.nox/tests/bin/python

self = <tests.unit.test_logger.TestLogger object at 0x7ad1b8336450>

    def test_logger(self):
        logger = pybamm.logger
>       assert logger.level == 30
E       assert 25 == 30
E        +  where 25 = <Logger pybamm.logger (NOTICE)>.level

tests/unit/test_logger.py:11: AssertionError

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
difficulty: hard Will take several weeks priority: low No existing plans to resolve
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Automate source wrap ArivoliR/PyBaMM
4 participants
@jaskiratsingh2000 @valentinsulzer @rtimms @ArivoliR