# Is OPENFOAM a good idea for direction/operator splitting in parabolic equation

I would like to know if I should switch to OPENFOAM for my task. I work only with Cartesian grids, right now in 2d, rectangular domains only. If $\mathbf{w} = (w_1, w_2)$ and suppose I want to solve

$\partial_t(\mathbf{w}) = \nabla(\mathrm{div} \mathbf{w}) = \mathbf{f}$

subject to Dirichlet boundary condition and an initial condition.

What I want to do is discretize in the following way and solve:

$\dfrac{w_1^{(n+1)} - w_1^n}{\tau} - \partial_{xx}w_1^{(n+1)} = f_1 + \partial_{xy}w_2^n \\ \dfrac{w_2^{(n+1)} - w_2^n}{\tau} - \partial_{yy}w_2^{(n+1)} = f_1 + \partial_{yx}w_1^{(n+1)}$.

1. Need 2 tridiaogonal solves ( so Ax = b is trivial)
2. I have made my choice of time stepping, and also the discretization.
3. I might modify my time stepping and discretization later.
4. Right now, I wrote a program myself in C

Basically I want to avoid doing the routine tasks like organize the MAC grid, store source values at the cell centre, compute $\partial_{xy}w_2^n$ given $w_2^n$ etc.

Q1) Is OPENFOAM overkill for my task and be more of a nuisance to actually implement the scheme I want to try ?

Q2) Am I better off doing what I am doing now, having a working code in C, which means though that each time I want to try something new, I have to change my program ?

Q3) Is it possible to use direct methods instead of iterative methods to solve my system in OPENFOAM ?

• I know it's totally outside of the scope of the question, but you might want to reconsider taking a Helmholtz decomposition approach to solving that system. Grad-div has a lot of nasty null space behaviours otherwise. – origimbo Mar 23 '16 at 6:34
• @origimbo: Thanks for the idea. How would that work? $w = \nabla \phi + u$ with $u$ having zero divergence ? I dont see how I will have enough information to solve for $\phi$ and $u$. – me10240 Mar 23 '16 at 6:46
• Set $\mathbf{w} = \nabla \phi + \nabla \times \mathbf{A}$ where, since you're in 2D, you have a sensible gauge condition for $\mathbf{A}$ (i.e. non-zero in 3rd component only). You can solve for the initial conditions (and justify the evolution equations) by taking divergence and curl of $\mathbf{w}$ and inverting the resultant Laplacians. – origimbo Mar 23 '16 at 7:06