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