FormSystemMatrix for NonlinearForm-types ?

[j-signorelli]

Hello,

I am not certain if there is any similar / preferred capability already, but if not: I think it would be useful to have an ability to assemble a Matrix with eliminated essential BCs from a given NonlinearForm, provided a solution vector $u$.

Just to explain my motivation for this: I am in the process of basically setting up Example 16p but w/ Newton iterations,

$$\mathbf{M}\dfrac{d\vec{u}}{dt}=-\mathbf{K}\vec{u} $$

where there is a nonlinear mass matrix $M_{ij}= \displaystyle \int_{\Omega}\rho C(u_h)\phi_i\phi_jd\vec{x}$ and stiffness matrix $K_{ij} = \displaystyle \int_{\Omega}k(u_h)\nabla\phi_{il}\nabla\phi_{lj}d\vec{x}$. I have written NonlinearFormIntegrator's for them that implement AssembleElementVector and AssembleElementGrad for Newton iterations.

For ConductionOperator::Mult, I can just call ParNonlinearForm::Mult for $\mathbf{K}$ to get the action of it on $\vec{u}$ to obtain a Vector for the RHS at the previous timestep, but unfortunately I have no easy way to get the LHS as a HypreParMatrix from the ParNonlinearForm for $\mathbf{M}$ at the previous timestep only.

I think this update would require moving AssembleElementMatrix from BilinearFormIntegrator up to its parent NonlinearFormIntegrator and then generalizing it to allow input for a solution vector, like:

NonlinearFormIntegrator::AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, DenseMatrix &elmat)

and then having a similar header for FormSystemMatrix that takes in a solution vector.

This just seemed more straightforward/cleaner to me than having an entirely separate ParBilinearForm that I have to update + reassemble at every timestep just for accomplishing this, but absolutely please feel free to let me know if there is any other/better approaches! I am still relatively new to MFEM/FEA.

[vladotomov]

@j-signorelli can't you get it from NonlinearForm::GetGradient(u), something like this?

[j-signorelli]

Hi @vladotomov

Thank you for your reply. Wouldn't that just return the Jacobian though? Unless I am missing something / not fully understanding. I'm looking for a way to assemble the operator in a matrix form for a given $u$.

[vladotomov]

Aha I see. No, we don't have such option. I'm not sure that moving AssembleElementMatrix into NonlinearFormIntegrator makes sense. What is this matrix for a general nonlinear operator F(u)?
Even if the function is in NonlinearFormIntegrator, won't you still have to reassemble the matrix at every step?

[j-signorelli]

@vladotomov I was thinking it would just be an evaluation of the nonlinear operator for a given $u$, but thinking more about it you are correct that it is generally not applicable, as for $\displaystyle\int_\Omega e^{\alpha u}vd\vec{x}$ for example. The nonlinear forms I am dealing with are just bilinear form types but with coefficients that make them not , which is why I had thought that. Apologies for the confusion and thank you for helping clear that up.

It would have just been for convenience/cleanliness in the code -- I would very easily be able to write these particular NonlinearFormIntegrator's to have an AssembleElementMatrix capability and then avoid having to have additional ParBilinearForm's allocating and assembling matrices just for ConductionOperator::Mult on top of the two ParNonlinearForm's I have for ConductionOperator::ImplicitSolve.

[tangqi]

You can checkout ex10p. I think you can assemble matrices when evaluating F(u). Handling bc there is also easy. I found it is more flexible than NonlinearForm.

[j-signorelli]

Hi @tangqi,

Thanks for your reply, yes ex10p has definitely been very helpful in understanding how to deal with a nonlinear term in a time-dependent problem with MFEM. For my ConductionOperator::Mult, I am doing similar with K(u) as is done for H(x) by just using ParNonlinearForm::Mult:

mfem/examples/ex10p.cpp

Lines 589 to 608 in 148047f

	void HyperelasticOperator::Mult(const Vector &vx, Vector &dvx_dt) const
	{
	// Create views to the sub-vectors v, x of vx, and dv_dt, dx_dt of dvx_dt
	int sc = height/2;
	Vector v(vx.GetData() + 0, sc);
	Vector x(vx.GetData() + sc, sc);
	Vector dv_dt(dvx_dt.GetData() + 0, sc);
	Vector dx_dt(dvx_dt.GetData() + sc, sc);
	H.Mult(x, z);
	if (viscosity != 0.0)
	{
	S.TrueAddMult(v, z);
	z.SetSubVector(ess_tdof_list, 0.0);
	}
	z.Neg(); // z = -z
	M_solver.Mult(z, dv_dt);
	dx_dt = v;
	}

However this presumes a constant mass matrix, while the one I have is nonlinear so it must be updated with the $u_n$ solution every timestep.

It seems that the best route is to just maintain two forms for M(u): one ParBilinearForm that I just update and reassemble every timestep for ConductionOperator::Mult and one ParNonlinearForm for ConductionOperator::ImplicitSolve.

[tangqi]

Yes, I would assemble M(u) inside this HyperelasticOperator::Mult when evaluating F(u) on the fly. I think you can define a reduced system operator inside HyperelasticOperator where it owns and assemble a separate mass matrix for Jacobian. This is likely the better way, because depending which nonlinear solver is called, you may need to keep a copy of both mass matrices.

[j-signorelli]

Thanks again. Although I'm starting to question if the above form of the problem even makes sense. I had previously wanted to say it should essentially be

$$\mathbf{M}(\vec{u})\dfrac{d\vec{u}}{dt}=-\mathbf{K}(\vec{u})\vec{u}$$

But this obviously doesn't make sense since it is wrong to think about a matrix for a NonlinearForm, as I had just recognized from @vladotomov. I wonder if it would be better to just divide the two coefficients $\rho$ and $C(u)$ to the other side before taken the inner product with a test function? But I also have a Neumann term (which I excluded from the discussion just for simplicity) which this would then complicate by making it nonlinear specifically because I want to be able to specify it as a heat flux. Then maybe it would look something like

$$\mathbf{M}\dfrac{d\vec{u}}{dt}=-\kappa(\vec{u}) + N(\vec{u}, t)$$

Just thinking out loud, I'll have to think more about it, thanks again very much for all your help.

[tangqi]

Both (making the solution dependent coefficient appearing on either side) can be done properly. However, when it appears on the left, your system becomes F(t, u, udot) = 0 because you have to differentiate M(u) too. You just need to be careful with your reduced system to solve for udot.