Accessing MaterialPropertyOld of neighboring elements from ElementUserObject

[chopsticks322]

Hi all,

I am implementing a nonlocal mean calculation of equivalent plastic strain (EQPS) in my application and would like to ask about the correct / supported way to access material values from neighboring elements.

My current approach is to implement an ElementUserObject (EqpsNonlocalMean) that:

• loops over neighboring elements (n-ring)
• reconstructs quadrature points using libMesh FEBase
• computes a spatial average of EQPS at each quadrature point
• uses EQPS values from the previous converged timestep

My core question is:

Is it supported to access MaterialPropertyOld values from different elements inside an ElementUserObject?
Or should a different MOOSE abstraction be used instead?

I chose ElementUserObject because:

• the nonlocal mean is evaluated once per timestep (EXEC_TIMESTEP_BEGIN)
• it does not contribute to residual or Jacobian
• conceptually it is a post-processing / averaging operation

Specifically, I found that
_eqps_old[qp]
always refers to the current element's quadrature point,
even when FEBase::reinit(elem) is called for a neighboring element.

This makes sense, but then I am not sure what the recommended way is to
query the EQPS value for (neighbor_elem, qp2) inside a UserObject.

The simplified programs are as follows:

myApp/include/userobjects/EqpsNonlocalMean.h

#include <map>
#include <utility>

/**
 * EqpsNonlocalMean
 * 
 * ElementUserObject that computes the nonlocal mean of EQPS at each quadrature point,
 * using values from the previous converged timestep.
 * 
 * The nonlocal domain is a sphere whose radius is defined by _radius.
 */

 class EqpsNonlocalMean : public ElementUserObject
 {
    public:
    static InputParameters validParams();

    EqpsNonlocalMean(const InputParameters & parameters);

    /***
     * Access the nonlocal mean EQPS for the given (element id, local qp index).
     */
    const std::vector<Real> & value(dof_id_type elem_id, unsigned int qp) const;

    /// Clear cached values at the beginning of a timestep
    virtual void initialize() override;

    /// Compute nonlocal means for the current qp on current element
    virtual void execute() override;

    /// Finalize 
    virtual void finalize() override;

    ///
    virtual void threadJoin(const UserObject & uo) override;

    protected:

    /// Geometry / neighborhood parameters
    const Real _radius;
    const unsigned int _n_neighbor;

    /// Material properties (old state)
    const MaterialProperty<Real> & _eqps_old;

    /// FE type (Necessary for reinit?)
    FEType _fe_type;

    /// Key type for identifying a quadrature point
    /// (element id, local qp index)
    using QPKey = std::pair<dof_id_type, unsigned int>
    
    /// Nonlocal averaged EQPS at each (elem, qp)
    std::map<QPKey, Real> _nonlocal_mean;
 };

myApp/src/userobjects/EqpsNonlocalMean.C

#include "EqpsNonlocalMean.h"

#include "MooseMesh.h"
#include "FEProblem.h"

#include "libmesh/fe.h"
#include "libmesh/quadrature.h"

#include <queue>
#include <unordered_set>

registerMooseObject("myApp", EqpsNonlocalMean);

InputParameters
EqpsNonlocalMean::validParams()
{
    InputParameters params = ElementUserObject::validParams();

    params.addRequiredParam<Real>("radius", "Radius of the sphere nonlocal domain");
    params.addRequiredParam<unsigned int>("n_neighbor", "Number of neighboring element layers (n-ring)");
    params.addRequiredParam<MaterialPropertyName>("eqps", "Equivalent plastic strain (old state is used)");
    params.addRequiredParam<std::string>("fe_type", "FE type string, e.g., TET4, TET10, HEX8, HEX20, TRI3, TRI6, QUAD4");

    params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_BEGIN;

    return params;
}

EqpsNonlocalMean::EqpsNonlocalMean(const Inputparameters & params)
    :   ElementUserObject(params),
        _radius(getParam<Real>("radius")),
        _n_neighbor(getParam<unsigned int>("n_neighbor")),
        _eqps_old(getMaterialPropertyOld<Real>(getParam<MaterialPropertyName>("eqps")))
{
    std::string fe_string = getParam<std::string>("fe_type");
    std::string fe_string = getParam<std::string>("fe_type");
    if (fe_string == "TET4")
        _fe_type = FEType(FIRST, LAGRANGE);
    else if (fe_string == "TET10")
        _fe_type = FEType(SECOND, LAGRANGE);
    else if (fe_string == "HEX8")
        _fe_type = FEType(FIRST, LAGRANGE);
    else if (fe_string == "HEX20")
        _fe_type = FEType(SECOND, LAGRANGE);
    else if (fe_string == "TRI3")
        _fe_type = FEType(FIRST, LAGRANGE);
    else if (fe_string == "TRI6") 
        _fe_type = FEType(SECOND, LAGRANGE);
    else if (fe_string == "QUAD4")
        _fe_type = FEType(FIRST, LAGRANGE);
    else
        mooseError("Unknown FE type string: " + fe_string);
}

void EqpsNonlocalMean::initialize()
{
    _nonlocal_mean.clear();
}

void EqpsNonlocalMean::finalize()
{
}

void EqpsNonlocalMean::threadJoin(const UserObject & uo)
{
}

static void
collectNRing(const Elem * root, unsigned int n_ring, std::unordered_set<const Elem *> & elems)
{
    std::queue<std::pair<const Elem *, unsigned int>> q;
    elems.insert(root);
    q.push({root, 0});

    // BFS loop
    while (!q.empty())
    {
        const Elem * e = q.front().first;
        unsigned int d = q.front().second;
        q.pop();

        if (d == n_ring)
            continue;
        
        for (const Elem * nbr : e->neighbor_ptr_range())
        {
            if (!nbr)
                continue;
            
            if (elems.insert(nbr).second)
                q.push({nbr, d+1});
        }
    }
}

void EqpsNonlocalMean::execute()
{
    const dof_id_type eid = _current_elem -> id();

    // n-ring elements
    std::unordered_set<const Elem *> nring_elems;
    collectNRing(_current_elem, _n_neighbor, nring_elems);

    const unsigned int n_qp = _qrule->n_points();

    std::unique_ptr<FEBase> fe(FEBase::build(_current_elem->dim(), _fe_type));

    std::unique_ptr<QBase> qrule_copy(_qrule->clone());
    fe->attach_quadrature_rule(qrule_copy.get());

    //It seems I need to call the functions I intend to use before `fe->reinit(elem)`.
    fe->get_xyz();
    fe->get_JxW();

    for (unsigned int qp = 0; qp < n_qp; ++qp)
    {
        Real numerator(0.0);
        Real denominator(0.0);

        const Point & x0 = _q_point[qp];
        
        for (const Elem *elem : nring_elems)
        {
            fe->reinit(elem);
            const auto & q_pts = fe->get_xyz();
            const auto & JxW = fe->get_JxW();

            const unsigned int n_qp2 = q_pts.size();
            for (unsigned int qp2 = 0; qp2 < n_qp2; ++qp2)
            {
                const Point r = q_pts[qp2] - x0;
                if (r.norm() <= _radius)
                {
                    //******This refers to current_elem, not elem
                    // How should I access EQPS_old at (elem, qp2)?
                    // numerator += _eqps_old[qp2] * JxW[qp2];
                    denominator += JxW[qp2];
                }
            }
        }

        Real mean(0.0);
        mean = numerator / denominator;
        _nonlocal_mean[{eid, qp}] = mean;
    }
}

const std::vector<Real> & EqpsNonlocalMean::value(dof_id_type elem_id, unsigned int qp) const
{
    return _nonlocal_mean.at({elem_id, qp});
}

[GiudGiud]

Hello

This is definitely a non-standard use of the ElementUO because this UO does not reinit element neighbors, only the current element.

On the neighbor elements you need:

• the quadrature & Jacobian JxW (you need to reinit on the neighbor elements one by one to obtain that.
• the old material properties (you can get this from storage, although it's not standard)
• anything else?

[chopsticks322]

Hi GiudGiud,

Thank you for your clarification.

Based on your suggestion, I switched to accessing the material properties through MaterialData instead of using MaterialPropertyOld directly.

In particular, inside the neighbor loop I now do:

• MaterialData::swap(elem) to switch context to the neighbor element
• retrieve the old state vid getProperty(..., 1, ...)
• use that value at each quadrature point
• then call swapBack(elem)

Modifications:

//EqpsNonlocalMean.h
#include "MaterialData.h"

class EqpsNonlocalMean : public ElementUserObject
{
...
protected:
MaterialData & _material_data;
...
}

//EqpsNonlocalMean.C

/* Constructor */
    _material_data(_fe_problem.getMaterialData(Moose::BLOCK_MATERIAL_DATA, _tid, this))

/* execute() */
void EqpsNonlocalMean::execute()
{
    const dof_id_type eid = _current_elem -> id();

    // n-ring elements
    std::unordered_set<const Elem *> nring_elems;
    collectNRing(_current_elem, _n_neighbor, nring_elems);

    std::unique_ptr<FEBase> fe(FEBase::build(_current_elem->dim(), _fe_type));
    std::unique_ptr<QBase> qrule_copy(_qrule->clone());
    fe->attach_quadrature_rule(qrule_copy.get());
    fe->get_xyz();
    fe->get_JxW();

    const unsigned int n_qp = _qrule->n_points();
    for (unsigned int qp = 0; qp < n_qp; ++qp)
    {
        Real numerator(0.0);
        Real denominator(0.0);

        const Point & x0 = _q_point[qp];
        
        for (const Elem *elem : nring_elems)
        {
            fe->reinit(elem); //reinit with the neighbor element
            const auto & q_pts = fe->get_xyz();
            const auto & JxW = fe->get_JxW();

            const unsigned int n_qp2 = q_pts.size();

            //MaterialData
            _material_data.swap(*elem);
            const auto & eqps_old_elem = _material_data.getProperty<Real>(
                getParam<MaterialPropertyName>("eqps"), 1, *this); // 1=old?

            for (unsigned int qp2 = 0; qp2 < n_qp2; ++qp2)
            {
                const Point r = q_pts[qp2] - x0;
                if (r.norm() <= _radius)
                {
                    // Now qp2 means qp2 in elem (neighbor element)?
                    numerator += eqps_old_elem[qp2] * JxW[qp2];
                    denominator += JxW[qp2];
                }
            }
            _material_data.swapBack(*elem);
        }

        Real mean(0.0);
        mean = numerator / denominator;
        _nonlocal_mean[{eid, qp}] = mean;
    }
}

This seems to correctly give access to the quadrature point values on the neighbor elements.
Could you please confirm if this is the appropriate way to access old material properties in neighbor elements through storage, or if there is a more recommended approach?

Thanks again for your help.

[GiudGiud]

hello

getParam("eqps"),

you'll want to do this outside of the loop, it's on the slower side since it has a string comparison

1, *this); // 1=old?

1=0ld yes and 2="older" which is two steps behind

[chopsticks322]

Thank you for your suggestion.

I have another question regarding the execution ordering of UserObject and Material.
Is it guaranteed that finalize() of such a UserObject will always be completed before or after any Material calculations are performed within the same time step?

If this ordering is guaranteed, I am considering using the results of UserObject inside material models, assuming they are fully finalized at the beginning of the time step. Would this be a valid assumption, or is there any possibility that material evaluation could occur before finalize() complete?

If this is not guaranteed, is there a recommended or standard pattern in MOOSE to enforce or achieve this kind of dependency between UserObjects and Material?

Thanks in advance for your patience and clarification.

[GiudGiud]

Hello

Is it guaranteed that finalize() of such a UserObject will always be completed before or after any Material calculations are performed within the same time step?

Material calculations should be on the fly. So if a UO uses materials, the material properties should be re-computed right before they are needed. An "old" matprop is precomputed from the previous time step so that should always be "computed" before UOs.

What could be happening is that UO does not know about material property dependencies if you are accessing them in a non-standard way. You could use "getMaterialProperty" in the UO constructor to signal the dependency if you are not.

If this ordering is guaranteed, I am considering using the results of UserObject inside material models, assuming they are fully finalized at the beginning of the time step. Would this be a valid assumption, or is there any possibility that material evaluation could occur before finalize() complete?

Now on the other direction, if a material uses a UO, you need to make sure the UO is executed before the material. The execution of UO (and the finalize() operation after that execution) is performed on the specified execute_on schedule. For example, execute_on = TIMESTEP_BEGIN will make the UOs execute before the kernels, BCs etc use the material properties that use results from UOs.

We don't handle dependency on UOs automatically because both options (before and after) make sense to different people depending on whether they want a lagged/explicit/constant over the step or implicitly treated term.

If this is not guaranteed, is there a recommended or standard pattern in MOOSE to enforce or achieve this kind of dependency between UserObjects and Material?

See first response for UOs depending on materials.
For materials depending on UOs see second response.

[chopsticks322]

Thank you for the detailed explanation.

I will make sure to control the execution timing of the UserObject via execute_on.
Also, I will explicitly declare material property dependencies using getMaterialProperty in the UserObject constructor.

Thank you again for your patience and help.