computeQpOffDiagJacobian in the MOOSE kernel

[Roma Gurung]

Hi,

1. What is the mathematical interpretation of computeQpOffDiagJacobian? In my perspective, i realize that it is the Jacobian w.r.t the coupled variable i,e, differential with respect to the second variable. For example , in phase field model with primary variable (c), the coupling to temperature (T) variable under thermal gradient highlights the need of computeQpOffDiagJacobian in SoretDiffusion.C code. So, is it partial_Residual/partial_T in SoretDiffusion kernel?

2. Recently, i am studying thermotransport in A-B alloys in accordance to http://www.sciencedirect.com/science/article/pii/S0022311511002960. The equation (1) described in the paper composes of two major parts in accordance to MOOSE phase field model. The first part consisting of the Cahn-Hilliard terms can be modeled by the  CoupledTimeDerivative, SplitCHParsed and SplitCHWRes kernels of the iron-chromium spinodal decomposition tutorial (https://github.com/idaholab/moose/blob/devel/modules/phase_field/tutorials/spinodal_decomposition/s4_mobility.i) . The term for thermal gradient

i.e. J_th = M_q * grad_T/ T where, M_q = c(1-c) (beta_1*Q_1 - beta_2 *Q*2) has been tried to be coded by modifying the SoretDiffusion kernel as in: https://github.com/anilkunwar/danphe/blob/master/src/kernel/MultiSoretDiffusion.C . In this part, I am pondering over the computations of the Jacobians. The code snippet from the above links are pasted below:

	MultiSoretDiffusion::computeQpResidual() {
	Real T_term = _Mq[_qp]*(- _c[_qp]*_c[_qp] + _c[_qp]) / ( _T[_qp] );
	return T_term * _grad_T[_qp] * _grad_test[_i][_qp];
	}
	Real
	MultiSoretDiffusion::computeQpJacobian()
	{
	if (_c_var == _var.number()) //Requires c jacobian
	return computeQpCJacobian();
	return 0.0;
	}
	Real
	MultiSoretDiffusion::computeQpOffDiagJacobian(unsigned int jvar)
	{
	if (_c_var == jvar) //Requires c jacobian
	return computeQpCJacobian();
	else if (_T_var == jvar) //Requires T jacobian
	//return _D[_qp] * _Q[_qp] * _c[_qp] * _grad_test[_i][_qp] *
	//(_grad_phi[_j][_qp]/(_kb * _T[_qp] * _T[_qp]) - 2.0 * _grad_T[_qp] * _phi[_j][_qp] / (_kb * _T[_qp] * _T[_qp] * _T[_qp]));
	return _Mq[_qp] * _c[_qp]* (1.0 - _c[_qp])* _grad_test[_i][_qp] *
	(_grad_phi[_j][_qp]/(_T[_qp]) - _grad_T[_qp] * _phi[_j][_qp]/(_T[_qp] * _T[_qp]));
	return 0.0;
	}
	Real
	MultiSoretDiffusion::computeQpCJacobian()
	{
	//Calculate the Jacobian for the c variable
	//return _D[_qp] * _Q[_qp] * _phi[_j][_qp] * _grad_T[_qp] / (_kb * _T[_qp] * _T[_qp]) * _grad_test[_i][_qp];
	return _Mq[_qp]*_grad_test[_i][_qp]*(_grad_T[_qp]/_T[_qp]- 2.0 *_phi[_j][_qp]*_grad_T[_qp]/_T[_qp]);
	}

Are the Jacobians written correctly? Notably, when  i don't write the _grad_test[_i][_qp] *  in the  computeQpOffDiagJacobian part, the following error arises:

/home/username/moose_projects/danphe/src/kernel/MultiSoretDiffusion.C:58:91: error: cannot convert 'libMesh::boostcopy::enable_if_c<true, libMesh::TypeVector<double> >::type {aka libMesh::TypeVector<double>}' to 'libMesh::Real {aka double}' in return
           (_grad_phi[_j][_qp]/(_T[_qp]) - _grad_T[_qp] * _phi[_j][_qp]/(_T[_qp] * _T[_qp]));

Can anyone guide me in the theoretical formulation of offdiagonal jacobian?

Yours Sincerely

Anil Kunwar

[SudiptaBiswas]

Hi Anil,

You can find the Jacobian Calculation method here, 

http://mooseframework.org/wiki/MooseTraining/FEM/NumericalImplementation/

For Off diagonal Jacobian you have to calculate the derivative of the residual w.r.t any coupled variable other than the one the kernel is acting on. 

Thanks!

Sudipta

[Roma Gurung]

Hi Sudipta,

Thank you very much for the link. Moreover, your concise reply on off diagonal jacobian has helped me understand the numerical coding more vividly.

Can you share your perspective on "what does _phi[_j][_qp] replace in the off diagonal jacobian? Does it replace the coupled variable w.r.t. which the residual is differentiated?"?

Also, do you find the off diagonal jacobian expression numerically consistent [meaning can the solver solver c(1-c) expression]?

	return _Mq[_qp] * _c[_qp]* (1.0 - _c[_qp])* _grad_test[_i][_qp] * (_grad_phi[_j][_qp]/(_T[_qp]) - _grad_T[_qp] * _phi[_j][_qp]/(_T[_qp] * _T[_qp]));

In this expression c is the phase field variable and T (temperature) is the coupled variable. 

Yours Sincerely,

Anil Kunwar

[SudiptaBiswas]

_phi[_j][_qp] represents the shape function and it does not change with jacobian terms, remains same. What do you mean by numeriacally consistent?  I did not quite understand your point.  Why do you think it would have difficulty in solving?

[Derek Gaston]

_phi[_j][_qp] will be the value of the shape function underlying the variable you're taking the derivative with respect to.

So if you have two finite element variables:

u = Sum(u_j * phi_j)

v = Sum(v_j * phi_j)

And your computeQpResidual for a Kernel *** operating on u *** looks like:

R = 2*u*v*(1-v)

What goes in computeQpJacobian() is dR/d(u_j):

2*phi_j*v*(1-v)

And what goes in computeQpOffDiagJacobian() for the case where MOOSE is looking for the derivative with respect to v is dR/d(v_j):

2*u*phi_j - 2*2*u*v*phi_j

Just so that is clear it came from:

d/d(v_j) (2*u*v - 2*u*v^2)

The first term is simple: 2*u*phi_j

With the second term you have to remember to do chain rule properly:

dR/d(v_j) = dR/dv * dv/d(v_j)

So, for the second term dR/dv = 2*2*u*v

(Note: I left the 2*2 there just to emphasize what's going on... obviously it's 4 ;-)

And then dv/d(v_j) = phi_j

So the whole second term looks like: 2*2*u*v*phi_j

Therefore:

dR/dv = 2*u*phi_j - 2*2*u*v*phi_j

Derek

--
You received this message because you are subscribed to the Google Groups "moose-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to moose-users...@googlegroups.com.
Visit this group at http://groups.google.com/group/moose-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/be0a20d1-8754-4209-a92e-8f771a7be3aa%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Roma Gurung]

Hi Sudipta and Derek,

Thank you both for your notes on finite element formulation. I have realized that i forgot to note the difference between u and u_j and this was confusing me a lot in Jacobian formulation. Now, things have become quite clear to me. Sudipta, my point on numerical consistency is "How will the solver deal with that expression to compute the things?" I have only dealt with expression with linearity in "c" variable and have not dealt with residuals with non-linearity in variables ? As Derek has presented residual R = 2*u*v(1-v) i.e. non-linearity in v, i assume that such expressions are numerically consistent within MOOSE framework. I need to study more about FEM: the MOOSE way.

Now, with the skills on Jacobian computations, i think i will be able to derive the correct expressions for computeQpJacobian()  and computeQpOffDiagJacobian(). Derek, your notes are really lucid. And Sudipta, you can state expressions concisely. Cheers.

Yours Sincerely,

Anil Kunwar

--
You received this message because you are subscribed to a topic in the Google Groups "moose-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/moose-users/mHHA6NDHkGk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to moose-users...@googlegroups.com.

Visit this group at http://groups.google.com/group/moose-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/CAFfxPjqVtTVaoGfUrx5KDStntzNCPuKTxYDzz9FM%3D%3DfCh76dZg%40mail.gmail.com.

[Roma Gurung]

Hi Derek and Sudipta,

I have tried to utilize your instructions in deriving the mathematical relation for computeQpResidual, computeQpJacobian and computeQpOffDiagJacobian for the MultiSorretDiffusion kernel as coded in https://github.com/anilkunwar/danphe/blob/master/src/kernel/MultiSoretDiffusion.C.

Moreover, I have come to discover that i had written the computeQpJacobian incorrectly for the given residual as shown in the above posts. Now, the kernel contains the corrected code.

The attached pdf file consists of the mathematical expressions.

Thank you very much to both of you for your valuable instructions.

Cheers.

Yours Sincerely

Anil Kunwar

To unsubscribe from this group and stop receiving emails from it, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at http://groups.google.com/group/moose-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/moose-users/be0a20d1-8754-4209-a92e-8f771a7be3aa%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "moose-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/moose-users/mHHA6NDHkGk/unsubscribe.

To unsubscribe from this group and all its topics, send an email to moose-users+unsubscribe@googlegroups.com.

Visit this group at http://groups.google.com/group/moose-users.