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
Evaluatormethodseval_residualandeval_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:
- 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.
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