Can using the FiPy library result in a different solution even if the equation and variables are the same?

[0pago]

I'm working on the electron distribution with the FiPy library. I will simply our problem for brevity. The governing equation has a form of the convective-diffusive equation with time-dependent coefficients.

Our distribution solver contains four main parts : (a) set up grids and initialize variables, (b) construct equation with eq = (TransientTerm == DiffusionTerm + PowerLawConvectionTerm), (c) update diffusion / convection coefficients and (d) solve equation.

I experienced that the solutions from two distribution solver differ from each other even though both have the same distribution function and convection / diffusion coefficients. I tested this by the following way.

1. Initialize the ordinary solver instance.

re = DistributionSolver( ... )

2. Initialize the distribution function (f is distribution function here)

re.f.setValue( ... )

3. Update plasma parameters (this means updating convective / diffusive coefficients)

re.update_parameters( ... )

4. Evolve the distribution

re.eq.solve( ... )

5. Initialize another solver instance and initialize the distribution function with the ordinary distribution function

re2 = DistributionSolver( ... )

re2.f.setValue(re.f.setValue)

6. Update plasma parameters in both

re.update_parameters( ... )

re2.update_parameters( ... )

7. Evolve the distribution in both

re.eq.solve( ... )

re2.eq.solve( ... )

8. Repeat 5-7.

If I compare distribution function in re and re2 at step 8), there is discrepancy. I want to understand the reason.

I tried the followings.

A) the distribution function, f, is the same before eq.solve.
B) convective / diffusive coefficients are the same before eq.solve.
C) our default solver is scipy.LinearLUSolver.
D) I tested the above-mentioned way using the hasOld attribute. But it doesn't resolve this issue.

I'm still confused. Is there any clue for discrepancy?

[guyer]

Can you show some real code? Too much of the pseudo-code doesn't make any sense.

[0pago]

Dear guyer,

Thank you for your reply. I will add more detail after simplifying the code and checking this issue.

[0pago]

Hello! This is my code simplified. Please check the below and I appreciate your help in advance.

import numpy as np
import fipy as fp

class DistributionSolver:
    def __init__(self, p_size=100, c_size=50, p_min=0., p_max=80., c_q=1.0, p_q=1.0):
        from scipy.special import erf
        #Generate mesh
        dp1 = (p_max-p_min) * (1-p_q)/(1-p_q**p_size)
        p_grid = p_min+np.cumsum([0]+[dp1 * p_q ** (n - 1) for n in range(1, p_size + 1)]); dp = np.diff(p_grid)
        c_grid = np.linspace(-1,1,c_size,endpoint=False);     dcos = c_grid[1]-c_grid[0]
        self.mesh = fp.Grid2D(dx=dp, dy=dcos, nx=p_size, ny=c_size) + [[p_grid[0]], [-1.0]]

        #Deffine sought-for function
        self.f = fp.CellVariable(name="Distribution function F",
                                 mesh=self.mesh,
                                 value=0.)

        #Deffine variable parameter
        self.E = fp.Variable()
        self.E.setValue(0.)

        self.Te = fp.Variable()
        self.Te.setValue(1000/511000)

        self.n_D = fp.Variable()
        self.n_D.setValue(0.)

        #Initialize convective/diffusive coefficients
        p        = self.mesh.x
        cos      = self.mesh.y
        self.p     = p; self.cos   = cos
        X_G        = p / fp.numerix.sqrt(2.0 * self.Te * (1.0 + p ** 2))
        self.G     = (erf(X_G) - X_G * 2.0 / np.sqrt(np.pi) * fp.numerix.exp(-X_G**2)) / (2.0 * X_G ** 2)

        self.convCoeff_p   = [-1.0, 0.0] * (\
                          self.E * cos\
                        - (self.G / self.Te - 2.0 * self.G * fp.numerix.sqrt(1.0 + p ** 2) / p ** 2))
        self.convCoeff_cos = [0.0, -1.0] * (self.E * (1.0 - cos ** 2) / p)

        self.Diff_cos      = 0.5 * (1. + (erf(X_G) - self.G + self.Te * p ** 2/(1.0 + p ** 2))) * fp.numerix.sqrt(1.0 + p ** 2) * (1.0 - cos ** 2) / p ** 3
        self.D_cos         = [[[0.0, 0.0], [0.0, 1.0]] * self.Diff_cos]

        self.Diff_p        = self.G * fp.numerix.sqrt(1.0 + p ** 2) / p
        self.D_p           = [[[1.0, 0.0], [0.0, 0.0]] * self.Diff_p]

        # The equation
        self.eq = (fp.TransientTerm() == (fp.DiffusionTerm(coeff=self.D_cos)+
                                          fp.DiffusionTerm(coeff=self.D_p)+
                                          fp.PowerLawConvectionTerm(coeff=self.convCoeff_p)+
                                          fp.PowerLawConvectionTerm(coeff=self.convCoeff_cos)))
        return

    def update_parameters(self, E_Ec = 0, Te = 1./511., n_D = 1.):
        self.E.setValue(E_Ec)
        self.Te.setValue(Te)
        self.n_D.setValue(n_D)

    def solve(self, dt):
        self.eq.solve(var=self.f, dt=dt)

    def Maxwell(self, p,Te):
        gamma = fp.numerix.sqrt(p*p + 1.0)
        beta  = p/gamma
        return 4.0 * np.pi / (2.0 * np.pi * Te) ** 1.5 * beta ** 2 * fp.numerix.exp(- 0.5 * beta ** 2 / Te) / gamma **3

def main():
    ne = np.array([1.0613087, 1.0635650, 1.0658213, 1.0680777]) * 1.e18
    Te = np.array([0.4632629, 0.4650261, 0.4667893, 0.4685526]) * 1.e-3
    EEc = np.array([30.855909, 30.820590, 30.785421, 30.750401])

    for i in range(len(ne)):
        print('# of time steps = {}'.format(i))
        if i == 0:
            # 1) Initialize the ordinary solve instance
            re = DistributionSolver(p_min= 0.01 * np.sqrt(Te[0]),\
                    p_max=100., p_size=320, c_size=60, c_q=1., p_q=1.03)

            # 2) Initialize the distribution function
            re.f.setValue(ne[i]*re.Maxwell(re.p, Te[i]) / 2)

            # 3) Update parameter (this means updating convective / diffusive coefficients)
            re.update_parameters(E_Ec = EEc[0], Te=Te[0], n_D=ne[i])
        if i > 0:
            # 4) Initialize another solver instance at every time step
            re2 = DistributionSolver(p_min= 0.01 * np.sqrt(Te[0]),\
                    p_max=100., p_size=320, c_size=60, c_q=1., p_q=1.03)

            # 4) Initialize another distibution function with the original one every time step
            re2.f.setValue(re.f.value)
            # 4') Estimate rel error of f before solve
            print('Rel err of f before solve = {}'.format(np.nan_to_num(((re.f.value-re2.f.value)**2/re.f.value**2)**0.5).sum()))

            # 5) Update both parameters
            re.update_parameters(E_Ec = EEc[i], Te=Te[i], n_D=ne[i])
            re2.update_parameters(E_Ec = EEc[i], Te=Te[i], n_D=ne[i])
            # 5') Estimate rel error of convective/diffusive coefficients before solve
            rel_err_conv_p = np.nan_to_num(((re.convCoeff_p.value - re2.convCoeff_p.value)**2/re.convCoeff_p.value**2)**0.5).sum()
            rel_err_conv_cos = np.nan_to_num(((re.convCoeff_cos.value - re2.convCoeff_cos.value)**2/re.convCoeff_p.value**2)**0.5).sum()
            rel_err_D_p = np.nan_to_num(((re.Diff_p.value - re2.Diff_p.value)**2/re.Diff_p.value**2)**0.5).sum()
            rel_err_D_cos = np.nan_to_num(((re.Diff_cos.value - re2.Diff_cos.value)**2/re.Diff_cos.value**2)**0.5).sum()
            print('Rel err of conv_p before solve = {}'.format(rel_err_conv_p))
            print('Rel err of conv_cos before solve = {}'.format(rel_err_conv_cos))
            print('Rel err of D_p before solve = {}'.format(rel_err_D_p))
            print('Rel err of D_cos before solve = {}'.format(rel_err_D_cos))

            # 6) Solve both equations
            re.solve(dt=(2*re.Te.value)**1.5)
            re2.solve(dt=(2*re.Te.value)**1.5)
            # 6') Estimate rel error of f after solve
            print('Rel err of f after solve = {}'.format(np.nan_to_num(((re.f.value - re2.f.value)**2/re.f.value**2)**0.5).sum()))

if __name__ == "__main__":
    main()

[guyer]

I see the different solutions you're talking about and I have no idea why it's happening, but I also don't really understand what you're trying to do.

You seem to have done a thorough job of isolating the two DistributionSolvers from each other, but I still suspect that something is shared between them. I recommend you simplify to try to isolate the difference.

• Does it still happen with a uniform grid?
• Does it still happen with a 1D grid?
• Does it still happen with simpler term coefficients?
• Does it still happen with fewer terms?
• Does it still happen if you completely duplicate and separate the code for the two instances?

One thing I found is that the difference goes away if I change to

            re.update_parameters(E_Ec = EEc[2], Te=Te[2], n_D=ne[2])
            re2.update_parameters(E_Ec = EEc[2], Te=Te[2], n_D=ne[2])

or

            re.update_parameters(E_Ec = EEc[3], Te=Te[3], n_D=ne[3])
            re2.update_parameters(E_Ec = EEc[3], Te=Te[3], n_D=ne[3])

I see no reason why this should be different from using index i, if i==2 or i==3.

[0pago]

I see the different solutions you're talking about and I have no idea why it's happening, but I also don't really understand what you're trying to do.
What I was doing is to simulate an evolution of the electron distribution function with dynamical plasma parameters (EEc, Te and ne). The plasma parameters were prescribed by other codes. In this task, I encountered a problem that a solution is different with the reference solution. I did many tests and eventually obtained a solution which makes sense when I reinitialize solver like the re2 instance.
From this experience, I suspected my implementation of code would be incorrect and checked if convective and diffusive coefficients are updated correctly and the distribution function too. However, those are updated correctly. The only difference was solutions after eq.solve. This code I uploaded here is the test code I used to debug.

You seem to have done a thorough job of isolating the two DistributionSolvers from each other, but I still suspect that something is shared between them. I recommend you simplify to try to isolate the difference.
This is a very good comment to me. I appreciate it. I accept your recommendation and will check a list of bullets below.

One thing I found is that the difference goes away if I change to
Here, do you mean i == 2 or i == 3 after for i in range(2) ? If I use i == 2 or i == 3 after for loop, there is still discrepancy.

[0pago]

Dear guyer,

Hello! I leave this message because I predicted I will take more time to test them (your check list). If you recommend, I close this issue now and raise this again after following the list.

[guyer]

Leaving the issue open is OK