Eikonal equation

AlexGubkin · Nov 23, 2022 · e9e516a · e9e516a
1 parent bb2dff1
commit e9e516a
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 137 deletions.
diff --git a/applications/solvers/eikonalFoam/createFields.H b/applications/solvers/eikonalFoam/createFields.H
@@ -13,18 +13,34 @@ wordList d2Types
     zeroGradientFvPatchField<scalar>::typeName
 );
 
-// wordList ndTypes
-// (
-//     mesh.boundary().size(),
-//     zeroGradientFvPatchField<vector>::typeName
-// );
-
 forAllConstIter(labelHashSet, patchIDs, iter)
 {
     dTypes[iter.key()] = fixedValueFvPatchField<scalar>::typeName;
-//     ndTypes[iter.key()] = fixedValueFvPatchField<vector>::typeName;
 }
 
+
+Info<< "Creation field alpha\n" << endl;
+volScalarField alpha
+(
+    IOobject
+    (
+        "alpha",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(dimless, 0)
+);
+
+surfaceScalarField deltaf
+(
+    "deltaf",
+    fvc::snGrad(alpha)
+);
+
+
 Info<< "Creation distance field d\n" << endl;
 volScalarField d
 (
@@ -40,6 +56,7 @@ volScalarField d
 );
 
 d = dimensionedScalar(dimLength, 1.0);
+// d = (mag(deltaf) > small) ? dimensionedScalar(dimLength, small) : dimensionedScalar(dimLength, 1.0);
 
 Info<< "Creation field gradd\n" << endl;
 volVectorField gradd
@@ -115,59 +132,3 @@ volVectorField graddPsi
     mesh,
     dimensionedVector(dimLength, Zero)
 );
-
-// Info<< "Creation field d2Psi\n" << endl;
-// volScalarField d2Psi
-// (
-//     IOobject
-//     (
-//         "d2Psi",
-//         runTime.timeName(),
-//         mesh
-//     ),
-//     mesh,
-//     dimensionedScalar(sqr(dimLength), 0),
-//     d2.boundaryFieldRef().types()
-// );
-//
-// Info<< "Creation field gradd2Psi\n" << endl;
-// volVectorField gradd2Psi
-// (
-//     IOobject
-//     (
-//         "gradd2Psi",
-//         runTime.timeName(),
-//         mesh,
-//         IOobject::NO_READ,
-//         IOobject::AUTO_WRITE
-//     ),
-//     mesh,
-//     dimensionedVector(dimLength, Zero)
-// );
-
-
-
-// Info<< "Creation field nd\n" << endl;
-// volVectorField nd
-// (
-//     IOobject
-//     (
-//         "nd",
-//         runTime.timeName(),
-//         mesh,
-//         IOobject::NO_READ,
-//         IOobject::AUTO_WRITE
-//     ),
-//     mesh,
-//     dimensionedVector(dimless, Zero),
-//     ndTypes
-// );
-
-// const fvPatchList& patches = mesh.boundary();
-// volVectorField::Boundary& ndbf = nd.boundaryFieldRef();
-// 
-// forAllConstIter(labelHashSet, patchIDs, iter)
-// {
-//     label patchi = iter.key();
-//     ndbf[patchi] == -patches[patchi].nf();
-// }
diff --git a/applications/solvers/eikonalFoam/eikonalFoam.C b/applications/solvers/eikonalFoam/eikonalFoam.C
@@ -87,9 +87,9 @@ int main(int argc, char *argv[])
         if (simple.finalNonOrthogonalIter())
         {
             graddPsi = fvc::grad(dPsi);
-//             volScalarField magGraddPsi(mag(graddPsi));
+            volScalarField magGraddPsi(mag(graddPsi));
 
-            d = sqrt(magSqr(graddPsi) + 2*dPsi) - mag(graddPsi);
+            d = sqrt(magGraddPsi*magGraddPsi + 2*dPsi) - magGraddPsi;
         }
     }
 
@@ -123,50 +123,7 @@ int main(int argc, char *argv[])
     d.write();
     gradd.write();
 
-//     // Non-orthogonal velocity potential corrector loop
-//     label d2PsiRefCell = 0;
-//     scalar d2PsiRefValue = 0.0;
-//
-//     while (simple.correctNonOrthogonal())
-//     {
-//         fvScalarMatrix d2PsiEqn
-//         (
-//             fvm::laplacian(d2Psi)
-//          ==
-//           - neg0(mag(graddPsi)/dimensionedScalar(dimLength, 1.0) - dimensionedScalar(dimless, 0.05))
-//         );
-//
-//         d2PsiEqn.relax();
-//         d2PsiEqn.setReference(d2PsiRefCell, d2PsiRefValue);
-//         d2PsiEqn.solve();
-//
-//         if (simple.finalNonOrthogonalIter())
-//         {
-// //             d2Psi -= dimensionedScalar(sqr(dimLength), gMin(d2Psi));
-//             gradd2Psi = fvc::grad(d2Psi);
-// //             volScalarField magGraddPsi(mag(graddPsi));
-//
-//             d2 = sqrt(magSqr(gradd2Psi) + 2*d2Psi) - mag(gradd2Psi);
-//         }
-//     }
-//
-//     // Write dPsi and graddPsi
-//     d2Psi.write();
-//     gradd2Psi.write();
-//
-//     // Non-orthogonal eikonal corrector loop
-// //     label d2RefCell = 0;
-// //     scalar d2RefValue = 0.0;
-//
-// //     d2 =
-// //         - dPsi/dimensionedScalar(dimLength, 1.0)
-// //         + dimensionedScalar(dimLength, gMax(dPsi))
-// //         + dimensionedScalar(dimLength, small);
-//
-//     d2 = d2Psi/dimensionedScalar(dimLength, 1.0);
-
-//     DynamicList<label>  removedCellsFromMatrix;
-//     DynamicList<scalar> valuesToImpose;
+
     List<label>  removedCellsFromMatrix;
     List<scalar> valuesToImpose;
 
@@ -200,8 +157,6 @@ int main(int argc, char *argv[])
           - epsilon*d2*fvm::laplacian(d2)
          ==
             dimensionedScalar(dimless, 1.0)
-//             pos(mag(gradd) - dimensionedScalar(dimless, 0.8))
-//           + dimensionedScalar(dimless, 1.0)*neg(mag(gradd) - dimensionedScalar(dimless, 0.8))
         );
 
         d2Eqn.setValues(removedCellsFromMatrix, valuesToImpose);
@@ -219,32 +174,6 @@ int main(int argc, char *argv[])
     // Write d and gradd
     d2.write();
     gradd2.write();
-
-//     int iter = 0;
-//     scalar initialResidual = 0;
-//     do
-//     {
-//         nd = fvc::grad(d);
-//         nd /= (mag(nd) + small);
-// 
-//         surfaceVectorField nf(fvc::interpolate(nd));
-//         nf /= (mag(nf) + small);
-// 
-//         surfaceScalarField dPhi("dPhi", nf & mesh.Sf());
-// 
-//         fvScalarMatrix dEqn
-//         (
-//             fvm::div(dPhi, d)
-//           - fvm::Sp(fvc::div(dPhi), d)
-//           - epsilon*d*fvm::laplacian(d)
-//         ==
-//             dimensionedScalar(dimless, 1.0)
-//         );
-// 
-//         dEqn.relax();
-//         initialResidual = dEqn.solve().initialResidual();
-//  
-//     } while (initialResidual > tolerance && ++iter < maxIter);
 }