Strang splitting for solving stiff source FiPy

Question

I have a problem with trying to use strang splitting with fipy to solve for stiff systems. Here is a sample code attached. At steady state I have balance issues. Would there be any problem with the implementation? Thanks.

import numpy
import scipy
import math
import fipy

def fun(t,u):
    return 20.

nx=50
dx=0.02
u_b=1.
dt=dx/u_b*0.5
mesh=fipy.Grid1D(nx=nx,dx=dx)
X=mesh.cellCenters

u=fipy.CellVariable(mesh=mesh,name="Var",value=0.,hasOld=True)
u.constrain(1.,where=mesh.facesLeft)
u.faceGrad.constrain((0.,),where=mesh.facesRight)



eq=fipy.TransientTerm(var=u) == -fipy.ExponentialConvectionTerm(coeff=(u_b,),var=u)
time_steps=120

for i in range(0,time_steps):
    
    
    eq.solve(dt=dt/2)
    u.updateOld()
    
    u_temp=u.value.copy()
    
    result=scipy.integrate.solve_ivp(fun, [0,dt], u_temp)
    
    for o in range(0,nx):
        u.setValue(result.y[o][-1],where=((X>=o*dx) & (X<=(o+1)*dx)))
    u.updateOld()
    
    eq.solve(dt=dt/2)
    
    u.updateOld()

FWIW, you can save about 1/3 of the runtime by doing `u.setValue(result.y[..., -1])` instead of looping `o` over `nx`. — jeguyer, Aug 01 '22 at 19:29
The amount added by the source term is not reflected in the change of flux across the domain at steady state. — Markiv, Aug 02 '22 at 10:27
The source side for example has 20*nx*dx=20. added to the domain, while the outlet - inlet flow is given by (with the u values at the end of the simulation) (20.74*u_b-1*u_b) = 19.74. — Markiv, Aug 02 '22 at 10:37
The problem hasn't reached steady state at the end of 120 steps. If I let it keep running and add `print(i, (u_b * u.faceValue[..., mesh.facesRight] - u_b * u.faceValue[..., mesh.facesLeft]))`, I get `119 [19.74703003]`, `169 [20.0498342]`, `217 [20.05]`, ... at which point it appears to be at steady state . The remaining discrepancy is confirmed to be `dx / 4` by varying the mesh resolution. If I had to guess, I suspect that is due to the integration range not matching the simulation range. — jeguyer, Aug 02 '22 at 12:49
I'm sorry I didn't understand your last point. For the time range, a full step of advection and the source integration is performed in a single time iteration. For the spatial range, as far as I understand, the strang splitting takes the integrates over each cell with initial cell value obtained from the first half step of advection. So in theory it covers the whole domain. — Markiv, Aug 02 '22 at 15:14
Like I said, it was a guess. It's something I would explore carefully if I were you. Another possibility is that, if the analytical solution is supposed to be a line, is that the FiPy solution with have the same value in the rightmost cell (`u[..., -1]`) and at the rightmost face (`u.faceValue[..., mesh.facesRight]`). This is what `u.faceGrad.constrain((0.,),where=mesh.facesRight)` means. It's possible this leads to a discretization error giving rise to the 0.05 excess. I don't know. I don't know anything about Strang splitting. — jeguyer, Aug 02 '22 at 15:46
[This notebook](https://gist.github.com/guyer/6b7699531f75b353f49a0cf36d683aa4) illustrates some of the counterintuitive things that can happen at boundaries. I'm not sure whether they're at play for this problem or not. — jeguyer, Aug 02 '22 at 16:00

Strang splitting for solving stiff source FiPy

0 Answers0