Hi Robert, Thank you for taking the time to review these bits of code.
In general, I am not too happy about LCBCs, as they work only for the nodal (Lagrange) basis. Remind me what kind of constraints do you use? Enforcing those by some penalty or Lagrange multipliers might be a more future-proof way to go.
Basically I have two domains with a different PDE defined on each and I want the solution to match at a frontier shared by the two domains (a bit like [1]). To do so, I use something like
def zero_shift(ts, coors, region):
return np.zeros_like(coors[:, 0], dtype=np.float64)
match_coors = Function('match_coors', per.match_coors)
shift_fun = Function('shift_fun', zero_shift)
lcbcs = Conditions([LinearCombinationBC('lc12', [gamma1, gamma2],
{'u2.all' : 'u1.all'}, match_coors, 'shifted_periodic',
arguments=(shift_fun,))])
I tend to find it quite convenient as Sfepy automatically builds the full stiffness matrix, with DOFs from both domains. Could this be easily achieved with Lagrange multipliers?
Yes, but it is not your case - you solve a different linear system in each iteration, right?
Yes indeed, so no need to pre-compute the factorization beforehand.
I have an additional question. Say I have assembled both the matrix and rhs
(in the reduced form,
just as detailed in the previous messages) and I want to modify the
essential boundary conditions
of my problem. According to [2], the stiffness matrix stays the same while
the rhs vector is modified.
Is there a way to update the rhs vector without having to make a full
ev.eval_residual() call?
Either way, you have answered my initial questions and my implementation seems to be working.
Cheers!
Hugo
[1] https://github.com/sfepy/sfepy/issues/689 [2] https://sfepy.org/doc-devel/ebcs_implementation.html?
Le lun. 2 mai 2022 à 20:01, Robert Cimrman <cimrman3@ntc.zcu.cz> a écrit :
Hi Hugo,
On 29/04/2022 15:57, Hugo Lévy wrote:
Hi Robert,
Thank you for your swift reply, very useful as always!
I have dug into the code (
Problem.solve,Newtonclass 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 linesvec_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.
Yes.
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)
In general, I am not too happy about LCBCs, as they work only for the nodal (Lagrange) basis. Remind me what kind of constraints do you use? Enforcing those by some penalty or Lagrange multipliers might be a more future-proof way to go.
- LCBCs modify the residual vector as well as the tangent matrix (class
Evaluatormethodseval_residualandeval_tangent_matrix), which would not be accounted for if I domtx = 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.
Yes, but it is not your case - you solve a different linear system in each iteration, right?
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=Trueoption does exactly. Should it be set to True when active_only is True?The
is_fulljust says whether the passed-in vector is full or reduced - in the latter case the evaluator function calls make_full_vec(). So yes, it interplays with active_only - when active_only is True, is_full in the Newton solver is False, as it works with the reduced linear system.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.
It looks OK at a glance.
r.
Thanks again,
Hugo
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