Rheological Flow around Confined Cylinder (Part 2, Postprocessing)

Results published: Kikker, Kummer, Oberlack: A fully coupled high-order discontinuous Galerkin solver for viscoelastic fluid flow, see also at tubiblio.

#r "BoSSSpad.dll"
using System;
using System.Collections.Generic;
using System.Linq;
using ilPSP;
using ilPSP.Utils;
using BoSSS.Platform;
using BoSSS.Platform.LinAlg;
using BoSSS.Foundation;
using BoSSS.Foundation.XDG;
using BoSSS.Foundation.Grid;
using BoSSS.Foundation.Grid.Classic;
using BoSSS.Foundation.Grid.RefElements;
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.AdvancedSolvers;
using BoSSS.Solution.Gnuplot;
using BoSSS.Application.BoSSSpad;
using BoSSS.Application.XNSE_Solver;
using BoSSS.Application.XNSFE_Solver;
using static BoSSS.Application.BoSSSpad.BoSSSshell;
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("ConfinedCylinder_ConvergenceStudy");

Project name is set to 'ConfinedCylinder_ConvergenceStudy'.

Default Execution queue is chosen for the database.

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

OpenOrCreateDatabase(@"\\fdygitrunner\ValidationTests\ConfinedCylinder_ConvergenceStudy");

Opening existing database '\\fdygitrunner\ValidationTests\ConfinedCylinder_ConvergenceStudy'.

databases

#0: { Session Count = 32; Grid Count = 4; Path = \\fdygitrunner\ValidationTests\databases\ConfinedCylinder_ConvergenceStudy }

#1: { Session Count = 0; Grid Count = 0; Path = \\fdygitrunner\ValidationTests\ConfinedCylinder_ConvergenceStudy }

using BoSSS.Application.Rheology;

// Remark: we need to instantiate one object out of the rheology solver assembly,
// so that the de-seialization in in below work.
var dummy = new RheologyControl();

Error Convergence Plot for Velocity

var ExpPlotS_velX = new List<Plot2Ddata>();
int[] pOrder = new int[] {1, 2, 3, 4};

foreach (int p in pOrder){ // loop over polynomial degrees
    // extract sessions with DG degree pDeg
    //Console.WriteLine("order " + p + "...");
    var pDegSessions = BoSSSshell.WorkflowMgm.Sessions.Where(
    // function which is true on all sessions we are interested in:
    Si => Convert.ToInt32(Si.KeysAndQueries["DGdegree:VelocityX"]) == p 
                           && Convert.ToDouble(Si.KeysAndQueries["Weissenberg"]) == 0.2
                           && Si.SuccessfulTermination).ToArray();
    //Console.WriteLine("got the sessions.");
    // now , create a plot from the selected sessions :
    Plot2Ddata pDegPlot = pDegSessions.ToEstimatedGridConvergenceData("VelocityX", 
      xAxis_Is_hOrDof:false, // false selects DOFs for x- axis
      normType:NormType.L2_approximate);
    Plot2Ddata pDegPlot_Y = pDegSessions.ToEstimatedGridConvergenceData("VelocityY", 
      xAxis_Is_hOrDof:false, // false selects DOFs for x- axis
      normType:NormType.L2_approximate);
    //Console.WriteLine("got the plots.");
    // Combine VelocityX and VelocityX-plot
    int N = 3;
    if(pDegPlot.dataGroups[0].Values.Length != N)
        throw new Exception();
    if(pDegPlot_Y.dataGroups[0].Values.Length != N)
        throw new Exception();
    
    for(int n = 0; n < N; n++) {
        double vx = pDegPlot.dataGroups[0].Values[n];
        double vy = pDegPlot_Y.dataGroups[0].Values[n];
        pDegPlot.dataGroups[0].Values[n] = Math.Sqrt(vx*vx + vy*vy);
    }  

    // remember the freshly created plot object in a list :
    ExpPlotS_velX.Add(pDegPlot);
    //Console.WriteLine("done with p = " + p);
}

ExpPlotS_velX[0].dataGroups[0].Format.PointType = PointTypes.Diamond ;
ExpPlotS_velX[1].dataGroups[0].Format.PointType = PointTypes.Box;
ExpPlotS_velX[2].dataGroups[0].Format.PointType = PointTypes.LowerTriangle ;
ExpPlotS_velX[3].dataGroups[0].Format.PointType = PointTypes.Asterisk ;
ExpPlotS_velX[0].dataGroups[0].Name = "k = 1";
ExpPlotS_velX[1].dataGroups[0].Name = "k = 2";
ExpPlotS_velX[2].dataGroups[0].Name = "k = 3";
ExpPlotS_velX[3].dataGroups[0].Name = "k = 4";
foreach (var p in ExpPlotS_velX ){
  p.dataGroups[0].Format.PointSize = 0.8;
  p.dataGroups[0].Format.DashType = DashTypes.Dotted ;
  p.dataGroups[0].Format.LineWidth = 2;
}

var ExpPlot_velX = ExpPlotS_velX[0]; // select 0-th object
foreach (var p in ExpPlotS_velX.Skip (1)) { // loop over other ( skip 0-th entry )
  ExpPlot_velX = ExpPlot_velX.Merge(p); // merge
}

ExpPlot_velX.Regression()

index	Key	Value
0	k = 1	-2.197317841967774
1	k = 2	-2.5758854097204806
2	k = 3	-3.6434120989046925
3	k = 4	-4.629626924528952

var velX = ExpPlot_velX.ToGnuplot().PlotCairolatex(xSize:14,ySize:12);
velX.SaveTo(@"vel_Convergence.tex");
ExpPlot_velX.PlotNow()

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!

This example is also part of the BoSSS testing framework, therefore it containes NUnit assertions to check the convergence slopes automatically:

using NUnit.Framework;

var reg = ExpPlot_velX.Regression();
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("1")).Value, -1.5, "Velocity convergence for k=1 failed");
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("2")).Value, -2.5, "Velocity convergence for k=2 failed");
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("3")).Value, -3.5, "Velocity convergence for k=3 failed");
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("4")).Value, -4.4, "Velocity convergence for k=4 failed");

Error Convergence Plot for Pressure

var ExpPlotS_pres = new List<Plot2Ddata>();
int[] pOrder2 = new int[] {0, 1, 2, 3};

foreach (int p in pOrder2){ // loop over polynomial degrees
    // extract sessions with DG degree pDeg
    var pDegSessions = BoSSSshell.WorkflowMgm.Sessions.Where(
    // function which is true on all sessions we are interested in:
    Si => Convert.ToInt32(Si.KeysAndQueries["DGdegree:Pressure"]) == p 
                           && Convert.ToDouble(Si.KeysAndQueries["Weissenberg"]) == 0.2
                           && Si.SuccessfulTermination).ToArray();

    // now , create a plot from the selected sessions :
    // (we could also do other things )
    Plot2Ddata pDegPlot = pDegSessions.ToEstimatedGridConvergenceData("Pressure",
      xAxis_Is_hOrDof:false, // false selects DOFs for x- axis
      normType:NormType.L2noMean_approximate);

    // remember the freshly created plot object in a list :
    ExpPlotS_pres.Add(pDegPlot);
}

ExpPlotS_pres[0].dataGroups[0].Format.PointType = PointTypes.Diamond ;
ExpPlotS_pres[1].dataGroups[0].Format.PointType = PointTypes.Box;
ExpPlotS_pres[2].dataGroups[0].Format.PointType = PointTypes.LowerTriangle ;
ExpPlotS_pres[3].dataGroups[0].Format.PointType = PointTypes.Asterisk ;
ExpPlotS_pres[0].dataGroups[0].Name = "k = 0";
ExpPlotS_pres[1].dataGroups[0].Name = "k = 1";
ExpPlotS_pres[2].dataGroups[0].Name = "k = 2";
ExpPlotS_pres[3].dataGroups[0].Name = "k = 3";
foreach (var p in ExpPlotS_pres ){
  p.dataGroups[0].Format.PointSize = 0.8;
  p.dataGroups[0].Format.DashType = DashTypes.Dotted ;
  p.dataGroups[0].Format.LineWidth = 2;
}

var ExpPlot_pres = ExpPlotS_pres[0]; // select 0-th object
foreach (var p in ExpPlotS_pres.Skip (1)) { // loop over other ( skip 0-th entry )
  ExpPlot_pres = ExpPlot_pres.Merge(p); // merge
}

ExpPlot_pres.Regression()

index	Key	Value
0	k = 0	-1.7571296430673014
1	k = 1	-2.26823681628724
2	k = 2	-2.7884756067494796
3	k = 3	-3.922533834210022

var pres = ExpPlot_pres.ToGnuplot().PlotCairolatex(xSize:14,ySize:12);
pres.SaveTo(@"pres_Convergence.tex");
ExpPlot_pres.PlotNow()

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!

This example is also part of the BoSSS testing framework, therefore it containes assertions to check the convergence slopes automatically:

var reg = ExpPlot_pres.Regression();
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("0")).Value, -1, "Pressure convergence for k=1 failed");
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("1")).Value, -2, "Pressure convergence for k=2 failed");
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("2")).Value, -2.5, "Pressure convergence for k=3 failed");
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("3")).Value, -3.5, "Pressure convergence for k=4 failed");

Error Convergence Plot for Stress

var ExpPlotS_strXX = new List<Plot2Ddata>();

foreach (int p in pOrder){ // loop over polynomial degrees
    // extract sessions with DG degree pDeg
    var pDegSessions = BoSSSshell.WorkflowMgm.Sessions.Where(
    // function which is true on all sessions we are interested in:
    Si => Convert.ToInt32(Si.KeysAndQueries["DGdegree:StressXX"]) == p 
                           && Convert.ToDouble(Si.KeysAndQueries["Weissenberg"]) == 0.2
                           && Si.SuccessfulTermination).ToArray();

    // now , create a plot from the selected sessions :
    // (we could also do other things )

    Plot2Ddata pDegPlot = pDegSessions.ToEstimatedGridConvergenceData("StressXX", xAxis_Is_hOrDof:false, // false selects DOFs for x- axis
      normType:NormType.L2_approximate);
    Plot2Ddata pDegPlot_XY = pDegSessions.ToEstimatedGridConvergenceData("StressXY", xAxis_Is_hOrDof:false, // false selects DOFs for x- axis
      normType:NormType.L2_approximate);
    Plot2Ddata pDegPlot_YY = pDegSessions.ToEstimatedGridConvergenceData("StressYY", xAxis_Is_hOrDof:false, // false selects DOFs for x- axis
      normType:NormType.L2_approximate);
      
    // Combine XX-, XY- and YY-Stress-plot
    int N = 3;
    if(pDegPlot.dataGroups[0].Values.Length != N)
        throw new Exception();
    if(pDegPlot_XY.dataGroups[0].Values.Length != N)
        throw new Exception();
    if(pDegPlot_YY.dataGroups[0].Values.Length != N)
        throw new Exception();    
    for(int n = 0; n < N; n++) {
        double sxx = pDegPlot.dataGroups[0].Values[n];
        double sxy = pDegPlot_XY.dataGroups[0].Values[n];
        double syy = pDegPlot_YY.dataGroups[0].Values[n];
        pDegPlot.dataGroups[0].Values[n] = Math.Sqrt(sxx*sxx + sxy*sxy + syy*syy);
    }    

    // remember the freshly created plot object in a list :
    ExpPlotS_strXX.Add(pDegPlot);
}

ExpPlotS_strXX[0].dataGroups[0].Format.PointType = PointTypes.Diamond ;
ExpPlotS_strXX[1].dataGroups[0].Format.PointType = PointTypes.Box;
ExpPlotS_strXX[2].dataGroups[0].Format.PointType = PointTypes.LowerTriangle ;
ExpPlotS_strXX[3].dataGroups[0].Format.PointType = PointTypes.Asterisk ;
ExpPlotS_strXX[0].dataGroups[0].Name = "k = 1";
ExpPlotS_strXX[1].dataGroups[0].Name = "k = 2";
ExpPlotS_strXX[2].dataGroups[0].Name = "k = 3";
ExpPlotS_strXX[3].dataGroups[0].Name = "k = 4";
foreach (var p in ExpPlotS_strXX ){
  p.dataGroups[0].Format.PointSize = 0.8;
  p.dataGroups[0].Format.DashType = DashTypes.Dotted ;
  p.dataGroups[0].Format.LineWidth = 2;
}

var ExpPlot_strXX = ExpPlotS_strXX[0]; // select 0-th object
foreach (var p in ExpPlotS_strXX.Skip (1)) { // loop over other ( skip 0-th entry )
  ExpPlot_strXX = ExpPlot_strXX.Merge(p); // merge
}

ExpPlot_strXX.Regression()

index	Key	Value
0	k = 1	-1.2753908296690686
1	k = 2	-1.8707295548812675
2	k = 3	-2.8416847752385626
3	k = 4	-4.138830919122578

var strXX = ExpPlot_strXX.ToGnuplot().PlotCairolatex(xSize:14,ySize:12);
strXX.SaveTo(@"str_Convergence.tex");
ExpPlot_strXX.PlotNow()

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!

This example is also part of the BoSSS testing framework, therefore it containes assertions to check the convergence slopes automatically:

var reg = ExpPlot_strXX.Regression();
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("1")).Value, -1, "Stress convergence for k=1 failed");
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("2")).Value, -1.5, "Stress convergence for k=2 failed");
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("3")).Value, -2.5, "Stress convergence for k=3 failed");
Assert.LessOrEqual(reg.Single(tt => tt.Key.Contains("4")).Value, -3.5, "Stress convergence for k=4 failed");