# Finite Volumes in Scilab

## Introduction

The goal of this page is to present the finite volume tools which are available in the Scilab environment.

## The PDE block in XCos

This block is an implementation of several numerical schemes (Finite Elements (1st and 2nd order), Finite Differences (1st and 2nd order), Finite Volumes (1st order)) to solve mono dimensional PDE (Partial Differential Equation) within SCICOS. The mathematical framwork was restricts in PDEs linear scalars with maximum order 2 in time and space. The goal is to provide engineers and physicists with an easy to use toolbox in SCICOS that will let them graphically describe the PDE to be solved. A decision system selects the most efficient numerical scheme depending on the type of the PDE and runs the solvers.

The block can manage the following discretization methods:

• Choice (check box) : is the choice for the manual or the automatic mode.
• Type : in the manual mode we can give the method type (Finite differences, Finite elements or Finite volumes).
• Degree : method degree (1 or 2 for the FD and FE methods, 1 for the FV method).
• Number of nodes : to give the number of the nodal points.

The boundary conditions are :

• Type : two types of the boundary conditions are possible : Dirichlet or Neumann.
• Expression : to give then boundary conditions expression.

The feature is based on a set of macros, within the PDE directory of the scicos_blocks module. This module manages sparse matrices.

The reference seems to be "Finite element, An introduction", Vol. 1 by E.Becker, G.Carey, and J.Oden.

The functions are:

• arbre_decision: type_meth=arbre_decision(delta): Cette fonction renvoie le type de la methode apres le parcours du fichier xml en testant le signe du discriminant
• assemb: [gk,gf]=assemb(gk,gf,ek,ef,nel,n,nodes) la fonction assemb assemble la matrice de regidite gk et le second membre gf en bouclant sur les nel- elements.
• assemblage_FE1D: [xi,w] = setint(): la fonction fournit les point d'integration x(i) et les poids pour la formule quadratique gaussienne.
• cb_IHM_EDP
• cformatlinedp
• CFORTREDP: cette fonction est pour la compilation et le link dynamique du bloc EDP avec Scilab voir CFORTR2.sci pour le bloc Scicos c_block2
• code_generation
• coef_FEM1d: Cette fonction renvoie les differentes matrices de discretisation des operateurs du/dx et dq/dx (q pour le changement de variable dans le cas de l'operateur d2u/dx2)
• Disc_diff_oper2: Cette fonction construit l'ecriture de la discritisation de l'operateur d2u/dt2 sous forme matricielle.
• Disc_diff_oper45: Cette fonction construit l'ecriture de la discritisation de l'operateur d2u/dtdx ou de l'operateur du/dx sous forme matricielle.
• Disc2_5_FVM1d: Cette fonction renvoie les differentes matrices de discretisation des operateurs du/dx et dq/dx (q pour le changement de variable dans le cas de l'operateur d2u/dx2)
• elemoper: la fonction elem evalue la matrice gk et le second memebre gf
• eval_pts_df: Cette fonction renvoie les equations de sorties correspondent aux points de mesures en approchant leus solution au noeud le plus poche
• eval_pts_EF: Cette fonction renvoie les equations de sorties correspondent aux points de mesures en utilisant l'interpolation en polynomme de Lagranges
• eval_pts_vf: Cette fonction renvoie les equations de sorties correspondent aux points de mesures en approchant leus solution par la moyenne
• extraction_infos
• formkf: la fonction formkf, construit le systeme discret de l'element finis en appelant la fonction elem pour avoir la matrices locales ek, ef et la fonction assemb pour les ajoutees aux matrices globales gk et le second membre gf.
• gen_code_FDM: Cette fonction est pour la generation des equations DAE du bloc
• gen_code_FEM: Cette fonction est pour la generation des equations DAE du bloc
• gen_code_FVM
• IHM_EDP
• lecture_xml
• maillage_FE1D: rentre les donnees des partitions (les elements)
• msprintfv
• mulf_string
• mulf3
• mulf3v
• mulfstring
• mulfv
• multMatVect
• multVectStr
• nombre_etats: cette fonction teste selon le type de discretisation differences finies et le type des conditions aux limites, elle renvoie la nouvelle taille du systemes (nombre d'etats) avec les noeuds fictifs.
• PDE: fonction graphic du bloc, elle permet le dessin et l'initialisation du bloc
• setint
• setvalue_IHM_EDP
• setvalue_IHM-PDE
• shape: evalue les valeurs des fonction de base et derivees en un point xi.
• subf_mat: fonction pour la soustraction ele/ele matrice (string) M .- N matrice(string)
• subfv
• translate: Cette fonction contient les différents algrithme de discretisation spaciale, ainsi que la génération du code du bloc EDP. Elle est appelée par la fonction graphic du bloc EDP.Sci
• unimesh1D: maillage pour les volumes finis 1D en incluant les noeuds aux limites.

## The Allaire - Coquel course

This course is devoted to hyperbolic systems of conservation laws, the most famous example of which is gaz dynamic (fully studied in details during the course). Both theoretical and numerical aspects are emphasized during each class of three hours. There are therefore a double table of contents which are followed in parallel.

• Sur la nécessité de l'hyperbolicité des systèmes de lois de conservation hyperb.sce avec le sous-programme shiftp.sci
• Divers schémas sur l'équation de transport linéaire en 1-D transp.sce avec les sous-programmes soltr.sci et shiftp.sci
• Divers schémas sur l'équation de Burgers en 1-D avec donnée initiale régulière et condition aux limites périodiques burgers0.sce avec le sous-programme shiftp.sci
• Divers schémas sur l'équation de Burgers en 1-D pour le problème de Riemann burgers.sce avec les sous-programmes solbu.sci et shiftn.sci
• Comparaison de schémas sur l'équation de Burgers en 1-D compar.sce avec le sous-programme shiftn.sci
• Convergence par raffinement de maillage de schémas sur l'équation de Burgers en 1-D rafin.sce avec le sous-programme shiftp.sci
• Résolution du problème de Riemann pour un flux non-convexe nonconv.sce avec le sous-programme shiftn.sci
• Résolution du problème de Riemann pour les équations d'Euler de la dynamique des gaz par divers schémas euler.sce avec les sous-programmes shiftn.sci flux.sci et fluvl.sci
• Quelques problèmes de Riemann pour les équations d'Euler de la dynamique des gaz riemann.sce avec le sous-programme shiftn.sci
• Schéma de Roe pour les équations d'Euler (sans et avec correction entropique) roe.sce avec le sous-programme shiftn.sci
• Divers schémas sur l'équation de transport linéaire en 2-D (comparaison de schémas vraiment 2-D avec des schémas par directions alternées) bidim.sce avec le sous-programme shiftp2d.sci

See also "Quelques programmes informatiques qui ont servi à illustrer le polycopié" (HTML).

## The Boyer codes

Un petit code de calcul 1D et 2D en scilab illustrant différents schémas de volumes finis pour les problèmes de type Darcy. Ce code a été développé avec S. Krell à des fins pédagogiques à l'occasion des journées numériques de Lille 2010.

## Rupture d'un barrage

• www.ann.jussieu.fr/~kray/Documents/RuptureBarrage_HKV.pdf

"Rupture d'un barrage", Hecketsweiler Adrien, Kray Marie, Vernier Claire, Avril 2007

On cherche à simuler un écoulement d'eau suite à la rupture d'un barrage dans le cas 1D. Il faut donc considérer la hauteur et la vitesse de l'eau. Pour cela, on utilise comme inconnues les fonctions :

```h(x, t) pour la hauteur
u(x, t) pour la vitesse```

où x est la variable d'espace et t est la variable du temps. De la même façon, on modélise la topographie du sol par la fonction a(x). On veut déterminer h et u à partir d'un système de deux équations physiques.

```//**********************************************************//
// Rupture d'un barrage : Cas de la Topographie Plate
//
//
en dimension 1
//
//**********************************************************//

// I]. Variables et paramètres :
// =========================

g=9.81;
t0=100;

// 1). Paramètres en espace :
// ----------------------
xmin=-1;
xmax=1;
N=100;               //nombre de mailles
p=(xmax-xmin)/N;     //pas en espace
x=(xmin+p*(1:N))/t0; //vecteur d'espace
B=int(N/2);          //position du barrage

// 2). Paramètres en temps :
// ---------------------
tmin=0;
tmax=40;
cfl=0.7;

// 3). Variables conservatives et flux numériques :
// --------------------------------------------
w=zeros(N,2); //w=(h,hu)
w1=zeros(N,2); //w Euler
h=zeros(N,1); //hauteur de l'eau
u=zeros(N,1); //vitesse de l'écoulement
f=zeros(N,2); //f(w)
f2=zeros(N,2); //f (i+1/2) à t=tn

// II]. Initialisation des vecteurs :
//=============================
for i=1:N;
if i<=B then
w(i,1)=2; w1(i,1)=w(i,1);
w(i,2)=0; w1(i,2)=w(i,2);
f(i,1)=w(i,2);
f(i,2)=(w(i,2)^2)/w(i,1)+g*(w(i,1)^2)/2;
else
w(i,1)=1; w1(i,1)=w(i,1);
w(i,2)=0; w1(i,2)=w(i,2);
f(i,1)=w(i,2);
f(i,2)=(w(i,2)^2)/w(i,1)+g*(w(i,1)^2)/2;
end;
end;

// III]. Résolution numérique par le schéma de Rusanov :
// ===============================================

t=tmin;
while(t < tmax )

// 1). Détermination du pas de temps :
// -------------------------------
vmax=0;
for i=2:N-1;
c=sqrt(g*w(i,1));
v1=abs(w(i,2)/w(i,1)-c);
v2=abs(w(i,2)/w(i,1)+c);

//condition de stabilité

vmax=max(vmax,v1);
vmax=max(vmax,v2);

end;

tau=cfl*p/vmax; // tau= pas de temps

// 2). Calcul du flux numérique f(i+1/2) :
// -----------------------------------
for i=1:N-1;
m1=abs(w(i,2)/w(i,1)-sqrt(g*w(i,1)));
m2=abs(w(i,2)/w(i,1)+sqrt(g*w(i,1)));
m3=abs(w(i+1,2)/w(i+1,1)-sqrt(g*w(i+1,1)));
m4=abs(w(i+1,2)/w(i+1,1)+sqrt(g*w(i+1,1)));
s=max(m1,m2,m3,m4);

f2(i,1) = 1/2*( f(i,1) + f(i+1,1) ) - s/2*( w(i+1,1) - w(i,1) );
f2(i,2) = 1/2*( f(i,2) + f(i+1,2) ) - s/2*( w(i+1,2) - w(i,2) );
end;

// 3). Schéma d'Euler :
// ----------------
for i=2:N-1;
w1(i,1)= w(i,1) -tau/p*(f2(i,1)-f2(i-1,1));
w1(i,2)= w(i,2) -tau/p*(f2(i,2)-f2(i-1,2));
end;

for i=1:N;
w(i,1)=w1(i,1);
w(i,2)=w1(i,2);
f(i,1)=w(i,2);
f(i,2)=(w(i,2)^2)/w(i,1)+g*(w(i,1)^2)/2;
h(i)=w(i,1);
u(i)=w(i,2)/w(i,1);
end;

t=t+0.05;

// 4). Tracé du résultat :
// -------------------

clf(0);
plot2d(x,h)
plot2d(x,u);
end;```

## Finite volumes in Matlab by Recktenwald

This page contains links to MATLAB codes used to demonstrate the finite difference and finite volume methods. These are simple codes applied to simple problems, hence the adjective "Toy". The toy finite volume codes can handle non-uniform meshes and non-uniform material properties. Not all sample problems exercise these features.

The codes are contained in Zip archives. The documentation is in PDF. The code were written from 2003 to 2008.

• Finite difference solutions to the one-dimensional heat equation
• One-dimensional, fully-developed flow in a round pipe
• One-dimensional upwind and central difference models of convection
• Two-dimensional finite-volume code with several test problems

The author is Gerald Recktenwald, Associate Professor and Chair Mechanical and Materials Engineering Dept, Portland State University

## Finite volumes in Matlab by Sanders

The author is Brett F. Sanders, Professor Civil and Environmental Engineering, University of California, Irvine.

On this page are a few basic codes that may be of interest to students of computational hydraulics. These codes have been developed for instructional and research purposes, and have been shared with students of my graduate courses in the past.

• 1D Wet-Bed Shallow-Water Solver : Here is a zip file containing a set of Matlab files that implement a Godunov-type finite volume scheme for solving the 1D shallow-water equations.
• 2D Wet-Bed Shallow-Water Solver : Here is a zip file containing a set of Matlab files that implement a numerical solution to the 2D shallow-water equations on a Cartesian grid. The numerical method is a first-order accurate Godunov-type finite volume scheme that utilizes Roe's approximate Riemann solver.
• 2D Wet-Bed Shallow-Water Solver Using Local Time Stepping : Here is a zip file containing a set of Matlab files very similar to the 2D solver above, except that the local time-stepping scheme described in Sanders (2008) is implemented.
• 1D Gradually Varied Flow Solver : Here is a zip file containing a pair of Matlab files that solve the steady state gradually varied flow equation using Matlab's built-in Runge Kutta solver, ode45 .
• Unstructured Grid Model for 2D Scalar Transport : Here is a zip file containing a Matlab program to solve the 2D advection equation on an unstructured grid. This program was developed to introduce students to unstructured grids, and those seeking an introduction to unstructured grids might find it worthwhile to run.
• Particle Tracking Model for 2D Diffusion : Here is a zip file containing a Matlab program to solve the 2D diffusion equation using a random-walk particle tracking method. The solution corresponds to an instantaneous load of particles at the origin at time zero.
• Particle Tracking Model for 2D Taylor Dispersion : Here is a script file taylor.m containing a Matlab program to solve the advection diffusion equation in a 2D channel flow with a parabolic velocity distribution (laminar flow). The solution corresponds to an instantaneous load of particles along an x=0 line at time zero.

In principle, it should be easy to translate them into Scilab.

## The Matlab codes by Randall Leveque

These Matlab codes are given for the book "Finite Difference Methods for Ordinary and Partial Differential Equations" :

In principle, it should be easy to translate them into Scilab.

## The CLAWPACK library

CLAWPACK is a software package designed to compute numerical solutions to hyperbolic partial differential equations using a wave propagation approach.

In principle, it should be feasible to make the link between Scilab and Python.

## The Matlab solvers by Shampine

Hyperbolic PDEs

We have developed software in MATLAB to solve initial-boundary value problems for first order systems of hyperbolic partial differential equations (PDEs) in one space variable x and time t. We allow PDEs of three general forms, viz.

• u_t = f(x,t,u,u_x)
• u_t = f(x,t,u)_x + s(x,t,u)
• u_t = f(u)_x

and we allow general boundary conditions. Solving PDEs of this generality is not routine and the success of our software is not assured. On the other hand, it is very easy to use and has performed well on a wide variety of problems. Provided here is a manuscript about this work and the solver itself. A revised version of the manuscript has appeared as Solving Hyperbolic PDEs in Matlab, Appl. Numer. Anal. & Comput. Math., 2 (2005) 346-358. Accompanying the solver are many examples.