Cylindrical diffusion solver with uncertainty quantification campaign

analysis.py

from easyvvuq.actions import CreateRunDirectory, Encode, Decode
from easyvvuq.actions import CleanUp, ExecuteLocal, Actions


params = {
    "max_time": {"type": "float", "default": 1.5 }, 
    "alpha": {"type": "float"}, 
    "outfile": {"type": "string", "default": "output.json"}
}

encoder = uq.encoders.GenericEncoder(
    template_fname='diffsolve.template', delimiter='$', 
    target_filename='input.json')

decoder = uq.decoders.JSONDecoder(
    target_filename='output.json', output_columns=['g1'])

execute = ExecuteLocal('{}/diffsolve.py input.json'.format(os.getcwd()))

actions = Actions(CreateRunDirectory('/tmp'), 
                  Encode(encoder), execute, Decode(decoder))

campaign = uq.Campaign(name='beam', params=params, actions=actions)


vary = {
    "alpha": cp.Uniform(0.1, 0.01),
}

campaign.set_sampler(uq.sampling.PCESampler(vary=vary, polynomial_order=1))

diffsolve.py

#!/usr/bin/env python3

import sys
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import json


def sphere_fdm(t,y,r,alpha,dr):
    """Naive FDM discretization for spherically-symmetric diffusion
       
       References:
        [1] Langtangen, H. P., & Linge, S. (2017). Finite difference computing 
            with PDEs: a modern software approach. Springer Nature. (pg. 252)
    """
    n = y.shape[0]
    dydt = np.empty_like(y)

    # Symmetry condition
    dydt[0] = 6*alpha*(y[1] - y[0])/dr**2

    for i in range(1,n-1):

        rl = 0.5*(r[i] + r[i-1])
        rr = 0.5*(r[i+1] + r[i])

        gl = rl**2 * alpha * (y[i] - y[i-1])
        gr = rr**2 * alpha * (y[i+1] - y[i])

        dydt[i] = (1./r[i])**2 * (gr - gl)/dr**2

    # Dirichlet boundary
    dydt[-1] = 0

    return dydt

def sphere_sol(t,r,alpha,nterms=40):
    """Series solution of heat equation with Dirichlet boundaries

       References:
        [1] Crank, J. (1979). The Mathematics of Diffusion. 
            Oxford University Press. (pg. 91)
    """
    a_ = 1
    C = np.zeros_like(r)
    csum = 0
    for n in range(1,nterms+1):
        expterm = np.exp(-alpha*n**2*np.pi**2*t/a_**2)
        C[1:] += (-1)**n/n * np.sin(n*np.pi*r[1:]/a_) * expterm
        csum += (-1)**n * expterm

    C[1:] *= 2*a_/(np.pi*r[1:])
    C[0] = 2*csum

    C += 1
    return 1 - C

def simulate(t,alpha,N=21):

    r, dr = np.linspace(0,1,N,retstep=True)

    method="BDF"

    y = np.ones(N)
    y[-1] = 0

    rhs = lambda t, y: sphere_fdm(t,y,r,alpha,dr)
    sol = solve_ivp(rhs,[0,t],y,method=method)

    fig, ax = plt.subplots()
    iax = ax.inset_axes([0.1,0.2,0.5,0.5])

    ax.plot(r,sol.y)
    ax.plot(r,sphere_sol(t,r,alpha),'o')

    ax.set_xlabel("r")
    ax.set_ylabel("Y(r,t)")

    iax.plot(r,sol.y)
    iax.set_xlim(0,3*dr)
    iax.set_ylim(0.98,1.02)
    iax.set_xlabel("r")

    plt.show()

    return sol.y[-1,0]


def diffsolve():

    json_input = sys.argv[1]
    with open(json_input, "r") as f:
        inputs = json.load(f)

    outfile = inputs['outfile']

    max_time = float(inputs['max_time'])
    alpha = float(inputs['alpha'])

    out = simulate(
            t=max_time,
            alpha=alpha)

    out_dict = {
        'conc': out 
        }

    with open(outfile,'w') as f:
        json.dump(out_dict,f)

if __name__ == '__main__':
  
    diffsolve()

diffsolve.template

{
    "outfile": "$outfile", 
    "max_time": $max_time, 
    "alpha": $alpha
}