1d boundary conditions and equations setup

[uliw]

Looking at the examples, I get the idea of how fastfd defines a model, but I struggle to transfer this to a 1D setup. I am trying to model something like this

$$\frac{d C}{d t} = D_x \frac{d^2 C}{d x^2} - \omega \frac{d C}{d x} -f(x)$$

how would I define this equation, and how would I add a Dirichlet boundary condition for x = 0, and a Neuman condition for x[-1] ?

Thanks!

[uliw]

never mind, I figured it out. Pretty awesome library!

import fastfd as ffd
from matrepr import mprint
import matplotlib.pyplot as plt
import numpy as np


def msr(z):
    """ Define microbial activity as a function of depth
    The parameters below are probably bogus
    """
    import numpy as np

    a = 808.911
    b = 0.0205
    c = 53.28
    f = a * np.exp(-b * (z)) + c
    return -f / 1e16


# use scipy sparse matrix
ffd.sparse_lib("scipy")
z = ffd.LinearAxis("z", start=0, stop=301, num=301)
f = msr(z.coords)
C = ffd.Scalar("C", [z], accuracy=4)

# define constants
D = 2.9e-11  #  Diffusion constant m*2/s
w = 0  # advection velocity, negative is upwelling

# Instantiate the model
model = ffd.FDModel([C])

# Set model governing equations
model.update_equations({"diffusion": (D * C.d("z", 2) - w * C.d("z", 1), -f)})

# Set model boundary conditions
model.update_bocos(
    {
        "C z=0": (C.i[0], C.i[0], 28),  # upper end concentration
        "C z=-1": (C.i[-1], C.d("z")[-1], 0),  # lower end gradient
    }
)

# Solve the model
results = model.solve()
fig, ax = plt.subplots()
ax.plot(z.coords, results["C"])
fig.tight_layout()
plt.show()

[stefanmeili]

I'm excited that somebody's using this thing! The 1D wave example is probably the best starting point. Glad you found what you needed.

[uliw]

I guess it depends on ones background, how easy it is to understands someone elses API ;-) But here is a thing I cannot figure out. I have a second process that precipitates a fraction of the solute C in the above equation. It is easy to adjust f for this process, butr I also need the resulting mass of the precipitate. Since it does not diffuse, one could write

$$\frac{d S}{d x} = - \omega \frac{d S}{d x} +f(x) \gamma$$

but this results in a singular matrix when I set it up as

 "pyrite": (-w * FES2.d("x", 1), +f * g),

Is there a way to define this within your framework?

[uliw]

apologies for the noise, that actually works, it was a typo elsewhere

[stefanmeili]

No worries at all - this library isn't exactly well documented. Singular matrices generally a result of insufficient constraints.

This framework is pretty flexible - I've used this framework to model a pipe reactor tracking multiple species, byproducts and temperature scalars (making mononitrobenzene), and had a working CFD example that was laughably slow. This library's a good solution for small engineering problems, but for anything more, I'd recommend COMSOL, ANSYS, or OpenFOAM, depending on your application. It's been quite a few years since I've worked in this space, so I'm not up to date on the tools that are available.

There are problems where you need to linearize a second derivative by approximating it with two fist derivatives. The first is solved simultaneously with the rest of the equations, while the second is iteratively updated with the last solution value. You'll probably need some under relaxation to keep things stable. I don't have access to a python environment for three weeks, but will try to find an example when I'm home.

If you produce a working chemical reactor example that you're able to share, I'd welcome the contribution to the documentation (such as it is).

Cheers,

Stefan

[stefanmeili]

Found the example I was looking for in my email:

solving a diffusion + chemical reaction problem:
$$\frac{\delta C}{\delta z} = \frac{\delta^2 C}{\delta y^2} - k\cdot C^2$$

import fastfd as ffd
ffd.sparse_lib('scipy')

import numpy as np



# Set up the model
y = ffd.LinearAxis('y', start = 0, stop = 1, num = 51)
z = ffd.LinearAxis('z', start = 0, stop = 1, num = 51)

C = ffd.Scalar('C', [y, z], accuracy = 2)

model = ffd.FDModel([C])

model.update_bocos({
    'Cy0=0': (C.i[0, :], C.i[0, :], 0),
    'Cz0=0': (C.i[:, 0], C.i[:, 0], np.linspace(0, 1, 51)),
})



k = 10 # reaction coefficient?
C_last = np.ones((51, 51)) * 1e-6 # initial guess for C must have same shape as y x z. Initialized to 1e-6.


def update_model(C_last):
    model.update_equations({
        'Dif + Rx': (C.d('z') - C.d('y', 2) + k * C.i * C_last, 0),
    })
    
    return model.solve()



# Solve the model
relax = 0.9 # relaxation coefficient that damps solution between iterations
results, residuals = [], []
for i in range(50): # arbitrary number of iterations
    result = update_model(C_last)['C']
    
    results.append(result)
    
    residual = np.sqrt(np.mean(np.square(result - C_last)))
    residuals.append(residual)
    
    print(residual)
    '''
    residual measures how much the solution changes between iterations. when it falls below some threshold, say
    1e-8, you could break the loop and return the most recent result.
    '''
    
    if residual < 1e-8: break
       
    C_last = result * relax + C_last * (1 - relax) # Damping between iterations. This can help with instability.
    
    
    
# Dump out results
results[-1]

[uliw]

hmmh, that is close to what I am trying to do. Do I understand from your example that there is no way to access C during intergration? My reaction term has a dependence on C

[stefanmeili]

For nonlinearities like this you can treat the reaction constant like the updated derivative. Assume a value and iterate in a loop, substituting the new value until the solution converges. Depending on how steep the gradients get, you may need fourth order approximations and / or very small time steps and fine spatial resolution.

…

On Sat, Jun 8, 2024, 12:03 PM Ulrich Wortmann ***@***.***> wrote: hmmh, that is close to what I am trying to do. Do I understand from your example that there is no way to access C during intergration? My reaction term has a dependence on C — Reply to this email directly, view it on GitHub <#1 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEBBKYBK37EZUCYH34NEPP3ZGNIRRAVCNFSM6AAAAABI7HS7R6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNJWGE2DENBWGE> . You are receiving this because you commented.Message ID: ***@***.***>

[uliw]

Interesting approach. Would it not be more straightforward to treat it as a non-steady state problem to solve each time step explicitly (i.e., invert the coefficient matrix and multiply with the constraints to get the new C)?

[stefanmeili]

I guess the answer depend on if you're building a transient or steady state model. If transient, ideally the temperature dependant coefficients are solved simultaneously with the mass and energy balance. If your coefficients are linear wrt temperature this may be possible inside the matrix. They likely aren't and if this is the case you'll you need a coefficient update loop inside your time loop. Without this your model predictions will not be independent of time step size or spatial resolution. If this is a steady state model then you can just let the coefficients lag a time step and it should converge to the the correct solution.

…

On Sun, Jun 9, 2024, 4:50 PM Ulrich Wortmann ***@***.***> wrote: Interesting approach. Would it not be more straightforward to treat it as a non-steady state problem to solve each time step explicitly (i.e., invert the coefficient matrix and multiply with the constraints to get the new C)? — Reply to this email directly, view it on GitHub <#1 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEBBKYG2722YUUFLSJTQIJDZGRTVDAVCNFSM6AAAAABI7HS7R6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNJWGYZTQNRYGA> . You are receiving this because you commented.Message ID: ***@***.***>

[stefanmeili]

The same goes for density, viscosity, surface tension etc. it's probably easier to update them in a loop. You likely have correlations for those that won't play nice with a linear model.

…

On Sun, Jun 9, 2024, 5:17 PM Stefan Meili ***@***.***> wrote: I guess the answer depend on if you're building a transient or steady state model. If transient, ideally the temperature dependant coefficients are solved simultaneously with the mass and energy balance. If your coefficients are linear wrt temperature this may be possible inside the matrix. They likely aren't and if this is the case you'll you need a coefficient update loop inside your time loop. Without this your model predictions will not be independent of time step size or spatial resolution. If this is a steady state model then you can just let the coefficients lag a time step and it should converge to the the correct solution. On Sun, Jun 9, 2024, 4:50 PM Ulrich Wortmann ***@***.***> wrote: > Interesting approach. Would it not be more straightforward to treat it as > a non-steady state problem to solve each time step explicitly (i.e., invert > the coefficient matrix and multiply with the constraints to get the new C)? > > — > Reply to this email directly, view it on GitHub > <#1 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AEBBKYG2722YUUFLSJTQIJDZGRTVDAVCNFSM6AAAAABI7HS7R6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNJWGYZTQNRYGA> > . > You are receiving this because you commented.Message ID: > ***@***.***> >

[uliw]

I see. I played with this a bit, and your approach is fast. I'll reopen this discussion so it is visible in case someone else finds it useful. If you rather keep it closed, feel free to close it again. Thanks again!

[uliw]

I noticed that specifying two models in the same code like this

so4_model = ffd.FDModel([SO4])
fes_model = ffd.FDModel([SO4, H2S, FES2])

so4_model.update_bocos
fes2_model._update_bocos
so4_model.update_equations
fes2_model._update_equations
so4_results = so4_model.solve()
fes2_results = fes2.results

Results in errors. Rearranging this sequentially works, though.
BTW, how would I specify a flux boundary condition with fastdf?