periodicity for concentration

fluid/fluid_problem_concentration.c

@@ -76,8 +76,8 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) {
     }
   }
   /*interface_term*/
-  for (int ii = 0; ii < mesh->n_interfaces; ++ii) {
-    const int *interface = &mesh->interfaces[ii*4];
+  for (int ii = 0; ii < mesh->n_interfaces_periodic; ++ii) {
+    const int *interface = &mesh->interfaces_periodic[ii*4];
     const int iel0 = interface[0];
     const int iel1 = interface[2];
     if (iel1 < 0) continue;
@@ -211,7 +211,6 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) {
   }
   static double I_a_old = 0;
   double I_a = fluid_problem_a_integ_volume(problem);
-  printf("After limiter: Volume = %.16g, Bnd = %.16g, diff=%.16g\n",I_a,I_bnd,-I_a_old+I_a-I_bnd);
   I_a_old = I_a;
   for (int iel = 0; iel < mesh->n_elements*N_N; ++iel) {
     if(fabs(problem->concentration[iel])<1e-16) problem->concentration[iel]=0;

fluid/mesh.c

@@ -66,7 +66,7 @@ static void sort_edge_nodes(int *n) {
 #endif
 }
 
-static void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *bnd_tags) {
+static void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *bnd_tags, int periodic, int *n_edgesp, int **edgesp) {
   int n_half_edges = m->n_elements*(DIMENSION+1)+n_bnd_elements;
   HalfEdge *halfedges = malloc(sizeof(HalfEdge)*n_half_edges);
   HalfEdge *he = halfedges; // current
@@ -91,6 +91,14 @@ static void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *
     he->pos = bnd_tags[i];
     he ++;
   }
+  if (periodic) {
+    for (int i = 0; i < n_half_edges; ++i){
+      he = &halfedges[i];
+      for (int k = 0; k < DIMENSION; ++k) {
+          he->n[k] = m->periodic_mapping[he->n[k]];
+      }
+    }
+  }
   qsort(halfedges,n_half_edges,sizeof(HalfEdge),hedgecmp);
   int n_edges = 0;
   for (int i = 0; i < n_half_edges; ++i) {
@@ -107,12 +115,12 @@ static void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *
       edges[(n_edges-1)*4+3] = -1;
     }
     else {
-      edges[(n_edges-1)*4+2] = halfedges[i].element;
-      edges[(n_edges-1)*4+3] = halfedges[i].pos;
+      if(edges[(n_edges-1)*4+2] == -1) {
+        edges[(n_edges-1)*4+2] = halfedges[i].element;
+        edges[(n_edges-1)*4+3] = halfedges[i].pos;
+      }
     }
   }
-  m->n_interfaces = n_edges;
-  m->interfaces = edges;
   free(halfedges);
   for (int i = 0; i < n_edges; ++i) {
     int *e = edges + i*4;
@@ -130,10 +138,17 @@ static void mesh_gen_edges(Mesh *m, int n_bnd_elements, int *bnd_elements, int *
       const int *e0 = &m->elements[e[0]*4];
       const int *e1 = &m->elements[e[2]*4];
       e[3] += 12;
-      while (e0[elbnd[e[1]][0]] != e1[elbnd[e[3]][k]]) e[3] += 4;
+      if (periodic) {
+        while (m->periodic_mapping[e0[elbnd[e[1]][0]]] != m->periodic_mapping[e1[elbnd[e[3]][k]]]) e[3] += 4;
+      }
+      else {
+        while (e0[elbnd[e[1]][0]] != e1[elbnd[e[3]][k]]) e[3] += 4;
+      }
 #endif
     }
   }
+  *n_edgesp = n_edges;
+  *edgesp = edges;
 }
 
 void mesh_free(Mesh *m)
@@ -170,7 +185,8 @@ Mesh *mesh_new_from_elements(int n_nodes, double *x, int n_elements, int *elemen
   for(int i = 0; i < m->n_nodes; ++i){
     m->periodic_mapping[i] = periodic ? periodic[i] : i;
   }
-  mesh_gen_edges(m, n_boundaries, boundaries, boundary_tags);
+  mesh_gen_edges(m, n_boundaries, boundaries, boundary_tags, 0, &m->n_interfaces, &m->interfaces);
+  mesh_gen_edges(m, n_boundaries, boundaries, boundary_tags, 1, &m->n_interfaces_periodic, &m->interfaces_periodic);
   return m;
 }
 

fluid/mesh.h

@@ -36,6 +36,8 @@ typedef struct {
   char **boundary_names;
   int n_interfaces;
   int *interfaces; // eid0, closure0, {(eid1, closure1),(-1,phys_id)}
+  int n_interfaces_periodic; // same as interfaces but taking into account the periodicity
+  int *interfaces_periodic; // eid0, closure0, {(eid1, closure1),(-1,phys_id)}
   int n_periodic;
   int *periodic_mapping; // n_nodes
 } Mesh;

testcases/periodicTest/2D/2fluid.py

+# MigFlow - Copyright (C) <2010-2020>
+# <Universite catholique de Louvain (UCL), Belgium
+#  Universite de Montpellier, France>
+#
+# List of the contributors to the development of MigFlow: see AUTHORS file.
+# Description and complete License: see LICENSE file.
+#
+# This program (MigFlow) is free software:
+# you can redistribute it and/or modify it under the terms of the GNU Lesser General
+# Public License as published by the Free Software Foundation, either version
+# 3 of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program (see COPYING and COPYING.LESSER files).  If not,
+# see <http://www.gnu.org/licenses/>.
+
+#!/usr/bin/env python
+from migflow import fluid as mbfluid
+from migflow import time_integration
+import numpy as np
+import os
+import subprocess
+import time
+import shutil
+dir_path = os.path.dirname(os.path.realpath(__file__))
+os.chdir(dir_path)
+# Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
+
+outputdir = "output-2fluid"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","0.5"])
+t = 0
+ii = 0
+
+
+#physical parameters
+g =  np.array([0.0,0])                                      # gravity
+rho = 1000                                      # fluid density
+nu = 1e-3                                   # kinematic viscosity
+mu = nu*rho                                     # dynamic viscosity
+tEnd = 1000                                     # final time
+
+#numerical parameters
+dt = 1                                         # time step
+outf = 10                                       # number of iterations between output files
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+
+#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
+fluid = mbfluid.FluidProblem(2,g,[nu*rho, nu*rho],[rho,rho],petsc_solver_type="-pc_type lu")
+fluid.load_msh("mesh.msh")
+fluid.set_wall_boundary("Bottom",velocity=[0.003,0],pressure=0)
+fluid.set_wall_boundary("Top",velocity=[0.003,0])
+fluid.solution()[:,0] = 0.003
+x = fluid.coordinates()
+fluid.set_concentration_cg(np.exp(-10*((x[:,0]-0.5)**2+(x[:,1]-0.5)**2)))
+
+
+ii = 0
+t = ii*dt
+
+#set initial_condition
+fluid.write_vtk(outputdir,ii,t)
+
+tic = time.time()
+while t < tEnd:
+    #Fluid solver
+    time_integration.iterate(fluid,None,dt)
+    t += dt
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        fluid.write_vtk(outputdir, ioutput, t)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))