roystgnr / MetaPhysicL

Metaprogramming and operator-overloaded classes for numerical simulations
Other
22 stars 12 forks source link

the (AD) derivative of the function **atan2** #83

Open saturn00000 opened 2 years ago

saturn00000 commented 2 years ago

Hi MetaPhysicL team,

is the (AD) derivative of the function atan2 not defined in "MetaPhysicL"?

If we use std::atan2(_grad_u_qp, _grad_u_qp) for the residual with an ADKernel in MOOSE , then we will have the error:

Linear solve did not converge due to DIVERGED_NANORINF iterations 0
Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
Solve Did NOT Converge!
Finished Solving
roystgnr commented 2 years ago

atan2() is defined for DualNumber, DualExpression, and all the other MetaPhysicL types.

What is the type of _grad_u_qp? git grep _grad_u_qp isn't finding anything in MOOSE for me ... there's a libMesh example where grad_u_qp is a Gradient variable, as I guess one might expect from grad in the name. But atan2(a,b) wouldn't make sense with vector (as opposed to array) arguments.

roystgnr commented 2 years ago

wouldn't make sense with vector (as opposed to array) arguments.

which I guess should qualify

all the other MetaPhysicL types.

e.g. in metaphysicl/numberarray.h we define atan2 term-by-term, but in numbervector.h we just mention it in a commented list under the note

// These functions are hard to consistently define on vectors, as
// opposed to arrays, so let's not define them for now.
saturn00000 commented 2 years ago

Hello @roystgnr

thanks for your reply!

P.S. std::atan2(_grad_u_qp , _grad_u_qp) grafik

saturn00000 commented 2 years ago

@roystgnr

we can use std::atan(_temperature[_qp]) in "MetaPhysicL"

but not for std::atan2(_grad_u[_qp] (1), _grad_u[_qp] (0)),

_grad_u[_qp] (1), _grad_u[_qp] (0) both are not vectors...

roystgnr commented 2 years ago

Ah, I see; the markdown formatting turned those gradient indices into hyperlinks before.

I'd still want to know the type of those arguments; the exact compiler error message you get would probably have the type in it.

saturn00000 commented 2 years ago

Ah, I see; the markdown formatting turned those gradient indices into hyperlinks before.

I'd still want to know the type of those arguments; the exact compiler error message you get would probably have the type in it.

As i mentioned before:

If we use std::atan2(_grad_u[_qp] (1), _grad_u[_qp] (0)) the residual with an ADKernel in MOOSE , then we will have the error:

Linear solve did not converge due to DIVERGED_NANORINF iterations 0
Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
Solve Did NOT Converge!
Finished Solving

Hi MetaPhysicL team,

is the (AD) derivative of the function atan2 not defined in "MetaPhysicL"?

If we use std::atan2(_grad_u_qp, _grad_u_qp) for the residual with an ADKernel in MOOSE , then we will have the error:

Linear solve did not converge due to DIVERGED_NANORINF iterations 0
Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
Solve Did NOT Converge!
Finished Solving

Cite: I'd still want to know the type of those arguments; the exact compiler error message you get would probably have the type in it.

Linear solve did not converge due to DIVERGED_NANORINF iterations 0
Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
Solve Did NOT Converge!
Finished Solving

.......

roystgnr commented 2 years ago

Oh, I misunderstood, then - I thought you'd had some difficulty compiling that line so you tried to get away without it (taking the raw_type, say). If that line compiled then this is the answer to your question:

atan2() is defined for DualNumber

Because in any case where a derivative exists, we never "half-define" a function on DualNumber; it'll either be undefined (e.g. all the new functions in https://en.cppreference.com/w/cpp/numeric/special_functions we haven't gotten to yet) or it'll be defined with an implementation that calculates derivatives too. The macro invocation for atan2 is around line 893 of dualnumber.h.

Forgive my suspicion here, but: what are your initial conditions? If you have a zero gradient at any quadrature point in the field, and if you try to use atan2 to get the gradient's direction, you're going to get NaNs for both the value and the derivatives that come out of that. If that comes up in the middle of a transient solve then depending on your solver options you might find PETSc or MOOSE able to back off and retry gracefully, but if it comes up from your initial conditions you may be stuck.