validation : check_norm

fluid/fluid_problem.c

@@ -199,14 +199,14 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
   }
   double h = problem->element_size[eid];
   if (wbnd->type == BND_WALL && pid >= 0) {
-      f0[P] = -(s[P]-data[pid])/h*problem->taup[eid];
+      f0[P] = -(p-data[pid])/h*problem->taup[eid];
       f00[P*n_fields+P] += -1/h*problem->taup[eid];
   }
   if(wbnd->type == BND_OPEN) {
+    f0[P] = unext;
     for (int i = 0; i < D; ++i) {
       f00[P*n_fields+U+i] += vid<0?n[i]:0;
     }
-    f0[P] = unext;
   }
   for (int id = 0; id < D; ++id) {
     f0[U+id] = c*(pid<0?0:data[pid]-p)*n[id];
@@ -216,7 +216,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
   double sigma = problem->ip_factor*(1+D)/(D*h)*mu*N_N;
   for (int iD = 0; iD < D; ++iD) {
     for (int jD = 0; jD < D; ++jD)
-      c_du_o_c[iD][jD] = ds[U+iD][jD] -u[iD]/c*dc[jD];
+      c_du_o_c[iD][jD] = ds[U+iD][jD]-u[iD]/c*dc[jD];
   }
   if(wbnd->compute_viscous_term == 1){
     for (int id = 0; id < D; ++id) {
@@ -371,15 +371,6 @@ static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, dou
   }
   *dfdu = gammastar/c;
   *dfddp = gammastar*dt/mass*vol;//-vol;
-  if(ip == 2){
-    for(int i = 0; i < D; ++i){
-      printf("u : %14.7f\n", u[i]);
-    }
-    printf("c : %14.7f\n", c);
-    printf("f : %14.7f\n", f[0]);
-    printf("dfdu : %14.7f\n", *dfdu);
-    printf("dfddp : %14.7f\n", *dfddp);
-  }
 }
 
 static void fluid_problem_f(const FluidProblem *problem, const double *sol, double dsol[][D], const double *solold, const double dsolold[][D], const double* mesh_velocity, const double c, const double *dc, const double *bf, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double f1[][D], double f00[], double f10[][D], double f01[][D], double f11[][D][D]) {

validation/cavity/cavity.py

@@ -28,7 +28,6 @@ import os
 import subprocess
 import time
 import shutil
-import random
 import unittest
 #Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
 
@@ -91,7 +90,7 @@ class Poiseuille(unittest.TestCase) :
         tic = time.time()
         while ii < 1000 : 
             #Fluid solver
-            time_integration.iterate(fluid,None,dt)
+            time_integration.iterate(fluid,None,dt, check_residual_norm=1e-4)
             t += dt
             #Output files writting
             if ii %outf == 0 :

validation/darcy/darcy.py

@@ -28,7 +28,6 @@ import numpy as np
 import os
 import subprocess
 import time
-import shutil
 import random
 import unittest
 
@@ -124,7 +123,7 @@ class Darcy(unittest.TestCase) :
             #Computation loop
             while t < tEnd :
                 #Fluid solver
-                time_integration.iterate(fluid,p,dt,fixed_grains=True)
+                time_integration.iterate(fluid,p,dt,fixed_grains=True,check_residual_norm=1e-4)
                 t += dt
                 #Output files writting
                 if ii %outf == 0 :

validation/dieSwell/dieSwell.py

@@ -27,8 +27,6 @@ import numpy as np
 import os
 import subprocess
 import time
-import shutil
-import random
 import unittest
 
 class DieSwell(unittest.TestCase):
@@ -92,8 +90,8 @@ class DieSwell(unittest.TestCase):
         # Set the mesh geometry for the fluid computation
         fluid.load_msh("mesh.msh")
         fluid.set_open_boundary("Left", velocity = [lambda x : 2*(1-x[:,1]**2),0])
-        fluid.set_strong_boundary("Left",0, lambda x : 2*(1-x[:,1]**2))
-        fluid.set_strong_boundary("Left",1,0)
+        # fluid.set_strong_boundary("Left",0, lambda x : 2*(1-x[:,1]**2))
+        # fluid.set_strong_boundary("Left",1,0)
         fluid.set_wall_boundary("TopLeft", velocity = [0,0])
         fluid.set_open_boundary("TopRight", pressure = apply_tension_weak, compute_viscous_term = 0)
         fluid.set_symmetry_boundary("BottomLeft")
@@ -112,7 +110,8 @@ class DieSwell(unittest.TestCase):
         #
         swelling = 0
         while t < tEnd :
-            time_integration.iterate(fluid,None,dt)
+            # time_integration.iterate(fluid,None,dt)
+            fluid.implicit_euler(dt, check_residual_norm=5e-4)
             if t >= 50.0:
                 fs.iterate(dt)
                 swelling = fluid.coordinates()[fs.fs_nodes[-1],1]/H

validation/poiseuille/poiseuille.py

@@ -27,7 +27,6 @@ import os
 import subprocess
 import time
 import shutil
-import random
 import unittest
 
 class Poiseuille(unittest.TestCase) :
@@ -87,7 +86,7 @@ class Poiseuille(unittest.TestCase) :
         tic = time.time()
         while ii < 100 : 
             #Fluid solver
-            time_integration.iterate(fluid,None,dt)
+            time_integration.iterate(fluid, None, dt, check_residual_norm=1e-4)
             t += dt
             #Output files writting
             if ii %outf == 0 :

validation/weight/weight.py

@@ -26,10 +26,6 @@ from migflow import time_integration
 
 import numpy as np
 import os
-import subprocess
-import time
-import shutil
-import random
 import unittest
 def genInitialPosition(filename, rmin, rmax, H, lx, ly, rhop) :
     """Set all the particles centre positions and create the particles objects to add in the computing structure
@@ -113,7 +109,6 @@ class Weight(unittest.TestCase) :
         # COMPUTATION LOOP
         #
         forces = g*p.mass()
-        print("coucou")
         def accumulate(bnd_forces,n_divide) :
           bnd_forces += total_boundary_force(p)/n_divide
         while t < tEnd :