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

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

- 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

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)

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

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)

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

...

...# 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