-
Notifications
You must be signed in to change notification settings - Fork 1k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
ADInterfaceKernel #15397
ADInterfaceKernel #15397
Conversation
Job Partial Modules Parallel Test (In progress) on e1789d9 : invalidated by @lindsayad huh? |
Job Documentation on f7aad99 wanted to post the following: View the site here This comment will be updated on new commits. |
be221a4
to
9ef963b
Compare
|
||
std::size_t ad_offset; | ||
if (type == Moose::ElementElement || type == Moose::NeighborElement) | ||
ad_offset = jvar * _sys.getMaxVarNDofsPerElem(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@cticenhour I believe the key change was here. I think you mostly copy-pasted this from ADDGKernel
maybe where this isn't a concern...but previously I believe you had _var.number()
. So we were getting the jacobian wrong when computing the Jacobians with respect to the neighbor residual because we were always using the "elemental" variable. So the change was to go from _var.number()
-> jvar
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah you're right - I originally had ivar
here. I actually didn't look at ADDGKernel
, but a combo of ADIntegratedBC
and ADKernel
and then built something up from what I saw there and my understanding of InterfaceKernel
. My core misunderstanding was the concept of an ad_offset
in the first place, it seems. Could you explain that a bit more?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have to remind myself about this sometimes because I forget...
Let's just think about continuous finite elements for a second. We use a variable major numbering scheme, such that each variables indices are in a contiguous block. Let's imagine we have two variables, u and v, and we're on a QUAD4. The AD indices will be ordered like this:
u0, u1, u2, u3, v0, v1, v2, v3
_sys.getMaxVarNDofsPerElem
should return for a QUAD4: 4. For a QUAD9, 9. HEX27, 27. Etc. For CFEM the offset will be simply that element number times the var number. So for u on a QUAD4: 4 * 0 = 0. For v: 4 * 1. So u indices start at index 0, v indices start at index 4.
With DFEM or interface kernels it's a little more complicated. We essentially already have an indices block that is num_vars * num_dofs_per_elem long, so in our two var, QUAD4 example: 4 * 2 = 8. So what we do is that if we are on a neighbor element, we do an additional offset by num_vars * num_dofs_per_elem. So now our derivative indices are ordered like this:
u0, u1, u2, u3, v0, v1, v2, v3, u0_neighbor, u1_neighbor, u2_neighbor, u3_neighbor, v0_neighbor, v1_neighbor, v2_neighbor, v3_neighbor
Does that make sense?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Very clear - thanks! :)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let's just think about continuous finite elements for a second. We use a variable major numbering scheme, such that each variables indices are in a contiguous block. Let's imagine we have two variables, u and v, and we're on a QUAD4. The AD indices will be ordered like this:
u0, u1, u2, u3, v0, v1, v2, v3
_sys.getMaxVarNDofsPerElem
should return for a QUAD4: 4. For a QUAD9, 9. HEX27, 27. Etc. For CFEM the offset will be simply that element number times the var number. So for u on a QUAD4: 4 * 0 = 0. For v: 4 * 1. So u indices start at index 0, v indices start at index 4.With DFEM or interface kernels it's a little more complicated. We essentially already have an indices block that is num_vars * num_dofs_per_elem long, so in our two var, QUAD4 example: 4 * 2 = 8. So what we do is that if we are on a neighbor element, we do an additional offset by num_vars * num_dofs_per_elem. So now our derivative indices are ordered like this:
u0, u1, u2, u3, v0, v1, v2, v3, u0_neighbor, u1_neighbor, u2_neighbor, u3_neighbor, v0_neighbor, v1_neighbor, v2_neighbor, v3_neighbor
I would really appreciate it if you could add these sentences as comments into code.
Job Mac Test on 9ef963b : invalidated by @cticenhour Unrelated test crash |
@cticenhour feel free to review this although we'll need someone else to give approval since you have an important commit here |
@rwcarlsen this has your |
|
||
/// The ad version of JxW | ||
const MooseArray<ADReal> & _ad_JxW; | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Interesting.
|
||
std::size_t ad_offset; | ||
if (type == Moose::ElementElement || type == Moose::NeighborElement) | ||
ad_offset = jvar * _sys.getMaxVarNDofsPerElem(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Let's just think about continuous finite elements for a second. We use a variable major numbering scheme, such that each variables indices are in a contiguous block. Let's imagine we have two variables, u and v, and we're on a QUAD4. The AD indices will be ordered like this:
u0, u1, u2, u3, v0, v1, v2, v3
_sys.getMaxVarNDofsPerElem
should return for a QUAD4: 4. For a QUAD9, 9. HEX27, 27. Etc. For CFEM the offset will be simply that element number times the var number. So for u on a QUAD4: 4 * 0 = 0. For v: 4 * 1. So u indices start at index 0, v indices start at index 4.With DFEM or interface kernels it's a little more complicated. We essentially already have an indices block that is num_vars * num_dofs_per_elem long, so in our two var, QUAD4 example: 4 * 2 = 8. So what we do is that if we are on a neighbor element, we do an additional offset by num_vars * num_dofs_per_elem. So now our derivative indices are ordered like this:
u0, u1, u2, u3, v0, v1, v2, v3, u0_neighbor, u1_neighbor, u2_neighbor, u3_neighbor, v0_neighbor, v1_neighbor, v2_neighbor, v3_neighbor
I would really appreciate it if you could add these sentences as comments into code.
ad_offset = jvar * _sys.getMaxVarNDofsPerElem(); | ||
else | ||
ad_offset = jvar * _sys.getMaxVarNDofsPerElem() + | ||
(_sys.system().n_vars() * _sys.getMaxVarNDofsPerElem()); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ad_offset
has to match what you are doing seed? If so, it may be a good idea to create an API somewhere, so we can retrieve the same value
Just minor things |
I agree with @fdkong - that discussion about indexing for derivatives would be very helpful to have as comments in the code. Especially looking forward to FV reconstruction stuff I will probably end up working on eventually - where we will have to add in derivs. w.r.t. even more than just two elements. |
9ef963b
to
258ecaf
Compare
@fdkong @rwcarlsen please see 258ecaf |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The adOffset func and doc comments are fantastic - a huge improvement. It looks like you missed a couple of replacements though - in FV flux kernels/bcs.
Ok hopefully I got the framework holdovers in the latest commit. @joshuahansel please checkout the API I introduced in 0b8daa6 for computing the map you want |
65f9efc
to
7727869
Compare
@lindsayad Thanks, might take me a day or so to test the API, but it looks like it'll do the trick. |
|
||
/** | ||
* A material that couples a material property | ||
*/ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not a very useful doc comment.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Changed to something a bit more meaningful
framework/src/utils/ADUtils.C
Outdated
std::unordered_map<dof_id_type, Real> | ||
globalDofIndexToDerivative(const ADReal & ad_real, | ||
const SystemBase & sys, | ||
const Elem * elem, | ||
const ElementType elem_type /*=ElementType::Element*/) | ||
{ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What do we anticipate this being used for in the future? Did you have anything in mind?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See #15450
dof_map.dof_indices(elem, global_indices, var_num); | ||
|
||
// determine the AD offset for the current var | ||
const auto ad_offset = adOffset(var_num, max_dofs_per_elem, elem_type, num_vars); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What happens if elem
isn't currently assembly's current elem or neighbor? Then we just get nonsense right? Maybe we should guard against that.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Or at least we need to know which dof indices were used when the dualnumbers that contributed to ad_real
were seeded/initialized.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Might as well just remove elem
from the argument list then and just take it from assembly
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I went ahead and removed elem
from the argument list, and am just using the assembly elem APIs to get the current elements
Also move ADOffset to be ADUtils Closes idaholab#15450
This is actually only really a test in debug mode
7727869
to
5f8a0e0
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I anticipate more work being done here eventually. I expect as Josh starts using this, and as we eventually start trying to to fv gradient reconstruction, more ideas will pop out about how to further improve this situation. Wouldn't it be nice if we just stored the global dof id's straight in the dual numbers?
I definitely agree with you. We could do this right now with the sparse container... |
Support for both standard and vector variables.
This PR also contains a lot of
const
additions toCoupleable/NeighborCoupleable/MooseVariable*
APIs since they should be such.Closes #15307