Permalink
Browse files

Merge pull request #12814 from lindsayad/ad-mat-prop-deps

Fix AD material property dependencies
  • Loading branch information...
lindsayad committed Feb 4, 2019
2 parents db9b647 + 3ef146f commit 099e66315bd0570e81df74fde3a43a2cab1db4de
@@ -54,9 +54,11 @@ class ComputeMaterialsObjectThread : public ThreadedElementLoop<ConstElemRange>
MaterialPropertyStorage & _bnd_material_props;
MaterialPropertyStorage & _neighbor_material_props;

/// Reference to the Material object warehouses
/// This is populated using _fe_problem.getResidualMaterialsWarehouse because it has the union
/// of traditional materials and the residual version of AD materials. We don't need the Jacobian
/// version of the ADMaterial for doing stateful stuff
const MaterialWarehouse & _materials;
const MaterialWarehouse & _ad_materials;

const MaterialWarehouse & _discrete_materials;

std::vector<std::unique_ptr<Assembly>> & _assembly;
@@ -223,6 +223,10 @@ MaterialData::resizeProps(unsigned int size, bool declared_ad)

if (_props[size] == nullptr)
_props[size] = new ADMaterialPropertyObject<T>(declared_ad);
// This branch is necessary in case the frist call to resizeProps for this property id was
// initiated through a getMaterialProperty call, which will have declared_ad = false
else if (declared_ad)
_props[size]->markAD(true);
if (_props_old[size] == nullptr)
_props_old[size] = new ADMaterialPropertyObject<T>;
if (_props_older[size] == nullptr)
@@ -637,7 +637,7 @@ class FEProblemBase : public SubProblem, public Restartable
virtual void addADJacobianMaterial(const std::string & kernel_name,
const std::string & name,
InputParameters parameters);
virtual void addMaterialHelper(MaterialWarehouse & warehouse,
virtual void addMaterialHelper(std::vector<MaterialWarehouse *> warehouse,
const std::string & kernel_name,
const std::string & name,
InputParameters parameters);
@@ -1378,7 +1378,6 @@ class FEProblemBase : public SubProblem, public Restartable
/*
* Return a reference to the material warehouse of Material objects to be computed.
*/
const MaterialWarehouse & getComputeMaterialWarehouse() const { return _materials; }
const MaterialWarehouse & getResidualMaterialsWarehouse() const { return _residual_materials; }
const MaterialWarehouse & getJacobianMaterialsWarehouse() const { return _jacobian_materials; }
const MaterialWarehouse & getDiscreteMaterialWarehouse() const { return _discrete_materials; }
@@ -1671,9 +1670,10 @@ class FEProblemBase : public SubProblem, public Restartable

///@{
// Material Warehouses
MaterialWarehouse _materials; // Traditional materials that MOOSE computes
MaterialWarehouse _residual_materials; // ADMaterials for computing residuals
MaterialWarehouse _jacobian_materials; // ADMaterials for computing jacobians
MaterialWarehouse _residual_materials; // Residual materials. This is the union of traditional
// materials and the residual copy of an ADMaterial
MaterialWarehouse _jacobian_materials; // Jacobian materials. This is the union of traditional
// materials and the Jacobian copy of an ADMaterial
MaterialWarehouse _discrete_materials; // Materials that the user must compute
MaterialWarehouse _all_materials; // All materials for error checking and MaterialData storage
///@}
@@ -39,8 +39,7 @@ ComputeMaterialsObjectThread::ComputeMaterialsObjectThread(
_material_props(material_props),
_bnd_material_props(bnd_material_props),
_neighbor_material_props(neighbor_material_props),
_materials(_fe_problem.getComputeMaterialWarehouse()),
_ad_materials(_fe_problem.getResidualMaterialsWarehouse()),
_materials(_fe_problem.getResidualMaterialsWarehouse()),
_discrete_materials(_fe_problem.getDiscreteMaterialWarehouse()),
_assembly(assembly),
_need_internal_side_material(false),
@@ -63,7 +62,6 @@ ComputeMaterialsObjectThread::ComputeMaterialsObjectThread(ComputeMaterialsObjec
_bnd_material_props(x._bnd_material_props),
_neighbor_material_props(x._neighbor_material_props),
_materials(x._materials),
_ad_materials(x._ad_materials),
_discrete_materials(x._discrete_materials),
_assembly(x._assembly),
_need_internal_side_material(x._need_internal_side_material),
@@ -81,16 +79,14 @@ ComputeMaterialsObjectThread::subdomainChanged()

std::set<MooseVariableFEBase *> needed_moose_vars;
_materials.updateVariableDependency(needed_moose_vars, _tid);
_ad_materials.updateVariableDependency(needed_moose_vars, _tid);
_fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
}

void
ComputeMaterialsObjectThread::onElement(const Elem * elem)
{
if (_materials.hasActiveBlockObjects(_subdomain, _tid) ||
_discrete_materials.hasActiveBlockObjects(_subdomain, _tid) ||
_ad_materials.hasActiveBlockObjects(_subdomain, _tid))
_discrete_materials.hasActiveBlockObjects(_subdomain, _tid))
{
_fe_problem.prepare(elem, _tid);
_fe_problem.reinitElem(elem, _tid);
@@ -111,11 +107,6 @@ ComputeMaterialsObjectThread::onElement(const Elem * elem)
_materials.getActiveBlockObjects(_subdomain, _tid),
n_points,
*elem);
if (_ad_materials.hasActiveBlockObjects(_subdomain, _tid))
_material_props.initStatefulProps(*_material_data[_tid],
_ad_materials.getActiveBlockObjects(_subdomain, _tid),
n_points,
*elem);
}
}
}
@@ -161,12 +152,6 @@ ComputeMaterialsObjectThread::onBoundary(const Elem * elem, unsigned int side, B
face_n_points,
*elem,
side);
if (_ad_materials.hasActiveBoundaryObjects(bnd_id, _tid))
_bnd_material_props.initStatefulProps(*_bnd_material_data[_tid],
_ad_materials.getActiveBoundaryObjects(bnd_id, _tid),
face_n_points,
*elem,
side);
}
}
}
@@ -197,13 +182,6 @@ ComputeMaterialsObjectThread::onInternalSide(const Elem * elem, unsigned int sid
face_n_points,
*elem,
side);
if (_ad_materials[Moose::FACE_MATERIAL_DATA].hasActiveBlockObjects(_subdomain, _tid))
_bnd_material_props.initStatefulProps(
*_bnd_material_data[_tid],
_ad_materials[Moose::FACE_MATERIAL_DATA].getActiveBlockObjects(_subdomain, _tid),
face_n_points,
*elem,
side);
}

const Elem * neighbor = elem->neighbor_ptr(side);
@@ -235,15 +213,6 @@ ComputeMaterialsObjectThread::onInternalSide(const Elem * elem, unsigned int sid
face_n_points,
*neighbor,
neighbor_side);
if (_ad_materials[Moose::NEIGHBOR_MATERIAL_DATA].hasActiveBlockObjects(
neighbor->subdomain_id(), _tid))
_neighbor_material_props.initStatefulProps(
*_neighbor_material_data[_tid],
_ad_materials[Moose::NEIGHBOR_MATERIAL_DATA].getActiveBlockObjects(
neighbor->subdomain_id(), _tid),
face_n_points,
*neighbor,
neighbor_side);
}
}
}
@@ -682,7 +682,8 @@ FEProblemBase::initialSetup()
for (THREAD_ID tid = 0; tid < n_threads; tid++)
{
// Sort the Material objects, these will be actually computed by MOOSE in reinit methods.
_materials.sort(tid);
_residual_materials.sort(tid);
_jacobian_materials.sort(tid);

// Call initialSetup on both Material and Material objects
_all_materials.initialSetup(tid);
@@ -2342,27 +2343,27 @@ FEProblemBase::addMaterial(const std::string & mat_name,
const std::string & name,
InputParameters parameters)
{
addMaterialHelper(_materials, mat_name, name, parameters);
addMaterialHelper({&_residual_materials, &_jacobian_materials}, mat_name, name, parameters);
}

void
FEProblemBase::addADResidualMaterial(const std::string & mat_name,
const std::string & name,
InputParameters parameters)
{
addMaterialHelper(_residual_materials, mat_name, name, parameters);
addMaterialHelper({&_residual_materials}, mat_name, name, parameters);
}

void
FEProblemBase::addADJacobianMaterial(const std::string & mat_name,
const std::string & name,
InputParameters parameters)
{
addMaterialHelper(_jacobian_materials, mat_name, name, parameters);
addMaterialHelper({&_jacobian_materials}, mat_name, name, parameters);
}

void
FEProblemBase::addMaterialHelper(MaterialWarehouse & warehouse,
FEProblemBase::addMaterialHelper(std::vector<MaterialWarehouse *> warehouses,
const std::string & mat_name,
const std::string & name,
InputParameters parameters)
@@ -2400,7 +2401,8 @@ FEProblemBase::addMaterialHelper(MaterialWarehouse & warehouse,
if (discrete)
_discrete_materials.addObject(material, tid);
else
warehouse.addObject(material, tid);
for (auto && warehouse : warehouses)
warehouse->addObject(material, tid);
}

// Non-boundary restricted require face and neighbor objects
@@ -2435,7 +2437,8 @@ FEProblemBase::addMaterialHelper(MaterialWarehouse & warehouse,
if (discrete)
_discrete_materials.addObjects(material, neighbor_material, face_material, tid);
else
warehouse.addObjects(material, neighbor_material, face_material, tid);
for (auto && warehouse : warehouses)
warehouse->addObjects(material, neighbor_material, face_material, tid);

// link parameters of face and neighbor materials
MooseObjectParameterName name(MooseObjectName("Material", material->name()), "*");
@@ -2464,8 +2467,6 @@ FEProblemBase::prepareMaterials(SubdomainID blk_id, THREAD_ID tid)
const std::set<BoundaryID> & ids = _mesh.getSubdomainBoundaryIds(blk_id);
for (const auto & id : ids)
{
_materials.updateBoundaryVariableDependency(id, needed_moose_vars, tid);
_materials.updateBoundaryMatPropDependency(id, needed_mat_props, tid);
if (_currently_computing_jacobian)
{
_jacobian_materials.updateBoundaryVariableDependency(id, needed_moose_vars, tid);
@@ -2508,9 +2509,6 @@ FEProblemBase::reinitMaterials(SubdomainID blk_id, THREAD_ID tid, bool swap_stat
if (_discrete_materials.hasActiveBlockObjects(blk_id, tid))
_material_data[tid]->reset(_discrete_materials.getActiveBlockObjects(blk_id, tid));

if (_materials.hasActiveBlockObjects(blk_id, tid))
_material_data[tid]->reinit(_materials.getActiveBlockObjects(blk_id, tid));

if (_jacobian_materials.hasActiveBlockObjects(blk_id, tid) && _currently_computing_jacobian)
_material_data[tid]->reinit(_jacobian_materials.getActiveBlockObjects(blk_id, tid));

@@ -2537,10 +2535,6 @@ FEProblemBase::reinitMaterialsFace(SubdomainID blk_id, THREAD_ID tid, bool swap_
_bnd_material_data[tid]->reset(
_discrete_materials[Moose::FACE_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));

if (_materials[Moose::FACE_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid))
_bnd_material_data[tid]->reinit(
_materials[Moose::FACE_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));

if (_jacobian_materials[Moose::FACE_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid) &&
_currently_computing_jacobian)
_bnd_material_data[tid]->reinit(
@@ -2572,10 +2566,6 @@ FEProblemBase::reinitMaterialsNeighbor(SubdomainID blk_id, THREAD_ID tid, bool s
_neighbor_material_data[tid]->reset(
_discrete_materials[Moose::NEIGHBOR_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));

if (_materials[Moose::NEIGHBOR_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid))
_neighbor_material_data[tid]->reinit(
_materials[Moose::NEIGHBOR_MATERIAL_DATA].getActiveBlockObjects(blk_id, tid));

if (_jacobian_materials[Moose::NEIGHBOR_MATERIAL_DATA].hasActiveBlockObjects(blk_id, tid) &&
_currently_computing_jacobian)
_neighbor_material_data[tid]->reinit(
@@ -2605,9 +2595,6 @@ FEProblemBase::reinitMaterialsBoundary(BoundaryID boundary_id, THREAD_ID tid, bo
_bnd_material_data[tid]->reset(
_discrete_materials.getActiveBoundaryObjects(boundary_id, tid));

if (_materials.hasActiveBoundaryObjects(boundary_id, tid))
_bnd_material_data[tid]->reinit(_materials.getActiveBoundaryObjects(boundary_id, tid));

if (_jacobian_materials.hasActiveBoundaryObjects(boundary_id, tid) &&
_currently_computing_jacobian)
_bnd_material_data[tid]->reinit(
@@ -3293,7 +3280,6 @@ FEProblemBase::updateActiveObjects()
_internal_side_indicators.updateActive(tid);
_markers.updateActive(tid);
_all_materials.updateActive(tid);
_materials.updateActive(tid);
_residual_materials.updateActive(tid);
_jacobian_materials.updateActive(tid);
_discrete_materials.updateActive(tid);
@@ -0,0 +1,106 @@
[Mesh]
type = GeneratedMesh
dim = 3
nx = 2
ny = 2
nz = 2
[]

[GlobalParams]
displacements = 'disp_x disp_y disp_z'
[]

[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[]

[ADKernels]
[./stress_x]
type = ADStressDivergenceTensors
component = 0
variable = disp_x
[../]
[./stress_y]
type = ADStressDivergenceTensors
component = 1
variable = disp_y
[../]
[./stress_z]
type = ADStressDivergenceTensors
component = 2
variable = disp_z
[../]
[]

[BCs]
[./symmy]
type = PresetBC
variable = disp_y
boundary = bottom
value = 0
[../]
[./symmx]
type = PresetBC
variable = disp_x
boundary = left
value = 0
[../]
[./symmz]
type = PresetBC
variable = disp_z
boundary = back
value = 0
[../]
[./tdisp]
type = PresetBC
variable = disp_z
boundary = front
value = 0.1
[../]
[]

[Materials]
[./elasticity]
type = ComputeIsotropicElasticityTensor
poissons_ratio = 0.3
youngs_modulus = 1e10
[../]
[]

[ADMaterials]
[./stress]
type = ADComputeLinearElasticStress
[../]
[./strain]
type = ADComputeSmallStrain
[../]
[]

[Preconditioning]
[./smp]
type = SMP
full = true
[../]
[]

[Executioner]
type = Transient
dt = 0.05
solve_type = 'NEWTON'

petsc_options_iname = -pc_hypre_type
petsc_options_value = boomeramg

dtmin = 0.05
num_steps = 1
[]

[Outputs]
exodus = true
file_base = "linear-out"
[]
Oops, something went wrong.

0 comments on commit 099e663

Please sign in to comment.