Capacitor2 example.

share/examples/capacitor2/capacitor.geo

+R = 0.01;   // disk radius
+d = 0.001;  // dielectric thickness
+t = 0.0003; // disk thickness
+
+SetFactory("OpenCASCADE");
+Mesh.Algorithm3D = 10;
+
+Cylinder(1) = {0, 0, d/2, 0, 0, t, R};
+Cylinder(2) = {0, 0, -d/2, 0, 0, d, R};
+Cylinder(3) = {0, 0, -d/2-t, 0, 0, t, R};
+Sphere(4) = {0, 0, 0, 2*R, -Pi/2, Pi/2, 2*Pi};
+Coherence;
+
+MeshSize{ PointsOf{ Volume{2}; } } = 0.0003;

share/examples/capacitor2/params_capacitor.lua

+sim.name = "capacitor"              -- simulation name
+sim.dt = 1e-14                      -- timestep size
+sim.timesteps = 1001                -- num of iterations
+sim.gmsh_model = "capacitor.geo"    -- gmsh model filename
+sim.use_gpu = 0                     -- 0: cpu, 1: gpu
+sim.approx_order = 1                -- approximation order
+sim.time_integrator = "leapfrog"    -- time integration method
+postpro.silo_output_rate = 10       -- rate at which to write silo files
+postpro.cycle_print_rate = 10       -- console print rate
+
+local diel_epsr = 10    -- Permittivity of the capacitor dielectric
+
+-- Aluminum plates material parameters
+local alu = {1, 3}
+for i,v in ipairs(alu) do
+    materials[v] = {}
+    materials[v].epsilon = 1
+    materials[v].mu = 1
+    materials[v].sigma = 3.69e7
+end
+
+-- Dielectric parameters
+local diel = { 2 }
+for i,v in ipairs(diel) do
+    materials[v] = {}
+    materials[v].epsilon = diel_epsr
+    materials[v].mu = 1
+    materials[v].sigma = 0
+end
+
+-- Air parameters
+local air = { 4 }
+for i,v in ipairs(air) do
+    materials[v] = {}
+    materials[v].epsilon = 1
+    materials[v].mu = 1
+    materials[v].sigma = 0
+end
+
+-- ** Boundary conditions **
+local absorbing_bcs = {1}
+for i,v in ipairs(absorbing_bcs) do
+    bndconds[v] = {}
+    bndconds[v].kind = "impedance"
+end
+
+local R   = 0.01    -- Capacitor disk radius (must match capacitor.geo)
+local d   = 0.001   -- Dielectric thickness (must match capacitor.geo)
+local tV  = 1       -- Target capacitor voltage
+
+-- Charging pulse area in time is 16*c_T*c_A/15 
+local c_t0 =    31 * sim.dt     -- Charging pulse center
+local c_T =     30 * sim.dt     -- Charging pulse half width
+local c_A =     (15*diel_epsr*const.eps0*tV)/(16*c_T*d)
+local d_t0 =    200 * sim.dt    -- Discharging pulse center
+local d_T =     50 * sim.dt     -- Discharging pulse half width
+local d_A =     c_A * c_T / d_T -- Discharging pulse amplitude
+
+
+local do_IO = (not parallel) or (parallel and parallel.comm_rank == 0)
+
+if ( do_IO ) then
+    local E = tV/d
+    print("\x1b[1mExpected capacitor field: " .. E .. " V/m")
+    print("Current density max amplitude: " .. c_A .. " A/m^2\x1b[0m")
+end
+
+-- Define the shape of charging current
+function source_modulation(t)
+    if (t >= c_t0-c_T) and (t <= c_t0+c_T) then
+        A = c_A
+        T = c_T
+        t0 = c_t0
+        dir = 1
+    elseif (t >= d_t0-d_T) and (t <= d_t0+d_T) then
+        A = d_A
+        T = d_T
+        t0 = d_t0
+        dir = -1
+    else
+        return 0
+    end
+
+    local t2 = (t - t0)*(t-t0)
+    local t4 = t2*t2
+    local T2 = T*T
+    local T4 = T2*T2
+
+    return dir * (A*t4/T4 - 2*A*t2/T2 + A)
+end
+
+
+-- Define the function evaluating the current source
+function current_source(tag, x, y, z, t)
+    local Ex = 0.0
+    local Ey = 0.0
+    local Ez = source_modulation(t)
+
+    return Ex, Ey, Ez
+end
+
+
+-- Apply the current source to entity 2
+sources[2] = current_source
+
+-- Dump charging current profile for debug
+if ( do_IO ) then
+    local file = io.open("current_profile.dat", "w")
+    for ii = 1,sim.timesteps do
+        local t = ii*sim.dt
+        local curr = source_modulation(t)
+        file:write(t .. " " .. curr .. "\n")
+    end
+    file:close()
+end
+