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

20 days 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)

a = dot(u, grad(v))*dx 
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
auto G_Test = std::make_shared<jGradVel::Form_a::TestSpace>(mesh);//scalar

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);
jGradVel::BilinearForm  a_Vel(G_Trial,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!


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 20 days 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 20 days ago by jwinkle  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »