Trying to understand dFEM qfunction derivative action

[lindsayad]

I'm trying to understand exactly where the action of the derivative happens. Does it happen within this call?

mfem/fem/dfem/doperator.hpp

Lines 810 to 812 in cd4e583

	call_qfunction_derivative_action<qf_param_ts>(
	qfunc, input_shmem, shadow_shmem, residual_shmem,
	da_size_on_qp, num_qp, q1d, dimension, use_sum_factorization);

If so, where exactly? I can't see the action happening within the apply* callees. I see the qfunction being called and with automatic differentiation this will return dual numbers (I'm not building with enzyme) and the Jacobian information will be contained within those duals. I see extraction of J (or at least I believe that is what is happening) in

mfem/fem/dfem/qfunction_transform.hpp

Lines 172 to 182 in cd4e583

	template <typename T, int n>
	MFEM_HOST_DEVICE inline
	void process_derivative_from_native_dual(
	DeviceTensor<1, T> &r,
	const tensor<dual<T, T>, n> &x)
	{
	for (size_t i = 0; i < n; i++)
	{
	r(i) = x(i).gradient;
	}
	}

.
But where is the multiplication of J by the input vector occurring?

[jandrej]

Does it happen within this call?

Yes

I see extraction of J (or at least I believe that is what is happening)

That's the extraction part, yes.

But where is the multiplication of J by the input vector occurring

If you're using the dual number approach, the derivative part of the dual (gradient part) is interpolated here and assigned here, by being called on the properly overloaded types here.

For example a matrix with it's dual entries is assigned with the direction here.

This is a rather long chain of invocations but if you want to plug in your favorite AD library, you can include a header that overloads these functions with the appropriate types and now dFEM is using your favorite backend.

Now dual::gradient holds the "direction" or "input" vector. Each consecutive operation on dual::gradient is multiplying the direction with the gradient chains.

Dataflow is similar for Enzyme backed operations.

[lindsayad]

Sorry, I'm still confused. In the dfem-minimal-surface example, dual_t is just dual<real_t, real_t> so the qfunction as I understand it will only support a single derivative since the type of gradient is just real_t. So when you talk about direction/input are you talking about specifying what partial derivative to take? E.g. specify $j$ for $\frac{\partial f}{\partial u_j}$ ? What I'm still not understanding is where the multiplication by $du$ is happening for the action $df = \frac{\partial f}{\partial u} du$

[jandrej]

The action (multiplication by du) is in each operation of the dual types gradient variable. This is the standard way to do it.

[lindsayad]

Ok I believe I've wrapped my head around it. I am so used to the CSR matrix world where we explicitly compute J. So we use our AD package (MetaPhysicL) and always seed our DualNumber derivative with $du = 1$, and then $J du$ is just $J$ which we historically need. But yes in the matrix-free world where what we care about is $J du$ for arbitrary $du$, the below makes sense

Now dual::gradient holds the "direction" or "input" vector. Each consecutive operation on dual::gradient is multiplying the direction with the gradient chains.

The initial seed corresponds to the input vector $du$, not just 1 which is what I'm used to. Thank you!

[jandrej]

Exactly as you described. If we compute J explicitly (e.g. for matrix assemble or Partial Assembly caching) then we do the same as you mentioned above, seed with 1 and write out the result.