defect and sort fo working mechanics

bin/Automatic_CV_holds.py

@@ -0,0 +1,130 @@
+import os
+import subprocess
+import itertools
+import configparser
+import numpy as np
+
+
+def run_MPET(cwd, config):
+    os.chdir(cwd)
+    subprocess.call(["python", os.path.join(cwd, "bin", "mpetrun.py"), config])
+
+
+def ensemble_definitions(parameters):
+    keys, vals = zip(*parameters)
+    return keys, vals
+
+k_B = 1.38064852e-23
+T = 298.15
+e = 1.60217662e-19
+
+def run_params_mpet(config_file, material_file,
+                    system_properties, material_properties, output_folder):
+    cwd = os.getcwd()
+    os.chdir(cwd)
+    if not os.path.exists(output_folder):
+        os.mkdir(output_folder)
+    os.chdir("configs")
+
+    keys_system, val_system = ensemble_definitions(system_properties)
+    cfg_sys = configparser.ConfigParser()
+    cfg_sys.optionxform = str
+    cfg_sys.read(config_file)
+    combinations_system = list(itertools.product(*val_system))
+    num_sys = len(combinations_system)
+
+    keys_mat, val_mat = ensemble_definitions(material_properties)
+    cfg_mat = configparser.ConfigParser()
+    cfg_mat.optionxform = str
+    cfg_mat.read(material_file)
+    combinations_material = list(itertools.product(*val_mat))
+    num_mat = len(combinations_material)
+    ind = 0
+
+    for comb_mat in combinations_material:
+        param_mat = dict(zip(keys_mat, comb_mat))
+        new_mat = cfg_mat
+        nicename_mat = []
+        for key, val in param_mat.items():
+            new_mat[key[0]][key[1]] = val
+            nicename_mat.append(f"{key[1]}={val}")
+        with open(material_file, "w") as f:
+            new_mat.write(f)
+
+        for combin_sys in combinations_system:
+            params_sys = dict(zip(keys_system, combin_sys))
+            new_sys = cfg_sys
+            nicename_sys = []
+            for key, val in params_sys.items():
+                new_sys[key[0]][key[1]] = val
+                if key[1] == "prevDir":
+                    continue
+                new_sys["Sim Params"]["profileType"] = "CVsegments"
+                nicename_sys.append(f"{key[1]}={val}")
+            with open(config_file, "w") as f:
+                new_sys.write(f)
+
+            os.chdir(cwd)
+            run_MPET(cwd, os.path.join(cwd, "configs", config_file))
+
+            os.chdir("history")
+            folders = os.listdir()
+            folders.sort()
+            last_folder = folders[0]
+            new_folder_name = "-".join(nicename_mat + nicename_sys)
+            new_folder_name += "-def10"
+            new_folder_path = os.path.join(cwd, output_folder, new_folder_name)
+            if os.path.exists(new_folder_path):
+                os.rename(new_folder_path, os.path.join(cwd, output_folder,
+                                                        (new_folder_name + 'old')))
+            os.rename(last_folder, new_folder_path)
+            os.chdir(cwd)
+            os.chdir("configs")
+            ind += 1
+            print(f"Simulation {ind} of {num_mat * num_sys} completed")
+
+# prev_dir = r"C:\Users\pierfrancescoo\Documents\Phase-field\mpet-LFMP\mpet\LFMP_dyn\pulses_y04\50percMn\base_50pMn"
+# prev_dir = r"C:\Users\pierfrancescoo\Documents\Phase-field\mpet-LFMP\mpet\LFP_CV\Iarchuk_1\base"
+# ocv = 3.9998
+ocv = 3.422
+etas = [1.5,2,4]
+
+holds = 3000 # sec
+holds = holds/60 # min
+etas = k_B*T/e*np.array(etas)
+segments = []
+for i in range(len(etas)):
+    Vp = str(ocv + etas[i])
+    Vm = str(ocv - etas[i])
+    holds = str(holds)
+    stringp = f"[(3.38,20),({Vp},{holds})]"
+    # stringm = f"[(3.5,3),({Vm},{holds})]"
+    # stringp = f"[({Vp},{holds})]"
+    # stringm = f"[({Vm},{holds})]"
+    segments.append(str(stringp))
+    # segments.append(str(stringm))
+
+
+system_properties = [
+    [("Sim Params","segments"), segments],
+    [("Particles","cs0_c"), ["0.5"]],
+    ]
+
+material_properties = [
+    # [("Reactions", 'k0'), ["3"]],
+    [("Reactions", 'surface_diffusion'), ["false"]],
+    [("Material", 'Bx'), ["0.1916e9"]],
+    # [("Material", 'By'), ["0.5e9"]],
+    # [("Material", 'D'), ["1e-13"]],
+    # [("Material", 'D_surf'), ["1e-17"]],
+    [("Material", 'dgammadc'), ["-177e-30","0"]],
+    # [("Material", 'kappa_x'), ["5e-10"]],
+    # [("Material", 'kappa_y'), ["15e-10"]],
+    ]
+
+
+output_folder = "LFP_CV_2D/defect_lognorm_charge"
+config_file = 'params_system_LFP_2D.cfg'
+material_file = 'params_LFP2D.cfg'
+
+run_params_mpet(config_file, material_file, system_properties, material_properties, output_folder)

configs/params_LFP2D.cfg

@@ -13,20 +13,20 @@ numnoise = 200
 Omega_a = 1.8560e-20
 kappa_x = 5e-10
 kappa_y = 15e-10
-Bx = 0.01e9
+Bx = 0.1916e9
 By = 0.5e9
-mechanics = false
+mechanics = true
 rho_s = 1.3793e28
 D = 1e-13
 E_D = 0
 Dfunc = lattice
-dgammadc = -250e-30
+dgammadc = -177e-30
 cwet = 0.98
-D_surf = 1e-18
+D_surf = 1e-15
 E_D_surf = 0
 
 [Reactions]
-surface_diffusion = true
+surface_diffusion = false
 rxnType = CIET
 k0 = 3
 E_A = 13000

configs/params_system_LFP_2D.cfg

@@ -1,16 +1,16 @@
 [Sim Params]
-profileType = CVsegments
+profileType = CV
 Crate = 1
 Vmax = 4.0
-Vmin = 2.5
-Vset = 0.12
+Vmin = 3.0
+Vset = 3.4
 power = 1
-segments = [(3.472,100)]
+segments = [(3.38,20),(3.524770281652968,50.0)]
 prevDir = false
 tend = 1.2e3
-tsteps = 2000
-relTol = 1e-7
-absTol = 1e-7
+tsteps = 1000
+relTol = 1e-6
+absTol = 1e-6
 T = 298
 nonisothermal = false
 randomSeed = true
@@ -20,10 +20,10 @@ Rser = 0.
 Nvol_c = 1
 Nvol_s = 0
 Nvol_a = 0
-Npart_c = 3
+Npart_c = 1
 Npart_a = 0
 agglomerate_c = false
-pAggl_c = 3
+pAggl_c = 1
 Naggl_c = 1
 agglomerate_a = false
 pAggl_a = 0
@@ -36,13 +36,13 @@ k0_foil = 8
 Rfilm_foil = 0
 
 [Particles]
-mean_c = 50e-9
-stddev_c = 10e-9
+mean_c = 60e-9
+stddev_c = 20e-9
 mean_a = 100e-9
 stddev_a = 1e-9
 specified_psd_c = False
 specified_psd_a = False
-cs0_c = 0.98
+cs0_c = 0.02
 cs0_a = 0.99
 
 [Conductivity]

mpet/electrode/materials/LiFePO42D.py

@@ -1,4 +1,5 @@
 import numpy as np
+from scipy.ndimage import gaussian_filter
 
 
 def LiFePO42D(self, c_mat, ybar, T, muR_ref):
@@ -16,6 +17,16 @@ def LiFePO42D(self, c_mat, ybar, T, muR_ref):
         beta_s = self.get_trode_param("beta_s")
         Ny = np.size(c_mat, 1)
         Nx = np.size(c_mat, 0)
+        # defect mat matrix Nx x Ny of random numbers between 0 and 0.06
+        # mean = 0.05
+        # stddevs = 0.03
+        # var = stddevs**2
+        # mu = np.log(mean**2/np.sqrt(var+mean**2))
+        # sigma = np.sqrt(np.log(var/mean**2 + 1))
+        # defect_mat = np.random.lognormal(mu, sigma, (Nx,Ny))
+        defect_mat = np.random.rand(Nx,Ny)*0
+        # gaussian smoothing
+        # defect_mat = gaussian_filter(defect_mat, sigma=0.5)
         dys = 1./(Ny-1)
         dxs = 1./Nx
         ywet = self.get_trode_param("cwet")*np.ones((1,Ny), dtype=object)
@@ -31,7 +42,7 @@ def LiFePO42D(self, c_mat, ybar, T, muR_ref):
         y_vert_avg = np.average(c_mat, axis=1)
         y_oriz_avg = np.average(c_mat, axis=0)
         # regular solution
-        muR_mat = T*np.log(c_mat/(1-c_mat)) + self.get_trode_param("Omega_a")*(1-2*c_mat)
+        muR_mat = T*np.log(c_mat/(1- defect_mat -c_mat)) + self.get_trode_param("Omega_a")*(1- defect_mat -2*c_mat)
         # non-homogeneous
         muR_mat += -self.get_trode_param("kappa_x")*curvx - self.get_trode_param("kappa_y")*curvy
         muR_mat += self.get_trode_param("Bx")*np.subtract(c_mat,y_oriz_avg)