### How to solve div(v) = alpha (a modified incompressibility condition for "cell growth")

119
views
0
12 weeks ago by
I am trying to solve for a (2D vector) velocity field that is governed by a modified incompressibility condition:
$\nabla \cdot v = \alpha$
where $\alpha$ is some scalar constant (a growth rate) and I will give a suitable Dirichlet BC at an interior region, and homogeneous Neumann on the perimeter.   I expect the velocity solution to be a linear function (in magnitude) of the distance from the interior fixed BC (e.g., in 1-D solution is v = $\alpha x + v_0$, where $v_0$ is the Dirichlet BC at some point $x_0$ )

My .ufl code attempt (using C++) is (in a file named "jGradVel.ufl"):

element 	= FiniteElement("Lagrange", triangle, 2)
velement 	= VectorElement("Lagrange", triangle, 2)

alpha 		= Constant(triangle)
v  			= TestFunction(element)
u 			= TrialFunction(velement)

L = inner(alpha,v)*dx
​

However, when I compile this I do not have defined a "FunctionSpace" and I am reverting to using:

auto G_Trial = std::make_shared<jGradVel::Form_a::TrialSpace>(mesh);//vector


This must already be a problem.  At the bottom of "jGradVel.h" header file formed after ffc compiling the .ufl file, there are only the following typedefs (usually is defined "FunctionSpace" here, but it is not generated):

...
// Class typedefs
typedef Form_a BilinearForm;
typedef MultiMeshForm_a MultiMeshBilinearForm;
typedef Form_a JacobianForm;
typedef MultiMeshForm_a MultiMeshJacobianForm;
typedef Form_L LinearForm;
typedef MultiMeshForm_L MultiMeshLinearForm;
typedef Form_L ResidualForm;
typedef MultiMeshForm_L MultiMeshResidualForm;

My goal is to solve for the vector velocity field that satisfies the modified compressibility condition, and I expect to use:

jGradVel::LinearForm    L_Vel(G_Test);
auto vel = std::make_shared<Function>(G_Trial);

solve(a_Vel == L_Vel, *vel);


Where 'vel" is my velocity field solution that I can then plot/write to file.

This can compile, but then on running issues an error:

*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 62 (Invalid argument).
*** Where:   This error was encountered inside /tmp/src/dolfin/dolfin/la/PETScKrylovSolver.cpp.

but I assume I must have a semantic error somewhere else in the code in setting up the problem correctly in variational form.

Any ideas in Fenics-Sphere how to help to properly set-up this simple problem?  Thanks for any help!

--JW

Community: FEniCS Project
Hi, please note that your system matrix is rectangular. Moreover, the continuous problem is ill-posed as for any v which is a curl the left hand side equals zero.
written 12 weeks ago by micro
@micro Yes, you are correct that the problem will be ill-posed.  Can this system be reduced to a 1-D problem on the rectangle by, say imposing $v = 0$ on the "south" boundary,  $v_x = 0$ on the east/west boundaries, and solve $v_y = \alpha$?  This should be a well-posed 1-D problem (but treating $v$ still as a 2D vector in fenics) with solution $v = v(y) = \alpha y$ (with $v = 0$ at $y = 0$ ).

I am wanting to make sure I understand how to properly form the problem in Fenics (and understanding the continuous problem requires special consideration of the BC to be well-posed).  In particular, in the previous problem I do not understand why a "FunctionSpace" was not typedef-defined in the .h file (after compiling via ffc the .ufl file) and whether the .ufl file per se  was ill-formed.
written 12 weeks ago by jwinkle