From 813b876ade8580e31b9b58f9b1baad1854d14ed4 Mon Sep 17 00:00:00 2001
From: Alex <alexshtil@gmail.com>
Date: Thu, 24 Nov 2022 18:37:34 +0500
Subject: [PATCH] Solver for eikonal equation

---
 .../solvers/eikonalFoam/createFields.H        | 28 +++++--
 .../solvers/eikonalFoam/eikonalFoam.C         | 27 ++++++-
 tests/eikonal/cylinders2D/0.orig/alpha        | 74 +++++++++++++++++++
 tests/eikonal/cylinders2D/Allclean            |  2 +-
 tests/eikonal/cylinders2D/Allrun              | 19 ++++-
 .../cylinders2D/system/decomposeParDict       |  2 +-
 .../eikonal/cylinders2D/system/setFieldsDict  | 46 ++++++++++++
 7 files changed, 182 insertions(+), 16 deletions(-)
 create mode 100644 tests/eikonal/cylinders2D/0.orig/alpha
 create mode 100644 tests/eikonal/cylinders2D/system/setFieldsDict

diff --git a/applications/solvers/eikonalFoam/createFields.H b/applications/solvers/eikonalFoam/createFields.H
index ae915ee..e2c5e7e 100644
--- a/applications/solvers/eikonalFoam/createFields.H
+++ b/applications/solvers/eikonalFoam/createFields.H
@@ -27,19 +27,27 @@ volScalarField alpha
         "alpha",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::MUST_READ,
         IOobject::AUTO_WRITE
     ),
-    mesh,
-    dimensionedScalar(dimless, 0)
+    mesh
 );
 
-surfaceScalarField deltaf
+volVectorField gradAlpha
 (
-    "deltaf",
-    fvc::snGrad(alpha)
+    IOobject
+    (
+        "gradAlpha",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    fvc::grad(alpha)
 );
 
+gradAlpha.write();
+
 
 Info<< "Creation distance field d\n" << endl;
 volScalarField d
@@ -56,7 +64,13 @@ volScalarField d
 );
 
 d = dimensionedScalar(dimLength, 1.0);
-// d = (mag(deltaf) > small) ? dimensionedScalar(dimLength, small) : dimensionedScalar(dimLength, 1.0);
+forAll(d, celli)
+{
+    if (mag(gradAlpha[celli]) > small)
+    {
+        d[celli] = small;
+    }
+}
 
 Info<< "Creation field gradd\n" << endl;
 volVectorField gradd
diff --git a/applications/solvers/eikonalFoam/eikonalFoam.C b/applications/solvers/eikonalFoam/eikonalFoam.C
index cddbb21..1caf269 100644
--- a/applications/solvers/eikonalFoam/eikonalFoam.C
+++ b/applications/solvers/eikonalFoam/eikonalFoam.C
@@ -97,6 +97,24 @@ int main(int argc, char *argv[])
     dPsi.write();
     graddPsi.write();
 
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    List<label>  removedCellsFromMatrix;
+    List<scalar> valuesToImpose;
+
+    forAll(d, celli)
+    {
+        if (mag(gradAlpha[celli]) > small)
+        {
+            removedCellsFromMatrix.append(celli);
+        }
+    }
+
+    forAll(removedCellsFromMatrix, celli)
+    {
+        valuesToImpose.append(small);
+    }
+
     // Non-orthogonal eikonal corrector loop
     while (eikonal.correctNonOrthogonal())
     {
@@ -115,6 +133,8 @@ int main(int argc, char *argv[])
             dimensionedScalar(dimless, 1.0)
         );
 
+        dEqn.setValues(removedCellsFromMatrix, valuesToImpose);
+
         dEqn.relax();
         dEqn.solve();
     }
@@ -123,9 +143,10 @@ int main(int argc, char *argv[])
     d.write();
     gradd.write();
 
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    List<label>  removedCellsFromMatrix;
-    List<scalar> valuesToImpose;
+    removedCellsFromMatrix.clear();
+    valuesToImpose.clear();
 
     forAll(d2, celli)
     {
@@ -142,8 +163,6 @@ int main(int argc, char *argv[])
 
     while (eikonal.correctNonOrthogonal())
     {
-//         d2 = d2*pos(mag(gradd) - dimensionedScalar(dimless, 0.8));
-
         gradd2 = fvc::grad(d2);
 
         surfaceVectorField gradd2f(fvc::interpolate(gradd2));
diff --git a/tests/eikonal/cylinders2D/0.orig/alpha b/tests/eikonal/cylinders2D/0.orig/alpha
new file mode 100644
index 0000000..c2cd187
--- /dev/null
+++ b/tests/eikonal/cylinders2D/0.orig/alpha
@@ -0,0 +1,74 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  9
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alpha;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    back
+    {
+        type            empty;
+    }
+    front
+    {
+        type            empty;
+    }
+    left
+    {
+        type            zeroGradient;
+    }
+    bottom
+    {
+        type            zeroGradient;
+    }
+    right
+    {
+        type            zeroGradient;
+    }
+    top
+    {
+        type            zeroGradient;
+    }
+
+    cylinder1
+    {
+        type            zeroGradient;
+    }
+    cylinder2
+    {
+        type            zeroGradient;
+    }
+    cylinder3
+    {
+        type            zeroGradient;
+    }
+    cylinder4
+    {
+        type            zeroGradient;
+    }
+    cylinder5
+    {
+        type            zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tests/eikonal/cylinders2D/Allclean b/tests/eikonal/cylinders2D/Allclean
index 7122bca..c35b10f 100755
--- a/tests/eikonal/cylinders2D/Allclean
+++ b/tests/eikonal/cylinders2D/Allclean
@@ -6,7 +6,7 @@ cd ${0%/*} || exit 1    # Run from this directory
 rm -rf processor*/
 rm -rf dynamicCode/
 rm -rf constant/
-rm -rf *[0-9]*
+rm -rf 0/
 rm -f log.*
 
 #------------------------------------------------------------------------------
diff --git a/tests/eikonal/cylinders2D/Allrun b/tests/eikonal/cylinders2D/Allrun
index 03e1829..398020a 100755
--- a/tests/eikonal/cylinders2D/Allrun
+++ b/tests/eikonal/cylinders2D/Allrun
@@ -6,19 +6,32 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 application=$(getApplication)
 
+cp -r\
+    0.orig\
+    0
+
 runApplication  blockMesh
 runApplication -a\
     foamDictionary  -entry numberOfSubdomains -set 4 system/decomposeParDict
 runApplication -a\
     foamDictionary  -entry method -set scotch system/decomposeParDict
-runApplication  decomposePar
+runApplication      decomposePar -copyZero
 runParallel         snappyHexMesh -overwrite
+# runParallel         patchSummary
+# runParallel         renumberMesh -noZero -overwrite
 runParallel         extrudeMesh
 runParallel         checkMesh -allTopology -allGeometry
 runParallel -a\
-    foamDictionary constant/polyMesh/boundary -entry entry0/back/type -set empty
+    foamDictionary  constant/polyMesh/boundary -entry entry0/back/type -set empty
 runParallel -a\
-    foamDictionary constant/polyMesh/boundary -entry entry0/front/type -set empty
+    foamDictionary  constant/polyMesh/boundary -entry entry0/front/type -set empty
+runApplication -a\
+    foamDictionary  constant/polyMesh/boundary -entry entry0/back/type -set empty
+runApplication -a\
+    foamDictionary  constant/polyMesh/boundary -entry entry0/front/type -set empty
+runApplication      setFields
+runParallel -a      setFields
+
 runParallel         $application
 
 #------------------------------------------------------------------------------
diff --git a/tests/eikonal/cylinders2D/system/decomposeParDict b/tests/eikonal/cylinders2D/system/decomposeParDict
index 082bab1..3d45f3c 100644
--- a/tests/eikonal/cylinders2D/system/decomposeParDict
+++ b/tests/eikonal/cylinders2D/system/decomposeParDict
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Version:  9
+    \\  /    A nd           | Version:  10
      \\/     M anipulation  |
 \*---------------------------------------------------------------------------*/
 FoamFile
diff --git a/tests/eikonal/cylinders2D/system/setFieldsDict b/tests/eikonal/cylinders2D/system/setFieldsDict
new file mode 100644
index 0000000..85eec88
--- /dev/null
+++ b/tests/eikonal/cylinders2D/system/setFieldsDict
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:                                        |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue alpha 0
+);
+
+regions
+(
+    boxToCell
+    {
+        box (-1 -1 -1) (1 1 1);
+        fieldValues
+        (
+            volScalarFieldValue alpha 1
+        );
+    }
+
+//     sphereToCell
+//     {
+//         centre (0 0 0);
+//         radius 2.5;
+//         fieldValues
+//         (
+//             volScalarFieldValue alpha.phase1 0
+//         );
+//     }
+);
+
+
+// ************************************************************************* //