Skip to content

Commit aacc1d6

Browse files
author
Lesur Geoffroy
committed
add a Python setup example
1 parent ce0a55a commit aacc1d6

9 files changed

Lines changed: 305 additions & 4 deletions

File tree

‎PythonSetup/problem/definitions.hpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
#define COMPONENTS 1
2+
#define DIMENSIONS 1
3+
#define ISOTHERMAL
4+
5+
#define GEOMETRY CARTESIAN

‎PythonSetup/problem/idefix.ini‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
[Grid]
2+
X1-grid 1 -1.0 1000 u 1.0
3+
4+
[TimeIntegrator]
5+
CFL 0.8
6+
tstop 10.0
7+
first_dt 1.e-4
8+
nstages 2
9+
10+
[Hydro]
11+
solver hll
12+
csiso constant 1.0
13+
14+
#[Python]
15+
#...
16+
#
17+
18+
[Setup]
19+
forcing_frequency 25.0
20+
forcing_amplitude 1e-3
21+
22+
[Boundary]
23+
X1-beg outflow
24+
X1-end outflow
25+
26+
[Output]
27+
#python 1.0 # we call the python output_function every 1.0
28+
vtk 1.0

‎PythonSetup/problem/myscript.py‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
from pydefix import *
2+
import numpy as np
3+
import matplotlib.pyplot as plt
4+
import inifix
5+
6+
# The output function
7+
# the only argument is dataBlockHost python object, wrapping a dataBlockHost Idefix object
8+
def output(data,grid,n):
9+
# GatherIdefixArray is provided by pydefix.
10+
# It combines all MPI subdomain into one array (not a problem here, as the full domain is reasonably small)
11+
# If MPI is not enabled, GatherIdefixArray merely does a local copy
12+
vx = GatherIdefixArray(data.Vc[VX1,:,:,:],data,broadcast=False,keepBoundaries=False)
13+
14+
# fetch frequency from idefix input file
15+
conf = inifix.load("idefix.ini")
16+
forcing_frequency = conf["Setup"]["forcing_frequency"]
17+
18+
# only process #0 performs the output
19+
if prank==0:
20+
vxf = np.fft.fft(vx[0,0,:]) # fourier transform of vx
21+
freqs = np.fft.fftfreq(grid.np_int[IDIR],grid.dx[IDIR][0]) # array of spatial frequencies
22+
23+
# TODO: find the peak of the fft
24+
#lambda_peak = ..
25+
#cs = ...
26+
#print("Pydefix(output): @t=%f, peak wavelength=%e, measured sound speed=%e"%(data.t,lambda_peak,cs))
27+
28+
29+
30+
31+
def initflow(data):
32+
# Initialize the flow
33+
data.Vc[RHO,:,:,:] = 1.0
34+
data.Vc[VX1,:,:,:] = 0
35+

‎PythonSetup/problem/setup.cpp‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
#include "idefix.hpp"
2+
#include "setup.hpp"
3+
4+
namespace SetupVariables
5+
{
6+
real forcingFrequency;
7+
real forcingAmplitude;
8+
}
9+
10+
// Add a sinusoidal forcing in the middle of the domain
11+
void MyForcingTerm(Hydro *hydro, const real t, const real dt) {
12+
const real forcing_amplitude = 1e-2;
13+
const real forcing_frequency = 25.0;
14+
15+
// Local copy of needed arrays
16+
IdefixArray4D<real> Uc = hydro->Uc; // Array of conservative variables to be updated
17+
IdefixArray1D<real> x1 = hydro->data->x[IDIR];
18+
19+
const real fAmplitude = SetupVariables::forcingAmplitude;
20+
const real fFrequency = SetupVariables::forcingFrequency;
21+
idefix_for("MySourceTerm",
22+
0, hydro->data->np_tot[KDIR],
23+
0, hydro->data->np_tot[JDIR],
24+
0, hydro->data->np_tot[IDIR],
25+
KOKKOS_LAMBDA (int k, int j, int i) {
26+
if(std::fabs(x1(i))<1e-2) {
27+
Uc(MX1,k,j,i) += fAmplitude*std::cos(2.*M_PI*fFrequency*t)*dt;
28+
}
29+
});
30+
}
31+
32+
// Default constructor
33+
// Can be used to allocate
34+
// Arrays or variables which are used later on
35+
Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) {
36+
// Tell idefix we want to add the source term "MyForcingTerm"
37+
data.hydro->EnrollUserSourceTerm(&MyForcingTerm);
38+
39+
// Fetch user-def variables
40+
SetupVariables::forcingFrequency = input.Get<real>("Setup","forcing_frequency",0);
41+
SetupVariables::forcingAmplitude = input.Get<real>("Setup","forcing_amplitude",0);
42+
}
43+
44+
// This routine initialize the flow
45+
// Note that data is on the device.
46+
// One can therefore define locally
47+
// a datahost and sync it, if needed
48+
void Setup::InitFlow(DataBlock &data) {
49+
// This is not called!
50+
IDEFIX_ERROR("You are using the Setup::Initflow method to initialise the flow. Try to use pydefix instead.");
51+
}

‎PythonSetup/solution/definitions.hpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
#define COMPONENTS 1
2+
#define DIMENSIONS 1
3+
#define ISOTHERMAL
4+
5+
#define GEOMETRY CARTESIAN

‎PythonSetup/solution/idefix.ini‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
[Grid]
2+
X1-grid 1 -1.0 1000 u 1.0
3+
4+
[TimeIntegrator]
5+
CFL 0.8
6+
tstop 10.0
7+
first_dt 1.e-4
8+
nstages 2
9+
10+
[Hydro]
11+
solver hll
12+
csiso constant 1.0
13+
14+
[Python]
15+
script myscript # always ommit the ".py"!
16+
output_function output
17+
initflow_function initflow
18+
19+
[Setup]
20+
forcing_frequency 25.0
21+
forcing_amplitude 1e-3
22+
23+
[Boundary]
24+
X1-beg outflow
25+
X1-end outflow
26+
27+
[Output]
28+
python 1.0 # we call the python output_function every 0.2
29+
vtk 1.0

‎PythonSetup/solution/myscript.py‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
from pydefix import *
2+
import numpy as np
3+
import matplotlib.pyplot as plt
4+
import inifix
5+
6+
# The output function
7+
# the only argument is dataBlockHost python object, wrapping a dataBlockHost Idefix object
8+
def output(data,grid,n):
9+
# GatherIdefixArray is provided by pydefix.
10+
# It combines all MPI subdomain into one array (not a problem here, as the full domain is reasonably small)
11+
# If MPI is not enabled, GatherIdefixArray merely does a local copy
12+
vx = GatherIdefixArray(data.Vc[VX1,:,:,:],data,broadcast=False,keepBoundaries=False)
13+
14+
# fetch frequency from idefix input file
15+
conf = inifix.load("idefix.ini")
16+
forcing_frequency = conf["Setup"]["forcing_frequency"]
17+
18+
# only process #0 performs the output
19+
if prank==0:
20+
vxf = np.fft.fft(vx[0,0,:]) # fourier transform of vx
21+
freqs = np.fft.fftfreq(grid.np_int[IDIR],grid.dx[IDIR][0]) # array of spatial frequencies
22+
23+
# find the peak of the fft
24+
ipeak = np.argmax(np.abs(vxf))
25+
lambda_peak = 1/np.fabs(freqs[ipeak])
26+
cs = forcing_frequency*lambda_peak
27+
print("Pydefix(output): @t=%f, peak wavelength=%e, measured sound speed=%e"%(data.t,lambda_peak,cs))
28+
29+
30+
31+
32+
def initflow(data):
33+
# Initialize the flow
34+
data.Vc[RHO,:,:,:] = 1.0
35+
data.Vc[VX1,:,:,:] = 0
36+

‎PythonSetup/solution/setup.cpp‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
#include "idefix.hpp"
2+
#include "setup.hpp"
3+
4+
namespace SetupVariables
5+
{
6+
real forcingFrequency;
7+
real forcingAmplitude;
8+
}
9+
10+
// Add a sinusoidal forcing in the middle of the domain
11+
void MyForcingTerm(Hydro *hydro, const real t, const real dt) {
12+
const real forcing_amplitude = 1e-2;
13+
const real forcing_frequency = 25.0;
14+
15+
// Local copy of needed arrays
16+
IdefixArray4D<real> Uc = hydro->Uc; // Array of conservative variables to be updated
17+
IdefixArray1D<real> x1 = hydro->data->x[IDIR];
18+
19+
const real fAmplitude = SetupVariables::forcingAmplitude;
20+
const real fFrequency = SetupVariables::forcingFrequency;
21+
idefix_for("MySourceTerm",
22+
0, hydro->data->np_tot[KDIR],
23+
0, hydro->data->np_tot[JDIR],
24+
0, hydro->data->np_tot[IDIR],
25+
KOKKOS_LAMBDA (int k, int j, int i) {
26+
if(std::fabs(x1(i))<1e-2) {
27+
Uc(MX1,k,j,i) += fAmplitude*std::cos(2.*M_PI*fFrequency*t)*dt;
28+
}
29+
});
30+
}
31+
32+
// Default constructor
33+
// Can be used to allocate
34+
// Arrays or variables which are used later on
35+
Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) {
36+
// Tell idefix we want to add the source term "MyForcingTerm"
37+
data.hydro->EnrollUserSourceTerm(&MyForcingTerm);
38+
39+
// Fetch user-def variables
40+
SetupVariables::forcingFrequency = input.Get<real>("Setup","forcing_frequency",0);
41+
SetupVariables::forcingAmplitude = input.Get<real>("Setup","forcing_amplitude",0);
42+
}
43+
44+
// This routine initialize the flow
45+
// Note that data is on the device.
46+
// One can therefore define locally
47+
// a datahost and sync it, if needed
48+
void Setup::InitFlow(DataBlock &data) {
49+
// This is not called!
50+
IDEFIX_ERROR("You are using the Setup::Initflow method to initialise the flow. Try to use pydefix instead.");
51+
}

‎start.ipynb‎

Lines changed: 65 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,9 @@
5050
{
5151
"cell_type": "markdown",
5252
"id": "2221ea27-c37b-41e4-b203-a716c742897d",
53-
"metadata": {},
53+
"metadata": {
54+
"jp-MarkdownHeadingCollapsed": true
55+
},
5456
"source": [
5557
"# 2- Basics: configuration, compilation, run on CPUs\n",
5658
"<a id=\"compilation\"></a>\n",
@@ -382,19 +384,78 @@
382384
"to see the effect on compressibility. You can also try to use [parallelism with MPI](#mpi) to speed up the computation\n"
383385
]
384386
},
387+
{
388+
"cell_type": "markdown",
389+
"id": "c8659a27-6da2-4bdf-85bb-a58e54940253",
390+
"metadata": {},
391+
"source": [
392+
"# 5- a python-based setup\n",
393+
"\n",
394+
"In this example, we're going to use `pydefix`, an Idefix module designed to allow Idefix to talk directly to a python interpreter, on-the-fly, without needing disk access.\n",
395+
"\n",
396+
"## The problem\n",
397+
"\n",
398+
"The problem we consider is a 1D isothermal problem. The flow is initially at rest, and we have already coded a forcing in the middle of the domain which forces a sinusoidal oscillation on VX1. This forcing is going to lead to a stationary sound wave with a sinusoidal wave pattern. Our goal is to measure the wavelength of this sound wave.\n",
399+
"\n",
400+
"To do this, we are going to fourier-transform the spatial domain on the fly. To avoid rewriting an fft routine, we will simply use the fft library from numpy.\n",
401+
"\n",
402+
"We first move to the problem directory\n",
403+
"\n",
404+
"```shell\n",
405+
"cd <path_to_idefix_tutorial>/PythonSetup/problem/\n",
406+
"```\n",
407+
"\n",
408+
"## Your tasks\n",
409+
"\n",
410+
"1 - configure the code (with cmake) to enable Python support and recompile (feel free to choose a CPU or a GPU target). The [pydefix documentation](https://idefix.readthedocs.io/latest/modules/pydefix.html) will probably be handy.\n",
411+
"\n",
412+
"2 - if you try to run the code straight away, you will see an error asking you to experiment with pydefix. This is because while idefix has been compiled with python support, pydefix has not yet been enabled. This is done in the `Python` block of idefix.ini. Check the [pydefix documentation](https://idefix.readthedocs.io/latest/modules/pydefix.html) and add a `Python` block in idefix.ini to use the script \"myscript\" and the function \"initflow\" (already defined) to make your initial conditions (constant density, no velocity). Have a look at the script \"myscript.py\" to see how this is done.\n",
413+
"\n",
414+
"3- if you now run the code, you can check that the code effectively generates a sound wave from the domain center. Let's open the vtk file series and plot VX1 to see that:"
415+
]
416+
},
417+
{
418+
"cell_type": "code",
419+
"execution_count": null,
420+
"id": "473343ff-eaa5-4e52-976b-895e3665a216",
421+
"metadata": {},
422+
"outputs": [],
423+
"source": [
424+
"# Load the VTK file series produced by Idefix\n",
425+
"plt.figure()\n",
426+
"for n in range(10):\n",
427+
" V=readVTK(\"./PythonSetup/solution/data.%04d.vtk\"%n)\n",
428+
" plt.plot(V.x,V.data[\"VX1\"][:,0,0])\n",
429+
"plt.show()"
430+
]
431+
},
432+
{
433+
"cell_type": "markdown",
434+
"id": "af4f6947-41ec-409f-afe5-4dfc61f57429",
435+
"metadata": {},
436+
"source": [
437+
"4- now the meat: let's make an output function in python. We have already defined an output function in \"myscript.py\", so let's tell idefix that we to call this output function in the `Python` block of idefix.ini (check the doc, again)\n",
438+
"\n",
439+
"5- in the `Output` block, we need to enable `Python` output and specify an output period. Let's use the same output period as the vtk output: 1.0\n",
440+
"\n",
441+
"6- now you need to finish coding the output function in python. Most of it is already done: you just need to find the location of the maximum of the spectrum, get its associated spatial frequency, and compute the associated sound speed (keep in mind $c_s=f \\lambda$)\n",
442+
"\n",
443+
"7- if you did your job well, the output function should now print directly the measured wavelength and sound speed. You can now vary the forcing frequency, amplitude, and sound speed directly in idefix.ini and check how this works!"
444+
]
445+
},
385446
{
386447
"cell_type": "markdown",
387448
"id": "16d63439-77e0-4c7b-852a-f12abd705834",
388449
"metadata": {
389450
"tags": []
390451
},
391452
"source": [
392-
"# 5- A more advanced setup\n",
453+
"# 6- A more advanced setup\n",
393454
"\n",
394455
"## Before we start\n",
395456
"\n",
396457
"\n",
397-
"In this tutorial we will introduce several important aspects hidden in the Simple Setup tutorial: Host and Device memory space, the `idefix_loop` construct and the tricks associated with it.\n",
458+
"In this part of the tutorial we will introduce several important aspects hidden in the Simple Setup tutorial: Host and Device memory space, the `idefix_loop` construct and the tricks associated with it.\n",
398459
"\n",
399460
"This tutorial is not intended to duplicate Idefix documentation. It is strongly recommended to read the introduction in the programming guide regarding [Host and device](https://idefix.readthedocs.io/latest/programmingguide.html#host-and-device), [Arrays](https://idefix.readthedocs.io/latest/programmingguide.html#arrays) and [Loops](https://idefix.readthedocs.io/latest/programmingguide.html#execution-space-and-loops).\n",
400461
"\n",
@@ -588,7 +649,7 @@
588649
"id": "e0db2593-d0c6-4779-9dfd-5da4ba886d75",
589650
"metadata": {},
590651
"source": [
591-
"# 6- Programming in Idefix: Debugging and profiling\n",
652+
"# 7- Programming in Idefix: Debugging and profiling\n",
592653
"\n",
593654
"## Problem1: a CPU segmentation fault\n",
594655
"\n",

0 commit comments

Comments
 (0)