Return to Question

The ratio of bisulfate in the concentrated tank  : dilute tank that I model agrees with the ratios they got, so I feel okay about my code. I think there may be an error in how I implemented the adaptive time step and I think I can speed it up. It takes 9 minutes to model 5 days of diffusion on my machine with no scaling constant (x_c = t_c = 1).

The ratio of bisulfate in the dilute tank:conc. tank that I model agrees with the ratios they got, so I feel okay about my code. I think there may be an error in how I implemented the adaptive time step and I think I can speed it up. It takes 9 minutes to model 5 days of diffusion on my machine with no scaling constant (x_c = t_c = 1).

I incorporated boundary conditions 1 and 2 (eqs 14d/e) by taking the time derivative of the boundary conditions so then $\frac{dC_m}{dxt} = K\cdot \frac{dC_c}{dxt}$ and then I can use eq 13b to update $C_m$ at x=0 and perform a similar operation for x = L and eq 13c

I incorporated boundary conditions 1 and 2 (eqs 14d/e) by taking the time derivative of the boundary conditions so then $\frac{dC_m}{dt} = K\cdot \frac{dC_c}{dt}$ and then I can use eq 13b to update $C_m$ at x=0 and perform a similar operation for x = L and eq 13c

Update

After continuing to work on this I have updated my code. I used the authors system parameters (V, A, C_0 etc.) to ensure that I can recreate their plots as a check.

I am now using Runge-Kutta with adaptive time steps. The adaptive timesteps seem to be key, because at low times, very far from membrane equilibration, dt must be very small because there are very large gradients in the membrane ($\frac{dC_m}{dx}$ very very large). But then dt must be allowed to increase as diffusion occurs (and $\frac{dC_m}{dx}$ decreases) so that modeling longer than a second does not take an hour.

I incorporated boundary conditions 1 and 2 (eqs 14d/e) by taking the time derivative of the boundary conditions so then $\frac{dC_m}{dxt} = K\cdot \frac{dC_c}{dxt}$ and then I can use eq 13b to update $C_m$ at x=0 and perform a similar operation for x = L and eq 13c

Cell setup using the author parameters:

D = (36e-3 * 10**-4) * np.exp(-26.3e3 / (8.314 * (273 + 25))) # m^2 / s, Diffusion Constant
A = .001 # m^2, area of the membrane
V = 0.001 # m^3, the volume in each bath
Cc_0 = 9.606e4 # ppm (unitless), the initial concentration in the cathode compartment
Ca_0 = 0 #ppm, initial concentration of the anode compartment

membrane_thickness = 183e-6 # m, the thickness of the nafion membrane in meters
Vc = Va = V # Setting volume of each tank, tank volumes are both 1 L

C_c = Cc_0 # Scaling constant for concentration in the concentrated tank set to inital concentration
Ca_c = C_c * (Vc / Va) # scaling constant for concentration in the dilute tank

k_partition = .23

Functions for solving the equations:

def f(mesh, dx_, x_c = 1, t_c = 1): #f = dCi / dt at each mesh point, performs calculations for membrane and both tanks
    # x_c and t_c are scaling constants to make equations dimensionless, default is 1 (not scaled)
    kappa_13b = (A * D * t_c)/(Vc * x_c) 
    kappa_13c = -(A * D * t_c)/(Va * x_c)
    kappa_13a = ( D * t_c ) / x_c**2
    
    dCm_dx = np.gradient(mesh[padding:-padding], dx_, axis = 0) # for eq 13b/c, only use the membrane region in the gradient
    
    dx0_ = dCm_dx[0] # for eq 13b. dCm/dx at x = 0. A scaler value
    dxl_ = dCm_dx[-1] #for eq 13c. dCm/dx at x = L. A scaler value
    
    # for eq 13a
    mesh_plus = np.roll(mesh, -1) # gives mesh value at (x + 1) for each mesh point. A vector 
    mesh_minus = np.roll(mesh, 1) # gives mesh value at (x - 1) for each mesh point. A vector
    # The first and last entry of mesh_plus and mesh_minus are invalid values. But are not used in eq 13a, so no error.
    
    second_deriv = ((mesh_plus+mesh_minus - (2*mesh)) / dx_**2)[1:-1] # for eq 13a. use entire mesh to calculate second derivative, keep only membrane region
    # The first and last values of the 2nd derivative within the membrane region are incorrect because they use values from the tanks during their calculation
    # However, these two values will be replaced with boundary conditions so no error
    
    #put calculations together to form the RHS of eqns 13a, 13b, and 13c. A single vector is returned
    fun = np.empty(mesh.shape) # the return mesh which will contain dCi/dt at each mesh point
     
    # f(t,x) for the cathode tank
    fun[:padding] = kappa_13b * dx0_ # eq 13b, f = dC/dt for the cathode tank (concentrated tank)
    
    #f(t,x) for the membrane region
    fun[padding:-padding] = kappa_13a * second_deriv # eq 13a, f = dC/dt to update interior of membrane
    
    #f(t,x) for the anode tank
    fun[-padding:] = kappa_13c * dxl_ # eq 13c, f = dC/dt for the anode tank (dilute tank)
    
    # Boundary Conditions
    # f(t,x) = dC/dt at x = 0 and x = L  set by the boundary conditions BC1 and BC2 from the paper
    fun[padding] = k_partition * kappa_13b * dx0_ # time derivative of BC1 (eq 14d) gives dCm/dt = k_partition * dCc/dt at x=0
    fun[-(padding+1)] = k_partition * kappa_13c * dxl_ # time derivative of BC2 (eq 14e) gives dCm/dt = k_partition * dCa/dt at x=L
    
    return fun

@jit
def rk4_step(f, mesh, dx_, dt_, x_c, t_c):
    
    k1 = dt_ * f(mesh, dx_, x_c, t_c)
    k2 = dt_ * f(mesh + .5*k1, dx_, x_c, t_c)
    k3 = dt_ * f(mesh + .5*k2, dx_, x_c, t_c)
    k4 = dt_ * f(mesh + k3, dx_, x_c, t_c)
    
    return (k1 + 2*k2 + 2*k3 + k4)/6

@jit
def solver(mesh, T, dt, dx, t_c=1, x_c=1, target_accuracy = 1e-5):
    cell = mesh.copy()
    
    # t_c and x_c are scaling constants
    # default scaling constant is 1 (no scaling)
    dt_ = dt / t_c
    dx_ = dx / x_c
    
    T_ = T / t_c
    
    t_ = 0
    cells = [(t_c*t_, t_c*dt_, cell)] # (current time, timestep for stepping to next point, cell state)
    while t_ < T_:
        # adaptive timestep method
        #     first test_step starts at t and steps twice using dt_ to get to test_step_1
        single_step  = cell + rk4_step(f, cell, dx_, dt_, x_c, t_c)
        test_step_1 = single_step + rk4_step(f, single_step, dx_, dt_, x_c, t_c)
        
        #     second test_step starts at t and steps once using 2*dt_ to get to test_step_2
        test_step_2 = cell + rk4_step(f, cell, dx_, 2*dt_, x_c, t_c)
        
        difference = np.abs(test_step_2 - test_step_1) # a vector of differences for each mesh point
        largest_diff = np.max(difference) # want to look at the worst discrepency between the two test steps
        
        if largest_diff != 0: # if largest_diff == 0, then dt_ does not need to be updated
            rho = (30 * dt_ * target_accuracy) / largest_diff # ratio of target accuracy and actual accuracy
            if rho > 1: # accuracy is better at every mesh point than required, keep data and update dt_ to a larger step size
                t_ = t_ + dt_
            dt_ = np.power(rho, .25) * dt_ # calculate correct time step based on desired accuracy using the worst rho value
            #dt_ = min(dt_calc, 2*dt_) # do no allow the new dt_ be more than twice the original dt_
            
            if rho < 1: # dt_ was too large for target accuracy, calculate step again using new dt_
                single_step = cell + rk4_step(f, cell, dx_, dt_, x_c, t_c)
                t_ = t_+dt_
        else:
            t_ = t_ + dt_
        
        new_cell = single_step
        cells.append((t_c*t_, t_c*dt_, new_cell))        
        cell = new_cell
    
    return cells

And here is the jupyter cell I use to implement and play around with the above code:

x_c = membrane_thickness# scaling constant for position
t_c = np.power(x_c, 2) / D # scaling constant for time

membrane_mesh_size = 15 # found to strongly affect the solving time

dx = membrane_thickness / membrane_mesh_size
dt = 1e-12 # starting timestep, the adaptive timestep in the solver will correct this to a more appropriate value

padding = 1 # the number of extra mesh points padding the start and end of the membrane mesh
            # holds values for the tanks

n_x =int( membrane_mesh_size + 2*padding)

fresh = np.zeros((n_x, 1))
fresh[:padding] = Cc_0 #/ C_c # pad mesh points set to tank starting concentration
fresh[-padding:]= Ca_0 #/ Ca_c# terminal pad mesh points set to anode tank starting concentration

fresh[padding] = k_partition * fresh[0] # starting membrane conditions at x=0
fresh[-(padding + 1)] = k_partition * fresh[-1] # starting membrane conditions at x=L

hours = 24 * 5
T = hours * 3600

t_c = x_c = 1

%time _ = solver(fresh, T, 1e-6, dx, t_c, x_c, target_accuracy=1e-6)

With this code, I can replicate Figure 4 from the publication.

The ratio of bisulfate in the concentrated tank : dilute tank that I model agrees with the ratios they got, so I feel okay about my code. I think there may be an error in how I implemented the adaptive time step and I think I can speed it up. It takes 9 minutes to model 5 days of diffusion on my machine with no scaling constant (x_c = t_c = 1).