Open 0pago opened 1 year ago
Can you show some real code? Too much of the pseudo-code doesn't make any sense.
Dear guyer,
Thank you for your reply. I will add more detail after simplifying the code and checking this issue.
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()
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.
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
.
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 usei == 2
or i == 3
after for loop, there is still discrepancy.
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.
Leaving the issue open is OK
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?