add implicit implementation of the scheme

Aalto-CFD · Jan 17, 2022 · 993f21b · 993f21b
1 parent 3bbf1dc
commit 993f21b
Show file tree

Hide file tree

Showing 8 changed files with 178 additions and 7 deletions.
diff --git a/...nvectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.C b/...nvectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.C
@@ -43,8 +43,7 @@ template<class Type>
 const surfaceInterpolationScheme<Type>&
 gaussScaleSelectiveConvectionScheme<Type>::interpScheme() const
 {
-    NotImplemented;
-    return NullObjectRef<surfaceInterpolationScheme<Type>>();
+    return thiResInterpScheme_();
 }
 
 
@@ -83,8 +82,56 @@ gaussScaleSelectiveConvectionScheme<Type>::fvmDiv
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
-    NotImplemented;
-    return NullObjectRef<tmp<fvMatrix<Type>>>();
+    tmp<surfaceScalarField> tweights = thiResInterpScheme_().weights(vf);
+    const surfaceScalarField& weights = tweights();
+
+    tmp<fvMatrix<Type>> tfvm
+    (
+        new fvMatrix<Type>
+        (
+            vf,
+            faceFlux.dimensions()*vf.dimensions()
+        )
+    );
+    fvMatrix<Type>& fvm = tfvm.ref();
+
+    fvm.lower() = -weights.primitiveField()*faceFlux.primitiveField();
+    fvm.upper() = fvm.lower() + faceFlux.primitiveField();
+    fvm.negSumDiag();
+
+    forAll(vf.boundaryField(), patchi)
+    {
+        const fvPatchField<Type>& psf = vf.boundaryField()[patchi];
+        const fvsPatchScalarField& patchFlux = faceFlux.boundaryField()[patchi];
+        const fvsPatchScalarField& pw = weights.boundaryField()[patchi];
+
+        fvm.internalCoeffs()[patchi] = patchFlux*psf.valueInternalCoeffs(pw);
+        fvm.boundaryCoeffs()[patchi] = -patchFlux*psf.valueBoundaryCoeffs(pw);
+    }
+
+    if (thiResInterpScheme_().corrected())
+    {
+        fvm += fvc::surfaceIntegrate(faceFlux*thiResInterpScheme_().correction(vf));
+    }
+
+    GeometricField<Type, fvPatchField, volMesh> vfPrime =
+        thiPassFilter_.ref().fvcLaplacian(DxDyDz_, vf);
+
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tConvection
+    (
+        fvc::surfaceIntegrate(faceFlux*(
+            thiResInterpScheme_().interpolate(vfPrime)
+          - tloResInterpScheme_().interpolate(vfPrime)))
+    );
+
+    tConvection.ref().rename
+    (
+        "convection(" + faceFlux.name() + ',' + vf.name() + ')'
+    );
+
+    fvm += tConvection;
+
+    return tfvm;
 }
 
 

diff --git a/...nvectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.H b/...nvectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.H
@@ -109,6 +109,7 @@ public:
             tloResInterpScheme_(loResScheme),
             thiResInterpScheme_(hiResScheme),
             DxDyDz_(
+                "DxDyDz_",
                 filterScope_*filterScope_
               * Foam::pow(mesh_.surfaceInterpolation::deltaCoeffs() ,-2.0)
               / (4.0*constant::mathematical::pi*constant::mathematical::pi)

diff --git a/tutorials/convection1d/Allplot_ddt.py b/tutorials/convection1d/Allplot_ddt.py
@@ -0,0 +1,69 @@
+import os
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+
+font_size = 15
+params = {
+    'image.origin': 'lower',
+    'image.interpolation': 'nearest',
+    'image.cmap': 'gray',
+    'axes.grid': False,
+    'savefig.dpi': 500,
+    'axes.labelsize': font_size,
+    'axes.titlesize': font_size,
+    'lines.markersize': 3,
+    'lines.linewidth': 1.4,
+    'font.size': font_size,
+    'legend.fontsize': 13,
+    'xtick.labelsize': font_size,
+    'text.usetex': True,
+    'ytick.labelsize': font_size,
+    'font.family': 'serif',
+}
+
+matplotlib.rcParams.update(params)
+
+sample_path = os.path.join('singleGraph', '3', 'line_T.xy')
+cases = [
+    dict(
+        id='explicit',
+        style=dict(label= 'expl. RK4', color='k' , linestyle='solid'),
+    ),
+    dict(
+        id='implicit-Euler',
+        style=dict(label='impl. Euler', color='tab:red', linestyle='dashed'),
+    ),
+    dict(
+        id='implicit-backward',
+        style=dict(label='impl. backward', color='tab:purple'  , linestyle='dashdot'),
+    ),
+    dict(
+        id='implicit-Crank-Nicolson',
+        style=dict(label='impl. C-N 0.9', color='tab:blue'  , linestyle='dotted'),
+    ),
+
+]
+
+plt.figure()
+
+# analytical solution
+L = 2*np.pi
+x = np.linspace(-L/2, L/2, 200)
+sigma = L/20
+x0 = 0
+T = lambda x: np.exp(-((x-x0)*(x-x0))/(2*sigma*sigma))
+plt.plot(x, T(x), color='tab:gray', linestyle='dashed', label='analytical')
+
+for case in cases:
+    data_path = os.path.join('postProcessing-' + case['id'], sample_path)
+    d = np.loadtxt(data_path, unpack=True)
+    plt.plot(d[0], d[1], **case['style'])
+
+
+plt.xlabel('x, m')
+plt.ylabel('T, [-]')
+plt.xticks(ticks=[-np.pi, -0.5*np.pi, 0, 0.5*np.pi, np.pi], labels=['$-\pi$', '$-\pi/2$', '0', '$\pi$/2', '$\pi$'])
+plt.legend()
+plt.grid()
+plt.savefig('convection1d_ddt.pdf')
diff --git a/tutorials/convection1d/Allrun b/tutorials/convection1d/Allrun
@@ -29,7 +29,7 @@ mv postProcessing postProcessing-BD
 foamListTimes -rm
 
 foamDictionary -entry "divSchemes/div(phi,T)" \
-    -set "GaussScaleSelective Gauss linear orthogonal 2.0 upwind linear" \
+    -set "GaussScaleSelective Gauss linear orthogonal 2.0 upwind blended 1.0" \
     system/fvSchemes
 explicitRK4scalarTransportFoam
 postProcess -func singleGraph -latestTime

diff --git a/tutorials/convection1d/Allrun_ddt b/tutorials/convection1d/Allrun_ddt
@@ -0,0 +1,41 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+foamListTimes -rm
+rm -r postProcessing*
+rm -r 0
+cp -r 0.orig 0
+blockMesh
+
+foamDictionary -entry "divSchemes/div(phi,T)" \
+    -set "GaussScaleSelective Gauss linear orthogonal 2.0 upwind blended 1.0" \
+    system/fvSchemes
+
+explicitRK4scalarTransportFoam
+postProcess -func singleGraph -latestTime
+mv postProcessing postProcessing-explicit
+foamListTimes -rm
+
+foamDictionary -entry "ddtSchemes/default" -set "Euler" \
+    system/fvSchemes
+scalarTransportFoam
+postProcess -func singleGraph -latestTime
+mv postProcessing postProcessing-implicit-Euler
+foamListTimes -rm
+
+foamDictionary -entry "ddtSchemes/default" -set "backward" \
+    system/fvSchemes
+scalarTransportFoam
+postProcess -func singleGraph -latestTime
+mv postProcessing postProcessing-implicit-backward
+foamListTimes -rm
+
+foamDictionary -entry "ddtSchemes/default" -set "CrankNicolson 0.9" \
+    system/fvSchemes
+scalarTransportFoam
+postProcess -func singleGraph -latestTime
+mv postProcessing postProcessing-implicit-Crank-Nicolson
+foamListTimes -rm
+
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/convection1d/system/controlDict b/tutorials/convection1d/system/controlDict
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     explicitRK4scalarTransportFoam;
+application     scalarTransportFoam;
 
 startFrom       startTime;
 

diff --git a/tutorials/convection1d/system/fvSchemes b/tutorials/convection1d/system/fvSchemes
@@ -17,6 +17,7 @@ FoamFile
 
 ddtSchemes
 {
+    default         backward;
 }
 
 gradSchemes
@@ -26,7 +27,7 @@ gradSchemes
 
 divSchemes
 {
-    div(phi,T)      GaussScaleSelective Gauss linear orthogonal 2 upwind linear;
+    div(phi,T)      GaussScaleSelective Gauss linear orthogonal 2.0 upwind blended 1;
 }
 
 laplacianSchemes

diff --git a/tutorials/convection1d/system/fvSolution b/tutorials/convection1d/system/fvSolution
@@ -15,8 +15,20 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+solvers
+{
+    T
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
 SIMPLE
 {
+    nNonOrthogonalCorrectors 0;
 }
 
 // ************************************************************************* //