Hi Hugo,
On 27/04/2022 09:29, hugo.levy@gmail.com wrote:
Hi Robert,
I am currently facing a new difficulty using Sfepy for solving nonlinear PDEs. I have coded my own Newton’s iterative method that simply consists in linearizing the PDE F(u)=0 around the previous iteration u_k. I have two questions:
- At each iteration, I create a
Probleminstance for the linearized PDE and solve it the usual way in Sfepy. However, only the “linearized-nonlinear” terms need to be recomputed from one iteration to another. For the sake of computational efficiency, I wish not to reconstruct all matrices / rhs vectors but only the ones that need to be updated. How would you achieve that? In other words, how can you store some assembled terms while re-assembling some others?
There is no built-in support for this, but because you implement your own solver, it should be possible to get without too much hassle. I suggest you to split your equations to linear and nonlinear parts and evaluate those separately. There are several ways how to do that, e.g. (pb is a Problem instance):
eqs = { 'K' : '...', }
... pb.time_update() pb.update_materials()
mtx = pb.mtx_a.copy() mtx_k = eqs['K'].evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx)
etc. See also assemble_matrices() in [1].
Another way would be to have two problem instance for the linear and nonlinear parts and evaluate those.
[1] https://sfepy.org/doc-devel/examples/linear_elasticity-dispersion_analysis.h...
- For some reason, I need to evaluate another function, say G (a residual), at each iteration in the weak form. Again, some terms in G are constant while others are iteration-dependent (through their material). How would you efficiently evaluate G(u_k) at each iteration? Currently, the only way I found is to do
dummy_ts = TimeStepper(0, 100) # ficticious time-stepper eqs = Equations([my_equation]) # equation associated with function G eqs.init_state() ebcs = Conditions(my_ebcs) epbcs = Conditions(my_epbcs) eqs.time_update(dummy_ts, ebcs=ebcs, epbcs=epbcs) res_vec = eqs.eval_residuals(sol) # sol is the solution at current iteration
which I guess this is not Sfepy’s philosophy…
It looks OK. Note that you probably do not need to use the dummy_ts. Maybe you could use Problem.evaluate(), or Problem.create_evaluable(), .eval_equations() combo?
I hope my questions are clear enough. Thank you very much in advance for your help.
As usual, send a self-contained small-size example if the above does not help or you get stuck.
Cheers, r.