Hi Robert,

Thank you for your swift reply, very useful as always!

I have dug into the code (`Problem.solve`, `Newton` class from solvers\nls.py) in order
to see what's happening with the residual vector and the tangent matrix. Especially, since
I only solve linear systems at each iteration of my own Newton algorithm, I could by-pass the
whole Sfepy's Newton call and keep only the lines

vec_dx = lin_solver(vec_r, x0=vec_x, eps_a=eps_a, eps_r=eps_r, mtx=mtx_a, status=ls_status)

where I would have pre-computed vec_r and mtx_a the way you proposed.

There are a few catches, though (at least for me):

- I am using LCBCs. As you pointed out in [1], LCBSs are not compatible with setting
'active_only' to False at the creation of the Problem instance. This means the matrices &
vectors I assemble "by hand" will be reduced in size (because of the constrained DOFs), and so
will be my solution vector. I will then recover the full solution with

variables.set_state(vec_dx, pb.active_only)

- LCBCs modify the residual vector as well as the tangent matrix (class `Evaluator` methods
`eval_residual` and `eval_tangent_matrix`), which would not be accounted for if I do

mtx = pb.mtx_a.copy()
mtx_k = eqs['K'].evaluate(mode='weak', dw_mode='matrix', asm_obj=mtx)

So I believe creating two Problem instances (each with the same BCs) as you suggested is better
as I wouldn't have to worry about it.

- When pb.is_linear() is True, there is a presolve step which computes a LU matrix factorization
from what I understand. Does this mean I have to call ls.presolve(mtx) before calling the lin_solver()
function? Just like in [2] for instance.

So far, the code could look something like this

...

pb_constant = Problem('pb_constant', equations=eqs_constant, active_only=True)
pb_constant.set_bcs(ebcs=my_ebcs, epbcs=my_epbcs, lcbcs=my_lcbcs)

# Compute matrix & vector of linear terms
pb_constant.set_solver(nls)
tss = pb_constant.get_solver()
pb_constant.equations.set_data(None, ignore_unknown=True)
variables = pb_constant.get_initial_state(vec=None)
pb_constant.time_update(tss.ts)
variables.apply_ebc(force_values=None)
pb_constant.update_materials()
ev_constant = pb_constant.get_evaluator()
mtx_constant = ev.eval_tangent_matrix(variables(), is_full=True)
rhs_constant = ev.eval_residual(variables(), is_full=True)

pb_update = Problem('pb_update', equations=eqs_update, active_only=True)
pb_update.set_bcs(ebcs=my_ebcs, epbcs=my_epbcs, lcbcs=my_lcbcs)

while 1: # my Newton iteration loop

# Compute matrix & vector of nonlinear terms
pb_update.set_solver(nls)
tss = pb_constant.get_solver()
pb_update.equations.set_data(None, ignore_unknown=True)
variables = pb_update.get_initial_state(vec=None)
pb_update.time_update(tss.ts)
variables.apply_ebc(force_values=None)
pb_update.update_materials()
ev_constant = pb_update.get_evaluator()
mtx_update = ev.eval_tangent_matrix(variables(), is_full=True)
rhs_update = ev.eval_residual(variables(), is_full=True)

# Construct the full matrix / rhs
mtx = mtx_update + mtx_constant
rhs = rhs_update + rhs_constant

# Solve
ls = nls.lin_solver
ls.presolve(mtx)
vec_x = ls(rhs, x0=None, eps_a=eps_a, eps_r=eps_r, mtx=mtx, status={})
sol = variables.set_state(vec_x, pb_update.active_only) # This might not be the correct syntax with the last release!

if stop_condition:
break

...
...

I do not quite get what the `is_full=True` option does exactly. Should it be set to True when
active_only is True?

I would be very glad if you could give me a quick feedback on the code / remarks. If you think this
looks fine (i.e. the computational cost is actually reduced by proceeding in this manner), I will implement
this solution for both of my initial questions.

Thanks again,

Hugo

[1] https://github.com/sfepy/sfepy/issues/785
[2] https://sfepy.org/doc-devel/examples/diffusion-time_poisson_interactive.html

Le jeu. 28 avr. 2022 à 11:08, Robert Cimrman <cimrman3@ntc.zcu.cz> a écrit :
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:
>
> 1) At each iteration, I create a `Problem` instance 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.html

> 2) 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.
_______________________________________________
SfePy mailing list -- sfepy@python.org
To unsubscribe send an email to sfepy-leave@python.org
https://mail.python.org/mailman3/lists/sfepy.python.org/
Member address: hugo.levy@gmail.com