# Using velocity profiles with DiffEqOperators.jl for 2D convection-diffusion problem

**URL:** <https://discourse.julialang.org/t/using-velocity-profiles-with-diffeqoperators-jl-for-2d-convection-diffusion-problem/49243>\
**Category:** Modelling & Simulations\
**Created:** [October 29, 2020, 3:26pm UTC](https://discourse.julialang.org/t/using-velocity-profiles-with-diffeqoperators-jl-for-2d-convection-diffusion-problem/49243 "2020-10-29T15:26:14Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![jarias9](https://avatars.discourse-cdn.com/v4/letter/j/cc9497/32.png) [@jarias9](https://discourse.julialang.org/u/jarias9)\
**Post date:** [October 29, 2020, 3:26pm UTC](https://discourse.julialang.org/t/using-velocity-profiles-with-diffeqoperators-jl-for-2d-convection-diffusion-problem/49243/1 "2020-10-29T15:26:14Z")

</div>

Hello all. I am trying to discretize and solve the following PDE using DiffEqOperators.jl

\frac{\partial C}{\partial t} = D\nabla^2 C - \vec{v} \cdot \nabla C + R

for solving non-steady convection diffusion of a generic chemical species with a generic velocity profile v\_x(y).

So far I have been able to generate and solve the problem using scalar values for the velocity field in each direction in both 2D and 3D using the following code:

```julia

# Loading Dependencies

using DifferentialEquations

using DiffEqOperators

using Plots

# Physical Constants

α = 0.05 # Diffusion Coefficient

v_x = 0.5; # Advection speed in x [m/s]

v_y = 0.25; # Advection speed in y [m/s]

γ=0.05 # Source Term

# Discretization Parameters

s = x, y = (-5:0.2:5, -5:0.2:5)

dx = dy = x[2] - x[1]

t0 = 0.0

t1 = 20.0

# Initial Distribution

u_low = convert(Int64, round(2 / 10 * length(x)))

u_high = convert(Int64, round(4 / 10 * length(x)))

u0 = fill(0.0, (length(x), length(y)));

u0[u_low:u_high,u_low:u_high] .= 5.0; # a small square in the middle of the domain

# Generate Coefficient Matrices

laplacianx = CenteredDifference{1}(2, 4, dx, length(x))

laplaciany = CenteredDifference{2}(2, 4, dy, length(y))

Δ = laplacianx + laplaciany

divx = CenteredDifference{1}(1, 4, dx, length(x))

divy = CenteredDifference{2}(1, 4, dy, length(y))

# Generate Boundary Conditions

Qx, Qy = Neumann0BC(Float64, (dx, dy), 4, (length(x), length(y)))

Q = compose(Qx, Qy)

# Define matrix system and solve

function f!(du, u, p, t)

    du .= α .* (Δ * Q * u) + (((-v_x * divx) + (-v_y * divy)) * Q * u) + γ .* u

### #Diffusion #Advection in x #Advection in y #Source

end

prob = ODEProblem(f!, u0, (t0, t1))

alg = Tsit5()

sol2D = solve(prob, alg)

```

My question is whether or not it is possible to use a velocity profile for the v\_x, v\_y, v\_w terms and some recommendations on how to go about it. Likewise I was wondering if there is a more compact form of the \vec{v} \cdot \nabla C expression that can be expressed using DiffEqOperators.jl. Thanks!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 5, 2020, 11:06pm UTC](https://discourse.julialang.org/t/using-velocity-profiles-with-diffeqoperators-jl-for-2d-convection-diffusion-problem/49243/2 "2020-11-05T23:06:02Z")

</div>

> [@jarias9](#):
>
> My question is whether or not it is possible to use a velocity profile for the vx,vy,vwv\_x, v\_y, v\_w terms and some recommendations on how to go about it.

I don’t quite get the question. Try it again?

---

<div class="post-metadata">

**Author:** ![jarias9](https://avatars.discourse-cdn.com/v4/letter/j/cc9497/32.png) [@jarias9](https://discourse.julialang.org/u/jarias9)\
**Post date:** [November 6, 2020, 4:24pm UTC](https://discourse.julialang.org/t/using-velocity-profiles-with-diffeqoperators-jl-for-2d-convection-diffusion-problem/49243/3 "2020-11-06T16:24:03Z")

</div>

Alright so the \vec{v} \cdot \nabla C portion of my governing equation models advection right. What this looks like expanded for an incompressible flow is v\_x \frac{\partial C}{\partial x} \vec{i} + v\_y \frac{\partial C}{\partial y} \vec{j} + v\_z \frac{\partial C}{\partial z} \vec{k}. Now when I say velocity profile, I mean using a function to give a varying velocity speed dependent on the position. So because they are not allowed to vary along its prinicipal axis, this would look like v\_x = v\_x(y,z), v\_y = v\_y(x,z), v\_z(y,z) etc.  
Now for the application, I was able to model a fixed scalar advection speed using the following code:

```julia
function f!(du, u, p, t)
    du .= α .* (Δ * Q * u) + (((-v_x * divx) + (-v_y * divy)) * Q * u) + γ .* u
### #Diffusion #Advection in x #Advection in y #Source
end

```

Where v\_x, v\_y are constant scalars and then `divx` is the divergence in x and so on. So what I would like to do is instead use a parabolic flow field given by v\_x(y) = v\_{max} \left[1 - \left( \frac{y - H/2}{H/2} \right) \right], for example but this could be any arbitrary function for the flow. So I’m thinking that if I can generate `v_x ` to be an array of the appropriate shape then I could do `.*` with the divergence operator `divx` for each direction and this will give me the behavior I want. I just don’t know what the ordering/shape would be like so that it corresponds to each line of y. Then lastly, I was wondering if there were any functions I could use to just do the dot multiply operation of \vec{v} \cdot \nabla C instead of having to expand it like this `(((-v_x * divx) + (-v_y * divy)) * Q * u)` (this one is more just wishful thinking). Hopefully this clears things up.

TLDR cross or dot product of a vector function by a DiffEqOperator.jl derivative operator (in this case `divx = CenteredDifference{1}(1, 4, dx, length(x))` in order to perform operations such as \vec{v} \cdot \nabla C which come up often in PDEs.

Thanks again.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 7, 2020, 10:03pm UTC](https://discourse.julialang.org/t/using-velocity-profiles-with-diffeqoperators-jl-for-2d-convection-diffusion-problem/49243/4 "2020-11-07T22:03:55Z")

</div>

Make `v_x` a vector and use `-v_x .* divx`

---

<div class="post-metadata">

**Author:** ![jarias9](https://avatars.discourse-cdn.com/v4/letter/j/cc9497/32.png) [@jarias9](https://discourse.julialang.org/u/jarias9)\
**Post date:** [November 13, 2020, 3:51pm UTC](https://discourse.julialang.org/t/using-velocity-profiles-with-diffeqoperators-jl-for-2d-convection-diffusion-problem/49243/5 "2020-11-13T15:51:13Z")

</div>

Hi Chris, I tried your suggestion and I got an error about a missing method regarding the `.*` when doing `v_x .* divx`. However I got a solution for what I’m trying to do that I want to document for other people trying to do dot products in `DiffEqOperators.jl`. The trick was actually changing the order of the operations for the governing equation. Essentially we went from

```julia

function f!(du, u, p, t)

    du .= α .* (Δ * Q * u) + (((-v_x .* divx) + (-v_y .* divy)) * Q * u) + γ .* u

end

```

to this

```julia

function f!(du, u, p, t)

    du .= α .* (Δ * Q * u) + ((-v_x' .* (divx * Q * u) + -v_y' .* (divy * Q * u)) ) + γ .* u

end

```

Following this change, v\_x and v\_y and even \alpha can all be matricies and I can functionally perform the dot product \vec{v} \cdot \nabla C since I can multiply each part of the gradient by the correct velocity component a la  
v\_x \frac{\partial C}{\partial x} \vec{i} + v\_y \frac{\partial C}{\partial y} \vec{j} + v\_z \frac{\partial C}{\partial z} \vec{k}.

So now then, the reason all of this was important to me was so that now I can have arbitrary flow fields (or more generally nonconstant properties that vary along the domain) being modeled which opens up the possibilities of the work I can do with this package. Thanks for the help.