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

Fixed the inaccurate gradients for solidDisplacementFoam. #568

Merged
merged 1 commit into from
Jan 20, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 4 additions & 5 deletions src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,6 @@ DAObjFuncVonMisesStressKS::DAObjFuncVonMisesStressKS(
objFuncDict_.readEntry<scalar>("scale", scale_);

objFuncDict_.readEntry<scalar>("coeffKS", coeffKS_);

}

/// calculate the value of objective function
Expand Down Expand Up @@ -87,16 +86,16 @@ void DAObjFuncVonMisesStressKS::calcObjFunc(
objFuncValue = 0.0;

const objectRegistry& db = mesh_.thisDb();
//const volVectorField& D = db.lookupObject<volVectorField>("D");
const volTensorField& gradD = db.lookupObject<volTensorField>("gradD");
const volScalarField& lambda = db.lookupObject<volScalarField>("solid:lambda");
const volScalarField& mu = db.lookupObject<volScalarField>("solid:mu");
const volScalarField& rho = db.lookupObject<volScalarField>("solid:rho");
const volTensorField& gradD = db.lookupObject<volTensorField>("gradD");

volSymmTensorField sigma = rho * (mu * twoSymm(gradD) + lambda * (I * tr(gradD)));
// volTensorField gradD(fvc::grad(D));
volSymmTensorField sigma = rho * (mu * twoSymm(gradD) + (lambda * I) * tr(gradD));

// NOTE: vonMises stress is scaled by scale_ provided in the objFunc dict
volScalarField vonMises = scale_* sqrt((3.0 / 2.0) * magSqr(dev(sigma)));
volScalarField vonMises = scale_ * sqrt((3.0 / 2.0) * magSqr(dev(sigma)));

scalar objValTmp = 0.0;
forAll(objFuncCellSources, idxI)
Expand Down
75 changes: 6 additions & 69 deletions src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.C
Original file line number Diff line number Diff line change
Expand Up @@ -28,47 +28,11 @@ DAResidualSolidDisplacementFoam::DAResidualSolidDisplacementFoam(
// these are intermediate variables or objects
gradD_(const_cast<volTensorField&>(
mesh.thisDb().lookupObject<volTensorField>("gradD"))),
sigmaD_(const_cast<volSymmTensorField&>(
mesh.thisDb().lookupObject<volSymmTensorField>("sigmaD"))),
divSigmaExp_(const_cast<volVectorField&>(
mesh.thisDb().lookupObject<volVectorField>("divSigmaExp"))),
lambda_(const_cast<volScalarField&>(
mesh.thisDb().lookupObject<volScalarField>("solid:lambda"))),
mu_(const_cast<volScalarField&>(
mesh.thisDb().lookupObject<volScalarField>("solid:mu")))
{

IOdictionary thermalProperties(
IOobject(
"thermalProperties",
mesh.time().constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE));

Switch thermalStress(thermalProperties.lookup("thermalStress"));
if (thermalStress)
{
FatalErrorIn("") << "thermalStress=true not supported" << abort(FatalError);
}

const dictionary& stressControl = mesh.solutionDict().subDict("stressAnalysis");

Switch compactNormalStress(stressControl.lookup("compactNormalStress"));

if (!compactNormalStress)
{
FatalErrorIn("") << "compactNormalStress=false not supported" << abort(FatalError);
}

isTractionDisplacementBC_ = 0;
forAll(D_.boundaryField(), patchI)
{
if (D_.boundaryField()[patchI].type() == "tractionDisplacement")
{
isTractionDisplacementBC_ = 1;
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Expand Down Expand Up @@ -99,54 +63,27 @@ void DAResidualSolidDisplacementFoam::calcResiduals(const dictionary& options)
URes_, pRes_, TRes_, phiRes_: residual field variables
*/

volSymmTensorField sigmaD = mu_ * twoSymm(gradD_) + (lambda_ * I) * tr(gradD_);

volVectorField divSigmaExp = fvc::div(sigmaD - (2 * mu_ + lambda_) * gradD_, "div(sigmaD)");

fvVectorMatrix DEqn(
fvm::d2dt2(D_)
== fvm::laplacian(2 * mu_ + lambda_, D_, "laplacian(DD,D)")
+ divSigmaExp_);
+ divSigmaExp);

DRes_ = DEqn & D_;
normalizeResiduals(DRes);
}

void DAResidualSolidDisplacementFoam::updateDAndGradD()
{
/*
Description:
Update D and gradD.

NOTE: we need to update D boundary conditions iteratively if tractionDisplacement BC
is used this is because tractionDisplacement BC is dependent on gradD, while gradD
is dependent on the D bc values.
*/

// this will be called after doing perturbStates
if (isTractionDisplacementBC_)
{
for (label i = 0; i < daOption_.getOption<label>("maxTractionBCIters"); i++)
{
D_.correctBoundaryConditions();
gradD_ = fvc::grad(D_);
}
}
else
{
D_.correctBoundaryConditions();
gradD_ = fvc::grad(D_);
}
}

void DAResidualSolidDisplacementFoam::updateIntermediateVariables()
{
/*
Description:
Update the intermediate variables that depend on the state variables
*/

this->updateDAndGradD();

sigmaD_ = mu_ * twoSymm(gradD_) + (lambda_ * I) * tr(gradD_);

divSigmaExp_ = fvc::div(sigmaD_ - (2 * mu_ + lambda_) * gradD_, "div(sigmaD)");
gradD_ = fvc::grad(D_);
}

void DAResidualSolidDisplacementFoam::correctBoundaryConditions()
Expand Down
8 changes: 1 addition & 7 deletions src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.H
Original file line number Diff line number Diff line change
Expand Up @@ -28,22 +28,16 @@ class DAResidualSolidDisplacementFoam
{

protected:

/// \name These are state variables, state residuals, and partial derivatives
//@{
volVectorField& D_;
volVectorField DRes_;
//@}

volTensorField& gradD_;
volSymmTensorField& sigmaD_;
volVectorField& divSigmaExp_;
volScalarField& lambda_;
volScalarField& mu_;

label isTractionDisplacementBC_;

void updateDAndGradD();

public:
TypeName("DASolidDisplacementFoam");
// Constructors
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -42,16 +42,13 @@ DASolidDisplacementFoam::DASolidDisplacementFoam(
char* argsAll,
PyObject* pyOptions)
: DASolver(argsAll, pyOptions),
mechanicalPropertiesPtr_(nullptr),
rhoPtr_(nullptr),
muPtr_(nullptr),
lambdaPtr_(nullptr),
EPtr_(nullptr),
nuPtr_(nullptr),
DPtr_(nullptr),
sigmaDPtr_(nullptr),
gradDPtr_(nullptr),
divSigmaExpPtr_(nullptr)
gradDPtr_(nullptr)
{
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Expand All @@ -66,7 +63,6 @@ void DASolidDisplacementFoam::initSolver()
Info << "Initializing fields for DASolidDisplacementFoam" << endl;
Time& runTime = runTimePtr_();
fvMesh& mesh = meshPtr_();
#include "createControlsSolidDisplacement.H"
#include "createFieldsSolidDisplacement.H"
#include "createAdjointSolid.H"
// initialize checkMesh
Expand Down Expand Up @@ -125,43 +121,21 @@ label DASolidDisplacementFoam::solvePrimal(
Info << "Iteration = " << runTime.value() << nl << endl;
}

label iCorr = 0;
scalar initialResidual = 0;
gradD = fvc::grad(D);
volSymmTensorField sigmaD = mu * twoSymm(gradD) + (lambda * I) * tr(gradD);

do
{
fvVectorMatrix DEqn(
fvm::d2dt2(D)
== fvm::laplacian(2 * mu + lambda, D, "laplacian(DD,D)")
+ divSigmaExp);

// get the solver performance info such as initial
// and final residuals
SolverPerformance<vector> solverU = DEqn.solve();
initialResidual = solverU.max().initialResidual();

this->primalResidualControl<vector>(solverU, printToScreen, printInterval, "U");

if (!compactNormalStress_)
{
divSigmaExp = fvc::div(DEqn.flux());
}

gradD = fvc::grad(D);
sigmaD = mu * twoSymm(gradD) + (lambda * I) * tr(gradD);

if (compactNormalStress_)
{
divSigmaExp = fvc::div(
sigmaD - (2 * mu + lambda) * gradD,
"div(sigmaD)");
}
else
{
divSigmaExp += fvc::div(sigmaD);
}

} while (initialResidual > convergenceTolerance_ && ++iCorr < nCorr_);
volVectorField divSigmaExp = fvc::div(sigmaD - (2 * mu + lambda) * gradD, "div(sigmaD)");

fvVectorMatrix DEqn(
fvm::d2dt2(D)
== fvm::laplacian(2 * mu + lambda, D, "laplacian(DD,D)")
+ divSigmaExp);

// get the solver performance info such as initial
// and final residuals
SolverPerformance<vector> solverD = DEqn.solve();

this->primalResidualControl<vector>(solverD, printToScreen, printInterval, "D");

if (this->validateStates())
{
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -52,18 +52,6 @@ class DASolidDisplacementFoam

protected:

/// whether to use compact normal stress
Switch compactNormalStress_;

/// number of correctors
label nCorr_;

/// convergence tolerance for D
scalar convergenceTolerance_;

/// mechanicalProperties pointer
autoPtr<IOdictionary> mechanicalPropertiesPtr_;

/// density field pointer
autoPtr<volScalarField> rhoPtr_;

Expand All @@ -81,17 +69,10 @@ protected:

/// displacement field pointer
autoPtr<volVectorField> DPtr_;

/// stress field pointer
autoPtr<volSymmTensorField> sigmaDPtr_;

/// displacement gradient pointer
autoPtr<volTensorField> gradDPtr_;

/// explicit part of div(sigma) pointer
autoPtr<volVectorField> divSigmaExpPtr_;


public:
TypeName("DASolidDisplacementFoam");
// Constructors
Expand All @@ -113,7 +94,6 @@ public:
virtual label solvePrimal(
const Vec xvVec,
Vec wVec);

};

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Expand Down
Original file line number Diff line number Diff line change
@@ -1,23 +1,25 @@
{
volSymmTensorField sigmaD = mu * twoSymm(gradD) + (lambda * I) * tr(gradD);

volSymmTensorField sigma(
IOobject(
"sigma",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE),
rho* sigmaD);
volSymmTensorField sigma(
IOobject(
"sigma",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE),
rho * sigmaD);

volScalarField sigmaEq(
IOobject(
"sigmaEq",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE),
sqrt((3.0 / 2.0) * magSqr(dev(sigma))));
volScalarField sigmaEq(
IOobject(
"sigmaEq",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE),
sqrt((3.0 / 2.0) * magSqr(dev(sigma))));

Info << "Max sigmaEq = " << max(sigmaEq).value()
<< endl;
Info << "Max sigmaEq = " << gMax(sigmaEq) << endl;

runTime.write();
runTime.write();
}

This file was deleted.

Original file line number Diff line number Diff line change
Expand Up @@ -13,20 +13,6 @@ DPtr_.reset(
mesh));
volVectorField& D = DPtr_();

Info << "Calculating stress field sigmaD\n"
<< endl;

sigmaDPtr_.reset(
new volSymmTensorField(
IOobject(
"sigmaD",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE),
muPtr_() * twoSymm(fvc::grad(D)) + lambdaPtr_() * (I * tr(fvc::grad(D)))));
volSymmTensorField& sigmaD = sigmaDPtr_();

// gradD is used in the tractionDisplacement BC
gradDPtr_.reset(
new volTensorField(
Expand All @@ -37,27 +23,3 @@ gradDPtr_.reset(
IOobject::NO_READ,
IOobject::NO_WRITE),
fvc::grad(D)));

Info << "Calculating explicit part of div(sigma) divSigmaExp\n"
<< endl;
divSigmaExpPtr_.reset(
new volVectorField(
IOobject(
"divSigmaExp",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE),
fvc::div(sigmaD)));
volVectorField& divSigmaExp = divSigmaExpPtr_();

if (compactNormalStress_)
{
divSigmaExp -= fvc::laplacian(2 * muPtr_() + lambdaPtr_(), D, "laplacian(DD,D)");
}
else
{
divSigmaExp -= fvc::div((2 * muPtr_() + lambdaPtr_()) * fvc::grad(D), "div(sigmaD)");
}

mesh.setFluxRequired(D.name());
Loading
Loading