Tutorial example demonstrating the basic functionalities of the StabFem Software

This tutorial script is designed as a support for chapter 2 of the StabFem manual available here

We show how to compute, plot and manipulate the solution of three simple problems formulated in a L-shaped domain, namely steady and unsteady heat conduction, and Navier-Stokes flow with bottom-driven boundary conditions.

Source code for this program : SCRIPT_Lshape.m

NB in this simple tutorial the DataBase management is disabled. After studying the present example you may have a look at the related example SCRIPT_Lshape_DataBase.m

Contents

Initialization

close all;
addpath([fileparts(fileparts(pwd)) '/SOURCES_MATLAB']);

SF_Start('verbosity',3);
% higest levels of verbosity results in more messages.
% recommended values are 2 (minimal notice) and 4 (notice + freefem output). 6 and more is debug mode.

SF_DataBase('disable');
% (optional :) this is to initialize the database management.
% NB you may use 'create' instead of 'open' to clean up


mkdir('FIGURES');
set(0,'defaulttextfontsize',16); % set font size for figures
set(0,'defaultaxesfontsize',14); % idem

NOTICE  - Initialization already done
NOTICE  - Verbosity set to 3 : Notice messages
NOTICE  - Database directory :./WORK/
NOTICE  - This directory already exists
NOTICE  - Please type "SF_Status" to see what is available in this dir
NOTICE  -     or type "SF_core_arborescence('cleanall')" to erase any previous data
WARNING - Disabling database management. Some advanced functionalities may not work properly
Warning: Directory already exists.

Generation of the mesh

The program Lshape_Mesh.edp is a simple mesh generator producing a mesh.msh file. We launch this program and import the mesh as a structure using the SF_Mesh.m driver :

Ndensity =40;
ffmesh = SF_Mesh('Lshape_Mesh.edp','Options',{'meshdensity',Ndensity})

NOTICE  -       ### INITIAL MESH CREATED WITH np = 1429 vertices

ffmesh = 

  struct with fields:

             points: [3x1429 double]
             bounds: [3x160 double]
                tri: [4x2696 double]
                 np: 1429
                nbe: 160
                 nt: 2696
             labels: [1 2]
        problemtype: 'unspecified'
           filename: './mesh.msh'
           meshtype: '2D'
           datatype: 'mesh'
    DataDescription: 'Mesh (generated by generic macro)'
     generationmode: 'initial'
    datastoragemode: 'unknown'
       TriangleSize: [0x1 double]
             status: 'loaded'
                np2: 5553
          Vh_P2P2P1: [40440x1 double]
        Vh_P2P2P2P1: [56616x1 double]
              Vh_P2: [16176x1 double]
     meshgeneration: 0
           symmetry: 'N'

NB this sequence is equivalent to executing the program with FreeFem++ in a shell terminal as follows :

> FreeFem++ Lshape_Mesh.edp -meshdensity 40

(but in addition the function returns a "handle' to the mesh which can be used for subsequent computations, plotted, etc...)

How to plot the mesh

figure();SF_Plot(ffmesh,'mesh','title','Mesh for L-shape body');
saveas(gca,'FIGURES/Lshape_Mesh','png'); % to save the figure
pause(0.1);

NB the SF_Plot.m function, which is built upon the ffpdeplot.m function developped by Markus Meisinger, allows a variety of plots, including mesh plots, data plots using color levels, contour plots, or quiver plots. A large number of options are available, see the function or type "help SF_Plot".

First problem : steady conduction

Here is how to launch a simple FreeFem++ program Lshape_Steady.edp and import the data as a structure using the SF_Launch.m driver. Have a look at these programs to understand how it works !

heatS=SF_Launch('Lshape_Steady.edp','Mesh',ffmesh)

heatS = 

  struct with fields:

               mesh: [1x1 struct]
           filename: 'Data.txt'
    DataDescription: 'Temperature in a L-Shape domain'
           datatype: 'HeatSteady'
    datastoragemode: 'ReP2'
    datadescriptors: 'T'
                  T: [5553x1 double]
             status: 'loaded'

NB this sequence is equivalent to executing the program with FreeFem++ in a shell terminal as follows :

> FreeFem++ Lshape_Steady.edp

(but in addition the data is imported as an object which can be used for subsequent computations, plotted, etc...)

Note that the data is imported as a structure witch includes plottable fields (here "T", "dxT", "dyT"), numerical fields (here "HeatFlux") along as a number of string-valued descriptor fields. All this data is read from the file Data.ff2m generated by the FreeFem program ; the .ff2m format is designed for easy customization.

Plot the temperature field

figure(1);SF_Plot(heatS,'T','xlim',[0 1],'ylim',[0 1],'title', 'Steady conduction on a L-shaped domain: Temperature field ');
% add one line on the figure
hold on; plot([0 1],[.25,.25],'r--')
saveas(gca,'FIGURES/Lshape_T0','png');
hold off;

Here is how to extract data along a 1D line and plot this data

Ycut = 0.25;Xcut = 0:.01:1;
Tcut = SF_ExtractData(heatS,'T',Xcut,Ycut);

figure(2);plot(Xcut,real(Tcut),'r-',Xcut,imag(Tcut),'b--');
title('Temperature along a line : T(x,0.25)');xlabel('x');ylabel('T');
saveas(gca,'FIGURES/Lshape_T0_Cut','png')
pause(0.1);

Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

Second problem : unsteady conduction

Assuming harmonic dependance as $T(x,y,t) = Re(\hat{T}(x,y) e^{-i \omega t})$, this leads to:

$- i \omega \hat{T} = \kappa \Delta \hat{T}$ on domain, $\hat{T} = 1$ along boundary

Have a look at the program Lshape_Unsteady.edp to see how this problem is solved by FreeFem++ !

The program computes the field $\hat{T}(x,y)$ as well as the complex heat flux such as $\phi(t) = Re(\hat{\phi} e^{-i \omega t})$ ; namely : $\hat{\phi} = \int_{\partial \Omega} \nabla \hat{T} \cdot {\bf n}$.

First we solve the problem for a single value of $\omega$:

omega = 100;
heatU=SF_Launch('Lshape_Unsteady.edp','Options',{'omega',omega,'kappa',1},'Mesh',ffmesh,'DataFile','Heat_Unsteady.ff2m')

heatU = 

  struct with fields:

               mesh: [1x1 struct]
           filename: 'Heat_Unsteady.txt'
    DataDescription: 'Temperature field in a L-shaped region ; Unsteady conduction with harmonic forcing at boundaries'
           datatype: 'HeatUnsteady'
    datastoragemode: 'CxP2'
    datadescriptors: 'T'
              omega: 100
           INDEXING: [1x1 struct]
           HeatFlux: 17.1650 -21.3108i
          normGradT: [1429x1 double]
                  T: [5553x1 double]
             status: 'loaded'

NB this sequence is equivalent to executing the program with FreeFem++ in a shell terminal as follows :

> FreeFem++ Lshape_Unsteady.edp -omega 40 -kappa 1

the parameters passed as "-omega 40 -kappa 1" are detected within the program using the "getARGV" Freefem function.

We plot the results

figure(); SF_Plot(heatU,'T.re','colormap','redblue','title',{'Unsteady heat condition problem ' , 'real(colors) and imaginary(levels) parts of T_1(x)'});
hold on;  SF_Plot(heatU,'T.im','contour','on','xystyle','off');
saveas(gca,'FIGURES/Lshape_Tc','png');

Note that resulting value 'heatU' has a field 'HeatFlux' corresponding to $\hat{\phi}$

heatU.HeatFlux

ans =

  17.1650 -21.3108i

Adapts the mesh to both steady and unsteady solutions

We now use the SF_Adapt.m function to adapt the mesh to the results. This function (which is a wrapper for the FreeFem++ program AdaptMesh.edp) allows to perform mesh adaptation using up to 8 fields with various data type (scalar P1 or P2 data ; composite data for instance [P2xP2xP1] ,...) Have a look at these programs to understand how it works !

ffmesh_adapted = SF_Adapt(ffmesh,heatS,heatU)

WARNING -  When using SF_Adapd it is advised to define a database folder.  
NOTICE  - SF_Adapt : created new mesh ; np = 3331 vertices 

ffmesh_adapted = 

  struct with fields:

             points: [3x3331 double]
             bounds: [3x236 double]
                tri: [4x6424 double]
                 np: 3331
                nbe: 236
                 nt: 6424
             labels: [1 2]
        problemtype: 'unspecified'
           filename: './mesh_adapt.msh'
           meshtype: '2D'
           datatype: 'mesh'
    DataDescription: 'Mesh (generated by generic macro)'
     generationmode: 'adapted'
    datastoragemode: 'unknown'
       TriangleSize: [0x1 double]
             status: 'loaded'
                np2: 13085
          Vh_P2P2P1: [96360x1 double]
        Vh_P2P2P2P1: [134904x1 double]
              Vh_P2: [38544x1 double]
     meshgeneration: NaN
           symmetry: 'N'

NB here the first argument is the mesh handle and the return object is a handle to the new mesh. We could also have used the alternative syntax: [heatS,heatU] = SF_Adapt(heatS,heatU) which, in addition, will project back the two fields on the new mesh and return a handle to them. It is not mandatory to specify the mesh as first argument since the mesh is also a sub-object of the fields.

figure;
subplot(1,2,1); SF_Plot(ffmesh,'mesh','title','Initial mesh');
subplot(1,2,2); SF_Plot(ffmesh_adapted,'mesh','title','Adapted mesh');

recomputes the solution on new mesh and plot the sum of both temperature fields

heatS=SF_Launch('Lshape_Steady.edp','Mesh',ffmesh_adapted);
heatU=SF_Launch('Lshape_Unsteady.edp','Options',{'omega',omega},'Mesh',ffmesh_adapted,'DataFile','Heat_Unsteady.ff2m');
heatSum = SF_Add({heatS,heatU},'coefs',[1 1])

WARNING - in SFcore_AddMESHFilenameToFF2M : File Addition.ff2m exists in root directory ! this may lead to incorrect behavior
WARNING -  When using SF_Add it is advised to define a database folder.  
NOTICE  -  in SF_Add : processsing data from txt files :  T

heatSum = 

  struct with fields:

               mesh: [1x1 struct]
           filename: 'Addition.txt'
    DataDescription: '(no description provided when creating file)'
           datatype: 'Addition'
    datastoragemode: 'ReP2'
    datadescriptors: 'T'
       meshfilename: './mesh_adapt.msh'
                  T: [13085x1 double]
             status: 'loaded'

plot the results

figure();SF_Plot(heatSum,'T','colormap','redblue','title','Temperature (sum of two problems)');
saveas(gca,'FIGURES/Lshape_Tc','png');

Third problem : Loop over a range of omega

The next program Lshape_Impedance.edp solves the unsteady heat conduction problem over a range of $\omega$ and returns the corresponding values of $\hat{\phi}$.

res = SF_Launch('Lshape_Impedance.edp','Mesh',ffmesh_adapted,'Options',{'kappa',1},'DataFile','Heat_Unsteady_loop.ff2m')

res = 

  struct with fields:

               mesh: [1x1 struct]
           filename: 'Heat_Unsteady_loop.txt'
    DataDescription: 'Solution of the unsteady thermal problem over a range of omega'
           datatype: 'EXAMPLE'
    datastoragemode: 'columns'
    datadescriptors: 'omega,Flux_r,Flux_i'
              kappa: 1
           INDEXING: [1x1 struct]
              omega: [50x1 double]
               Flux: [50x1 double]
             status: 'loaded'

We plot the 'thermal impedance' defined as $Z_{\theta} = 1/\hat{\phi}$. We compare with the asymptotic solutions in the low-frequency and high-frequency limits.

figure;loglog(res.omega,real(1./res.Flux),'r','LineWidth',2);
hold on;loglog(res.omega,imag(1./res.Flux),'b','LineWidth',2);
rangeomegaLOW = logspace(-1,0,10);
loglog(rangeomegaLOW,2./rangeomegaLOW,'k--','LineWidth',2);
rangeomegaHIGH = logspace(2,3,10);
loglog(rangeomegaHIGH,.25.*rangeomegaHIGH.^(-.5),'k:','LineWidth',2);
xlabel('\omega');ylabel('Re(Z_\theta), Im(Z_\theta)');
title('Thermal impedance Z_\theta(\omega)')
legend('$Re(Z_\theta$)', '$Im(Z_\theta)$','Low-$\omega$ asymptotics','High-$\omega$ asymptotics','Interpreter','latex');
saveas(gca,'FIGURES/Z_omega','png');

Fourth problem : Navier-Stokes flow

Eventually we call a third program Lshape_NavierStokes.edp which computes the flow in the L-shape domain (with bottom-driven boundary conditions). by solving the Navier-Stokes Flow using Newton iteration.

Ndensity =30;
ffmesh = SF_Mesh('Lshape_Mesh.edp','Options',{'meshdensity',Ndensity});
viscosity = 1;Umax = 1;
Flow = SF_Launch('Lshape_NavierStokes.edp','Mesh',ffmesh,'Options',{'Re',10})

NOTICE  -       ### INITIAL MESH CREATED WITH np = 826 vertices

Flow = 

  struct with fields:

               mesh: [1x1 struct]
           filename: 'Data.txt'
    DataDescription: 'Navier-Stokes flow in a L-shaped bottom-driven cavity'
           datatype: 'StokesFlow'
    datastoragemode: 'ReP2P2P1.1'
    datadescriptors: 'ux,uy,p,Re'
                 Re: 10
           INDEXING: [1x1 struct]
                 Fx: -1.7808
          vorticity: [826x1 double]
                 ux: [9180x1 double]
                 uy: [9180x1 double]
                  p: [4590x1 double]
             status: 'loaded'

Remarks : - Here one option is specified. The FreeFem program is launched as follows: > FreeFem++ Lshape_NavierStokes.edpc -Re 10

- In this case the data is a vectorial dataset (ux,uy,p) defined using (P2,P2,P1) mixed elements by FreeFem. Note that the SF_Launch driver is able to recognize this and import the fields (ux,uy,p) along with the vorticity.

We adapt the mesh and recompute, this time for Re=20

ffmesh = SF_Adapt(ffmesh,Flow);
Flow = SF_Launch('Lshape_NavierStokes.edp','Mesh',ffmesh,'Options',{'Re',20})

WARNING -  When using SF_Adapd it is advised to define a database folder.  
NOTICE  - SF_Adapt : created new mesh ; np = 1933 vertices 

Flow = 

  struct with fields:

               mesh: [1x1 struct]
           filename: 'Data.txt'
    DataDescription: 'Navier-Stokes flow in a L-shaped bottom-driven cavity'
           datatype: 'StokesFlow'
    datastoragemode: 'ReP2P2P1.1'
    datadescriptors: 'ux,uy,p,Re'
                 Re: 20
           INDEXING: [1x1 struct]
                 Fx: -1.3219
          vorticity: [1933x1 double]
                 ux: [21984x1 double]
                 uy: [21984x1 double]
                  p: [10992x1 double]
             status: 'loaded'

We plot the vorticity fields with isocontours; We then add the velocity field using quiver plot

figure;SF_Plot(Flow,'vorticity','colorrange',[-20 20]);
hold on;SF_Plot(Flow,{'ux','uy'},'ylim',[0 1],'xlim',[0 1],'title','Navier-Stokes flow (vorticity and velocity fields)');
saveas(gca,'FIGURES/StokesFlow','png');

Final remark :

To understand what the FreeFem++ programs actually do behind the interface, you may have a look to the file ".stabfem_log.bash" generated at the end of this script, which regroups all the calls to FreeFem programs and file manipulations !

% [[PUBLISH]] -> this tag is here to tell the git runner to automatically generate a web page from this script


%
fclose(fopen('SCRIPT_Lshape_NoDB.success','w+'));


Published with MATLAB® R2022b