From 993f21ba2c78cbd19655976e7b3f2818c341e980 Mon Sep 17 00:00:00 2001 From: Morev Ilya Date: Mon, 17 Jan 2022 04:21:15 +0200 Subject: [PATCH] add implicit implementation of the scheme --- .../gaussScaleSelectiveConvectionScheme.C | 55 +++++++++++++-- .../gaussScaleSelectiveConvectionScheme.H | 1 + tutorials/convection1d/Allplot_ddt.py | 69 +++++++++++++++++++ tutorials/convection1d/Allrun | 2 +- tutorials/convection1d/Allrun_ddt | 41 +++++++++++ tutorials/convection1d/system/controlDict | 2 +- tutorials/convection1d/system/fvSchemes | 3 +- tutorials/convection1d/system/fvSolution | 12 ++++ 8 files changed, 178 insertions(+), 7 deletions(-) create mode 100644 tutorials/convection1d/Allplot_ddt.py create mode 100755 tutorials/convection1d/Allrun_ddt diff --git a/src/finiteVolume/finiteVolume/convectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.C b/src/finiteVolume/finiteVolume/convectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.C index b9bdc94..1da395b 100644 --- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.C +++ b/src/finiteVolume/finiteVolume/convectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.C @@ -43,8 +43,7 @@ template const surfaceInterpolationScheme& gaussScaleSelectiveConvectionScheme::interpScheme() const { - NotImplemented; - return NullObjectRef>(); + return thiResInterpScheme_(); } @@ -83,8 +82,56 @@ gaussScaleSelectiveConvectionScheme::fvmDiv const GeometricField& vf ) const { - NotImplemented; - return NullObjectRef>>(); + tmp tweights = thiResInterpScheme_().weights(vf); + const surfaceScalarField& weights = tweights(); + + tmp> tfvm + ( + new fvMatrix + ( + vf, + faceFlux.dimensions()*vf.dimensions() + ) + ); + fvMatrix& fvm = tfvm.ref(); + + fvm.lower() = -weights.primitiveField()*faceFlux.primitiveField(); + fvm.upper() = fvm.lower() + faceFlux.primitiveField(); + fvm.negSumDiag(); + + forAll(vf.boundaryField(), patchi) + { + const fvPatchField& 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 vfPrime = + thiPassFilter_.ref().fvcLaplacian(DxDyDz_, vf); + + tmp> 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/src/finiteVolume/finiteVolume/convectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.H b/src/finiteVolume/finiteVolume/convectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.H index b7f65f5..441145f 100644 --- a/src/finiteVolume/finiteVolume/convectionSchemes/gaussScaleSelectiveConvectionScheme/gaussScaleSelectiveConvectionScheme.H +++ b/src/finiteVolume/finiteVolume/convectionSchemes/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 new file mode 100644 index 0000000..9ae2380 --- /dev/null +++ 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 index 8cd6ddc..779851c 100755 --- 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 new file mode 100755 index 0000000..144c6b4 --- /dev/null +++ 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 index 0f836e0..c2d5b16 100644 --- 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 index 83a5798..29698a1 100644 --- 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 index 4969bdf..800cf97 100644 --- 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; } // ************************************************************************* //