started testing, two still failing

docs/src/tutorials.md

@@ -271,6 +271,340 @@ end
 
 Given the finite slippage we do not observe *large* deviations from the $\alpha=1/10$ powerlaw in the long time limit.
 
+
+## Active film
+
+In [Richter et all](https://arxiv.org/abs/2402.14635) we introduced a thin film model for chemical reactions inside a liquid film aimed at Supported liquid phase catalysts. For the details of the model look at the article. Here we give a minimal example of how to run a simulation, starting from a homogenous film, leading to film rupture. 
+
+Lets get straight in it 
+```julia 
+using Swalbe, DelimitedFiles
+function simulation(;
+    D=5e-6,			# Diffusion coefficient
+    Tmax=1e3, 			# Simulation time
+    tdump=Int(floor(Tmax/100)), # dump time
+    h= 1,  			# initial film height
+    rho=1,			# initial catalyst concentration
+    ϵ=0.001,			# initial noise
+    γ_0=0.001,			# reference surface tension in the absence of chemical compounds
+    Γ = 0.1, 			# Surface tension effect product rho_A
+    GammaB=0.1,			# Surface tension effect reactant rhoB
+    data = "data",		# location to save data
+    L=2^12, 			# system size
+    alpha=0.0, 			# catalyst vertical distribution parameter 
+    theta_0=1/9,		# reference contact angle in fractions of pi
+    n=9, m=3, hmin=0.1, hcrit=0.05, hcrit2=0.05,	# disjoining pressure parameters
+    delta=1.0,			# slip length
+    production_rate=0.1,	# production rate
+    sigma_A_up=0.1, 		# evaporation rate product rho_A
+    sigma_B_up=0.1, 		# evaporation rate reactant rho_B
+    rho_BRes_sigma_B_down=0.1,	# sorption rate reactant rho_B
+    rho_0B=rho_BRes_sigma_B_down*h/(production_rate*rho + sigma_B_up),	# initial reactant concentration 
+    rho_0A=production_rate*rho*rho_0B/sigma_A_up,			# initial product cocentration 
+    D_Product=D*20, D_B=D*20,	# Diffusion coefficients Product rho_A and reactant rho_B
+    rho_crit=0.05,		# precursor length catalyst
+    mu=1/6,			# viscosity
+)
+	# Set up sys consts
+    	    GammafacA=rho_0A/(h*γ_0)
+    	    GammafacB=rho_0B/(h*γ_0)
+	    sys = Swalbe.SysConstActive_1D{Float64}(
+        	L =  L,
+        	D=D,
+        	Γ = Γ/GammafacA,
+        	GammaB=GammaB/GammafacB,
+        	Tmax=Tmax,
+        	tdump=tdump,
+        	γ_0 = γ_0,
+        	δ=delta,
+        	alpha=alpha,
+        	θ_0 =theta_0,
+        	n=n,m=m,hmin=hmin,hcrit=hcrit, hcrit2=hcrit2,
+        	production_rate=production_rate,
+        	sigma_A_up=sigma_A_up,
+        	sigma_B_up=sigma_B_up,
+        	rho_BRes_sigma_B_down=rho_BRes_sigma_B_down,
+        	D_Product=D_Product,
+        	D_B=D_B,
+        	rho_crit=rho_crit,
+        	b_A=rho_crit,
+        	b_B=rho_crit,
+        	μ=mu
+	)
+	# set up simulation fields 
+	state = Swalbe.Sys(sys)
+	# Set initial values
+	state.rho.=rho
+        state.rho_A .= rho_0A
+        state.rho_B .= rho_0B
+	# This set up the film height as constant plus a small noise without high wavenumber frequencies
+        Swalbe.browniannoise!(state.height, h, ϵ, L/2-L/10)
+	i = 0
+	# The LBM time loop
+	for t in 0:sys.Tmax
+		# dump data
+		if t % sys.tdump == 0 
+			writedlm("$(data)/$(Swalbe.to_4_digits(i))_height.csv", state.height)
+			writedlm("$(data)/$(Swalbe.to_4_digits(i))_rho.csv", state.rho)
+			writedlm("$(data)/$(Swalbe.to_4_digits(i))_rho_A.csv", state.rho_A)
+			writedlm("$(data)/$(Swalbe.to_4_digits(i))_rho_B.csv", state.rho_B)
+			i += 1
+			println("Time step t=$t")
+		end
+		# system update, mostly as in a standard TFE simulation, the parts that are different are commented
+		Swalbe.surface_tension_4_fields!(state, sys)
+		Swalbe.filmpressure_fast!(state, sys)
+		Swalbe.h∇p!(state)
+            	Swalbe.slippage!(state, sys)
+		# claculate the Marangoni stress
+		Swalbe.∇γ!(state, sys)
+		# Include the Marangoni stress into the forcing
+            	state.F .= -state.h∇p .- state.slip .+ state.stress
+		# Do the LBM loops for catalys state.rho, product state.rho_A, and reactant state.rho_B
+		Swalbe.update_rho_LBM!(state,sys)
+		Swalbe.update_rho_A_LBM_precursor!(state,sys)
+		Swalbe.update_rho_B_LBM_precursor!(state,sys)
+		Swalbe.equilibrium!(state)
+		Swalbe.BGKandStream!(state, sys)
+		Swalbe.moments!(state)	
+	end
+end
+
+# run the simulation 
+simulation()
+```
+
+
+That's it
+
+## Thin films on curved substrate
+
+It is possible to include a curved substrate into the model using the substrate profile `b` into the pressure
+```math
+p_{\text{cap}} = -(1 + 0.5 \nabla s^2) \gamma (1 - \cos(\theta)) \phi'(h) - \gamma \nabla^2 (h+b)
+```
+
+Here is an example on how to do this 
+
+```julia 
+using Swalbe, DelimitedFiles, Dates
+function simulation(;
+	Tmax=1e8, 			# Simulation time
+	tdump=Int(floor(Tmax/100)),	# dump time
+	h= 1,				# initial film height  
+	eps=0.001, 			# initial noise
+	gamma=1/6^2,			# surface tension
+	data = "data",			# folder for saving data
+ 	L=2^12,				# system size
+    	theta=1/9,			# contact angle
+    	n=9, m=3, hmin=0.1, 		# parameters for disjoining pressure
+	delta=1.0,			# slip length, in combination with heigher hmin i
+	omega=3,			# number of periods of the underlying substrate
+)
+	# set up system
+    	sys = Swalbe.SysConst_1D{Float64}(L =  L ,Tmax=Tmax, tdump=tdump, γ = gamma, δ=delta,θ =theta, n=n,m=m,hmin=hmin)
+     	try
+        	mkdir(data)
+            	mkdir("$(data)/figs")
+        catch
+        end
+    	state = Swalbe.Sys(sys, kind="curved")
+        # Initial data
+    	state.height.= h .+ eps .* rand(sys.L)
+	# set up system profile
+    	state.substrate .=[1/4*sin(i*omega*2*pi/sys.L) for i in 1:sys.L]
+	# We need the gradient of the substrate profile
+    	Swalbe.∇f!(state.grad_substrate, state.substrate, state.dgrad)
+    	println("Starting Lattice Boltzmann time loop for active thin film")
+    	i=0
+    	for t in 0:sys.Tmax
+		# save data
+        	if t % sys.tdump == 0
+                	writedlm("$(data)/$(Swalbe.to_4_digits(i))_height.csv", state.height)
+            		i+=1
+                end
+		# system update, using the pressure for curved substrates
+        	Swalbe.filmpressure_curved!(state, sys)
+        	Swalbe.h∇p!(state)
+        	Swalbe.slippage!(state, sys)
+        	state.F .= -state.h∇p .- state.slip
+        	Swalbe.equilibrium!(state)
+        	Swalbe.BGKandStream!(state, sys)
+        	Swalbe.moments!(state)
+    	end
+    println("Finished timeloop")
+    return nothing
+end
+simulation()
+```
+
+
+## Miscible liquids
+
+Here I show to simulate the colaescence of two droplets of miscible liquids with different surface tensions, lets say water and ethanol. 
+
+```julia
+using Swalbe, DelimitedFiles,  Dates
+
+
+function run_Miscible(;
+    eps=0.01,                                   # initial noise
+    Tmax=Int(1e7),                              # Simulation time
+    dumps=100,                                  # number of dumps
+    tdump=Int(floor(Tmax/dumps)),               # dump time
+    L=2^9,                                      # system size
+    delta=1.0,                                  # slip length
+    hmin=0.1, hcrit=0.05, n=9, m=3,             # disjoining pressure parameters
+    gamma= 0.01 .* ones(3,2),                   # surface tensions
+    data="path/to/data",                        # save path
+    D=1e-2,                                     # diffusion coeficient
+    r=50,                                       # Initial droplet radiusj
+    center_weight=0.5                           # weight for surface tension smoothing
+    )
+    #Make folders
+    datathis= "$(data)/Tmax=$Tmax,L=$L,delta=$delta,gamma=$(replace(replace("$gamma", ";"=>","), " " => "_")),D=$D,r=$r,hmin=$hmin,hcrit=$hcrit"
+        try
+        mkdir(dataall)
+        catch
+        end
+        try
+        mkdir(datathis)
+        catch
+        end
+        try
+            mkdir("$(datathis)/figs")
+        catch
+        end
+    #Setup system
+    sys = Swalbe.SysConstMiscible_1D(L=L, Tmax=Tmax, tdump=tdump, gamma=gamma, delta=delta, n=n, m=m, D=D, hmin=hmin, hcrit=hcrit)
+    state = Swalbe.Sys(sys)
+    # Initial data
+    theta_1=acos((sys.gamma[3,1]-sys.gamma[2,1])/sys.gamma[1,1])
+    theta_2=acos((sys.gamma[3,2]-sys.gamma[2,2])/sys.gamma[1,2])
+    state.height[:,1] .= Swalbe.droplet_base(state.height[:,1], r, theta_1/pi, Int(sys.L/2 - r/2), precursor=(sys.hmin-2*sys.hcrit)/2)
+    state.height[:,2] .= Swalbe.droplet_base(state.height[:,2], r, theta_2/pi, Int(sys.L/2 + r/2), precursor=(sys.hmin-2*sys.hcrit)/2)
+    bridge_height=[]
+    i=0
+    println("Starting LBM time loop")
+    for t in 0:sys.Tmax
+        if t % sys.tdump == 0
+            #Store data
+            writedlm("$(datathis)/$(Swalbe.to_4_digits(i))_h_t=$t",state.height)
+            println("time step t=$t")
+            i+=1
+        end
+        #actual simulations
+            # surface tension calculation from the height fields
+            # Swalbe.surface_tension!(state,sys)
+            # surface tension calculation from the height fields with moving averadge smoothing applied
+            Swalbe.surface_tension_smooth!(state,sys, center_weight=center_weight)
+            Swalbe.filmpressure!(state, sys)
+            Swalbe.h∇p!(state,sys)
+            Swalbe.slippage!(state,sys)
+            # calculate Marangoni shear stress
+            Swalbe.∇gamma!(state, sys)
+            state.F .= -state.h∇p  .-  state.slip .+ state.stress
+            Swalbe.equilibrium!(state)
+            Swalbe.BGKandStream!(state, sys)
+            Swalbe.moments!(state, sys)
+    end
+end
+
+
+dataall = "data/$(Dates.today())Miscible_coal_test_00"
+gamma=ones(3,2)
+gamma_0=0.002
+# fixed contact angles
+Delta=0.2
+gamma[1,1]=gamma_0*(1+Delta/2)
+gamma[1,2]=gamma_0*(1-Delta/2)
+# solid vapour
+gamma[3,1]=2*gamma_0
+gamma[3,2]=2*gamma_0
+theta=1/9
+# fixed contact angle
+gamma[2,1]=gamma[3,1]-gamma[1,1]*cospi(theta)
+gamma[2,2]=gamma[3,2]-gamma[1,2]*cospi(theta)
+run_Miscible(Tmax=Int(1e7),gamma=gamma, D=1.5e-4, r=200, delta=1.0, hmin=0.2, L=2^10, data=dataall)
+```
+
+Have fun with it
+
+## Multilayer thin films
+
+This shows how to simulate a layer film resting on a layer of water. This we call a two layer multilayer system. To read more about the theory see [Richter et al.](https://arxiv.org/abs/2409.16659). 
+
+This is the code on how to do this: 
+```julia 
+using Pkg
+
+Pkg.activate(".")
+
+
+using Swalbe, DelimitedFiles, CUDA
+
+
+function run_multilayer(;
+    h_1=2,		# initial height lower layer
+    h_2=1,		# initial height upper layer
+    eps=[0.001,0.001],	# initial noise
+    Tmax=Int(1e7),	# simulation time
+    dumps=100,		# number of dumps
+    tdump=Int(floor(Tmax/dumps)), 	# dump time
+    Lx=2^8, Ly=2^8,	# system size
+    delta=1.0,		# slip lower layer
+    delta_2=1.0,	# slip upper layer
+    n=9, m=3,		# disjoining pressure exponents
+    gamma= [0 0.001 0.001 0.001 ; 0 0 0.001 0.001 ; 0 0 0 0.001; 0 0 0 0],	# surface tensions gamma[i,j] is the interaction between layer i-1 and j-1 having 0 for the solid substrat and 3 the atmosphere
+    data="data", 	# path to save simulation dumps
+    hmin=0.1,		# precursor layer
+    mu=[1/6 1/6],	# viscosities 
+    )
+    #Setup system
+    sys = Swalbe.SysConstMultiLayer(Lx=Lx, Ly=Ly, Tmax=Tmax, tdump=tdump, gamma=gamma, delta=[delta,delta_2], n=n, m=m, hmin=hmin, mu=mu)
+    state = Swalbe.Sys(sys, device=device)
+    # this is needed to send back and force in between GPU and CPU
+    host_height=zeros(Lx,Ly,2)
+    #initial conditions
+    host_height[:,:,1] .= Swalbe.browniannoise2D!(host_height[:,:,1], h_1, eps[1], Lx/2-Lx/10)
+    host_height[:,:,2] .= Swalbe.browniannoise2D!(host_height[:,:,2], h_2, eps[2], Lx/2-Lx/10)
+    CUDA.copyto!(state.height, host_height)
+    i=0
+    println("Starting LBM time loop")
+    for t in 0:sys.Tmax
+        if t % sys.tdump == 0 
+            #Store data
+                CUDA.copyto!(host_height,state.height)
+                writedlm("$(datathis)/$(Swalbe.to_4_digits(i))_h1_t=$t",host_height[:,:,1])
+                writedlm("$(datathis)/$(Swalbe.to_4_digits(i))_h2_t=$t",host_height[:,:,2])
+            #status report
+            println("Time step t=$t)
+            i+=1
+        end
+        #actual simulations
+        if run 
+            Swalbe.filmpressure!(state, sys)
+	    Swalbe.h∇p!(state,sys)
+            Swalbe.slippage!(state, sys)
+	    state.Fx .= -state.h∇px  .+  state.slipx
+	    state.Fy .= -state.h∇py  .+  state.slipy
+            Swalbe.equilibrium!(state)
+            Swalbe.BGKandStream!(state)
+            Swalbe.moments!(state, sys)
+        end
+    end 
+end
+
+
+dataall = "path\to\save\dumps"
+gamma_0=0.01
+gamma=gamma_0*ones(4,4)
+gamma[2,4]=gamma_0*(1+cospi(2/9))
+gamma[1,3]=gamma_0*(2.1)
+gamma[1,4]=gamma[2,4]+gamma[1,3]-gamma[2,3]
+run_multilayer(Tmax=Int(1e6),gamma=gamma, h_2=1.0, h_1=6.0, delta=1.0, delta_2=1.0, Lx=2^8,Ly=2^8, mu=[0.1 0.1], hmin=0.2, devicenumber=2)
+```
+
 ## Further tutorials
 
 More tutorials will follow in the future.
@@ -279,4 +613,4 @@ So be sure to check out the docs every now and then.
 
 The next tutorial will be about switchable substrates. 
 In this case the wettability can not only addressed locally but also with a time dependency.
-Here is what happens if the time frequency is [high](https://gist.github.com/Zitzeronion/116d87978ece82c8ae64a3f7edb9dbb3#gistcomment-3811414) and this happens if we update with a [lower](https://gist.github.com/Zitzeronion/116d87978ece82c8ae64a3f7edb9dbb3#gistcomment-3811414) frequency.
+Here is what happens if the time frequency is [high](https://gist.github.com/Zitzeronion/116d87978ece82c8ae64a3f7edb9dbb3#gistcomment-3811414) and this happens if we update with a [lower](https://gist.github.com/Zitzeronion/116d87978ece82c8ae64a3f7edb9dbb3#gistcomment-3811414) frequency.

scripts/init.jl

@@ -0,0 +1,5 @@
+using Pkg
+
+Pkg.activate(".")
+
+using Swalbe

src/differences.jl

@@ -229,6 +229,20 @@ function ∇f!(output, f::Vector, dgrad)
     return nothing
 end
 
+
+function ∇f_Miscible!(output, f::Array, dgrad)
+    fip, fim = viewneighborsMiscible_1D(dgrad)
+    # One dim case, central differences
+    circshift!(fip, f, 1)
+    circshift!(fim, f, -1)
+    # In the end it is just a weighted sum...
+    output .= -0.5 .* (fip .- fim)
+    return nothing
+end
+
+
+
+
 """
     viewneighbors(f)