3434
3535===== Mathematical model =====
3636
37+ idx{`wave1D_dn_vc.py`}
38+
3739Let $u_t$, $u_{tt}$, $u_x$, $u_{xx}$ denote derivatives of $u$ with
3840respect to the subscript, i.e., $u_{tt}$ is a second-order time
3941derivative and $u_x$ is a first-order space derivative. The
@@ -119,6 +121,8 @@ label{softeng2:wave1D:filestorage:savez}
119121
120122=== Storing individual arrays ===
121123
124+ idx{`savez`} idx{`storez`} idx{`load`} idx{zip archive}
125+
122126The `numpy.storez` function can store a set of arrays to a named
123127file in a zip archive. An associated function
124128`numpy.load` can be used to read the file later.
@@ -179,6 +183,8 @@ for array_name in array_names:
179183===== Using `joblib` to store arrays in files =====
180184label{softeng2:wave1D:filestorage:joblib}
181185
186+ idx{`joblib`} idx{memoize function} idx{hash}
187+
182188The Python package `joblib` has nice functionality for efficient storage
183189of arrays on disk. The following class applies this functionality so that
184190one can save an array, or in fact any Python data structure (e.g., a
@@ -409,14 +415,16 @@ cite{Langtangen_deqbook_wave}]).
409415
410416=== Convergence rates ===
411417
418+ idx{verification! convergence rates}
419+
412420A more general method for verification is to check the convergence rates.
413421We must introduce one discretization parameter $h$ and assume an error
414422model $E=Ch^r$, where $C$ and $r$ are constants to be determine (i.e.,
415423$r$ is the rate that we are interested in). Given two experiments with
416424different resolutions $h_i$ and $h_i{-1}$, we can estimate $r$ by
417425
418426!bt
419- \[ r = \frac{\ln(E_{i}/E_{i-1})}{\ln(h_{i}h_{i-1}},tp \]
427+ \[ r = \frac{\ln(E_{i}/E_{i-1})}{\ln(h_{i}/ h_{i-1})}, \]
420428!et
421429where $E_i$ is the error corresponding to $h_i$ and $E_{i-1}$ corresponds to
422430$h_{i-1}$. ref[Section ref{wave:pde2:fd:standing:waves}][ in cite{Langtangen_deqbook_wave}][The
@@ -427,7 +435,7 @@ for some constant $\hat c$. To compute the error, we had to rely on
427435a global variable in the user action function. Below is an implementation
428436where we have a more elegant solution in terms of a class: the `error`
429437variable is not a class attribute and there is no need for a global
430- error (which is always considered as an advantage).
438+ error (which is always considered an advantage).
431439
432440@@@CODE src-wave/wave1D/wave1D_dn_vc.py def convergence_rates@def test_convrate_sincos
433441The returned sequence `r` should converge to 2 since the error
@@ -445,6 +453,8 @@ cite{Langtangen_deqbook_wave}], see the file "`wave1D_dn_vc.py`":
445453
446454======= Programming the solver with classes =======
447455
456+ idx{`wave1D_oo.py`}
457+
448458Many who know about class programming prefer to organize their software
449459in terms of classes. This gives a richer application programming interface
450460(API) since a function solver must have all its input data in terms
@@ -506,7 +516,7 @@ documentation and a doc test should now be self-explanatory:
506516@@@CODE src-softeng2/UniformFDMesh.py fromto: import numpy@class Function
507517
508518!bnotice We rely on attribute access - not get/set functions!
509- Java programmers in particular are used to get/set functions in
519+ Java programmers, in particular, are used to get/set functions in
510520classes to access internal data. In Python, we usually apply direct
511521access of the attribute, such as `m.N[i]` if `m` is a `Mesh` object.
512522A widely used convention is to do this as long as access to
@@ -557,7 +567,7 @@ is reproduced (within machine precision).
557567======= Migrating loops to Cython =======
558568label{wave2D3D:impl:Cython}
559569
560- idx{Cython}
570+ idx{Cython} idx{`wave2D_u0.py`} idx{`wave2D_u0_adv.py`}
561571
562572We now consider the "`wave2D_u0.py`": "${src_wave}/wave2D_u0/wave2D_u0.py"
563573code for solving the 2D linear wave equation with constant wave
@@ -594,6 +604,8 @@ in one or more separate files with extension `.pyx`.
594604
595605===== Declaring variables and annotating the code =====
596606
607+ idx{`wave2D_u0_loop_cy.pyx`}
608+
597609Our starting point is the plain `advance_scalar` function for a scalar
598610implementation of the updating algorithm for new values
599611$u^{n+1}_{i,j}$:
@@ -786,14 +798,16 @@ the Cython version.
786798
787799======= Migrating loops to Fortran =======
788800
801+ idx{Fortran 77}
802+
789803Instead of relying on Cython's (excellent) ability to translate Python to C,
790804we can invoke a compiled language directly and write the loops ourselves.
791805Let us start with Fortran 77, because this is a language with more
792806convenient array handling than C (or plain C++), because
793807we can use the same multi-dimensional indices
794808in the Fortran code as in the `numpy`
795809arrays in the Python code, while in C these arrays are
796- one-dimensional and requires us to reduce multi-dimensional indices
810+ one-dimensional and require us to reduce multi-dimensional indices
797811to a single index.
798812
799813#Fortran compilers
@@ -804,6 +818,7 @@ to a single index.
804818
805819idx{wrapper code}
806820idx{Fortran subroutine}
821+ idx{`wave2D_u0_loop_f77.f`}
807822
808823We write a Fortran subroutine `advance` in a file
809824"`wave2D_u0_loop_f77.f`": "${src_wave}/wave2D_u0/wave2D_u0_loop_f77.f"
@@ -840,6 +855,8 @@ advance(u, u_n, u_nm1, f, Cx2, Cy2, dt2)
840855
841856===== Building the Fortran module with f2py =====
842857
858+ idx{Fortran 90}
859+
843860The nice feature of writing loops in Fortran is that, without
844861much effort, the tool `f2py`
845862can produce a C extension module such that
@@ -1094,10 +1111,10 @@ easier to debug.
10941111!bwarning Be careful with macro definitions
10951112Macros just perform simple text substitutions:
10961113`idx(hello,world)` is expanded to `(hello)*(Ny+1) + world`.
1097- The parenthesis in `(i)` are essential - using the natural mathematical
1114+ The parentheses in `(i)` are essential - using the natural mathematical
10981115formula `i*(Ny+1) + j` in the macro definition,
10991116`idx(i-1,j)` would expand to `i-1*(Ny+1) + j`, which is the wrong
1100- formula. Macros are handy, but requires careful use.
1117+ formula. Macros are handy, but require careful use.
11011118In C++, inline functions are safer and replace the need for macros.
11021119!ewarning
11031120
@@ -1109,6 +1126,8 @@ The C version of our function `advance` can be coded as follows.
11091126
11101127===== The Cython interface file =====
11111128
1129+ idx{`wave2D_u0_loop_c.c`} idx{`wave2D_u0_loop_c.h`}
1130+
11121131All the code above appears in a file "`wave2D_u0_loop_c.c`":
11131132"${src_wave}//wave2D_u0/wave2D_u0_loop_c.c".
11141133We need to compile this file together with C wrapper code such that
@@ -1178,6 +1197,8 @@ words, that we work with small meshes.
11781197
11791198======= Migrating loops to C via f2py =======
11801199
1200+ idx{`wave2D_u0_loop_c_f2py_signature.f`}
1201+
11811202An alternative to using Cython for interfacing C code is to apply
11821203`f2py`. The C code is the same, just the details of specifying how
11831204it is to be called from Python differ. The `f2py` tool requires
0 commit comments