Defining boundary conditions through bc_ops for 'n' coupled pde's

[martin-phys]

Hello,

I am working on a project to numerically evaluate a system of n coupled pde's in only one spatial dimension. The system is such that n1 of the equations are to have boundary conditions of a constant value (but unique to each equation) on the left side and a constant value for the derivative (here, set to 0 for all equations) on the right side of the model. n2 equations are the opposite, where n1+n2=n.

I am attempting to set these conditions via:

bc_ops1 = [({"value":n0L2[i]}, {"derivative":0}) for i in range(n1)] bc_ops2 = [({"derivative":0}, {"value":n0R2[i+n1]}) for i in range(n2)] bc_ops = [] bc_ops.append(bc_ops1) bc_ops.append(bc_ops2)

where n0L2 and n0R2 are lists of the constant values I'd like the equations to have on the left and right ends respectively.

The code runs without error but the boundary conditions do not seem to be enforced in the simulation. Am I doing this wrong?

Also, my code runs without specifying any boundary conditions at all, which seems odd to me. Is there a default that py-pde uses for bc's when they are not specified?

I hope my questions make sense, I can add details, more code snippets, or outputs if needed.

Thanks in advance, Martin

[david-zwicker]

I assume that you use the PDE class to enforce boundary conditions. If so, the argument bc defines boundary conditions common to all operators and fields. It defaults to auto_periodic_neumann, implying Neumann conditions for non-periodic boundaries.

Since you want to set specific boundary conditions for individual fields, you need to (additionally) use the argument bc_ops, which is explained in the documentation of the PDE class. This argument must be a dictionary (which wasn't really enforced so far). Modifying your code minimally, I would try the following:

bc_ops = {}
for i in range(n1):
    bc_ops["{variable[i]}:*"] = ({"value": n0L2[i]}, {"derivative": 0})
for i in range(n1, n2):
    bc_ops["{variable[i]}:*"] = ({"derivative": 0}, {"value": n0R2[i]})

where variable is a sequence of the names of the fields in the PDE.

[david-zwicker]

Closed by #542

[martin-phys]

Thank you for the fast help! I've tried your suggestion and it appears I am getting the same results (no errors but boundary conditions not properly enforced). I've given more of my code below in case you can spot anything else I might be doing wrong:

storage = MemoryStorage()

class SpeciesPDE(PDEBase):
    def __init__(self, D, A, r, bc_ops):
        self.D = D
       self.A = A
        self.r = r
        self.bc = bc_ops
        super().__init__()

    def evolution_rate(self, state, t=0):
        print("Time: ", t)
        u = state.data
        dudt = np.zeros_like(u)
        du_dx = np.gradient(u, dx, axis=1)
        d2u_dx2 = np.gradient(du_dx, dx, axis=1)
        interaction_terms = np.dot(self.A, u)  

        for i in range(n):
            dudt[i] += self.r[i] * u[i]  # Reproduction term
            dudt[i] += u[i] * interaction_terms[i]  # Interaction terms
            dudt[i] += self.D[i]*d2u_dx2[i]
        scalar_fields = [ScalarField(grid=state.grid, data=dudt[i, :]) for i in range(n)]
        return FieldCollection(scalar_fields)

# Create the grid
grid = CartesianGrid([[0, 150]], Nx)

bc_ops = {}
for i in range(n1):
    bc_ops["{u[i]}:*"] = ({"value": n0L2[i]}, {"derivative": 0})
for i in range(n1, n2):
    bc_ops["{u[i]}:*"] = ({"derivative": 0}, {"value": n0R2[i]})
# Create the initial state
initial_state = FieldCollection([ScalarField(grid, data=u0[i, :]) for i in range(n)])
#print(bc_ops)
# Create the PDE
eq = SpeciesPDE(D=D, A=Anew, r=r, bc_ops = bc_ops)

trackers = [
    "progress",  # show progress bar during simulation
    "steady_state",  # abort when steady state is reached
     storage.tracker(1) 
]

eq.solve(initial_state, t_range=(0, T), dt=Dt, tracker=trackers, adaptive = True)

[david-zwicker]

You're using the numpy gradient function and not the gradients of the fields in your evolution_rate. Try replacing d2u_dx2[i] by state[i].laplace(self.bc[i]) or something like that.

[martin-phys]

Thanks again for the help, I've tried your suggestion, as well as a few different other things but I cannot get the boundary conditions to work properly. Your suggestion, replacing the NumPy gradient functions with the laplace function gives me an error that is hard to interpret. The error output is simply "KeyError: 0" for the lined2u_dx2[i] = state[i].laplace(self.bc[i]) in:

# PDE solver

# Time points
t_span = (0, T)

storage = MemoryStorage()

# Define the PDE system
class SpeciesPDE(PDEBase):
    def __init__(self, D, A, r, bc_ops):
        self.D = D
        self.A = A
        self.r = r
        self.bc = bc_ops
        super().__init__()

    def evolution_rate(self, state, t=0):
        #print("Time: ", t)
        u = state.data
        dudt = np.zeros_like(u)
        interaction_terms = np.dot(self.A, u)  # Matrix multiplication to compute vector of interaction terms

        # Compute the Laplacian for each scalar field in the FieldCollection
        #d2u_dx2 = [field.laplace(bc=self.bc) for field in state]

        for i in range(n):
            d2u_dx2[i] = state[i].laplace(self.bc[i])
            dudt[i] += r[i] * u[i]  # Reproduction term
            dudt[i] += u[i] * interaction_terms[i]  # Interaction terms
            dudt[i] += self.D[i]*d2u_dx2[i]
            dudt[i,0] = 0
            dudt[i,-1] = 0
        scalar_fields = [ScalarField(grid=state.grid, data=dudt[i, :]) for i in range(n)]

        return FieldCollection(scalar_fields)

# Create the grid
grid = CartesianGrid([[0, L]], Nx)

bc_ops = {}
for i in range(n1):
    bc_ops["{u[i]}:*"] = ({"value": n0L2[i]}, {"derivative": 0})
for i in range(n1, n2):
    bc_ops["{u[i]}:*"] = ({"derivative": 0}, {"value": n0R2[i]})

# Create the initial state
initial_state = FieldCollection([ScalarField(grid, data=u0[i, :]) for i in range(n)])

# Create the PDE
A = Anew(ep)
eq = SpeciesPDE(D=D, A=A, r=r, bc_ops=bc_ops)

trackers = [
     storage.tracker(interval=Dt) 
]

eq.solve(initial_state, t_range=(0, T), dt=Dt, tracker=trackers, adaptive = True)

When I instead comment out the line d2u_dx2[i] = state[i].laplace(self.bc[i]) in the loop and uncomment the line d2u_dx2 = [field.laplace(bc=self.bc) for field in state] above the loop, I get a more interpretable error message:

    raise BCDataError(

BCDataError: Boundary condition `{u[i]}:*` not defined. Possible types of boundary conditions are 'user', 'virtual_point', 'value_expression', 'derivative_expression', 'mixed_expression', 'value', 'derivative', 'mixed', 'curvature', 'normal_value', 'normal_derivative', 'normal_mixed', 'normal_curvature'. Values can be set using {'type': TYPE, 'value': VALUE}. Here, VALUE can be a scalar number, a vector for tensorial boundary conditions, or a string, which can be interpreted as a sympy expression. In the latter case, the names of the axes not associated with this boundary can be used as variables to describe inhomogeneous boundary conditions.

But I am unable to figure this out. Am I doing anything else wrong that is obvious to you? here, u is an n by Nx array which contains the values of each n quantities at each Nx spatial position.

[david-zwicker]

The problem is likely in your code where bc_ops is constructed as a dictionary, but you later try to use it more like a list. This has probably nothing to do with py-pde.

[martin-phys]

Ah, yes I was making a mistake with referencing the dictionary values. I think I have fixed this, but I am now getting an incompatibility error that does seem to be coming from py-pde. The error is:

  File ~\anaconda3\Lib\site-packages\pde\fields\base.py:383 in assert_field_compatible
    raise TypeError(f"Fields {self} and {other} are incompatible")

TypeError: Fields ScalarField(grid=CartesianGrid(bounds=((0.0, 150.0),), shape=(500,), periodic=[False]), data=Array(500,)) and [ 0.00000000e+00 -6.19103585e-05 -1.23820717e-04 -1.85731075e-04 ...] are incompatible

for this updated code:

# Time points
t_span = (0, T)

storage = MemoryStorage()

# Define the PDE system
class SpeciesPDE(PDEBase):
    def __init__(self, D, A, r, bc_ops):
        self.D = D
        self.A = A
        self.r = r
        self.bc = bc_ops
        super().__init__()

    def evolution_rate(self, state, t=0):
        print("Time: ", t)
        u=state.data
        dudt = np.zeros_like(u)
        interaction_terms = np.dot(self.A, u)  # Matrix multiplication to compute vector of interaction terms
        d2u_dx2 = np.zeros((n,Nx))
        for i in range(n):
            d2u_dx2 = state[i].laplace(bc=self.bc['{u[%i]}:*' %i])
            dudt[i] += r[i] * u[i]  # Reproduction term
            dudt[i] += u[i] * interaction_terms[i]  # Interaction terms
            dudt[i] += self.D[i]*d2u_dx2
        scalar_fields = [ScalarField(grid=state.grid, data=dudt[i, :]) for i in range(n)]

        return FieldCollection(scalar_fields)

# Create the grid
grid = CartesianGrid([[0, L]], Nx)

bc_ops = {}
for i in range(n1):
    bc_ops["{u[%i]}:*" %i] = ({"value": n0L2[i]}, {"derivative": 0})
for i in range(n1, n):
    bc_ops["{u[%i]}:*" %i] = ({"derivative": 0}, {"value": n0R2[i]})

# Create the initial state
initial_state = FieldCollection([ScalarField(grid, data=u0[i, :]) for i in range(n)])

# Create the PDE
A = Anew(ep)
eq = SpeciesPDE(D=D, A=A, r=r, bc_ops=bc_ops)

trackers = [
    #"progress",  # show progress bar during simulation
    #"steady_state",  # abort when steady state is reached
     storage.tracker(interval=Dt)
]

eq.solve(initial_state, t_range=(0, T), dt=Dt, tracker=trackers, adaptive = True)

I am looking into what might be wrong but if you can help it would be appreciated.

Martin

[david-zwicker]

Sorry, but I don't have the bandwidth to debug your code for you. If you have a specific issue that is really related to py-pde I can help, but such general questions are better asked over at stackoverflow.com or similar webpages.

[martin-phys]

Sorry for bothering you, I understand you're very busy.

After some tinkering, I was able to figure out the issue with my code. I've attached a working version below in case anyone is interested:

##################################################################################
# PDE solver

# Time points
t_span = (0, T)

storage = MemoryStorage()

# Define the PDE system
class SpeciesPDE(PDEBase):
    def __init__(self, D, A, r, bc):
        self.D = D
        self.A = A
        self.r = r
        self.bc = bc_ops
        super().__init__()

    def evolution_rate(self, state, t=0):
        print("Time: ", t)
        u=state.data
        dudt = np.zeros_like(u)
        interaction_terms = np.dot(self.A, u)  # Matrix multiplication to compute vector of interaction terms
        d2u_dx2 = np.zeros((n,Nx))
        for i in range(n):
            d2u_dx2[i] = state[i].laplace(bc=self.bc['{u[%i]}:*' % i]).data
            dudt[i] = self.D[i]*d2u_dx2[i]
            dudt[i] += r[i] * u[i]  # Reproduction term
            dudt[i] += u[i] * interaction_terms[i]  # Interaction terms

        scalar_fields = [ScalarField(grid=state.grid, data=dudt[i, :]) for i in range(n)]

        return FieldCollection(scalar_fields)

# Create the grid
grid = CartesianGrid([[0, L]], Nx)

bc_ops = {}
for i in range(n1):
    bc_ops["{u[%i]}:*" %i] = ({"value": n0L2[i]}, {"derivative": 0})
for i in range(n1, n):
    bc_ops["{u[%i]}:*" %i] = ({"derivative": 0}, {"value": n0R2[i]})

print(bc_ops)
print(bc_ops['{u[0]}:*'])

# Create the initial state
initial_state = FieldCollection([ScalarField(grid, data=u0[i, :]) for i in range(n)])

# Create the PDE
A = Anew(ep)
eq = SpeciesPDE(D=D, A=A, r=r, bc= bc_ops)

trackers = [
    #"progress",  # show progress bar during simulation
    #"steady_state",  # abort when steady state is reached
     storage.tracker(interval=Dt)
]

eq.solve(initial_state, t_range=(0, T), dt=Dt, tracker=trackers, adaptive = True)

Thanks for the help, Martin