Merge branch 'master' into FOM_dev

analysis/morphology.cpp

@@ -702,10 +702,12 @@ double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
 
 		if (rank == 0) printf("     delta=%f, growth=%f, max. displacement = %f \n",morph_delta, GrowthEstimate, MAX_DISPLACEMENT);
 		// Now adjust morph_delta
+		if (fabs(GrowthEstimate - GrowthPrevious) > 0.0) {
 			double step_size = (TargetGrowth - GrowthEstimate)*(morph_delta - morph_delta_previous) / (GrowthEstimate - GrowthPrevious);
 			GrowthPrevious = GrowthEstimate;
 			morph_delta_previous = morph_delta;
 			morph_delta += step_size;
+		}
 		if (morph_delta / morph_delta_previous > 2.0 ) morph_delta = morph_delta_previous*2.0;
 
 		//MAX_DISPLACEMENT *= max(TargetGrowth/GrowthEstimate,1.25);

analysis/runAnalysis.cpp

@@ -706,6 +706,139 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> input_db, const RankInfoStru
     }
 
 
+    // Initialize the comms
+    for ( int i = 0; i < 1024; i++ )
+        d_comm_used[i] = false;
+    // Initialize the threads
+    int N_threads = db->getWithDefault<int>( "N_threads", 4 );
+    auto method   = db->getWithDefault<std::string>( "load_balance", "default" );
+    createThreads( method, N_threads );
+}
+
+runAnalysis::runAnalysis(  ScaLBL_ColorModel &ColorModel)
+/*	std::shared_ptr<Database> input_db, const RankInfoStruct &rank_info,
+    std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr<Domain> Dm, int Np,
+    bool Regular, IntArray Map )
+    : d_Np( Np ),
+      d_regular( Regular ),
+      d_rank_info( rank_info ),
+      d_Map( Map ),
+      d_comm( Dm->Comm.dup() ),
+      d_ScaLBL_Comm( ScaLBL_Comm )*/
+{
+	
+	d_comm = ColorModel.Dm->Comm.dup();
+	d_Np = ColorModel.Np;
+	bool Regular = false;
+	
+	auto input_db = ColorModel.db;
+    auto db     = input_db->getDatabase( "Analysis" );
+    auto vis_db = input_db->getDatabase( "Visualization" );
+
+    // Ids of work items to use for dependencies
+    ThreadPool::thread_id_t d_wait_blobID;
+    ThreadPool::thread_id_t d_wait_analysis;
+    ThreadPool::thread_id_t d_wait_vis;
+    ThreadPool::thread_id_t d_wait_restart;
+    ThreadPool::thread_id_t d_wait_subphase;
+
+    char rankString[20];
+    sprintf( rankString, "%05d", ColorModel.Dm->rank() );
+    d_n[0] = ColorModel.Dm->Nx - 2;
+    d_n[1] = ColorModel.Dm->Ny - 2;
+    d_n[2] = ColorModel.Dm->Nz - 2;
+    d_N[0] = ColorModel.Dm->Nx;
+    d_N[1] = ColorModel.Dm->Ny;
+    d_N[2] = ColorModel.Dm->Nz;
+
+    d_restart_interval           = db->getScalar<int>( "restart_interval" );
+    d_analysis_interval          = db->getScalar<int>( "analysis_interval" );
+    d_subphase_analysis_interval = INT_MAX;
+    d_visualization_interval     = INT_MAX;
+    d_blobid_interval            = INT_MAX;
+    if ( db->keyExists( "blobid_interval" ) ) {
+        d_blobid_interval = db->getScalar<int>( "blobid_interval" );
+    }
+    if ( db->keyExists( "visualization_interval" ) ) {
+        d_visualization_interval = db->getScalar<int>( "visualization_interval" );
+    }
+    if ( db->keyExists( "subphase_analysis_interval" ) ) {
+        d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
+    }
+
+    auto restart_file = db->getScalar<std::string>( "restart_file" );
+    d_restartFile     = restart_file + "." + rankString;
+
+
+    d_rank = d_comm.getRank();
+    writeIDMap( ID_map_struct(), 0, id_map_filename );
+    // Initialize IO for silo
+    IO::initialize( "", "silo", "false" );
+    // Create the MeshDataStruct
+    d_meshData.resize( 1 );
+
+    d_meshData[0].meshName = "domain";
+    d_meshData[0].mesh     = std::make_shared<IO::DomainMesh>(
+        d_rank_info, d_n[0], d_n[1], d_n[2], ColorModel.Dm->Lx, ColorModel.Dm->Ly, ColorModel.Dm->Lz );
+    auto PhaseVar    = std::make_shared<IO::Variable>();
+    auto PressVar    = std::make_shared<IO::Variable>();
+    auto VxVar       = std::make_shared<IO::Variable>();
+    auto VyVar       = std::make_shared<IO::Variable>();
+    auto VzVar       = std::make_shared<IO::Variable>();
+    auto SignDistVar = std::make_shared<IO::Variable>();
+    auto BlobIDVar   = std::make_shared<IO::Variable>();
+
+    if ( vis_db->getWithDefault<bool>( "save_phase_field", true ) ) {
+        PhaseVar->name = "phase";
+        PhaseVar->type = IO::VariableType::VolumeVariable;
+        PhaseVar->dim  = 1;
+        PhaseVar->data.resize( d_n[0], d_n[1], d_n[2] );
+        d_meshData[0].vars.push_back( PhaseVar );
+    }
+
+    if ( vis_db->getWithDefault<bool>( "save_pressure", false ) ) {
+        PressVar->name = "Pressure";
+        PressVar->type = IO::VariableType::VolumeVariable;
+        PressVar->dim  = 1;
+        PressVar->data.resize( d_n[0], d_n[1], d_n[2] );
+        d_meshData[0].vars.push_back( PressVar );
+    }
+
+    if ( vis_db->getWithDefault<bool>( "save_velocity", false ) ) {
+        VxVar->name = "Velocity_x";
+        VxVar->type = IO::VariableType::VolumeVariable;
+        VxVar->dim  = 1;
+        VxVar->data.resize( d_n[0], d_n[1], d_n[2] );
+        d_meshData[0].vars.push_back( VxVar );
+        VyVar->name = "Velocity_y";
+        VyVar->type = IO::VariableType::VolumeVariable;
+        VyVar->dim  = 1;
+        VyVar->data.resize( d_n[0], d_n[1], d_n[2] );
+        d_meshData[0].vars.push_back( VyVar );
+        VzVar->name = "Velocity_z";
+        VzVar->type = IO::VariableType::VolumeVariable;
+        VzVar->dim  = 1;
+        VzVar->data.resize( d_n[0], d_n[1], d_n[2] );
+        d_meshData[0].vars.push_back( VzVar );
+    }
+
+    if ( vis_db->getWithDefault<bool>( "save_distance", false ) ) {
+        SignDistVar->name = "SignDist";
+        SignDistVar->type = IO::VariableType::VolumeVariable;
+        SignDistVar->dim  = 1;
+        SignDistVar->data.resize( d_n[0], d_n[1], d_n[2] );
+        d_meshData[0].vars.push_back( SignDistVar );
+    }
+
+    if ( vis_db->getWithDefault<bool>( "save_connected_components", false ) ) {
+        BlobIDVar->name = "BlobID";
+        BlobIDVar->type = IO::VariableType::VolumeVariable;
+        BlobIDVar->dim  = 1;
+        BlobIDVar->data.resize( d_n[0], d_n[1], d_n[2] );
+        d_meshData[0].vars.push_back( BlobIDVar );
+    }
+
+
     // Initialize the comms
     for ( int i = 0; i < 1024; i++ )
         d_comm_used[i] = false;

analysis/runAnalysis.h

@@ -7,6 +7,7 @@
 #include "common/Communication.h"
 #include "common/ScaLBL.h"
 #include "threadpool/thread_pool.h"
+#include "models/ColorModel.h"
 #include <limits.h>
 
 
@@ -32,6 +33,8 @@ public:
         std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr<Domain> dm, int Np,
         bool Regular, IntArray Map );
     
+    runAnalysis( ScaLBL_ColorModel &ColorModel);
+
     //! Destructor
     ~runAnalysis();
 

cpu/MixedGradient.cpp

@@ -7,7 +7,7 @@ extern "C" void ScaLBL_D3Q19_MixedGradient(int *Map, double *Phi, double *Gradie
 			{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1},
 			{0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}};
 	
-	int i,j,k,n,N;
+	int i,j,k,n;
 	int np,np2,nm; // neighbors
 	double v,vp,vp2,vm; // values at neighbors
 	double grad;

example/Bubble/input.db

@@ -12,6 +12,22 @@ Color {
     ComponentAffinity = -1.0, 1.0, -1.0
 }
 
+FreeLee {
+    tauA = 1.0;  
+    tauB = 1.0;  
+    tauM = 1.0;//relaxation parameter for the phase field
+    rhoA   = 1.0;  
+    rhoB   = 1.0;  
+    gamma = 1.0e-4;//surface tension parameter in Lee model
+    W = 3.0; //theoretical interfacial thickness in Lee model; unit:[voxel]
+    F = 0, 0, 0
+    Restart = false
+    timestepMax = 1000
+    flux = 0.0
+    ComponentLabels = 0
+    ComponentAffinity = -1.0
+}
+
 Domain {
     nproc = 1, 1, 1     // Number of processors (Npx,Npy,Npz)
     n = 80, 80, 80      // Size of local domain (Nx,Ny,Nz)

models/ColorModel.cpp

@@ -531,6 +531,121 @@ void ScaLBL_ColorModel::Initialize(){
 	ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double));
 }
 
+double ScaLBL_ColorModel::Run(int returntime){
+	int nprocs=nprocx*nprocy*nprocz;
+
+	//************ MAIN ITERATION LOOP ***************************************/
+	comm.barrier();
+	PROFILE_START("Loop");
+	//std::shared_ptr<Database> analysis_db;
+	bool Regular = false;
+	auto current_db = db->cloneDatabase();
+	auto t1 = std::chrono::system_clock::now();
+	int START_TIMESTEP = timestep;
+	int EXIT_TIMESTEP = min(timestepMax,returntime);
+	while (timestep < EXIT_TIMESTEP ) {
+		//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
+		PROFILE_START("Update");
+		// *************ODD TIMESTEP*************
+		timestep++;
+		// Compute the Phase indicator field
+		// Read for Aq, Bq happens in this routine (requires communication)
+		ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
+		ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
+		ScaLBL_Comm->Barrier();
+		ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
+
+		// Perform the collision operation
+		ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
+		if (BoundaryCondition > 0 && BoundaryCondition < 5){
+			ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
+			ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
+		}
+		// Halo exchange for phase field
+		ScaLBL_Comm_Regular->SendHalo(Phi);
+
+		ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
+				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm_Regular->RecvHalo(Phi);
+		ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+		ScaLBL_Comm->Barrier();
+		// Set BCs
+		if (BoundaryCondition == 3){
+			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+		if (BoundaryCondition == 4){
+			din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+		else if (BoundaryCondition == 5){
+			ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
+			ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
+		}
+		ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
+				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
+		ScaLBL_Comm->Barrier(); 
+
+		// *************EVEN TIMESTEP*************
+		timestep++;
+		// Compute the Phase indicator field
+		ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
+		ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
+		ScaLBL_Comm->Barrier();
+		ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
+
+		// Perform the collision operation
+		ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
+		// Halo exchange for phase field
+		if (BoundaryCondition > 0 && BoundaryCondition < 5){
+			ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
+			ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
+		}
+		ScaLBL_Comm_Regular->SendHalo(Phi);
+		ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
+				alpha, beta, Fx, Fy, Fz,  Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
+		ScaLBL_Comm_Regular->RecvHalo(Phi);
+		ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+		ScaLBL_Comm->Barrier();
+		// Set boundary conditions
+		if (BoundaryCondition == 3){
+			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+		else if (BoundaryCondition == 4){
+			din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+		else if (BoundaryCondition == 5){
+			ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
+			ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
+		}
+		ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
+				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
+		ScaLBL_Comm->Barrier(); 
+		//************************************************************************
+	}
+	PROFILE_STOP("Update");
+
+	PROFILE_STOP("Loop");
+	PROFILE_SAVE("lbpm_color_simulator",1);
+	//************************************************************************
+	// Compute the walltime per timestep
+	auto t2 = std::chrono::system_clock::now();
+	double cputime = std::chrono::duration<double>( t2 - t1 ).count() / (timestep - START_TIMESTEP);
+	// Performance obtained from each node
+	double MLUPS = double(Np)/cputime/1000000;
+
+	if (rank==0) printf("********************************************************\n");
+	if (rank==0) printf("CPU time = %f \n", cputime);
+	if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
+	return(MLUPS);
+	MLUPS *= nprocs;
+
+}
+
 void ScaLBL_ColorModel::Run(){
 	int nprocs=nprocx*nprocy*nprocz;
 	const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
@@ -580,7 +695,6 @@ void ScaLBL_ColorModel::Run(){
 	if (color_db->keyExists( "krA_morph_factor" )){
 		KRA_MORPH_FACTOR  = color_db->getScalar<double>( "krA_morph_factor" );
 	}
-	
 	/* defaults for simulation protocols */
 	auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
 	if (protocol == "image sequence"){
@@ -625,7 +739,7 @@ void ScaLBL_ColorModel::Run(){
 	if (analysis_db->keyExists( "seed_water" )){
 		seed_water = analysis_db->getScalar<double>( "seed_water" );
 		if (rank == 0) printf("Seed water in oil %f (seed_water) \n",seed_water);
-		USE_SEED = true;
+		ASSERT(protocol == "seed water");
 	}
 	if (analysis_db->keyExists( "morph_delta" )){
 		morph_delta = analysis_db->getScalar<double>( "morph_delta" );
@@ -656,7 +770,6 @@ void ScaLBL_ColorModel::Run(){
 		MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
 	}
 
-
 	if (rank==0){
 		printf("********************************************************\n");
 		if (protocol == "image sequence"){
@@ -1320,7 +1433,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 	double vF = 0.f;
 	double vS = 0.f;
 	double delta_volume;
-	double WallFactor = 0.0;
+	double WallFactor = 1.0;
 	bool USE_CONNECTED_NWP = false;
 
 	DoubleArray phase(Nx,Ny,Nz);
@@ -1343,6 +1456,11 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 		}
 	}
 	double volume_initial = Dm->Comm.sumReduce(  count);
+	double PoreVolume = Dm->Volume*Dm->Porosity();
+	/*ensure target isn't an absurdly small fraction of pore volume */
+	if (volume_initial < target_delta_volume*PoreVolume){
+		volume_initial = target_delta_volume*PoreVolume;
+	}
 	/*
 	sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank);
 	FILE *INPUT = fopen(LocalRankFilename,"wb");

models/ColorModel.h

@@ -16,6 +16,10 @@ Implementation of color lattice boltzmann model
 #include "ProfilerApp.h"
 #include "threadpool/thread_pool.h"
 
+
+#ifndef ScaLBL_ColorModel_INC
+#define ScaLBL_ColorModel_INC
+
 class ScaLBL_ColorModel{
 public:
 	ScaLBL_ColorModel(int RANK, int NP, const Utilities::MPI& COMM);
@@ -29,6 +33,7 @@ public:
 	void Create();
 	void Initialize();
 	void Run();
+	double Run(int returntime);
 	void WriteDebug();
 	void getPhaseField(DoubleArray &f);
 	
@@ -99,4 +104,5 @@ private:
 	int timestep;
 	int timestep_previous;
 };
+#endif
 

tests/CMakeLists.txt

@@ -62,7 +62,7 @@ ADD_LBPM_TEST( TestMap )
 ADD_LBPM_TEST( TestWideHalo )
 ADD_LBPM_TEST( TestColorGradDFH )
 ADD_LBPM_TEST( TestBubbleDFH ../example/Bubble/input.db)
-ADD_LBPM_TEST( testGlobalMassFreeLee ../example/Bubble/input.db)
+#ADD_LBPM_TEST( testGlobalMassFreeLee ../example/Bubble/input.db)
 #ADD_LBPM_TEST( TestColorMassBounceback ../example/Bubble/input.db)
 ADD_LBPM_TEST( TestPressVel  ../example/Bubble/input.db)
 ADD_LBPM_TEST( TestPoiseuille ../example/Piston/poiseuille.db)

tests/lbpm_color_simulator.cpp

@@ -27,19 +27,24 @@ int main( int argc, char **argv )
     // Initialize
     Utilities::startup( argc, argv );
 
-    // Load the input database
-    auto db = std::make_shared<Database>( argv[1] );
-
     { // Limit scope so variables that contain communicators will free before MPI_Finialize
 
         Utilities::MPI comm( MPI_COMM_WORLD );
         int rank   = comm.getRank();
         int nprocs = comm.getSize();
+        std::string SimulationMode = "production";
+        // Load the input database
+        auto db = std::make_shared<Database>( argv[1] );
+        if (argc > 2) {
+        	SimulationMode = "development";
+        }
 
         if ( rank == 0 ) {
             printf( "********************************************************\n" );
             printf( "Running Color LBM	\n" );
             printf( "********************************************************\n" );
+            if (SimulationMode == "development")
+            	printf("**** DEVELOPMENT MODE ENABLED *************\n");
         }
         // Initialize compute device
         int device = ScaLBL_SetDevice( rank );
@@ -62,8 +67,29 @@ int main( int argc, char **argv )
         ColorModel.Create();     // creating the model will create data structure to match the pore
                                  // structure and allocate variables
         ColorModel.Initialize(); // initializing the model will set initial conditions for variables
+       
+        if (SimulationMode == "development"){
+            double MLUPS=0.0;
+            int timestep = 0;
+            int analysis_interval = ColorModel.timestepMax;
+        	if (ColorModel.analysis_db->keyExists( "" )){
+        		analysis_interval = ColorModel.analysis_db->getScalar<int>( "analysis_interval" );
+        	}
+        	FlowAdaptor Adapt(ColorModel);
+        	runAnalysis analysis(ColorModel);
+            while (ColorModel.timestep < ColorModel.timestepMax){
+            	timestep += analysis_interval;
+            	MLUPS = ColorModel.Run(timestep);
+            	if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS);
+            	
+            	Adapt.MoveInterface(ColorModel);
+            }
+        }            	//Analysis.WriteVis(LeeModel,LeeModel.db, timestep);
+
+        else
         	ColorModel.Run();
-        // ColorModel.WriteDebug();
+        
+        ColorModel.WriteDebug();
 
         PROFILE_STOP( "Main" );
         auto file  = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );

tests/lbpm_freelee_simulator.cpp

@@ -62,8 +62,8 @@ int main( int argc, char **argv )
         double MLUPS=0.0;
         int timestep = 0;
         int visualization_time = LeeModel.timestepMax;
-    	if (LeeModel.vis_db->keyExists( "visualizataion_interval" )){
+    	if (LeeModel.vis_db->keyExists( "visualization_interval" )){
-    		visualization_time = LeeModel.vis_db->getScalar<int>( "visualizataion_interval" );
+    		visualization_time = LeeModel.vis_db->getScalar<int>( "visualization_interval" );
     		timestep += visualization_time;
     	}
         while (LeeModel.timestep < LeeModel.timestepMax){