stupid mistake

fluid/fluid_problem.c

@@ -366,7 +366,7 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
   *dfddp = gammastar*dt/mass*vol;//-vol;
 }
 
-static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double *bf, const double dcdt, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {//rempalcer cold par dcdt
+static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double *bf, const double dcdt, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
   double p = sol[P];
   double taup = problem->taup[iel];
   double tauc = problem->tauc[iel];
@@ -810,7 +810,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
         f11[i] = 0;
       }
       double c = 0, a = 0;
-      double dcdt = 0, aold = 0; //cold par dcdt
+      double dcdt = 0, aold = 0;
       for (int i = 0; i < N_SF; ++i) {
         c += problem->porosity[el[i]]*phi[i];
         dcdt += problem->dcdt[el[i]]*phi[i];
@@ -833,7 +833,7 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
       const double jw = QW[iqp]*detj;
       problem->kinetic_energy += rho*sqrt(s[U]*s[U]+s[U+1]*s[U+1])*sqrt(s[U]*s[U]+s[U+1]*s[U+1])/2*jw;
 
-      fluid_problem_f(problem,s,ds,sold,dsold,c,dc,bf,dcdt,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);//dcdt a la palce de cold
+      fluid_problem_f(problem,s,ds,sold,dsold,c,dc,bf,dcdt,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
       for (int ifield = 0; ifield < n_fields; ++ifield) {
         for (int iphi = 0; iphi < N_SF; ++iphi){
           local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
@@ -1603,7 +1603,7 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
   problem->dcdt = (double*)malloc(mesh->n_nodes*sizeof(double));
   for(int i = 0; i < mesh->n_nodes; ++i){
     problem->porosity[i] = 1.;
-    problem->dcdt[i] = 1.;
+    problem->dcdt[i] = 0.;
   }
 
   free(problem->boundary_force);