# What's wrong with initial conditions?

I want to solve the nonlinear PDE for the anisotropic fluid flow:

$-\rho*\partial{v_i}/\partial{t} + (C_{ijkl}v_{k,l})_{,j}-p_{,j} = 0$

(The nonlinear term $-v_jv_{i,j}$ will be added later)

Here is the part of the code describing PDE (differential operator, boundary and initial conditions):

Ui[t, x, y, z] = {u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]};
coord = {x, y, z};

(*Operator*)
op = {
-rho*D[Ui[t, x, y, z][[1]], t]+Sum[D[Sum[Sum[Cijkl[[1, j, k, l]]*D[Ui[t, x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[t, x, y, z], x],
-rho*D[Ui[t, x, y, z][[2]], t]+Sum[D[Sum[Sum[Cijkl[[2, j, k, l]]*D[Ui[t, x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[t, x, y, z], y],
-rho*D[Ui[t, x, y, z][[3]], t]+Sum[D[Sum[Sum[Cijkl[[3, j, k, l]]*D[Ui[t, x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[t, x, y, z], z],
Sum[D[Ui[t, x, y, z][[k]], coord[[k]]], {k, 3}]
};

(*Boundary conditions*)
bcs = {
DirichletCondition[p[t, x, y, z] == 100, x == 0],
DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, w[t, x, y, z] == 0}, 0 < x < 8],
DirichletCondition[p[t, x, y, z] == 0, x == 8]
};

(*Initial conditions*)
ics = {
u[0, x, y, z] == 0,
v[0, x, y, z] == 0,
w[0, x, y, z] == 0
};

(*Solve*)
{xVel, yVel, zVel, pressure} = NDSolveValue[{op == {0, 0, 0, 0}, ics, bcs}, {u, v, w, p}, {t, 0, 1}, {x, y, z} \[Element] mesh, Method -> {"MethodOfLines","SpatialDiscretization" -> {"FiniteElement","InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}}}];


When I try to solve it, I get the following error message:

NDSolveValue::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable.

I cant understand what's wrong with the conditions? May there be a syntax error somewhere?

I suppose that boundary conditions are ok, because the stationary version of this equation $(C_{ijkl}v_{k,l})_{,j}-p_{,j} = 0$ is well solved with the following code:

Ui[x, y, z] = {u[x, y, z], v[x, y, z], w[x, y, z]};
coord = {x, y, z};

op = {
Sum[D[Sum[Sum[Cijkl[[1, j, k, l]]*D[Ui[x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[x, y, z], x](*-Sum[Ui[x,y,z][[j]]*D[Ui[x,y,z][[1]],coord[[j]]], {j,3}]*),
Sum[D[Sum[Sum[Cijkl[[2, j, k, l]]*D[Ui[x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[x, y, z], y](*-Sum[Ui[x,y,z][[j]]*D[Ui[x,y,z][[2]],coord[[j]]], {j,3}]*),
Sum[D[Sum[Sum[Cijkl[[3, j, k, l]]*D[Ui[x, y, z][[k]], coord[[l]]], {l,3}], {k, 3}], coord[[j]]], {j, 3}] - D[p[x, y, z], z](*-Sum[Ui[x,y,z][[j]]*D[Ui[x,y,z][[3]],coord[[j]]], {j,3}]*),
Sum[D[Ui[x, y, z][[k]], coord[[k]]], {k, 3}]
};

bcs = {
DirichletCondition[p[x, y, z] == 100, x == 0],
DirichletCondition[{u[x, y, z] == 0, v[x, y, z] == 0, w[x, y, z] == 0}, 0 < x < 8],
DirichletCondition[p[x, y, z] == 0, x == 8]
};

{xVel, yVel, zVel, pressure} = NDSolveValue[{op == {0,0,0,0}, bcs}, {u,v,w,p}, {x,y,z} \[Element] mesh, Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,w-> 2,p->1}}];


So, is there any mistake in the ICs or solve section?

P.S.: there are 2 versions of whole file for your inspection and execution:

1) Steady flow (it works, evaluate whole notebook): https://1drv.ms/u/s!AvmYXV0MDC0Jh4ZFYTSCBZ9vspyjKw

2) Unsteady flow (doesn't work): https://1drv.ms/u/s!AvmYXV0MDC0Jh4ZGjmCorU_DymHq3A

UPD: ok, there is the whole code with the constants and mesh definition:

Needs["NDSolveFEM"]

a2 = { {1, 0, 0}, {0, 0, 0}, {0, 0, 0} };

a4 = Table[a2[[i,j]]*a2[[k,l]],{i,3},{j,3},{k,3},{l,3}];
Np = 0;
Ns = 0;
mu = 1;
rho = 1;

Cijklt = Table[2mu*(KroneckerDelta[i,k]*KroneckerDelta[j,l]+Np*a4[[i,j,k,l]]+Ns*(a2[[l,j]]*KroneckerDelta[i,k]+a2[[i,k]]*KroneckerDelta[j,l])), {i,3},{j,3},{k,3},{l,3}];
Delta = Table[0.5*(KroneckerDelta[i,k]*KroneckerDelta[j,l]+KroneckerDelta[i,l]*KroneckerDelta[j,k]),{i,3},{j,3},{k,3},{l,3}];

Cijkl = Table[Sum[Sum[Cijklt[[i,j,m,n]]*Delta[[m,n,k,l]],{n,3}],{m,3}],{i,3},{j,3},{k,3},{l,3}];

omega = RegionUnion[Cylinder[{{0,0,0},{4,0,0}},2],Cylinder[{{4,0,0},{8,0,0}},1]];
mesh = ToElementMesh[omega, MaxCellMeasure-> 0.005];
RegionPlot3D[omega, PlotPoints->50];
mesh["Wireframe"];

Ui[t,x,y,z] = {u[t,x,y,z], v[t,x,y,z],w[t,x,y,z]};
coord = {x,y,z};

op = {
-rho*D[Ui[t,x,y,z][[1]], t]+Sum[D[Sum[Sum[Cijkl[[1,j,k,l]]*D[Ui[t,x,y,z][[k]],coord[[l]]],{l,3}],{k,3}], coord[[j]]],{j,3}]-D[p[t,x,y,z],x](*-Sum[Ui[t,x,y,z][[j]]*D[Ui[t,x,y,z][[1]], coord[[j]]], {j,3}]*),
-rho*D[Ui[t,x,y,z][[2]], t]+Sum[D[Sum[Sum[Cijkl[[2,j,k,l]]*D[Ui[t,x,y,z][[k]],coord[[l]]],{l,3}],{k,3}], coord[[j]]],{j,3}]-D[p[t,x,y,z],y](*-Sum[Ui[t,x,y,z][[j]]*D[Ui[t,x,y,z][[2]], coord[[j]]], {j,3}]*),
-rho*D[Ui[t,x,y,z][[3]], t]+Sum[D[Sum[Sum[Cijkl[[3,j,k,l]]*D[Ui[t,x,y,z][[k]],coord[[l]]],{l,3}],{k,3}], coord[[j]]],{j,3}]-D[p[t,x,y,z],z](*-Sum[Ui[t,x,y,z][[j]]*D[Ui[t,x,y,z][[3]], coord[[j]]], {j,3}]*),
Sum[D[Ui[t,x,y,z][[k]],coord[[k]]],{k,3}]
};

bcs={
DirichletCondition[p[t,x,y,z]==100,x== 0 ],
DirichletCondition[{u[t,x,y,z]==0,v[t,x,y,z]==0,w[t,x,y,z]==0},0<x<8],
DirichletCondition[p[t,x,y,z]==0,x== 8]
};

ics = {
u[0,x,y,z]==0,
v[0,x,y,z]==0,
w[0,x,y,z]==0
};

{xVel,yVel, zVel,pressure}=NDSolveValue[{op=={0,0,0,0}, ics,bcs},{u,v,w,p},{t, 0,1},{x,y,z}\[Element]mesh,Method->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","InterpolationOrder"->{u->2,v->2,w-> 2,p->1}}}];

• you should give Cijkl if you want folks to be able to try this. Commented May 7, 2018 at 15:04
• They are just numbers. You can find them in attached files. The whole code is too long to copy it into the post. The problem does not realy depend on which numbers they are. Commented May 7, 2018 at 15:48
• The thing is that it's better if you put dummy values by yourself in the original question rather than asking every user who is trying to help you to choose dummy values themselves. In short: help people help you :) Commented May 7, 2018 at 20:05
• Is it better to type the whole huge code in the question body or just attach a fite to press a 1 button to execute? =) Commented May 8, 2018 at 9:08
• UPD: ok, the whole code has been added, as you asked Commented May 8, 2018 at 10:31

It is a little bit unexpected, but NDSolveValue requires initial conditions for all unknown functions of PDE. Even those they have no time derivative.

There is the pressure function p[t,x,y,z] in equations. Usualy, the initial conditions for Navier-Stokes-like equations are the velocity field but Mathematica requires the pressure initial condition too.

For example, there is the working version of initial conditions that allows you to avoid error message:

ics = {
u[0, x, y, z] == 0,
v[0, x, y, z] == 0,
w[0, x, y, z] == 0,
p[0, x, y, z] == 0
};