I’m solving equations describing density of some physical tissues so obviously negative values have no sense in real world however I get some during calculation (and it’s getting worse with each time step).
Have you faced such issue? What kind of clever techniques can be used not to influence resolution too much. I’ve initial and boundary condition, already decreased time step and no further clues.
Hello, can you post some code?
You can check the sign of your linear/bilinear term. Also check that you’re writing your boundary condition correctly (the sign of the normal may have an influence).
Where should I check the sign? In problem definition?
Don’t want to place the whole code (it’s long). I have system of 3 equations, below 2 most connected. All functions, macros, constants are defined. Below problem is solved in 5-10 steps but etaa become negative.
problem Prob([etai,etaa],[o1,o3],solver=LU)=
//1st
int2d(Th)(etai/dto1) -int2d(Th)(etaiold/dto1)
-int2d(Th)(Detai*(dx(etai)*dx(o1) + dy(etai)dy(o1)))
-int2d(Th)(wektoru(uv1, uv2)'grad(etai)o1)
+int2d(Th)(phi(x,y)etaio1)
//3rd
+int2d(Th)(etaa/dto3) -int2d(Th)(etaaold/dto3)
-int2d(Th)(Detaa(dx(etaa)dx(o3) + dy(etaa)dy(o3)))
-int2d(Th)(wektoru(uv1, uv2)'grad(etaa)o3)
-int2d(Th)(phi(x,y)etaio3)
+int2d(Th)(alphaetaao3)
+int2d(Th)(omegaaetaao3)
//boundary conditions
+on(C04, etai = w2, etaa = 0)
Can you upload the complete code? Probably the best thing to do is
work through the discretezation eqns and see how the problem
occurs with coarse time or space steps. Sign mistakes are common
too ( esp for me doing int by parts in head lol ). An alterantive is
to work with logs but that is not easy either.
Slightly simplified code below but the main parts are there. Maybe I’m making some silly mistake
border C01(t=0.06, 2.08){x=t; y=0.002;}
border C02(t=0.002, 1.002){x=2.08; y=t;}
border C03(t=2.08, 0.06){x=t; y=1.002;}
border C04(t=1.002, 0.002){x=0.06; y=t;}
int n = 100;
mesh Th = buildmesh(C01(n) + C02(n/2) + C03(n) + C04(n/2));
fespace Vh1(Th,P2);
// Parameters
real Dgamma = 1000;
real Detai = 20;
real Detaa = 40;
real gammahypo = 10;
real fi0 = 0.9;
real w1=60;
real w2=50 ;
real alpha=0.5;
real beta=0.85;
real gamma0=0.85;
real omegaa=0.03010299956;
real t=0, T=5, dt=1;
Vh1 etai, etaa, gamma,utest,
etaiold, etaaold, gammaold,
o1,o2,o3,
uv1 = 2,
uv2 = 2;
//initial conditions
etai = 0;
gamma = gamma0;
etaa = 0;
//functions
macro grad(ff) [dx(ff),dy(ff)] //EOM
macro wektoru(w1,w2) [w1,w2] //EOM
func real fi(real xx, real yy)
{if (gamma(xx,yy) < gammahypo)
{return fi0;}
else {return 0; }
}
problem Prob([etai,gamma, etaa],[o1,o2,o3],solver=LU)=
int2d(Th)(etai/dto1) -int2d(Th)(etaiold/dto1)
-int2d(Th)(Detai*(dx(etai)dx(o1) + dy(etai)dy(o1)))
-int2d(Th)(wektoru(uv1, uv2)'grad(etai)o1)
+int2d(Th)(fi(x,y)etaio1)
//2nd
+int2d(Th)(gamma/dto2) -int2d(Th)(gammaold/dto2)
-int2d(Th)(Dgamma(dx(gamma)dx(o2) + dy(gamma)dy(o2)))
-int2d(Th)(wektoru(uv1, uv2)'grad(gamma)o2)
+int2d(Th)(betagammao2)
//3rd
+int2d(Th)(etaa/dto3) -int2d(Th)(etaaold/dto3)
-int2d(Th)(Detaa(dx(etaa)dx(o3) + dy(etaa)dy(o3)))
-int2d(Th)(wektoru(uv1, uv2)'grad(etaa)o3)
-int2d(Th)(fi(x,y)etaio3)
+int2d(Th)(alphaetaao3)
+int2d(Th)(omegaaetaao3)
//boundary conditions
+on(C04, etai = w2, gamma = w1, etaa = 0)
;
for(real t=t;t<=T;t+=dt)
{
cout<<"===================================Solution in step t: "<<t<<endl;
Prob;
//plot(etai,value=1,fill=1,wait=0, cmm="solution etai at time: "+t);
//plot(gamma,value=1,fill=1,wait=0, cmm="solution gamma at time: "+t);
//plot(etaa,value=1,fill=1,wait=0, cmm="solution etaa at time: "+t);
etaiold=etai; etaaold=etaa; gammaold=gamma;
}
It is not so simple to find a numerical method such that the solution is positive.with finite element method.
- for P2 finite element no method exist to day.
- for P1 finite element it is generally ok with mass lump quadrature if the mesh is acute but it is hard to build a mesh with only acute triangle, but in general if the no acute triangles are not in bad place it works.
for mass lumping for P1 add
change
fespace Vh1(Th,P2);
...
int2d(Th)(etai/dto1) -int2d(Th)(etaiold/dto1)
in
fespace Vh1(Th,P1);
....
int2d(Th,qft=qf1pTlump)(etai/dto1) -int2d(Th,qft=qf1pTlump)(etaiold/dto1)
Thanks for answering.
I introduced suggested modifications, results changed but still getting negative values.