No Convergence with PCE while MC converges

[hamzabuw]

Hi Jonathan,

I am trying to compute standard deviation of the output Mass loss Rate (MLR) from a CFD Pyrolysis Model. I am observing that my PCE doesn't converge even though Monte Carlo(MC) has already converged. Moreover the convergence behaviour below shows that it's quite far off for PCE and MC.

I used Sobol samples for three input parameters to the CFD model. The model includes an exponential term as well to define reaction behaviour based on the temperature.

m is mass here, T is temperature, Ru is universal gas constant, E and A are reaction parameters. The temperature is evaluated by energy and mass balance.

Could you kindly give your opinion as to why that could be so?

Best regards, Hamza

[jonathf]

I don't know. You mind making a minimal code example that illustrate how you used chaospy to estimate the standard deviations?

[hamzabuw]

num_samples = 200
num_cells = 1500
E2_actual = 142792.40627
rho_actual = 1236.5102342
h_actual =  813931.17434
E2 = chaospy.Uniform(0.8*E2_actual, 1.2*E2_actual)
rho = chaospy.Uniform(0.8*rho_actual, 1.2*rho_actual)
h = chaospy.Uniform(0.8*h_actual, 1.2*h_actual)
joint = chaospy.J(E2, rho, h)
samples = joint.sample(num_samples, rule="sobol")
expansion = chaospy.generate_expansion(6, joint,normed=True)

# ------------------------------------------------------------------------------
dataRoot = "../"
from numpy import genfromtxt
for cell_ID in range(num_cells):
    samples_cell = samples
    for timeStep in range(50,451,50):
        print('Cell %d at timestep %d' % (cell_ID, timeStep))
        data=dataRoot + "ddt_rho_all_csv_"+str(timeStep)+".csv"
        df_evals = pd.read_csv(data,header=None)
        eval_all_data = genfromtxt(data, delimiter=',')
        evals = eval_all_data[:,cell_ID]
        # cell volume and surfaceArea is constant below
        evals = -evals*cell_vol/ cell_surfaceArea * 1000.0 
        if debug:
            print(df_evals)
            print(df_evals.columns)
            print(evals)
        # ------------------------------------------------------------------------------
        model_approx = chaospy.fit_regression(expansion, samples_cell, evals)
        # ------------------------------------------------------------------------------
        E = chaospy.E(model_approx, joint)
        std = chaospy.Std(model_approx, joint)
        var = chaospy.Var(model_approx, joint)
        loggerPCE.info("%d,%d,%f,%f,%f" % ((timeStep,cell_ID,E,std,var)))

The CSV files contain simulated data output for each time step for the samples above.

[jonathf]

I can't be sure, since I don't have the full picture, but I see a potential problem with h is almost 100 times as big as rho. I highly recommend switching to generalized polynomial chaos. Use

joint = cp.J(cp.Uniform(-1, 1), cp.Uniform(-1, 1), cp.Uniform(-1, 1))

as your R distribution. Don't normalize.

[hamzabuw]

Thank you for the input. I did implement Rosenblatt Transformation as you suggested. This does improve convergence of PCE. However the PCE convergence is still slower than that of Monte Carlo (MC) as shown below. Is that an expected behaviour in your opinion?

I am currently running for higher PCE orders as well.

[jonathf]

Not all problems are a good fit for PCE. The fact the your expansions converges quickly with samples, but not with order, might indicate that your problem is an example of that unfortunatly.

[hamzabuw]

Thank you for your feedback