Low-Mach flow in a heated cavity (Rayleigh number sweep)

Part 1, Simulation Setup

Results published: Gutiérrez-Jorquera, Kummer: A fully coupled high-order discontinuous Galerkin method for diffusion flames in a low-Mach number framework

Note

This example can be found in the source code repository as as HeatedCavity_RaSweep.ipynb. One can directly load this into Jupyter to interactively work with the following code examples.

Note: First, BoSSS has to be loaded into the Jupyter kernel. Note: In the following line, the reference to BoSSSpad.dll is required. One must either set #r "BoSSSpad.dll" to something which is appropirate for the current computer (e.g. C:\Program Files (x86)\FDY\BoSSS\bin\Release\net5.0\BoSSSpad.dll if working with the binary distribution), or, if one is working with the source code, one must compile BoSSSpad and put it side-by-side to this worksheet file (from the original location in the repository, one can use the scripts getbossspad.sh, resp. getbossspad.bat).

// #r "C:\BoSSS2\experimental\public\src\L4-application\BoSSSpad\bin\Release\net5.0\bossspad.dll"
#r "BoSSSpad.dll"
using System;
using System.Collections.Generic;
using System.Linq;
using System.IO;
using System.Data;
using System.Globalization;
using System.Threading;
using ilPSP;
using ilPSP.Utils;
using BoSSS.Platform;
using BoSSS.Foundation;
using BoSSS.Foundation.Grid;
using BoSSS.Foundation.Grid.Classic;
using BoSSS.Foundation.IO;
using BoSSS.Solution;
using BoSSS.Solution.Control;
using BoSSS.Solution.GridImport;
using BoSSS.Solution.Statistic;
using BoSSS.Solution.Utils;
using BoSSS.Solution.Gnuplot;
using BoSSS.Application.BoSSSpad;
using BoSSS.Application.XNSE_Solver;
using static BoSSS.Application.BoSSSpad.BoSSSshell;
using BoSSS.Foundation.Grid.RefElements;
using BoSSS.Platform.LinAlg;
using BoSSS.Solution.NSECommon;
using BoSSS.Application.XNSEC;
Init();

The below script needs to be able to find the current output cell; this is an easy method to get it.

BoSSSshell.WorkflowMgm.Init("HeatedCavity_RayleighSweepStudy");

Project name is set to 'HeatedCavity_RayleighSweepStudy'.

Default Execution queue is chosen for the database.

Opening existing database '\\fdygitrunner\ValidationTests\databases\HeatedCavity_RayleighSweepStudy'.

Batch System to use:

GetDefaultQueue();

static var myDb = BoSSSshell.WorkflowMgm.DefaultDatabase;// = myDb;

BoSSSshell.WorkflowMgm.SetNameBasedSessionJobControlCorrelation()

Study of the heat cavity configuration for increasing Rayleigh number.

A comparison with the benchmark results from Vierendeels (DOI 10.1108/09615530310501957) for Ra = [1e2,1e3,1e4,1e5,1e6,1e7] is done

Create grid

int[] Resolutions = new int[]{5};
int[] DGdegree = new int[]{ 3 };
int[] nCells  =  Resolutions.Select(r => (int)(Math.Pow(2,r+1))).ToArray();
double[] Rayleighs = new double[] {1e2,1e3,1e4,1e5,1e6,1e7};

public static class GridFactory {
    public static double[] GetXNodes(int Res) { 
        var xNodes = GenericBlas.SinLinSpacing(-0.5, 0.5, 0.0, Res + 1);
        return xNodes;
    }
 
    static double[] GetYNodes(int Res) {
        double[] yNodes =  GenericBlas.SinLinSpacing(-0.5, 0.5,0.0, Res + 1);
        return yNodes;
    }
 
    public static Grid2D GenerateGrid(int Res) { 
        var xNodes = GetXNodes(Res);
        var yNodes = GetYNodes(Res);
        var grd    = Grid2D.Cartesian2DGrid(xNodes, yNodes);
        grd.EdgeTagNames.Add(1, "NoSlipNeumann");
        grd.EdgeTagNames.Add(2, "wall_tempfixed_left");
        grd.EdgeTagNames.Add(3, "wall_tempfixed_right");
        grd.DefineEdgeTags( delegate (double[] X) {
            double x = X[0];
            double y = X[1];

            //Edge tags
            //1: Adiabatic no slip wall
            //2: Temperature fixed no slip wall

            //right cold wall 
            if (Math.Abs(x - 0.5) < 1e-8)
                return 3;

            //bottom adiabatic Wall
            if (Math.Abs(y - 0.5 ) < 1e-8)
                return 1;

            // left hot wall
            if (Math.Abs(x + 0.5) < 1e-8)
                return 2;

            //top adiabatic Wall
            if (Math.Abs(y + 0.5 ) < 1e-8)
                return 1;

            else throw new ArgumentOutOfRangeException();
        });
        
        myDb.SaveGrid(ref grd);
        
        return grd;
     }
 
 }

public static class BoundaryValueFactory { 

    public static string GetPrefixCode(double Th, double Tc, double Froude) {
        using(var stw = new System.IO.StringWriter()) {
           
           stw.WriteLine("static class BoundaryValues {");
           stw.WriteLine("  static public double VelX(double[] X) {");
           stw.WriteLine("    return 0.0;");
           stw.WriteLine("  }");

           stw.WriteLine("  static public double VelY(double[] X) {");
           stw.WriteLine("    return 0.0;");
           stw.WriteLine("  }");

           stw.WriteLine("  static public double TemperatureHot(double[] X) {");
           stw.WriteLine("    return 1.6;");
           stw.WriteLine("  }");

           stw.WriteLine("  static public double TemperatureCold(double[] X) {");
           stw.WriteLine("    return 0.4;");
           stw.WriteLine("  }");

           stw.WriteLine("  static public double One(double[] X) {");
           stw.WriteLine("    return  1.0;");
           stw.WriteLine("  }");

           stw.WriteLine("  static public double Zero(double[] X) {");
           stw.WriteLine("    return  0.0;");
           stw.WriteLine("  }");

           stw.WriteLine(" static public double InitialPressure(double[] X) { ");
            stw.WriteLine("         return (-1)* X[1] / ("+Froude * Froude +") ;");
            stw.WriteLine("    }");
            stw.WriteLine("}"); 
           return stw.ToString();
        }
    }
   
    static public Formula Get_VelX(double Th, double Tc , double Froude) {
        return new Formula("BoundaryValues.VelX", AdditionalPrefixCode:GetPrefixCode(Th,Tc,Froude));
    }

    static public Formula Get_VelY(double Th, double Tc, double Froude){
        return new Formula("BoundaryValues.VelY", AdditionalPrefixCode:GetPrefixCode(Th,Tc,Froude));
    }
    static public Formula Get_TemperatureHot(double Th, double Tc, double Froude){
        return new Formula("BoundaryValues.TemperatureHot", AdditionalPrefixCode:GetPrefixCode(Th,Tc,Froude));
    }

     static public Formula Get_TemperatureCold(double Th, double Tc, double Froude){
        return new Formula("BoundaryValues.TemperatureCold", AdditionalPrefixCode:GetPrefixCode(Th,Tc,Froude));
    }

    static public Formula Get_One(double Th, double Tc, double Froude){
        return new Formula("BoundaryValues.One", AdditionalPrefixCode:GetPrefixCode(Th,Tc,Froude));
    }
    static public Formula Get_Zero(double Th, double Tc, double Froude){
        return new Formula("BoundaryValues.Zero", AdditionalPrefixCode:GetPrefixCode(Th,Tc,Froude));
    }
    static public Formula Get_InitialPressure(double Th, double Tc, double Froude){
        return new Formula("BoundaryValues.InitialPressure", AdditionalPrefixCode:GetPrefixCode(Th,Tc,Froude));
    }
}

Send and run jobs

var controls = new List<XNSEC_Control>();

var controls = new List<BoSSS.Application.XNSEC.XNSEC_Control>();
double Th = 1.6; double Tc = 0.4;
foreach(double Ra in Rayleighs){
foreach(int dg in DGdegree){
foreach(int Res in Resolutions) {

    var C = new BoSSS.Application.XNSEC.XNSEC_Control();
    C.SetDGdegree(dg);
    var nCells = (int) Math.Pow(2,Res+1);
    C.SetGrid(GridFactory.GenerateGrid(nCells));
    C.Paramstudy_CaseIdentification.Add(new Tuple<string, object>("Res", Res));
    C.Paramstudy_CaseIdentification.Add(new Tuple<string, object>("Dgdegree", dg));
    C.Paramstudy_CaseIdentification.Add(new Tuple<string, object>("Rayleigh", Ra));
    C.SessionName = "NaturalConvection_k" + Res + "_DG" + dg+ "_Ra"+Ra;

    C.EnableMassFractions = false;
    C.NumberOfChemicalSpecies = 1;
    C.ChemicalReactionActive = false;
    C.MatParamsMode = MaterialParamsMode.Sutherland;
    C.physicsMode = PhysicsMode.Combustion;
    C.TimesteppingMode = AppControl._TimesteppingMode.Steady;
    C.LinearSolver = LinearSolverCode.direct_pardiso.GetConfig();
    C.NonLinearSolver.verbose = true;
    C.NonLinearSolver.MaxSolverIterations = 100;

    C.PenaltyViscMomentum = 1.0 * 1;
    C.PenaltyHeatConduction = 1.0 * 1;
    C.PhysicalParameters.IncludeConvection = true;

    C.UseSelfMadeTemporalOperator = false;
    C.timeDerivativeEnergyp0_OK = false;
    C.timeDerivativeConti_OK = false;

 
    C.EdgeTagsNusselt = new string[] { "wall_tempfixed_left", "wall_tempfixed_right", "NoSlipNeumann" };

    C.Rayleigh = Ra;
    C.Reynolds = Math.Sqrt(Ra);
    C.Prandtl = 0.71;
    double Fr =Math.Sqrt(2 * C.Prandtl * (1.6 - 0.4) / (1.6 + 0.4)); 
    C.Froude = Fr;
    C.HeatCapacityRatio = 1.4;
    C.T_ref_Sutherland = 600;
    C.ThermodynamicPressureMode = ThermodynamicPressureMode.MassDetermined; // Because its a closed system, i.e. p0 = p0(time)
    C.PhysicalParameters.IncludeConvection = true;
    C.Timestepper_LevelSetHandling         = BoSSS.Solution.XdgTimestepping.LevelSetHandling.None;
    
     
  
    if(Ra > 1e5){ // For Rayleigh numbers greater than 1e5 the newton-dogleg algorithm doesnt find a solution within ~ 70 iterations => Use homotopy
        C.HomotopyVariable = XNSEC_Control.HomotopyVariableEnum.Reynolds;
        C.homotopieAimedValue = Math.Sqrt(Ra);
        C.StartingHomotopyValue = Math.Sqrt(1e5); // Suficiently easy to find solution
        C.HomotopyApproach = XNSEC_Control.HomotopyType.Automatic;
        C.NonLinearSolver.HomotopyStepLongFail = 100;
    }
    


    C.AddBoundaryValue("NoSlipNeumann", VariableNames.VelocityX,BoundaryValueFactory.Get_VelX(Th, Tc,Fr));
    C.AddBoundaryValue("NoSlipNeumann", VariableNames.VelocityY,BoundaryValueFactory.Get_VelY(Th, Tc,Fr));

    C.AddBoundaryValue("wall_tempfixed_left", VariableNames.Temperature,BoundaryValueFactory.Get_TemperatureHot(Th, Tc,Fr));
    C.AddBoundaryValue("wall_tempfixed_right", VariableNames.Temperature, BoundaryValueFactory.Get_TemperatureCold(Th, Tc,Fr));

    C.AddBoundaryValue("wall_tempfixed_left", VariableNames.MassFraction0, BoundaryValueFactory.Get_One(Th, Tc,Fr));
    C.AddBoundaryValue("wall_tempfixed_right", VariableNames.MassFraction0, BoundaryValueFactory.Get_One(Th, Tc,Fr));


    C.AddInitialValue(VariableNames.VelocityX,  BoundaryValueFactory.Get_Zero(Th, Tc,Fr));
    C.AddInitialValue(VariableNames.VelocityY,BoundaryValueFactory.Get_Zero(Th, Tc,Fr));
    C.AddInitialValue(VariableNames.Pressure,BoundaryValueFactory.Get_InitialPressure(Th, Tc,Fr));

    C.AddInitialValue(VariableNames.Temperature,BoundaryValueFactory.Get_One(Th, Tc,Fr));
    C.AddInitialValue(VariableNames.MassFraction0, BoundaryValueFactory.Get_One(Th, Tc,Fr));
    C.AddInitialValue(VariableNames.ThermodynamicPressure,BoundaryValueFactory.Get_One(Th, Tc,Fr));

    
    // Note (fk 09oct23): The following two settings seem to influence the value of the thermodynamic pressure.
    // The reason is not clear.
    // In recent merge requests, the default values for these two properties changed from false to true,
    // causing the follow-up nutebook (`HeatedCavity_RaSweepPostProc.ipynb`) to fail when checking the 
    // value of the thermodynamic pressure. All other checked properties seem un-affected.
    // This should be invesigated closer if any future works build upon this solver.
    C.DynamicLoadBalancing_On = false;
    C.DynamicLoadBalancing_RedistributeAtStartup = false;
    
    

    controls.Add(C);
}
}
}

Grid Edge Tags changed.

An equivalent grid (801ecb15-ca90-4b50-b1c7-7268059ebd23) is already present in the database -- the grid will not be saved.

Grid Edge Tags changed.

An equivalent grid (801ecb15-ca90-4b50-b1c7-7268059ebd23) is already present in the database -- the grid will not be saved.

Grid Edge Tags changed.

An equivalent grid (801ecb15-ca90-4b50-b1c7-7268059ebd23) is already present in the database -- the grid will not be saved.

Grid Edge Tags changed.

An equivalent grid (801ecb15-ca90-4b50-b1c7-7268059ebd23) is already present in the database -- the grid will not be saved.

Grid Edge Tags changed.

An equivalent grid (801ecb15-ca90-4b50-b1c7-7268059ebd23) is already present in the database -- the grid will not be saved.

Grid Edge Tags changed.

An equivalent grid (801ecb15-ca90-4b50-b1c7-7268059ebd23) is already present in the database -- the grid will not be saved.

Run Simulations

BoSSSshell.ExecutionQueues.ForEach(q => Console.WriteLine(q))

MS HPC client  MSHPC-Gitrunner-HighPrio @DC2, @\\fdygitrunner\BoSSStests

MS HPC client  MSHPC-Gitrunner-DefaultTest @DC2, @\\fdygitrunner\BoSSStests

MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy

MS HPC client  MSHPC-AllNodes-test @DC2, @\\fdygitrunner\BoSSStests

MS HPC client  MSHPC-FastNodes @DC2, @\\fdygitrunner\ValidationTests\deploy

//var myBatch = BoSSSshell.ExecutionQueues[0];

//myBatch.AllowedDatabasesPaths.Add(new AllowedDatabasesPair(myDb.Path,""));

Type solver = typeof(BoSSS.Application.XNSEC.XNSEC);
foreach(var c in controls) {
    string jobName       = c.SessionName;    
    var oneJob           = new Job(jobName, solver);
    oneJob.NumberOfMPIProcs = 6;
    oneJob.SetControlObject(c);
    Console.WriteLine("Executing " + jobName   );
    oneJob.Activate(); 
}

Executing NaturalConvection_k5_DG3_Ra100

Deployments so far (1): (Job token: 820429, FinishedSuccessful 'HeatedCavity_RayleighSweepStudy-XNSEC2023Nov09_131724.596787' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "HeatedCavity_RayleighSweepStudy	NaturalConvection_k5_DG3_Ra100	11/09/2023 13:17:51	082ec5d1..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

Executing NaturalConvection_k5_DG3_Ra1000

Deployments so far (1): (Job token: 820431, FinishedSuccessful 'HeatedCavity_RayleighSweepStudy-XNSEC2023Nov09_131741.409230' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "HeatedCavity_RayleighSweepStudy	NaturalConvection_k5_DG3_Ra1000	11/09/2023 13:18:08	cf5599f9..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

Executing NaturalConvection_k5_DG3_Ra10000

Deployments so far (1): (Job token: 820432, FinishedSuccessful 'HeatedCavity_RayleighSweepStudy-XNSEC2023Nov09_131800.303207' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "HeatedCavity_RayleighSweepStudy	NaturalConvection_k5_DG3_Ra10000	11/09/2023 13:18:27	b7690393..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

Executing NaturalConvection_k5_DG3_Ra100000

Deployments so far (1): (Job token: 820434, FinishedSuccessful 'HeatedCavity_RayleighSweepStudy-XNSEC2023Nov09_131819.452861' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "HeatedCavity_RayleighSweepStudy	NaturalConvection_k5_DG3_Ra100000	11/09/2023 13:18:47	543c6f35..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

Executing NaturalConvection_k5_DG3_Ra1000000

Deployments so far (1): (Job token: 820436, FinishedSuccessful 'HeatedCavity_RayleighSweepStudy-XNSEC2023Nov09_131838.650298' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "HeatedCavity_RayleighSweepStudy	NaturalConvection_k5_DG3_Ra1000000	11/09/2023 13:19:06	9ceb8fb7..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

Executing NaturalConvection_k5_DG3_Ra10000000

Deployments so far (1): (Job token: 820437, FinishedSuccessful 'HeatedCavity_RayleighSweepStudy-XNSEC2023Nov09_131857.832509' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "HeatedCavity_RayleighSweepStudy	NaturalConvection_k5_DG3_Ra10000000	11/09/2023 13:19:26	a6cc631b..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

// // wait for all jobs to finish (up to 1 day, check every 1 minutes)
BoSSSshell.WorkflowMgm.BlockUntilAllJobsTerminate(TimeOutSeconds:(3600*24*1), PollingIntervallSeconds:(60*1));

All jobs finished.

Test for failed Jobs:

// detect failed Jobs in the job management
var suspects = BoSSSshell.WorkflowMgm.AllJobs.Select(kv => kv.Value)
    .Where(job => job.LatestSession.Tags.Contains(SessionInfo.NOT_TERMINATED_TAG)
                  || job.LatestSession.Tags.Contains(SessionInfo.SOLVER_ERROR)).ToArray();
suspects

NUnit.Framework.Assert.IsTrue(suspects.Count() <= 0, $"{suspects.Count()} Failed Jobs of {BoSSSshell.WorkflowMgm.AllJobs.Count()} in total.");