Question about usage of model.generate()

[ottarg]

Hi all,

I am just getting familiar with PyBaMM. What I am trying to do is get a C or C++ version of a cell model, for use on an embedded system.

I see that PyBaMM models have methods for calling the CasADi code generation functionality.
https://docs.pybamm.org/en/stable/_modules/pybamm/models/base_model.html#BaseModel.generate
https://web.casadi.org/docs/#generating-c-code

and it seems as though the intent of this CasADi functionality is indeed to support embedded systems, but it isn't immediately clear to me if I am understand this feature correctly. Is the purpose of
model.generate(filename, variable_names, input_parameter_order=None, cg_options=None)
to be able to obtain C/C++ code that can run on a step-by-step basis? i.e. essentially to run in the same fashion as what is described in the BaseSolver.step() method? (https://docs.pybamm.org/en/stable/_modules/pybamm/solvers/base_solver.html#BaseSolver.step)
where I would then call this generated function each timestep with inputs being the previous system state and any further inputs like current and timestep?

If so, is there any example showing this workflow?

Thanks!

[valentinsulzer]

Once you have generated the code, you have to compile it yourself and then you can call functions from the generated file. See

PyBaMM/tests/unit/test_models/test_base_model.py

Lines 813 to 834 in 44d946f

	model.generate("test.c", ["a+b"], input_parameter_order=["p", "q"])
	# Compile
	subprocess.run(["gcc", "-fPIC", "-shared", "-o", "test.so", "test.c"]) # nosec
	# Read the generated functions
	x0_fn = casadi.external("x0", "./test.so")
	z0_fn = casadi.external("z0", "./test.so")
	rhs_fn = casadi.external("rhs_", "./test.so")
	alg_fn = casadi.external("alg_", "./test.so")
	jac_rhs_fn = casadi.external("jac_rhs", "./test.so")
	jac_alg_fn = casadi.external("jac_alg", "./test.so")
	var_fn = casadi.external("variables", "./test.so")
	# Test that function values are as expected
	self.assertEqual(x0_fn([2, 5]), 5)
	self.assertEqual(z0_fn([0, 0]), 1)
	self.assertEqual(rhs_fn(0, 3, 2, [7, 2]), -21)
	self.assertEqual(alg_fn(0, 3, 2, [7, 2]), 1)
	np.testing.assert_array_equal(np.array(jac_rhs_fn(5, 6, 7, [8, 9])), [[-8, 0]])
	np.testing.assert_array_equal(np.array(jac_alg_fn(5, 6, 7, [8, 9])), [[1, -1]])
	self.assertEqual(var_fn(6, 3, 2, [7, 2]), -1)

You are then responsible for writing your own solver and post-processing. There are no workflows for doing this in PyBaMM.

[valentinsulzer]

The inputs required are:

• t, a scalar
• x, the values of the ODE states, which can be a scalar or a vector. In the example I linked, it's a scalar (a)
• z, the values of the algebraic states, which can be a scalar or a vector. In the example I linked, it's a scalar (b)
• the inputs parameters, which can be a scalar or a vector. In the example I linked, it's a vector ([p, q])

For example, in a model with 3 ode states, 2 algebraic states, and 4 input parameters, this would be a valid function call:
rhs_fun(0, [1,2,3], [4,5], [6,7,8,9])

[ottarg]

Ah okay thank you! That makes sense

[ottarg]

Hello again @valentinsulzer, sorry to bother you further about this, but I have reached a somewhat confusing blocker and I am not sure whether I have done something wrong, or whether I am misunderstanding something.

So far, I have succeeded in generating code for the built in model
model = pybamm.lithium_ion.DFN()
I obtain the code, compile it, and I can run the model in matlab using ode15s and everything seems to work quite well. However, when I try to run the model with a solver of my own, I have a strange issue. (important to note: I have very little experience with such things and may be making a foolish mistake)

I am trying to put together a simple Newton-Raphson implicit solver, which looks something like this:

for i = 2:numTimeSteps
    yPrev = y(i-1, :);
    % Newton-Raphson iteration to find yCurrent
    yGuess = yPrev; % Initial guess
    for j = 1:maxIterations% Could also use a convergence-based stopping criterion
        % Evaluate the combined system residual at the current guess
        F = daeResiduals(t(i), yGuess, yPrev, dt, p);
        
        % Calculate Jacobian and solve the linear system for updates
        J = jacobianFunction(t(i), yGuess, p);
        delta = -J \ F;
        yGuess = yGuess + delta;

        % Check convergence
        if norm(delta, 'inf') < 1e-6
            break;
        end
    end
    y(i, :) = yGuess;
end

Where the method daeResiduals(t(i), yGuess, yPrev, dt, p) uses the generated code for

codeGenDFN("rhs_",t,x,z,p) and codeGenDFN("alg_",t,x,z,p) to construct the residuals, and the jacobianFunction(t(i), yGuess, p)
calculates the Jacobian using codeGenDFN("jac_rhs",t,x,z,p) and codeGenDFN("jac_alg",t,x,z,p).

My problem:

I am construction the Jacobian J as:

Jx = codeGenDFN("jac_rhs",t,x,z,p);
Jz = codeGenDFN("jac_alg",t,x,z,p);
J = [Jx;Jz];

The built-in DFN model has 860 ODE states, and 100 algebraic states and so the dimensions are: Jx = [860x960] and Jz = [100x960]

However, when I then execute the step delta = -J \ F;, I have the problem that my J is rank-deficient, and so the calculation goes crazy. The ranks of the matrices are:

rank(full(Jx) = 859
rank(full(Jz) = 100
rank(full(J)) = 957

Now, what I am unsure about is: Is this expected? I am not familiar enough with the DFN model to know the answer. What really confuses me though is that this seems to be very clearly the reason why my solution is not working, and yet when I use ode15s with the same Jacobian function, it works perfectly fine. Is there something obvious I am missing here, like an influence of boundary conditions etc which resolve this?

[valentinsulzer]

I'm not totally surprised that the DFN is rank deficient. ode15s probably does some smart approximations or reformulations to work even with rank deficient matrices. I can't think of an easy solution unfortunately. You'll have to either figure out how to make DFN not be rank deficient, or write a DAE solver that works with rank deficient matrices. Maybe others in the pybamm community will have ideas, I'm not an expert in numerical algorithms

[ottarg]

Okay! Thanks for the help, I'm at least glad to know that there is not an obvious thing I am missing.