Hi,
i try to write a FF++ code to resolve attached problem
I write the following code but it give me a false value (the values of solution must be positives and with my code, the solutions are negatives.
Thank you in advance for the help
There is my FF++ code:
macro model 2
macro typeFE() P1
verbosity=0;
bool delay=true;
//load "UMFPACK64"
// Parameters
real L=10.;
real T=50.;
real a1=1.e-4;
real a2=1.e-4;
real a3=1.e-4;
real sigma1=0.01;
real sigma2=0.01;
real sigma3=1.;
real sigma4=1.;
real D=1.e-3;
real b=1.e6;
real k1=1.;
real k2=2.;
real d=1.0;
int Nbx=1e4;
real dt=.01;
real e0=1.;
real f0=0.;
real m0=1.;
real n0=0.;
real s0=0.;
real v0=100.;
real q0=0.;
real x0=0.1;
real tau=0.;
int Nbt=1e2;
real Dx=L/Nbx; // calcul de dx
string viralloadtxt="viralload.txt";
string JNtxt="JN.txt";
string JFtxt="JF.txt";
ofstream foutJN(JNtxt);//créer un fichier log
ofstream foutJF(JFtxt);//créer un fichier log
load "msh3"
meshL Th=segment(Nbx,[x*L,0.]);
fespace Vh2(Th,typeFE);
Vh2 Id = (x<=x0) ? v0 : 0.;
IFMACRO(model,2)
Vh2 V,Q,F,E,M,N,VH,QH,FH,EH,MH,NH,V0=Id,F0=f0,E0=e0,N0=n0,M0=m0,Q0=q0;
macro SYSTEM
{
solve systemM(M,MH,solver=LU)=
int1d(Th)(M*MH/dt + a2*M*V0*MH)
- int1d(Th)(M0*MH/dt)
;
solve systemN(N,NH,solver=LU)=
int1d(Th)(N*NH/dt +sigma2*N*NH)
- int1d(Th)(N0*NH/dt + a2*M0*V0*NH)
;
solve systemE(E,EH,solver=LU)=
int1d(Th)(E*EH/dt + a1*E*V0*EH)
- int1d(Th)(E0*EH/dt)
;
solve systemF(F,FH,solver=LU)=
int1d(Th)(F*FH/dt + k1*s0*F*FH+sigma1*F*FH)
- int1d(Th)(F0*FH/dt + a1*E0*V0*FH)
;
solve systemQ(Q,QH,solver=LU)=
int1d(Th)(tau*(Q*QH/dt)+Q*QH)+int1d(Th)(D*dx(V0)*QH)
- int1d(Th)(tau*(Q0*QH/dt))
;
solve systemV(V,VH,solver=LU)=
int1d(Th)(V*VH/dt+sigma3*V*VH)+int1d(Th)(dx(Q0)*VH)
- int1d(Th)(V0*VH/dt+b*F0*VH)
;
V0=V; F0=F; E0=E; N0=N; M0=M; Q0=Q;
real JF=int1d(Th)(F0);
real JN=int1d(Th)(N0);
s0=(a3*JF+d*JN)/(k2*JF+sigma4);
}//
ENDIFMACRO
int op=0;
real ViralLoad;
real JN;
real JF;
int it,itsave;
int Thnv = Th.nv;
fespace Vh(Th,P1);
Vh Vs=x,Fs=x,Es=x, Ms=x, Ns=x;
real[int,int] Esave(Thnv,3),Vsave(Thnv,3),Fsave(Thnv,3), Nsave(Thnv,3),Msave(Thnv,3);
Nsave(:,itsave)=Ns[];
Msave(:,itsave)=Ms[];
Esave(:,itsave)=Es[];
Fsave(:,itsave)=Fs[];
Vsave(:,itsave)=Vs[];
for (real t=0;t<T;t+=dt){
SYSTEM;
if(!(it%(max(1.,1./dt/100.)))){
ViralLoad = int1d(Th)(V);
ofstream foutviralload(viralloadtxt,append);// append: pour ajouter des lignes
foutviralload << t << " " <<ViralLoad << endl;
cout << t << " " << ViralLoad<< endl;
}
//if(!(it%(max(1.,1./dt/100.)))){
// JN = int1d(Th)(N);
// ofstream foutJN(JNtxt,append);// append: pour ajouter des lignes
// foutJN << t << " " <<JN << endl;
//cout << t << " " << JN<< endl;
// }
//if(!(it%(max(1.,1./dt/100.)))){
// JF = int1d(Th)(F);
// ofstream foutJF(JFtxt,append);// append: pour ajouter des lignes
// foutJF << t << " " <<JF << endl;
//cout << t << " " << JF<< endl;
// }
if(!(it%max(1.,1./dt))){
Vs=V;
Fs=F;
Es=E;
Ns=N;
Ms=M;
itsave++;
Nsave(:,itsave)=Ns[];
Msave(:,itsave)=Ms[];
Esave(:,itsave)=Es[];
Fsave(:,itsave)=Fs[];
Vsave(:,itsave)=Vs[];
Nsave.resize(Thnv,itsave+3);
Msave.resize(Thnv,itsave+3);
Esave.resize(Thnv,itsave+3);
Fsave.resize(Thnv,itsave+3);
Vsave.resize(Thnv,itsave+3);
}
op=1;
it++;
}
//ofstream foutE("E.txt");
//ofstream foutM("M.txt");
//ofstream foutN("N.txt");
//ofstream foutF("F.txt");
//ofstream foutV("v.txt");
//for(int m=0;m<Th.nv;m++) {
// for(int n=0;n<itsave;n++) {
// foutV << Vsave(m,n) << " ";
// foutF << Fsave(m,n) << " ";
// foutM << Msave(m,n) << " ";
// foutE << Esave(m,n) << " ";
// foutN << Nsave(m,n) << " ";
//}
//foutV<<endl;
//foutF<<endl;
//foutN<<endl;
//foutE<<endl;
//foutM<<endl;
//}