Updated the solid solver. (#568)

src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C

@@ -46,7 +46,6 @@ DAObjFuncVonMisesStressKS::DAObjFuncVonMisesStressKS(
     objFuncDict_.readEntry<scalar>("scale", scale_);
 
     objFuncDict_.readEntry<scalar>("coeffKS", coeffKS_);
-
 }
 
 /// calculate the value of objective function
@@ -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)

src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.C

@@ -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;
-        }
-    }
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -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()

src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.H

@@ -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

src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.C

@@ -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)
 {
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -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
@@ -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())
         {

src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.H

@@ -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_;
 
@@ -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
@@ -113,7 +94,6 @@ public:
     virtual label solvePrimal(
         const Vec xvVec,
         Vec wVec);
-
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DASolver/DASolidDisplacementFoam/calculateStressSolidDisplacement.H

@@ -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();
+}

src/adjoint/DASolver/DASolidDisplacementFoam/createFieldsSolidDisplacement.H

@@ -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(
@@ -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());