Using boundary_value and setting boundary conditions with it

[CSchwet]

Hi, I do not quite understand what pybamm.boundary_value() returns. I have attached an example code below. It outputs at the end for the variable phi (an algebraic variable) at two different times:

1. the boundary condition as defined in model.boundary_conditions
2. the edge value of the variable
3. the value pybamm.boundary_value obtains.

pybamm.boundary_value seems to hit neither of these values. What is it calculating instead?

In the end I want to implement boundary conditions which depend on either a boundary value or gradient:
If I use the commented out verion, i.e replace lines 62 & 63 with

gradmu1 = pybamm.grad(phi + mu1) - gamma1 * rho * pybamm.grad(phi)
gradmu2 = - pybamm.grad(phi - mu2) - gamma2 * rho * pybamm.grad(phi)

and lines 90 & 94 with

"left": (gamma1 * pybamm.boundary_value(rho, "left") * pybamm.boundary_gradient(phi, "left"), "Neumann"),

"left": (-gamma2 * pybamm.boundary_value(rho, "left") * pybamm.boundary_gradient(phi, "left"), "Neumann"),

it fails at about t=7.6. Just before, the bc. I wanted to set, that gradmu1=gradmu2=0 on the left hand side also fails, most likely triggering the instability.

How would I enforce this boundary condition, if not using boundary_value, or what is the correct syntax?

Thank you!

import numpy as np
import pybamm

model = pybamm.BaseModel()

# Domain  variables, all normalised
c3 = pybamm.Variable("Concentration 3", domain="electrolyte")
dilution = pybamm.InputParameter("dilution")
rho = pybamm.Variable("Charge density", domain="electrolyte")
phi = pybamm.Variable("Potential", domain="electrolyte")
phi0 = pybamm.InputParameter("phi0")
V_ext = pybamm.InputParameter("V_ext")

# Species Parameters
z1 = pybamm.InputParameter("z1") # z2 has to take value -z1, z3 = 0
t1 = pybamm.InputParameter("t1") # t2 = 1 - t1
gamma1 = pybamm.InputParameter("gamma1") # molar volume 1 / molar volume salt/IL
gamma2 = pybamm.InputParameter("gamma2")
gamma3 = pybamm.InputParameter("gamma3")

# Global parameters
D = pybamm.InputParameter("D")

def Phi_boundarycond(times,rampingfraction):
    length = (int) (np.size(times)/rampingfraction)
    Phibc = np.concatenate((np.linspace(0,1,length),np.ones(np.size(times)-length)))
    Phibc = pybamm.Interpolant(times,Phibc,pybamm.t)
    return Phibc



times = np.linspace(0,7.5,10000) #np.concatenate((np.linspace(0,0.1,1000),np.linspace(0.1001,10,10000)))
Phibc = Phi_boundarycond(times,3) * V_ext

# Equations
def mu_fct(species, rho, c3, gamma1, gamma2, gamma3):
    c1 = 1 + gamma2 * rho - c3
    c2 = 1 - gamma1 * rho - c3
    c = c1 + c2 + c3 / gamma3

    rhophi = c1 * pybamm.log(c1 / c) + c2 * pybamm.log(c2 / c) + c3 / gamma3 * pybamm.log(c3 / gamma3 / c)

    mu1 = rhophi.diff(c1)
    mu2 = rhophi.diff(c2)

    if species == 1:
        return mu1
    else:
        if species == 2:
            return mu2
        else:
            raise AssertionError("Species must be 1 or 2.")

mu1 = mu_fct(1, rho, c3, gamma1, gamma2, gamma3)
mu2 = mu_fct(2, rho, c3, gamma1, gamma2, gamma3)


c1 = 1 + gamma2 * rho - c3
c2 = 1 - gamma1 * rho - c3


gradmu1 = pybamm.grad(phi + mu1) #- gamma1 * rho * pybamm.grad(phi)
gradmu2 = - pybamm.grad(phi - mu2) #- gamma2 * rho * pybamm.grad(phi)

zeta1 = gradmu1 + gamma1 / c3 * (c1 * gradmu1 + c2 * gradmu2)
zeta2 = gradmu2 + gamma2 / c3 * (c1 * gradmu1 + c2 * gradmu2)

J = - t1 / z1 * zeta1 + (1 - t1) / z1 * zeta2 # charge flux density
N3 = D * (zeta2 + zeta1)

drhodt = - pybamm.div(J)
dc3dt = - pybamm.div(N3)

poisson = rho + pybamm.div(pybamm.grad(phi))

phi_boundary = pybamm.boundary_value(phi, "left")

model.rhs = {rho: drhodt, c3: dc3dt}
model.algebraic = {phi: poisson}#, lambda1: l1fct, lambda2: l2fct}

# boundary conditions
c30 = 1 - dilution
model.boundary_conditions = { 
    rho: {"right": (0, "Dirichlet")}, 
    c3: {"right": (c30, "Dirichlet")},
    J: {"left": (0, "Dirichlet")},
    N3: {"left": (0, "Dirichlet")},
    phi: {"left": (Phibc, "Dirichlet"), "right": (phi0, "Dirichlet")},
    phi + mu1: {
        "left": (0, "Neumann"), #(gamma1 * pybamm.boundary_value(rho, "left") * pybamm.boundary_gradient(phi, "left"), "Neumann"), 
        "right": (phi0 + mu_fct(1, 0, c30, gamma1, gamma2, gamma3), "Dirichlet")
        },
    phi - mu2: {
        "left": (0, "Neumann"), #(-gamma2 * pybamm.boundary_value(rho, "left") * pybamm.boundary_gradient(phi, "left"), "Neumann"), 
        "right": (phi0 - mu_fct(2, 0, c30, gamma1, gamma2, gamma3), "Dirichlet")
        },
}
mu_bdr = mu_fct(2, 0, c30, gamma1, gamma2, gamma3)

# initial conditions
model.initial_conditions = {c3: c30, rho: 0, phi: phi0}

model.variables = {
    "Concentration 3": c3,
    "Charge density": rho,
    "Potential": phi,
    "Concentration 1": 1 + gamma2 * rho - c3,
    "Concentration 2": 1 - gamma1 * rho - c3,
    "Current density": J,
    "Solvent flux": N3,
    "Force 1": mu1,
    "Force 2": mu2,
    "Mu-bdr": mu_bdr,
    "Phi boundary": phi_boundary,
}

x = pybamm.SpatialVariable("x", domain=["electrolyte"], coord_sys="cartesian")

geometry = {"electrolyte": {x: {"min": 0, "max": 100}}}

spatial_methods = {"electrolyte": pybamm.FiniteVolume()}
submesh_types = {"electrolyte": pybamm.Uniform1DSubMesh}
var_pts = {x: 400}
# create a mesh of our geometry, using a uniform grid with 20 volumes
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

parameter_values = pybamm.ParameterValues({
    "Dilution": dilution,
    "Charge number 1": z1,
    "Initial potential difference": phi0,
    "Applied external potential": V_ext,
    "Transport number 1": t1,
    "Relative molar volume 1": gamma1,
    "Relative molar volume 2": gamma2,
    "Relative molar volume 3": gamma3,
    "Diffusion constant": D,
})


solver = pybamm.CasadiSolver()

sim = pybamm.Simulation(
        model,
        geometry=geometry,
        parameter_values=parameter_values,
        submesh_types=submesh_types,
        var_pts=var_pts,
        spatial_methods=spatial_methods,
        solver=solver,
)
 
solution = sim.solve(t_eval = times, inputs = {
    "z1": 1, 
    "t1": 0.5, 
    "dilution": 0.33, 
    "phi0": 0, 
    "V_ext": 10, 
    "gamma1": 0.5, 
    "gamma2": 0.5, 
    "gamma3": 0.5,
    "D": 10, 
})


print("Boundary values at timestep 1")
print("BC: ", 0.003)
print("Edge value: ", solution["Potential"].entries[0,1])
print("Boundary_value: ", solution["Phi boundary"].entries[1])
print("Boundary values at timestep 999 (end)")
print("BC: ", 10)
print("Edge value: ", solution["Phi boundary"].entries[9999])
print("Boundary_value: ", solution["Potential"].entries[0,9999])

[brosaplanella]

I am not sure I follow the boundary condition you are trying to implement, but if it is a Robin boundary condition the best way is to implement it as a Neumann boundary condition, where the RHS depends on the boundary value. Also note that the keys of the boundary condition dictionary should be the variables, so defining the boundary condition for phi + mu1 actually doesn't get used. Your current model will only use the BCs for rho, c3 and phi (the ones appearing as keys of model.rhs and model.algebraic).