Merge branch 'p2p1' of https://git.immc.ucl.ac.be/fluidparticles/migflow into p2p1

fluid/fluid_problem.c

@@ -151,7 +151,7 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
       f1[P][d] += pspg*f[d];
       for (int e = 0; e < D; ++e) {
         double supg = drag_in_stab?sold[U+e]*problem->taup[iel]:0;
-        f1[U+d][d] += supg*f[d];
+        f1[U+d][e] += supg*f[d];
       }
     }
     fe_fields_add_to_local_vector(problem->fields, f0, f1, sf, dsf, 1, local_vector);
@@ -471,6 +471,7 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
       f11[(U+i)*n_fields+(U+j)][i][j] += lsic;
     }
     f1[P][i] = -u[i];
+    
     f10[P*n_fields+U+i][i] += -1;
     // PSPG
     double pspg = taup/rho;
@@ -525,12 +526,12 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
         double s[fields->n], sold[fields->n], ds[fields->n][D], dsold[fields->n][D];
         double sf[fields->local_size], dsf[fields->local_size][D];
         fe_fields_f(fields, xiel, sf);
-        fe_fields_df(fields, QP[iqp], dxidx, dsf);
+        fe_fields_df(fields, xiel, dxidx, dsf);
         fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, problem->solution, s, ds);
         fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, solution_old, sold, dsold);
         double dc[D], c, sfporosity[problem->field_porosity->local_size], dsfporosity[problem->field_porosity->local_size][D] ;
-        fe_fields_f(problem->field_porosity, QP[iqp], sfporosity);
-        fe_fields_df(problem->field_porosity, QP[iqp], dxidx, dsfporosity);
+        fe_fields_f(problem->field_porosity, xiel, sfporosity);
+        fe_fields_df(problem->field_porosity, xiel, dxidx, dsfporosity);
         fe_fields_eval_grad_sf(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
         double f0[fields->n], f00[fields->n*fields->n],f01[fields->n*fields->n][D];
         for (int i = 0; i < fields->n; ++i) {
@@ -547,12 +548,12 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
         double a = 0;
         if (problem->n_fluids==2) {
           double sfconcentration[problem->field_concentration->local_size];
-          fe_fields_f(problem->field_concentration, QP[iqp], sfconcentration);
+          fe_fields_f(problem->field_concentration, xiel, sfconcentration);
           fe_fields_eval_sf(problem->field_concentration, mesh, iel, sfconcentration, problem->concentration, &a);
         }
         const FEElement *mesh_element = problem->mesh->element;
         double meshf[mesh_element->nlocal];
-        mesh_element->f(QP[iqp], meshf);
+        mesh_element->f(xiel, meshf);
         double umesh[D] = {0};
         for (int i = 0; i < mesh_element->nlocal; ++i) {
           for (int j=0; j<D; ++j) {
@@ -751,19 +752,20 @@ static void fluid_problem_volume(FluidProblem *problem, const double *solution_o
       fe_fields_df(problem->field_porosity, QP[iqp], dxidx, dsfporosity);
       double c, dc[D], cold;
       fe_fields_eval_grad_sf(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
-      fe_fields_eval_sf(problem->field_porosity, mesh, iel, sfporosity, problem->porosity, &cold);
+      fe_fields_eval_sf(problem->field_porosity, mesh, iel, sfporosity, problem->oldporosity, &cold);
       double a = 0;
       if (problem->n_fluids==2) {
         double sfconcentration[problem->field_concentration->local_size];
         fe_fields_f(problem->field_concentration, QP[iqp], sfconcentration);
         fe_fields_eval_sf(problem->field_concentration, mesh, iel, sfconcentration, problem->concentration, &a);
       }
-      double bf[D];
+      double bf[D]={0};
       const FEElement *mesh_element = problem->mesh->element;
       double meshf[mesh_element->nlocal];
       mesh_element->f(QP[iqp], meshf);
-      for (int i = 0; i < mesh_element->nlocal; ++i) {
-        for (int j=0; j<D; ++j) {
+      for (int j=0; j<D; ++j) {
+        mesh_velocity[j] = 0;
+        for (int i = 0; i < mesh_element->nlocal; ++i) {
           mesh_velocity[j] += meshf[i]*problem->mesh_velocity[el[i]*D+j];
           bf[j] += meshf[i]*problem->bulk_force[el[i]*D+j];
         }

python/fluid.py

@@ -476,7 +476,6 @@ class FluidProblem :
             norm0 = np.linalg.norm(rhs)
         self.solution()[:] -= sys.solve(rhs).reshape([-1,self.n_fields()])[periodic_map]
         if check_residual_norm > 0 :
-            sys.zero_matrix()
             rhs[:] = 0
             self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
             sys.local_to_global(localv,localm,rhs)