Passive scalar implementation

examples/passive_scalar/benchmark/cylinder2d/args.lua

@@ -0,0 +1,5 @@
+nelem = 40
+length_usr = 128
+level_usr = math.log(length_usr)/math.log(2)
+t_total = 10
+interval = t_total / 10

examples/passive_scalar/benchmark/cylinder2d/func.lua

@@ -0,0 +1,12 @@
+require 'args'
+
+x_origin = 0.
+y_origin = 0.
+
+function gauss_pulse_real(x, y, z, t)
+  D = (tau - 0.5) / 3
+  c = c0*sigma0^2/(sigma0^2+2.*D*t)*math.exp(
+    -0.5*(( x - x_origin - u_field*t )^2+( y - y_origin )^2) / (sigma0^2 + 2*D*t)
+  ) + c_add
+  return c
+end

examples/passive_scalar/benchmark/cylinder2d/harvester.lua

@@ -0,0 +1,24 @@
+require 'args'
+require 'func'
+
+variable = {
+  {
+    name = 'spcie_real',
+    ncomponents = 1,
+    vartype = 'st_fun',
+    st_fun = gauss_pulse_real,
+  }
+}
+restart   = {read = 'restart/simulation_lastHeader.lua'}
+NOtracking  = {
+  label     = 'main',
+  folder    = 'harvest/',
+  output    = { format = 'vtk' },
+  variable  = {'spc1_density','spcie_real'},
+  -- variable  = {'velocity_phy'},
+  shape     = {kind='all'},
+}
+output = { 
+  folder = 'harvest/', 
+  {format = 'VTU'}
+} 

examples/passive_scalar/benchmark/cylinder2d/hmesh.lua

@@ -0,0 +1,8 @@
+mesh      = 'mesh/'
+
+tracking  = { 
+  label     = 'all', 
+  output    = { format = 'vtk'}, 
+  shape     = { kind='all'  },
+  folder    = 'harvest/',
+}

examples/passive_scalar/benchmark/cylinder2d/index.md

@@ -0,0 +1,37 @@
+title: Diffusion from cylinder without flow
+@warning WORK IN PROGRESS @endwarning
+
+# Diffusion from cylinder without flow # {#eg_GPP}
+
+In this example, we will investigate a system consisting of a cylindrical 
+cavity of radius a filled with a stationary liquid.
+The control equation is
+$$
+      \frac{\partial}{\partial t} C(r, t) = D \frac{1}{r} \frac{\partial}{\partial r}
+       r \frac{\partial}{\partial r} C(r, t)
+$$
+
+  with the initial and boundary condition $ C(r < a, t = 0) = C_0 $, 
+  $ C(r = a, t) = C_c $. The analytical solution is 
+$$
+      \frac{C(r, t)-C_c}{C_0-C_c} = \sum_{n=1}^{\inf} \frac{2}{\mu_n J_1(\mu_n)} 
+      exp(-\mu_n^2 \frac{Dt}{a^2}) J_0(\mu_n \frac{r}{a})
+$$
+
+  $ J_0(x) $ and $ J_1(x) $ are 0th and 1st order Bessel functions. $ \mu_n $ is the 
+  n-th root of $ J_0(x) $. The first 5 roots are 2.2048, 5.5201, 8.6537, 11.7915, 
+  14.9309 respectively. By taking initial value $ C_0 = 0 $, the analytical solution becomes
+$$
+      C(r, t) = C_c (1 - \sum_{n=1}^{\inf} \frac{2}{\mu_n J_1(\mu_n)} 
+      exp(-\mu_n^2 \frac{Dt}{a^2}) J_0(\mu_n \frac{r}{a}))
+$$
+
+In our simulation, the cylinder diameter $a=40$ is used as a default value. The D2Q9 layout is used. 
+The pressure anti-bounce back boundary condition is used to compute and compare with the 
+analytical solutions.
+
+The objectives of this example is to introduce how to:
+* Simulate time evolution of the diffution process inside a 2D cylinder
+* Testify the correctness of the anti-bounce back dirichlet boundary condition
+* Examine the solver stability with different initial conditions and collision schemes
+

examples/passive_scalar/benchmark/cylinder2d/musubi.lua

@@ -0,0 +1,149 @@
+-- Diffusion process inside a 2D cylinder. The concentration profile is compared 
+-- to the analytical solution. In this simulation, anti-bounce back boundary 
+-- condition is used.
+
+require 'args'
+-------------------------------------------------------------------------------
+mesh               = './mesh/'  
+cs2 = 1./3.
+c_init = 1.
+maxN = 150.0
+tau = 0.516
+diff = (tau - 0.5) / 3
+tIni = 0
+
+function factorial(n)
+  if n < 0.1 and n > -0.1 then
+      return 1
+  else
+      return n * factorial(n - 1)
+  end
+end
+
+-- 0th and 1st Bessel function is computed for the concentration profile
+function bessel(param, order)
+  local sum = 0.
+  for iN = 0., maxN - 1 do
+    sum = sum + (-1)^iN / (factorial(iN) * factorial(iN + order)) * (
+      param / 2) ^ (2 * iN + order)
+  end
+  return sum
+end
+
+-- Analytical solution of the time evolution of the concentration profile
+function concentration_real(x, y, z, t)
+  local mu = {2.4048, 5.5201, 8.6537, 11.7915, 14.9309}
+  local sum = 0
+  local r = math.sqrt((x-nelem)^2 + (y-nelem)^2)
+  -- print("x = " .. x .. ", y = " .. y .. ", r = " .. r .. ", sum = ".. sum)
+  if (r > nelem) then
+    return cs2
+  end
+
+  for i = 1, #mu do
+    sum = sum + 2 / (mu[i] * bessel(mu[i], 1)) * math.exp(-mu[i] ^ 2
+      * diff * t / nelem ^ 2) * bessel(mu[i] * r / nelem, 0)
+    -- if x == 13 and y == 11 then  
+    --   print("mu[i] = " .. mu[i])
+    --   print("bessel(mu[i], 1) = " .. bessel(mu[i], 1))
+    -- end
+  end
+
+  local c = c_init * (1 - sum)
+  if c < 0.001 then
+    return 0.001 * cs2
+  end 
+  return c*cs2
+end
+
+-- Initial boundary condition can be assigned with analytical solution
+-- at a certain time step.
+function initial_concentration(x, y, z)
+  return concentration_real(x, y, z, tIni)
+end
+
+function concentration_afterIni(x, y, z, t)
+  return concentration_real(x, y, z, t + tIni) / cs2
+end
+
+-------------------------------------------------------------------------------
+sim_control        = {
+  time_control     = {
+    min      = { iter = 0        },
+    max      = { iter = t_total     },
+    interval = { iter = interval }
+  }
+}
+-------------------------------------------------------------------------------
+identify  = {
+  label = 'species', 
+  kind = 'passive_scalar', 
+  relaxation='trt', 
+  layout='d2q9',
+  order = 'second'
+}
+
+transport_velocity = 'velocity_fluid'
+
+variable = {
+  {
+    name = 'velocity_fluid',
+    ncomponents = 3,
+    vartype = 'st_fun',
+    st_fun = {0., 0.0, 0.0},
+  },
+  {
+    name = 'concentration_real',
+    ncomponents = 1,
+    vartype = 'st_fun',
+    st_fun = concentration_afterIni,
+  },
+  {
+    name = 'c_diff',
+    ncomponents = 1,
+    vartype = 'operation',
+    operation = {
+      kind = 'difference',
+      input_varname = {
+        'concentration_real', 'spc1_density'
+      },
+    },
+  }
+}
+tracking  = { 
+  { label     = 'spc1',
+    variable  = {'spc1_density', 'concentration_real'},
+    shape = {
+      -- kind = 'all'
+      kind  = 'canoND',
+      object  = {
+        origin = {nelem, nelem, 5},
+        vec = { {nelem, 0., 0.0} }
+        -- origin = {-50, -nelem-1, 0},
+        -- vec = { {0., 2.0*(nelem+1), 0.0} }
+      }
+    },
+    folder    = 'tracking/',
+    output    = {format = 'asciispatial'},  
+    time_control     = { 
+      min = { iter = 0 }, max = { iter = t_total }, interval = { iter = interval } }
+  }
+}
+
+field = { 
+  label   = 'spc1',
+  species = {diff_coeff = {(tau-0.5)/3} },
+  initial_condition = { pressure  = initial_concentration,
+                        velocityX = 0.0,
+                        velocityY = 0.0,
+                        velocityZ = 0.0 
+                      },
+  boundary_condition = { 
+                        { 
+                          label = 'circle', 
+                          -- anti-bounce back boundary condition is assigned
+                          kind = 'pressure_antiBounceBack_pasScal',
+                          pressure = c_init * cs2,
+                        }
+                       }
+}

examples/passive_scalar/benchmark/cylinder2d/postHarvest.lua

@@ -0,0 +1,59 @@
+require 'musubi'
+-- require 'params'
+
+restart   = {read = 'restart/'..simulation_name..'_lastHeader.lua'}
+variable = {
+  {
+     name = 'vel_x',
+     ncomponents = 1,
+     vartype = "operation",
+     operation = {kind='extract',
+                  input_varname={'velocity_phy'},
+                  input_varindex = {1}
+                 }
+  },
+  {
+      name = 'uxval',
+      ncomponents = 1,
+      vartype = 'st_fun',
+      st_fun = uxanaval,
+  },
+  {
+      name = 'pval',
+      ncomponents = 1,
+      vartype = 'st_fun',
+      st_fun = panaval,
+  },
+  {
+    name = 'uxdiff',
+    ncomponents = 1,
+    vartype = "operation",
+    operation = {
+      kind = 'difference',
+      input_varname = {
+        'uxval', 'vel_x'
+      },
+    },
+  },
+  {
+    name = 'pdiff',
+    ncomponents = 1,
+    vartype = "operation",
+    operation = {
+      kind = 'difference',
+      input_varname = {
+        'pval', 'pressure_phy'
+      },
+    },
+  },
+}
+
+tracking  = {
+  label     = 'main',
+  folder    = 'harvest/',
+  output    = { format = 'ascii' },
+  variable  = {'uxdiff', 'pdiff'},
+  reduction = {'l2norm', 'l2norm'},
+  -- variable  = {'velocity_phy'},
+  shape     = {kind='all'},
+}

examples/passive_scalar/benchmark/cylinder2d/resource/simulation_spc1_p00000_t17.500E+03.res

@@ -0,0 +1,42 @@
+# Rank of the process:       0
+#                   coordX                   coordY                   coordZ             spc1_density       concentration_real
+  0.4000000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7260586238429767E-001  0.2617440914218760E-001
+  0.4100000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7291150038562380E-001  0.2647236323318802E-001
+  0.4200000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7383999166467575E-001  0.2737153493488076E-001
+  0.4300000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7540044186601844E-001  0.2888781711036603E-001
+  0.4400000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7762565898098267E-001  0.3104757759350918E-001
+  0.4500000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.8054431255538549E-001  0.3388747121806701E-001
+  0.4600000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.8420949648148093E-001  0.3745416195138285E-001
+  0.4700000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.8866822625876189E-001  0.4180394139817190E-001
+  0.4800000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.9399147365576555E-001  0.4700222738106397E-001
+  0.4900000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.1002428395099528E+000  0.5312292472605962E-001
+  0.5000000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.1075069085888297E+000  0.6024762994975873E-001
+  0.5100000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.1158620282731709E+000  0.6846466239447835E-001
+  0.5200000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.1253991906242743E+000  0.7786790657898324E-001
+  0.5300000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.1362101647598444E+000  0.8855545416870392E-001
+  0.5400000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.1483827555909914E+000  0.1006280390050120E+000
+  0.5500000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.1620221155454586E+000  0.1141872649896066E+000
+  0.5600000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.1771996961748398E+000  0.1293336341502475E+000
+  0.5700000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.1940312492064198E+000  0.1461643907003999E+000
+  0.5800000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2125687908672545E+000  0.1647712060610986E+000
+  0.5900000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2329029803611678E+000  0.1852377392877257E+000
+  0.6000000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2551038172345962E+000  0.2076371167320342E+000
+  0.6100000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2791739837784324E+000  0.2320293836244058E+000
+  0.6200000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3052114067004179E+000  0.2584589881139420E+000
+  0.6300000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3331897985912062E+000  0.2869523646842796E+000
+  0.6400000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3630409586672670E+000  0.3175156883238681E+000
+  0.6500000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3949882730434113E+000  0.3501328729727197E+000
+  0.6600000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.4284807809458508E+000  0.3847638872701625E+000
+  0.6700000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.4641725675529820E+000  0.4213434572657949E+000
+  0.6800000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.5011422681668777E+000  0.4597802194245117E+000
+  0.6900000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.5399979837687606E+000  0.4999563779877385E+000
+  0.7000000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.5803131846426295E+000  0.5417279087232363E+000
+  0.7100000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.6215210517188060E+000  0.5849253366277714E+000
+  0.7200000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.6644090443088924E+000  0.6293550987025180E+000
+  0.7300000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7075663205290451E+000  0.6748014850860612E+000
+  0.7400000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7521083524988368E+000  0.7210291332898047E+000
+  0.7500000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7967147646166601E+000  0.7677860317900870E+000
+  0.7600000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.8419859312520177E+000  0.8148069715771766E+000
+  0.7700000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.8872412141562385E+000  0.8618173682226790E+000
+  0.7800000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.9326842612659030E+000  0.9085373633336671E+000
+  0.7900000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.9780624982873943E+000  0.9546861035547979E+000

examples/passive_scalar/benchmark/cylinder2d/resource/simulation_spc1_p00000_t35.000E+03.res

@@ -0,0 +1,42 @@
+# Rank of the process:       0
+#                   coordX                   coordY                   coordZ             spc1_density       concentration_real
+  0.4000000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2592159728710007E+000  0.2143884789743639E+000
+  0.4100000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2597920528281829E+000  0.2149824383129274E+000
+  0.4200000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2615240053362219E+000  0.2167643272476560E+000
+  0.4300000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2644075616869356E+000  0.2197341540019323E+000
+  0.4400000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2684462324409884E+000  0.2238918509395649E+000
+  0.4500000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2736345961848853E+000  0.2292371540769341E+000
+  0.4600000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2799755408523446E+000  0.2357694374296940E+000
+  0.4700000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2874620032162892E+000  0.2434875051769042E+000
+  0.4800000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.2960951773210355E+000  0.2523893454109081E+000
+  0.4900000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3058658511545372E+000  0.2624718499768011E+000
+  0.5000000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3167717348136568E+000  0.2737305055801280E+000
+  0.5100000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3288014105970458E+000  0.2861590619449489E+000
+  0.5200000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3419461388406246E+000  0.2997491833264820E+000
+  0.5300000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3561932370797801E+000  0.3144900901135100E+000
+  0.5400000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3715238862046530E+000  0.3303681975867612E+000
+  0.5500000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.3879266838962664E+000  0.3473667591225512E+000
+  0.5600000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.4053678349220342E+000  0.3654655212393408E+000
+  0.5700000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.4238390011061315E+000  0.3846403978730544E+000
+  0.5800000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.4432934315928271E+000  0.4048631711313234E+000
+  0.5900000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.4637086741038612E+000  0.4261012255152836E+000
+  0.6000000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.4850520865846121E+000  0.4483173222103884E+000
+  0.6100000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.5072541363157004E+000  0.4714694195372860E+000
+  0.6200000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.5303072571282453E+000  0.4955105450249587E+000
+  0.6300000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.5541352724620666E+000  0.5203887238283329E+000
+  0.6400000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.5786439920515989E+000  0.5460469673711568E+000
+  0.6500000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.6039227168382065E+000  0.5724233251643778E+000
+  0.6600000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.6296045776088262E+000  0.5994510017450291E+000
+  0.6700000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.6560569764293500E+000  0.6270585396174502E+000
+  0.6800000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.6827108861946439E+000  0.6551700679760450E+000
+  0.6900000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7099277039196251E+000  0.6837056158667838E+000
+  0.7000000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7374576144165081E+000  0.7125814873244668E+000
+  0.7100000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7649778577983826E+000  0.7417106949262501E+000
+  0.7200000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.7929895953618926E+000  0.7710034471510576E+000
+  0.7300000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.8206830967084426E+000  0.8003676839509747E+000
+  0.7400000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.8487753634847598E+000  0.8297096540452301E+000
+  0.7500000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.8765274205180027E+000  0.8589345266592576E+000
+  0.7600000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.9043615017642893E+000  0.8879470297679085E+000
+  0.7700000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.9319338831109123E+000  0.9166521063780749E+000
+  0.7800000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.9594381908599210E+000  0.9449555800138355E+000
+  0.7900000000000000E+002  0.4000000000000000E+002  0.5000000000000000E+001  0.9867932153506809E+000  0.9727648203555681E+000