Minor changes, maybe look at it during the weekend

scripts/sliding_analysis.jl

@@ -16,16 +16,13 @@ using DataFrames, CSV, CodecZlib, Mmap, Images
 # Define some constants
 lx = 512
 ly = 256
-maxtime = 50000
+maxtime = 200000
 θₑ = 1/9    # Contact angle on the substrate, expect on the patch
-δ₀ = 4.0    # Slip length, larger better for geometry 
+δ₀ = 2.0    # Slip length, larger better for geometry 
 T = Float64
 # The center of the lattice
 center = (lx÷2, ly÷2)
-# Things to be measured
-measureframe = DataFrame()
-dropletframe = DataFrame()
-contactlineframe = DataFrame()
+
 # Geometric parameters for the patterns
 boxside = 62
 triangleside = 94.3
@@ -41,6 +38,9 @@ lat = JuThinFilm.lattice{T}()
 for R₀ in [55.8 69.75 83.7]         
     # Loop over the four different patch shapes
     for shape in ["box" "triangle" "ellipse" "circle"]
+        measureframe = DataFrame()
+        contactlineframe = DataFrame()
+        dropletframe = DataFrame()
         # Loop over five different contact angle contrasts, last one is less wettable than the surrounding substrate
         for δᵪ in [1/36 1/18 1/12 1/9 -1/18]
             vol = π/3*R₀^3*(2+cospi(θₑ)*(1-cospi(θₑ))^2)
@@ -50,11 +50,10 @@ for R₀ in [55.8 69.75 83.7]
             diffsurf = []
             energy = []
             contactline = []
-            positions = [0, 0, 0]
+            positions = []
             MaxVel = []
             θₙ = fill(θₑ,lx,ly)
             angle_switch = true
-            # And the defect with the correct shape and contrast
             
             # Forces needed for the simulation
             Fpressure = default_force(T, lx, ly)
@@ -68,6 +67,7 @@ for R₀ in [55.8 69.75 83.7]
             copy!(dists.ftemp,dists.feq)
             # Start of the main time loop
             for t in 1:parameters.mt
+                # isnan.(mom.height) && break
                 # Only change the substrate after the droplet has started sliding!
                 if findmax(mom.height, dims=1)[2][:1][1][1] > lx÷2+150 && t > 1000 && angle_switch   
                     if shape == "box" 
@@ -83,13 +83,13 @@ for R₀ in [55.8 69.75 83.7]
                     angle_switch = false
                 end
                 # Performe the measurements
-                if t%10 == 0
+                if t%500 == 0
                     # Measurements of surface area and three phase contact line lenght
                     JuThinFilm.measuresurfaceandenergy!(diffsurf, energy, mom, θₙ)
                     JuThinFilm.measurecontactlinelengthandposition!(contactline, positions, t, mom)
-                    JuThinFilm.measuremaxvelocity(MaxVel, mom)
+                    JuThinFilm.measuremaxvelocity(MaxVel, mom) 
                 end
-                if t%5000  == 0
+                if t%10000  == 0
                     # IO for progression
                     mass = round(sum(mom.height),digits=3)
                     println("time step $t and mass $mass")
@@ -121,20 +121,22 @@ for R₀ in [55.8 69.75 83.7]
             dropletframe["vely_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = vec(mom.velocity.y)
             dropletframe["pres_shape_$(shape)_contrast_$(angle_diff)_vol_$(R₀)"] = vec(mom.pressure)
         end
+        # Save the collected data for any given shape to the disk
+        savepath = "data/sims/dynamic_different_shape_angle/"
+        results_measures = "measurements_shape_$(shape)_vol_$(R₀)_set1.csv" 
+        results_drop = "final_hdvars_shape_$(shape)_vol_$(R₀)_set1.csv"
+        results_TPC = "contactline_shape_$(shape)_vol_$(R₀))set1.csv"
+        CSV.write(savepath * results_measures, measureframe)
+        CSV.write(savepath * results_drop, dropletframe)
+        CSV.write(savepath * results_TPC, contactlineframe)
+        zip1 = `gzip data/sims/dynamic_different_shape_angle/measurements_shape_$(shape)_vol_$(R₀)_set1.csv`
+        zip2 = `gzip data/sims/dynamic_different_shape_angle/final_hdvars_shape_$(shape)_vol_$(R₀)_set1.csv`
+        zip3 = `gzip data/sims/dynamic_different_shape_angle/contactline_shape_$(shape)_vol_$(R₀)_set1.csv`
+        run(zip1)
+        run(zip2)
+        run(zip3)
+        println("Simulations done results are saved in $savepath")
     end
 end
-savepath = "data/sims/dynamic_different_shape_angle/"
-results_measures = "measurements_0.csv" 
-results_drop = "final_hdvars_0.csv"
-results_TPC = "contactline_0.csv"
-CSV.write(savepath * results_measures, measureframe)
-CSV.write(savepath * results_drop, dropletframe)
-CSV.write(savepath * results_TPC, contactlineframe)
-zip1 = `gzip data/sims/dynamic_different_shape_angle/measurements_0.csv`
-zip2 = `gzip data/sims/dynamic_different_shape_angle/final_hdvars_0.csv`
-zip3 = `gzip data/sims/dynamic_different_shape_angle/contactline_0.csv`
-run(zip1)
-run(zip2)
-run(zip3)
-println("Simulations done results are saved in $savepath")
+
 # TODO: Add a send email flag, to be notified if the simulation has finished
\ No newline at end of file

src/initialstates.jl

@@ -250,7 +250,7 @@ function boxpattern(lx, ly, stats::JuThinFilm.pressurestats; center=(lx÷2, ly÷
     vertex₃ = [xcoords[3], ycoords[3]]
     vertex₄ = [xcoords[4], ycoords[4]]
     P = VPolygon([vertex₁, vertex₂, vertex₃, vertex₄])
-    for i in 1:lx, j in 1:ly
+    for i in 1:ly, j in 1:lx
         if [Float64(i), Float64(j)] ∈ P
             θₙ[j,i] = stats.θ + δₐ
         else
@@ -307,7 +307,7 @@ function ellipsepattern(lx, ly, stats::JuThinFilm.pressurestats; center=(lx÷2,
     shape₁ = a^(2)
     shape₂ = b^(2)
     P = Ellipsoid(mid, [shape₁ 0.0; 0.0 shape₂])
-    for i=1:lx, j=1:ly
+    for i=1:ly, j=1:lx
         if [Float64(i), Float64(j)] ∈ P
             θₙ[j,i] = stats.θ + δₐ
         else
@@ -368,7 +368,7 @@ function trianglepattern(lx, ly, stats::JuThinFilm.pressurestats; center=(lx÷2,
     vertex₂ = [xcoords[2], ycoords[2]]
     vertex₃ = [xcoords[3], ycoords[3]]
     P = VPolygon([vertex₁, vertex₂, vertex₃])
-    for i in 1:lx, j in 1:ly
+    for i in 1:ly, j in 1:lx
         if [Float64(i), Float64(j)] ∈ P
             θₙ[j,i] = stats.θ + δₐ
         else

src/measurements.jl

@@ -192,6 +192,7 @@ end
 
 function measurecontactlinelengthandposition!(perimeter, position, t, mom)
     len, wid = size(mom.height)
+    ll = []
     height = copy(mom.height)
     height[height .<= 0.06] .= 0
     sol = 0
@@ -199,11 +200,14 @@ function measurecontactlinelengthandposition!(perimeter, position, t, mom)
         for j in 2:wid-1
             # Basically if one or more of my neighbors is dry I have to be contactline!
             if (height[i+1, j] == 0 && height[i, j] > 0) || (height[i-1, j] == 0 && height[i, j] > 0) ||(height[i, j+1] == 0 && height[i, j] > 0) ||(height[i, j-1] == 0 && height[i, j] > 0)
-                append!(position, (i, j, t))
+                push!(ll, (i, j, t))
                 sol += 1
+            else
+                push!(ll, missing)
             end
         end
     end 
+    push!(position, ll)
     push!(perimeter, sol)
 end