reload

src/fluid_problem_io.c

@@ -28,6 +28,8 @@
 
 int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
 {
+  const int n_fields = fluid_problem_n_fields(problem);
+
   MbXmlReader *r = mb_xml_reader_create(filename);
   MbXmlElement fileelement,gridelement,pieceelement;
   mb_xml_read_element(r,&fileelement);
@@ -76,9 +78,11 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
       }
     }
     else if (!strcasecmp(e.name,"PointData")) {
-      if (problem->solution)
+      if (problem->solution){
         free(problem->solution);
-      problem->solution = (double*)malloc(sizeof(double)*(DIMENSION+1)*m->n_nodes);
+        problem->solution = (double*)malloc(sizeof(double)*n_fields*m->n_nodes);
+    }
+      
       while(mb_read_vtk_data_array(r,&a) != 0){
         if(!strcasecmp(a.name,"Porosity")) {
           free(problem->porosity);
@@ -87,15 +91,20 @@ int fluid_problem_import_vtk(FluidProblem *problem, const char *filename)
         else if(!strcasecmp(a.name,"Pressure")) {
           double *p = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
           for (int i = 0; i < m->n_nodes; ++i)
-            problem->solution[(DIMENSION+1)*i+DIMENSION] = p[i];
+            problem->solution[n_fields*i+DIMENSION] = p[i];
           free(p);
         }
         else if(!strcasecmp(a.name,"Velocity")) {
           double *v = mb_vtk_data_array_to_double(&a,m->n_nodes,3);
           for (int i = 0; i < m->n_nodes; ++i)
             for (int d = 0; d < DIMENSION; ++d)
-              problem->solution[(DIMENSION+1)*d+DIMENSION] = v[i*DIMENSION+d];
+              problem->solution[n_fields*i+d] = v[i*DIMENSION+d];
           free(v);
+        }else if(!strcasecmp(a.name,"Concentration")) {
+          double *c = mb_vtk_data_array_to_double(&a,m->n_nodes,1);
+          for (int i = 0; i < m->n_nodes; ++i)
+            problem->solution[n_fields*i+DIMENSION+1] = c[i];
+          free(c);
         }
         mb_vtk_data_array_destroy(&a);
       }

testcases/injectionWGrains/inject2fGrainsReload.py

+# Marblesbag - Copyright (C) <2010-2018>
+# <Universite catholique de Louvain (UCL), Belgium
+#  Universite de Montpellier, France>
+# 	
+# List of the contributors to the development of Marblesbag: see AUTHORS file.
+# Description and complete License: see LICENSE file.
+# 	
+# This program (Marblesbag) 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 marblesbag import fluid as fluid
+from marblesbag import scontact2
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+
+
+
+outputdir = "outputquiplante"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+
+
+
+#physical parameters
+alpha = 0*np.pi/4.
+g =  -9.81*np.cos(alpha)                                      # gravity
+print(g)
+rho0 = 1.117*10                                      # fluid density
+rho1 = 1000
+nu1 = 1.57e-5                                   # kinematic viscosity
+nu0 = 1e-6
+tEnd = 10                                     # final time
+r = 5e-4
+N = 10000
+lx = .222
+ly = .14
+rhop = 2500
+#numerical parameters
+lcmin = .1                                  # mesh size
+dt = .001                                       # time step
+iReload = 30
+
+
+shutil.copy("mesh.msh", outputdir +"/mesh.msh")
+
+p = scontact2.ParticleProblem()
+p.read_vtk(outputdir,iReload)
+outf = 1                                       # number of iterations between output files
+ii = 0
+t = 0
+
+
+def outerBndV(x) :
+    return max(0.1*t,0.1)
+
+strong_boundaries = [("Top",2,2,0.),("Injection",1,1,outerBndV),("Injection",3,3,1),("Injection",0,0,0)
+#        ("Top",3,3,1),("Lateral",3,3,1),("Bottom",3,3,1)
+        ]#,("Lateral",0,0,0),("Bottom",1,1,0)]
+fluid = fluid.fluid_problem("mesh.msh",g,[nu0*rho0,nu1*rho1],[rho0,rho1],strong_boundaries,1)
+fluid.set_weak_boundary("Bottom","Wall")
+fluid.set_weak_boundary("Lateral","Wall")
+fluid.set_weak_boundary("Top","Outflow")
+fluid.set_weak_boundary("Injection","Inflow")
+
+fluid.import_vtk(outputdir+"/fluid_%05d.vtu"%iReload)
+
+fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+
+#set initial_condition
+#fluid.solution()[:,3] = 1
+
+
+ii = iReload*outf+1
+t = ii*dt
+tic = time.clock()
+while t < tEnd : 
+#Adaptation of the mesh. Args are minimal mesh radius, maximal mesh radius and number of elements
+    if (ii%5==0 and ii != 0):
+       #fluid.adapt_mesh(1e-3,1e-4,50000)
+       fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    else :
+       fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
+    #Fluid solver
+    fluid.implicit_euler(dt)
+    forces = fluid.compute_node_force(dt)
+    #Computation of the new velocities
+    vn = p.velocity() + forces * dt / p.mass()
+    vmax = np.max(np.hypot(vn[:, 0], vn[:, 1]))
+    #number of sub time step
+    nsub = max(1, int(np.ceil((vmax * dt * 4)/min(p.r()))))
+    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(p.r()))
+    #Contact solver
+    for i in range(nsub) :
+        p.iterate(dt/nsub, forces)  
+    t += dt
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        p.write_vtk(outputdir, ioutput, t)
+        fluid.export_vtk(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.clock() - tic))