This demo example shows ...



In [ ]:

    
%pylab inline

Se også notebook for den eksplisitte løsningen. Vi ønsker nå å løse varmelikningen implisitt,

$\frac{\delta u}{\delta t} = \kappa \frac{\delta^2u}{\delta x^2}$.

Diskretisert implisitt med hensyn på tid får vi

$\frac{1}{\Delta t}\left(u_i^{k} - u_i^{k-1}\right) = \frac{\kappa}{\Delta x^2} \left(u^k_{i-1} - 2u_i^k + u_{i+1}^k \right)$

og vår stensil

$u_i^{k-1} = -su^k_{i-1} +(1+2s)u_i^k - su_{i+1}^k , \quad s=\frac{\kappa\Delta t}{\Delta x^2}$.

Om vi lar ("flux-koeffisienter") $F_i = (L_i, R_i)$ være s.a. flux fra celle $i$ til celle $i+1$ er gitt ved $L_i\cdot u_i + R_i\cdot u_{i+1}$, kan vi skrive dette som

$-s\cdot F^n_{i-1} + u_i^n + s\cdot F^n_i = u_i^{n-1}$,

som illustrerer sammenhengen mellom matrisa for likningssystemet og fluksene samt verdiene $u_i^n$ for løsningen.



In [ ]:

    
import numpy as np
import time
import matplotlib.pyplot as plt
from IPython.display import clear_output

# Grid spacing
dx = 0.5

# Heat diffusion constant
kappa = 0.3

# Maximum time step size (constrained by CFL)
dt = dx*dx/(2*kappa)

# Hmm... That was the dt from the explicit case... Setting it a bit larger now:
dt = 5.0;

# Number of grid cells
n = 6

# Boundary conditions
r0 = 0.5
r1 = 1.5

s = kappa * dt / (dx*dx)

# Coefficients for fluxes, one flux for each face, i.e., n+1 of these, and two coefficients for each flux, 
# in essence a delta for face indices, resulting in cell indices... ugh.
f = np.zeros( (n+1, 2) )
for i in range(0, n+1):
    f[i] = np.array( [1, -1] )
f[0] = np.array( [nan, nan] ) # should not be used
f[n] = np.array( [nan, nan] ) # should not be used



In [ ]:

    
# Plot or not
pl = False

if (pl):
    fig, ax = plt.subplots()



In [ ]:

    
# Initial temperatures
u = np.linspace(0.0, 0.0, n)
u[0] = r0
u[n-1] = r1

A = np.zeros( (n, n) )

# Assemble the matrix one cell at a time. (This loops over interior cells, for indices 1, ..., n-2 inclusive.)
for i in range(1, n-1):
    # The equation system is: -s(u_{i-1}^n - u_i^n) + u_i^n +s(u_i^n - u_{i+1}^n) = u_i^{n-1}, for i=0, ..., n-1.
    # With f[i] = [1, -1], we get the following.
    A[i, i-1:i+1] += -s*f[i]
    A[i, i  ]     += 1
    A[i, i:i+2]   += s*f[i+1]

# Apply boundary conditions,  # u_i^n = u_i^n-1, for i=0 and i=n-1.
A[0, 0]     = 1
A[n-1, n-1] = 1



In [ ]:

    
if (not(pl)):
    print "A =\n", A
    print -1, "u =", u
    
for k in range(0, 20):
    # Perform one timestep
    
    # Solve
    v = np.linalg.solve(A, u);

    # Swap, for the next round
    u, v = v, u
    
    if (pl):
        # Plot the solution
        x = np.linspace(0.0, n-1, n) + 0.5
        plot(x, u, 'r.-', label='$u^{k+1}$')
        pylab.ylim([r0-0.5,r1+0.5])
        #time.sleep(0.2)
        clear_output()
        display(fig)
        ax.cla()
        title("Time = " + str(dt*k))
        legend()
    else:
        print k, "u =", u

if (pl):
    plt.close()

Ok, what can we imagine a DSL-program for this would look like?



In [ ]:

    
%%file 1D_heat_implicit.txt

# The grid details should be input to the DSL program
grid.type = RegularCartesian1D
grid.x0 = 0.0
grid.x1 = 1.0
grid.n = 6

# Physics that requires specification
physics.diffusion_constant = 0.3 # Heat diffusion constant. Or should this rather be an input parameter to the generated simulator?

boundary.condition_type = Dirichlet
boundary.val0 = 0.5
boundary.val1 = 1.5

initial.value = 0 # Or should this rather be an input parameter to the generated simulator?
initial.boundary = BoundaryConditions # Overriding the initial value specified above, for the boundaries

# What kind of simulator should be generated?
simulator.type = Implicit # Explicit or Implicit
simulator.basis = Constant # Temperatures averaged over cells

dt = MaxStableTimeStep # "Special" value for a variable ("Local variables" in the DSL program...)
dx = grid.getSpacing() # (x1-x0)/n for 1D grid
s = physics.diffusion_constant * dt / (dx*dx) # Think "Python name" more than "variable". The stencil-specification below makes use of this 's'...

simulator.stencil = (-s, 1+2s, -s) # -s(u_{i-1}^n - u_i^n) + u_i^n +s(u_i^n - u_{i+1}^n) = u_i^{n-1},
# Hmm. Does it make more sense not to specify this, but instead let this be generated from "flux considerations"?
# If the extremas consist of either telling the generator "solve for heat" or actually feeding it detailed C++-code, where in between should we define our sweet spot?!

simulator.generate("generated_1D_heat_implicit.cpp")

And what should the generator produce? Maybe something like the following. First a small routine for the EB3: (And no, I don't remember what EB3 stands for... EBBB = Equelle Basic Building Blocks?)



In [ ]:

    
%%file EB3_sample_linSystemSolver.cpp

#include <iostream>
#include <cmath>
#include <utility>
#include <cstring>

// Solving the n x n system a x = b, upper triangular matrix ends up in uppertri,
// b_scratch must have length n, solution ends up in x, returns true if non-singular,
// false otherwise. Matrices in row-major order. 
bool GaussElimIfNonSingular(const int n,
			    const double * const a, double * const uppertri,
			    const double * const b, double * const b_scratch,
			    double * const x)
{
    memcpy(uppertri, a, n*n*sizeof(double));
    memcpy(b_scratch, b, n*sizeof(double));
    memset(x, 0, n*sizeof(double));
    for (int j=0; j<n; j++) {
        int maxrow=j;
        for (int i=j+1; i<n; i++)
            if (fabs(uppertri[i*n+j]) > fabs(uppertri[maxrow*n+j]))
                maxrow = i;
        for (int k=j; k<n; k++)
            std::swap(uppertri[j*n+k], uppertri[maxrow*n+k]);
        std::swap(b_scratch[j], b_scratch[maxrow]);
        if (fabs(uppertri[j*n+j]) < 1e-12) // Singular?
            return false;
        for (int i=j+1; i<n; i++) {
            const double tmp = uppertri[i*n+j] / uppertri[j*n+j];
            b_scratch[i] -= b_scratch[j] * tmp;
            for (int k=n-1; k>=j; k--)
                uppertri[i*n+k] -= uppertri[j*n+k] * tmp;
        }
    }
    x[n-1] = b_scratch[n-1] / uppertri[(n-1)*n+(n-1)];
    for (int i=n-2; i>=0; i--) {
        double tmp = uppertri[i*n+i+1] * x[i+1];
        for (int k=i+2; k<n; k++)
            tmp += uppertri[i*n+k] * x[k];
        x[i] = (b_scratch[i] - tmp) / uppertri[i*n+i];
    }
    return true;
}



In [ ]:

    
!g++ -Wall -O2 -c -o EB3_sample_linSystemSolver.o EB3_sample_linSystemSolver.cpp

And then the "generated" simulator:



In [ ]:

    
%%file generated_1D_heat_implicit.cpp

// E3B-headers to be included here
bool GaussElimIfNonSingular(const int n,
                            const double * const a, double * const uppertri,
                            const double * const b, double * const b_scratch,
                            double * const x);

double u[6]; // Length from "grid.n = 6", type and dimension from "grid.type = RegularCartesian1D"
double A[6][6]; // From "simulator.type = Implicit", "grid.n = 6" and "simulator.basis = Constant"

void solution_init(void)
{
    for (int i=0; i<6; i++)
        u[i] = 0.0; // From "initial.value = 0"
    u[0] = 0.5; // From "boundary.val0 = 0.5"
    u[5] = 1.5; // From "boundary.val1 = 1.5"
}

void system_init(const double s)
{
    const int stencil_length = 3;
    const double stencil[stencil_length] = {-s, 1.0+2.0*s, -s}; // From "simulator.stencil = (-s, 1+2s, -s)"
    for (int i=0; i<6; i++) {
        for (int j=0; j<6; j++)
            A[i][j] = 0.0;
        // A[i][i] = 1.0; // got this in the stencil already
        if ( (i>0) && (i<5) )
            for (int j=0; j<3; j++)
                A[i][i-stencil_length/2+j] += stencil[j];
    }
    A[0][0] = 1.0; // From the "boundary.*" entries:
    A[5][5] = 1.0;
}

void init(void)
{
    // const double dx = 1.0/6.0; // From "dx = grid.getSpacing()"
    const double dx = 0.5; // Overriding with the number from the iPython notebook, so that we can compare the results!
    // const double dt = dx*dx/(2*0.3); // From "dt = MaxStableTimeStep"
    // Hmm... This was the max dt for the *explicit* stencil...
    const double dt = 5.0; // Trying a larger value (Note: Use the same as for the iPython notebook, so that we can compare results...)
    system_init( 0.3 * dt / (dx*dx) ); // From "physics.diffusion_constant = 0.3" and "s = physics.diffusion_constant * dt / (dx*dx)"
    solution_init();
}

void performOneTimeStep(void)
{
    double v[6], Ascratch[36], uScratch[6];

    const bool nonSingular = GaussElimIfNonSingular(6, (double *)A, Ascratch, u, uScratch, v);
    for (int i=0; i<6; i++)
        u[i] = v[i];
}



In [ ]:

    
!g++ -c -o generated_1D_heat_implicit.o generated_1D_heat_implicit.cpp

Let us test this in a small app



In [ ]:

    
%%file simtest.cpp

#include <iostream>

// Contained in the "generated" C++-code:
extern double u[6];
void init(void);
void performOneTimeStep(void);

int main(int argc, char *argv[])
{
    init();
    for (int i=0; i<20; i++) {
        std::cout << i << ": \t";
        for (int j=0; j<6; j++)
             std::cout << u[j] << " ";
        std::cout << std::endl;
        performOneTimeStep();
    }
    return 0;
}



In [ ]:

    
!g++ -o simtest simtest.cpp generated_1D_heat_implicit.o EB3_sample_linSystemSolver.o



In [ ]:

    
!simtest

Hmm... A bit strange that this converges faster than the iPython version?! Better linear solver? iPython using float instead of double internally?



In [ ]: