Skip to content
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

Add grad_zero for phi and test function #15205

Merged
merged 10 commits into from
May 19, 2020

Conversation

arovinelli
Copy link
Contributor

Just add two constants gradient variable

@lindsayad as per mailing list discussion

ref #15204

Comment on lines 1607 to 1608
_grad_phi_zero[tid].resize(max_qpts, {0., 0., 0.});
_grad_test_zero[tid].resize(max_qpts, {0., 0., 0.});
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lindsayad this init here is ugly but I can't figure out how to init with RealVectorValue...
suggestions?

@lindsayad
Copy link
Member

I need to see a test to believe this is working. For one thing the outer size should be according to the number of shape functions, not the number of q-points. Perhaps you are getting away with this because for first-order Lagrange those numbers are the same in the element volume? Also my current reading of the {0,0,0} initializer list is that its going to create three RealVectorValues with each component initialized to zero. That length of three doesn't seem like it's going to match an arbitrary number of quadrature points. What are you testing this on?

@moosebuild
Copy link
Contributor

moosebuild commented May 2, 2020

Job Documentation on e63ae8d wanted to post the following:

View the site here

This comment will be updated on new commits.

@arovinelli
Copy link
Contributor Author

I need to see a test to believe this is working. For one thing the outer size should be according to the number of shape functions, not the number of q-points. Perhaps you are getting away with this because for first-order Lagrange those numbers are the same in the element volume? Also my current reading of the {0,0,0} initializer list is that its going to create three RealVectorValues with each component initialized to zero. That length of three doesn't seem like it's going to match an arbitrary number of quadrature points. What are you testing this on?

wait, that is not how the already available _grad_zero works. If I understand the avaialble code correctly, it is simply a _grad_zero is simply a RealVectorValue(0) with a size equal to the number of _qp such that you can do something like in this example
So you just need _grad_phi_zero to be a constant reference to a RealVector.Am I wrong?

I'm not testing it right now, try to make a kernel or material for this

@lindsayad
Copy link
Member

A VariableGradient is a MooseArray<RealVectorValue>. A VariablePhiGradient is a MooseArray<std::vector<RealVectorValue>. So the latter takes two indexes, the first for dofs, the second for quadrature points. The former only takes a single index for quadrature points.

@arovinelli
Copy link
Contributor Author

arovinelli commented May 2, 2020

A VariableGradient is a MooseArray<RealVectorValue>. A VariablePhiGradient is a MooseArray<std::vector<RealVectorValue>. So the latter takes two indexes, the first for dofs, the second for quadrature points. The former only takes a single index for quadrature points.

This makes a lot of sense and explains why I was not able to find the right initialization for it.
So it should something like

    _grad_test_zero[tid].resize(getMaxShapeFunctions(),
                                std::vector<RealGradient>(max_qpts, RealGradient(0.)));

@lindsayad
Copy link
Member

Correct. You might actually want to look at the code around this line: https://github.com/idaholab/moose/blob/next/framework/src/problems/FEProblemBase.C#L627. We actually compute exactly what you're looking for there

@moosebuild
Copy link
Contributor

Job Modules debug PETSc alt on 2d2630c : invalidated by @arovinelli

uhm... this should be related to this pr

@arovinelli
Copy link
Contributor Author

@lindsay I'm also trying to construct a VariapliPhiValue _phi_zero. What type is and what size should a VariapliPhiValue ?

@lindsayad
Copy link
Member

A VariablePhiValue is a MooseArray<std::vector<Real>.

@arovinelli
Copy link
Contributor Author

@lindsay it isn't as straightforward as I thought.
first, the number of DOF might be different between nodes and elements _nl->getMaxVarNDofsPerElem _nl->getMaxVarNDofsPerNode so I have to choose which of this to use. Also when called in FEProbelmBase::createQRules they both return zero, therefore I assume FEProbelmBase::intialSetup has not been called yet.

I tried to use the code you indicated , inside FEProbelmBase::createQRules but it gives me a segfault

    MaxVarNDofsPerElem mvndpe(*this, *_nl);
    Threads::parallel_reduce(*_mesh.getActiveLocalElementRange(), mvndpe);
    auto max_var_n_dofs_per_elem = mvndpe.max();
    _communicator.max(max_var_n_dofs_per_elem);

here is the stack trace

#0  0x00007fffe8596a22 in __gnu_debug::_Error_formatter::_Parameter::_M_print_field(__gnu_debug::_Error_formatter const*, char const*) const () from /lib64/libstdc++.so.6
#1  0x00007fffe8596cfb in __gnu_debug::_Error_formatter::_Parameter::_M_print_description(__gnu_debug::_Error_formatter const*) const () from /lib64/libstdc++.so.6
#2  0x00007fffe859738f in __gnu_debug::_Error_formatter::_M_error() const () from /lib64/libstdc++.so.6
#3  0x00007ffff743ca0b in std::__debug::vector<unsigned int, std::allocator<unsigned int> >::operator[] (this=0x892f10, __n=0)
    at /opt/rh/devtoolset-7/root/usr/include/c++/7/debug/vector:424
#4  0x00007fffeef8a9e5 in libMesh::DofMap::dof_indices (this=0x892e80, elem=0x85b650, 
    di=std::__debug::vector of length 0, capacity 0, vn=0, p_level=0)
    at /home/arovinelli/projects/NEW_stuff/moose/scripts/../libmesh/src/base/dof_map.C:2093
#5  0x00007ffff5599308 in MaxVarNDofsPerElem::onElement (this=0x7fffffffb930, elem=0x85b650)
    at /home/arovinelli/projects/NEW_stuff/moose/framework/src/base/MaxVarNDofsPerElem.C:39
#6  0x00007ffff47ba86b in ThreadedElementLoopBase<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> >::operator() (this=0x7fffffffb930, range=..., bypass_threading=false)
    at /home/arovinelli/projects/NEW_stuff/moose/framework/build/header_symlinks/ThreadedElementLoopBase.h:209
#7  0x00007ffff4d4b073 in libMesh::Threads::parallel_reduce<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*>, MaxVarNDofsPerElem> (range=..., body=...)
    at /home/arovinelli/projects/NEW_stuff/moose/scripts/../libmesh/installed/include/libmesh/threads_pthread.h:380
#8  0x00007ffff4d21b23 in FEProblemBase::createQRules (this=0x8685a0, type=libMesh::QGAUSS, order=libMesh::THIRD, 
    volume_order=libMesh::THIRD, face_order=libMesh::THIRD)
    at /home/arovinelli/projects/NEW_stuff/moose/framework/src/problems/FEProblemBase.C:4580
#9  0x00007ffff4c676a2 in SetupQuadratureAction::act (this=0x81fc90)
    at /home/arovinelli/projects/NEW_stuff/moose/framework/src/actions/SetupQuadratureAction.C:51

so I can't compute either the maxNdof per elem or node

@lindsayad
Copy link
Member

lindsayad commented May 13, 2020

It might be helpful for you to look at the order of actions in Moose.C, starting around line 245. Quadrature setup happens well before the FEProblemBase is initialized (the init_problem task calls FEProblemBase::init), which is where the libMesh EquationSystems and DofMap objects are initialized, which is where dof indices and things will get setup.

Why don't you wait to do your zeroing and setup until FEProblemBase::initialSetup which is where we already do the MaxVarNDofsPerElem calculation? You can choose to make that calculation non-optional if you want (we currently only do it if there are AD objects in the problem)

@arovinelli
Copy link
Contributor Author

@lindsay, thanks a lot for the suggestion I've embraced it and I'm making good progress.
At this point, I believe the size of _second_phi_zero wrong and also initialized in the wrong place, could you check (see link)

_second_phi_zero[tid].resize(
max_qpts, std::vector<RealTensor>(getMaxShapeFunctions(), RealTensor(0.)));

@lindsayad
Copy link
Member

lindsayad commented May 13, 2020 via email

@moosebuild
Copy link
Contributor

Job Test debug on 2a4ef51 : invalidated by @arovinelli

these failures should not be related to this PR

@arovinelli
Copy link
Contributor Author

arovinelli commented May 14, 2020

So this

_second_phi_zero[tid].resize(
max_qpts, std::vector<RealTensor>(getMaxShapeFunctions(), RealTensor(0.)));

is a preexisting bug.
This PR should fix it.

Looking at MaxQpsThread::operator(), getMaxShapeFunctions will definitely return the wrong value if there are second order bases.

}

PhiZeroKernel::PhiZeroKernel(const InputParameters & parameters)
: NullKernel(parameters), _second_phi(_assembly.secondPhi(_var))
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lindsayad Why I can't get _second_phi.size() > 0? In my tests, I use second-order element and second-order lagrangian function

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just calling secondPhi in Assembly is not enough to get MOOSE to do second derivatives. You need to actually request something like _var.secondSln() in order to get it to happen. A hacky way that may also work is to directly call _assembly.feSecondPhi<Real>(FEType(LAGRANGE,FIRST))

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks, now it works!

@moosebuild
Copy link
Contributor

Job Create Linux Stack on 667552f : invalidated by @arovinelli

stucked

@arovinelli
Copy link
Contributor Author

@lindsayad when you have time, this is ready for review

Copy link
Member

@lindsayad lindsayad left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A nice addition. thank you @arovinelli . I'm curious what your use case is?

Mostly fixup the test specs file

if (haveADObjects() || (_displaced_problem && _displaced_problem->haveADObjects()))
{
CONSOLE_TIMED_PRINT("Computing max dofs per elem/node");
// alwasy esecute to get the max number of DoF per elment and node needed to initialized phi_zero
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

always execute


#include "NullKernel.h"

class PhiZeroKernel;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need for this forward declaration

Comment on lines 16 to 18
/**
*
*/
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Either add a doxygen comment or just delete this

Comment on lines 35 to 36
"_phi.size() " + std::to_string(_phi.size()) + "> _phi_zero.size() " +
std::to_string(_phi_zero.size()));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I understand that you were thinking about cases in which the "zero" members were of zero size, but your message should probably just say "!=" as opposed to ">". Same for all the following assertions

input = 'simple_transient_diffusion.i'
expect_out = ' Solve Converged!'
issues = '#15204'
design = ' '
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ruh roh. You should either write some documentation or find some existing documentation that roughly fits what you're working on

issues = '#15204'
design = ' '
requirement = 'The system shall be able to produce consestent size _phi_zero and _grad_phi_zero objects'
method = 'dbg'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You might as well remove this as well. I understand that your assertions won't be tested, but the problem should still run for all methods

cli_args = 'Mesh/elem_type=QUAD8 Variables/dummy/order=SECOND Variables/u/order=SECOND Kernels/phi_zero/second_order=true'
requirement = 'The system shall be able to produce consestent size _phi_zero and _grad_phi_zero'
'objects for quadratic shape functions'
expect_out = ' Solve Converged!'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again can remove

'objects for quadratic shape functions'
expect_out = ' Solve Converged!'
issues = '#15204'
design = ' '
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

need doc

type = 'RunApp'
input = 'simple_transient_diffusion.i'
cli_args = 'Mesh/elem_type=QUAD8 Variables/dummy/order=SECOND Variables/u/order=SECOND Kernels/phi_zero/second_order=true'
requirement = 'The system shall be able to produce consestent size _phi_zero and _grad_phi_zero'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could use a similar requirement to the one I proposed above

expect_out = ' Solve Converged!'
issues = '#15204'
design = ' '
method = 'dbg'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove

@arovinelli
Copy link
Contributor Author

A nice addition. thank you @arovinelli . I'm curious what your use case is?

@lindsay This is a preparation step for implementing large deformation cohesive zone modeling. The main reason I needed this is to ease the formulation of the jacobian of the traction w.r.t. the discreet displacements. It might not be optimal in terms of efficiency but make the code much more understandable. It will also allow using exactly the same formulation for 1D 2D and 3D problems, thus avoiding a bunch of if statements that are detrimental when performing maintenance.

Mostly fixup the test specs file

The last PR should address all your comments. Check the design and requirement parts to see if they comply with your specs.

Alex, thanks again for helping with this!!

@moosebuild
Copy link
Contributor

Job Private App tests PETSc alt on bd98388 : invalidated by @lindsayad

@moosebuild
Copy link
Contributor

Job Documentation on bd98388 : invalidated by @arovinelli

something is funny, try rerun

@arovinelli
Copy link
Contributor Author

I can't put the .md file in the test doc directory, as it will not be found. Suggestions?

@arovinelli
Copy link
Contributor Author

@lindsay this has passed all tests, please give it another review and let me know if there is something left to do.

Copy link
Member

@lindsayad lindsayad left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I really like how this looks. Just the two tiny doc requests, and then I'll approve!

# PhiZero

MOOSE has helper zero objects for shape functions and shape function gradients that can be accessed by all MOOSE objects. The names of the zero objects are: `_phi_zero`, `_grad_phi_zero` and `_second_phi_zero`.
These object first two dimensions are,: (i) the maximum number of degree of freedom per variable per element, and (ii) the maximum number of quadrature points per element.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove comma in are,


MOOSE has helper zero objects for shape functions and shape function gradients that can be accessed by all MOOSE objects. The names of the zero objects are: `_phi_zero`, `_grad_phi_zero` and `_second_phi_zero`.
These object first two dimensions are,: (i) the maximum number of degree of freedom per variable per element, and (ii) the maximum number of quadrature points per element.
The `PhiZeroKernel` checks that the first two dimensions of `_phi_zero`, `_grad_phi_zero` and `_second_phi_zero` are
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably say the "PhiZeroKernel test object"

@lindsayad lindsayad merged commit cec65dc into idaholab:next May 19, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants