Boundary correction

applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/heatSourceModel/heatSourceModel.C

@@ -87,6 +87,8 @@ Foam::heatSourceModel::heatSourceModel
     heatSourceModelCoeffs_(sourceDict_.optionalSubDict(type + "Coeffs")),
 
     mesh_(mesh),
+    scanPatchName_(""),
+    scanPatchPoints_(),
     absorptionModel_(nullptr),
     movingBeam_(nullptr)
 {
@@ -95,18 +97,145 @@ Foam::heatSourceModel::heatSourceModel
 
     dimensions_ = heatSourceModelCoeffs_.lookup<vector>("dimensions");
     staticDimensions_ = dimensions_;
-
+
+    normalize_ = sourceDict_.lookupOrDefault<Switch>("normalize", true);
     transient_ = heatSourceModelCoeffs_.lookupOrDefault<Switch>("transient", false);
 
     if (transient_)
     {
         isoValue_ = heatSourceModelCoeffs_.lookup<scalar>("isoValue");
     }
+
+    //- Set up scanPatch boundary point list if normalize is enabled
+    if (normalize_)
+    {
+        //- Get name of patch to scan over
+        scanPatchName_ = sourceDict_.lookup<word>("scanPatchName");
+
+        //- Get instance of polyPatch for scanPatchName
+        label scanPatchID = mesh_.boundaryMesh().findPatchID(scanPatchName_);
+        const polyPatch& scanPatch = mesh_.boundaryMesh()[scanPatchID];
+
+        //- Create list of LOCALLY ordered boundary points
+        const labelList& localScanPatchPoints = scanPatch.boundaryPoints();
+        const labelList& gMeshPoints = scanPatch.meshPoints();
+
+        //- Create list of points on myProc
+        List<point> myProcScanPatchPoints;
+
+        //- Append each LOCAL boundary point to a list of PROCESSOR boundary points
+        forAll (localScanPatchPoints, pointi)
+        {
+            //- Get list of GLOBAL faces from GLOBAL point label
+            const labelList& pointFaces = mesh_.pointFaces()[gMeshPoints[localScanPatchPoints[pointi]]];
+
+            bool isProcFace = false;
+
+            //- Check to see if any face connected to the point is coupled
+            forAll(pointFaces, facei)
+            {              
+                //- Get patch of current face
+                const label& currPatchID = mesh_.boundaryMesh().whichPatch(pointFaces[facei]);
+
+                //- If a valid ID is obtained (i.e. the face is not internal), test if coupled
+                if (currPatchID >= 0)
+                {
+                    const polyPatch& currPatch = mesh_.boundaryMesh()[currPatchID];
+
+                    if (currPatch.coupled())
+                    {
+                        isProcFace = true;
+                        break;
+                    }
+                }
+            }
+
+            //- If the point is not connected to a coupled face, add the coordinates of the point
+            //  to the list of processor scan patch points
+            if (!isProcFace)
+            {
+                myProcScanPatchPoints.append
+                (
+                    mesh_.points()[gMeshPoints[localScanPatchPoints[pointi]]]
+                );
+            }
+        }
+
+        //- Create list of lists of the scan patch points
+        List<List<point>> gatheredScanPatchPoints(Pstream::nProcs());
+
+        //- Populate and gather list onto master processor
+        gatheredScanPatchPoints[Pstream::myProcNo()] = myProcScanPatchPoints;
+        Pstream::gatherList(gatheredScanPatchPoints);
+
+        //- Combine list from all processors into single list
+        for (int i = 0; i < Pstream::nProcs(); i++)
+        {
+            for (int j = 0; j < gatheredScanPatchPoints[i].size(); j++)
+            {
+                scanPatchPoints_.append(gatheredScanPatchPoints[i][j]);
+            }
+        }
+    }      
 }
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+void Foam::heatSourceModel::correctPower
+(
+    volScalarField& qDot_,
+    const scalar& eta_
+)
+{
+    //- Find total resolved power in domain
+    scalar resPower = (fvc::domainIntegrate(qDot_)).value() / eta_;
+
+    //- Only normalize if resolved power is > 0
+    if (resPower > SMALL)
+    {
+        //- Initialize boundary test flag
+        bool isNearBoundary = false;
+
+        //- Get beam position from movingBeam
+        vector position = movingBeam_->position();
+
+        //- Loop over points to determine if beam center is within 2sigma of
+        //  any patch boundary points
+        forAll (scanPatchPoints_, pointi)
+        {
+            scalar bDist = mag(position - scanPatchPoints_[pointi]);
+
+            //- Set flag and break loop if position is near boundary
+            if (bDist < D4sigma()/2.0)
+            {
+                Info << "Beam center within 2sigma of boundary. Skipping normalization." << endl;
+                isNearBoundary = true;
+                break;
+            }
+        }
+
+        //- Normalize if beam isn't near boundary
+        if (!isNearBoundary)
+        {
+            Info << "Beam center not within 2sigma of boundary. Applying normalization." << endl;
+
+            //- Get expected power from movingBeam class
+            scalar expPower = movingBeam_->power();
+
+            //- Normalize qDot
+            qDot_ *= expPower / resPower;
+            qDot_.correctBoundaryConditions();
+        }
+        //- Code hangs without the operations in this else statement, not sure why
+        else
+        {
+            qDot_ *= 1.0;
+            qDot_.correctBoundaryConditions();
+        }
+    }
+}
+
 bool Foam::heatSourceModel::read()
 {
     if (regIOobject::read())

applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/heatSourceModel/heatSourceModel.H

@@ -82,6 +82,15 @@ protected:
 
         //- Pointer to mesh info
         const fvMesh& mesh_;
+
+        //- Flag to enable/disable power correction (default TRUE)
+        Switch normalize_;
+
+        //- Name of patch that beam scans over
+        word scanPatchName_;
+
+        //- List of scan patch bounding points
+        List<point> scanPatchPoints_;
 
         //- Flag for transient or static heat source dimensions
         Switch transient_;
@@ -187,7 +196,17 @@ public:
         }
 
         //- Return the volumetric heating field
-        virtual tmp<volScalarField> qDot() const = 0;
+        virtual tmp<volScalarField> qDot() = 0;
+
+        //- Return D4sigma value of beam shape
+        virtual scalar D4sigma() const = 0;
+
+        //- Correct integrated power
+        virtual void correctPower
+        (
+            volScalarField& qDot_,
+            const scalar& eta_
+        );
 
         //- Read the heat source dictionary
         virtual bool read() = 0;

applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/modifiedSuperGaussian/modifiedSuperGaussian.C

@@ -61,7 +61,7 @@ Foam::heatSourceModels::modifiedSuperGaussian::modifiedSuperGaussian
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::heatSourceModels::modifiedSuperGaussian::qDot() const
+Foam::heatSourceModels::modifiedSuperGaussian::qDot()
 {
     tmp<volScalarField> tqDot
     (
@@ -128,7 +128,7 @@ Foam::heatSourceModels::modifiedSuperGaussian::qDot() const
         {
             point d = cmptMag(mesh_.Cf()[facei] - movingBeam_->position());
 
-            if (d.z() < s.z())
+            if ((mag(d) < D4sigma()/2.0) && (d.z() < s.z()))
             {                
                 scalar g = Foam::pow(1.0 - Foam::pow(d.z() / s.z(), m_), 1.0/m_);
 
@@ -148,7 +148,7 @@ Foam::heatSourceModels::modifiedSuperGaussian::qDot() const
             {
                 point d = cmptMag(faceCentres[facei] - movingBeam_->position());
 
-                if (d.z() < s.z())
+                if ((mag(d) < D4sigma()/2.0) && (d.z() < s.z()))
                 {
                     scalar g =
                         Foam::pow
@@ -165,23 +165,28 @@ Foam::heatSourceModels::modifiedSuperGaussian::qDot() const
                 }
             }
         }
-
+
+        //- Assemble qDot from normalized power and face weights
         qDot_ = I0 * fvc::average(factor);
-
-        //- Rescale the total integrated power to the applied power
-        scalar integratedPower_ = fvc::domainIntegrate(qDot_).value() / eta_;
-
-        if (integratedPower_ > small)
+        qDot_.correctBoundaryConditions();
+
+        //- Perform power correction for coarse meshes
+        if (normalize_)
         {
-            qDot_ *= power_ / integratedPower_;
+            correctPower(qDot_, eta_);
         }
-
-        qDot_.correctBoundaryConditions();
     }
 
     return tqDot;
 }
 
+
+scalar Foam::heatSourceModels::modifiedSuperGaussian::D4sigma() const
+{
+    return 2.0 * max(dimensions_.x(), dimensions_.y());
+}
+
+
 bool Foam::heatSourceModels::modifiedSuperGaussian::read()
 {
     if (heatSourceModel::read())

applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/modifiedSuperGaussian/modifiedSuperGaussian.H

@@ -94,7 +94,10 @@ public:
     // Member Functions
 
         //- Return the source volumetric heating
-        virtual tmp<volScalarField> qDot() const;
+        virtual tmp<volScalarField> qDot();
+
+        //- Return the D4sigma value of the modifiedSuperGaussian
+        virtual scalar D4sigma() const;
 
         //- Read the heatSourceProperties dictionary
         virtual bool read();

applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/superGaussian/superGaussian.C

@@ -61,7 +61,6 @@ Foam::heatSourceModels::superGaussian::superGaussian
 
 Foam::tmp<Foam::volScalarField>
 Foam::heatSourceModels::superGaussian::qDot()
-const
 {
     tmp<volScalarField> tqDot
     (
@@ -84,7 +83,7 @@ const
     const scalar power_ = movingBeam_->power();
 
     if (power_ > small)
-    {        
+    {
         const scalar AR = 
             dimensions_.z() / min(dimensions_.x(), dimensions_.y());
 
@@ -119,10 +118,13 @@ const
         for (label facei=0; facei < mesh_.nInternalFaces(); facei++)
         {
             point d = cmptMag(mesh_.Cf()[facei] - movingBeam_->position());
+
+            if (mag(d) < D4sigma()/2.0)
+            {
+                scalar f = magSqr(cmptDivide(d, s));
 
-            scalar f = magSqr(cmptDivide(d, s));
-
-            factor[facei] = Foam::exp(-Foam::pow(f, k_/2.0));
+                factor[facei] = Foam::exp(-Foam::pow(f, k_/2.0));
+            }
         }
 
         surfaceScalarField::Boundary& bFactor = factor.boundaryFieldRef();
@@ -132,29 +134,37 @@ const
             forAll(faceCentres, facei)
             {
                 point d = cmptMag(faceCentres[facei] - movingBeam_->position());
+
+                if (mag(d) < D4sigma()/2.0)
+                {
+                    scalar f = magSqr(cmptDivide(d, s));
 
-                scalar f = magSqr(cmptDivide(d, s));
-
-                bFactor[patchi][facei] = Foam::exp(-Foam::pow(f, k_/2.0));
+                    bFactor[patchi][facei] = Foam::exp(-Foam::pow(f, k_/2.0));
+                }
             }
         }
 
+        //- Assemble qDot from normalized power and face weights
         qDot_ = I0 * fvc::average(factor);
-
-        //- Rescale the total integrated power to the applied power
-        scalar integratedPower_ = fvc::domainIntegrate(qDot_).value() / eta_;
-
-        if (integratedPower_ > small)
+        qDot_.correctBoundaryConditions();
+
+        //- Perform power correction for coarse meshes
+        if (normalize_)
         {
-            qDot_ *= power_ / integratedPower_;
+            correctPower(qDot_, eta_);
         }
-
-        qDot_.correctBoundaryConditions();
     }
 
     return tqDot;
 }
 
+
+scalar Foam::heatSourceModels::superGaussian::D4sigma() const
+{
+    return 2.0 * max(dimensions_.x(), dimensions_.y());
+}
+
+
 bool Foam::heatSourceModels::superGaussian::read()
 {
     if (heatSourceModel::read())

applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/superGaussian/superGaussian.H

@@ -87,7 +87,10 @@ public:
     // Member Functions
 
         //- Return the source volumetric heating
-        virtual tmp<volScalarField> qDot() const;
+        virtual tmp<volScalarField> qDot();
+
+        //- Return the D4sigma value of the superGaussian
+        virtual scalar D4sigma() const;
 
         //- Read the heatSourceProperties dictionary
         virtual bool read();

tutorials/AMB2018-02-B/constant/heatSourceDict

@@ -19,7 +19,11 @@ sources (beam);
 
 beam
 {
-    pathName            scanPath;    
+    pathName            scanPath;  
+
+    normalize           true;
+
+    scanPatchName       top;  
 
     absorptionModel     constant;
 

tutorials/multiBeam/constant/heatSourceDict

@@ -21,6 +21,10 @@ sources (beam1 beam2);
 beam1
 {
     pathName            scanPath_1;
+
+    normalize           true;
+
+    scanPatchName       top;
 
     absorptionModel     constant;    
     constantCoeffs
@@ -41,7 +45,11 @@ beam1
 beam2
 {
     pathName            scanPath_2;
-
+
+    normalize           true;
+
+    scanPatchName       top;
+
     absorptionModel     Kelly;    
     KellyCoeffs
     {