diff --git a/__main__.py b/__main__.py
index 66ccadf76742dc0c49302a9e9acdd5d54485cbc7..a8b6ce06052c526bf24f0ab4aa2ad18209c5d2fa 100755
--- a/__main__.py
+++ b/__main__.py
@@ -37,6 +37,11 @@ parser_field.add_argument('-v', '--verbose', help='VERBOSE level',
action='count', default=0)
parser_field.add_argument('-cv', help='CV simulation', action="store_true")
parser_field.add_argument("-wf", help="WeightField Simulation", action="store_true")
+parser_field.add_argument("-step", help="Voltage step-by-step simulation", action="store_true")
+parser_field.add_argument("-loop", help="Voltage step-by-step simulation", action="store_true")
+parser_field.add_argument("-v_current", help="Current voltage for step-by-step simulation", type=float)
+parser_field.add_argument("-noise", help="Detector Noise simulation", action="store_true")
+
parser_fpga = subparsers.add_parser('fpga', help='FPGA design')
parser_fpga.add_argument('label', help='LABEL to identify FPGA design')
@@ -63,6 +68,7 @@ parser_tct.add_argument('-vol', '--voltage', type=str, help='bias voltage')
parser_tct.add_argument('-amp', '--amplifier', type=str, help='amplifier')
parser_tct.add_argument('-s', '--scan', type=int, help='instance number for scan mode')
+
args = parser.parse_args()
if len(sys.argv) == 1:
@@ -103,8 +109,8 @@ elif kwargs['shell'] == False: # not in shell
+ BINDPATH + " " \
+ IMGFILE + " " \
+ "python3 raser"
- subprocess.run([raser_shell+' '+command], shell=True, executable='/bin/bash')
+ subprocess.run([raser_shell+' '+command], shell=True, executable='/bin/bash')
else: # in shell
submodule = importlib.import_module(submodule)
submodule.main(kwargs)
diff --git a/field/__init__.py b/field/__init__.py
index 2e68a4572ccbd3546bbb6b2bf1216151b44bdaa4..8cf12639ba27d24bc3daac9409d11feba1a79a46 100755
--- a/field/__init__.py
+++ b/field/__init__.py
@@ -21,6 +21,7 @@ def main(kwargs):
elif label == "3d_ringcontact_Elefield":
from . import test4hsic
test4hsic.main("3d_ringcontact")
+
else:
- from . import devsim_solve
- devsim_solve.main(kwargs)
+ from . import solver_section
+ solver_section.main(kwargs)
diff --git a/field/build_device.py b/field/build_device.py
index 6a6dcb51e0acb34faa7ceec20649f781f84a1c0c..56276a598538419e7ab048c7405c8202a03892f6 100644
--- a/field/build_device.py
+++ b/field/build_device.py
@@ -27,7 +27,7 @@ class Detector:
---------
2023/12/03
"""
- def __init__(self, device_name, devsim_solve_paras=None):
+ def __init__(self, device_name, devsim_solve_paras):
self.det_name = device_name
self.device = device_name
self.region = device_name
@@ -81,7 +81,10 @@ class Detector:
if self.dimension == 1:
self.create1DMesh()
elif self.dimension == 2:
- self.create2DMesh()
+ if self.device_dict.get("mesh", {}).get("gmsh_mesh", {}):
+ self.createGmshMesh()
+ else:
+ self.create2DMesh()
elif self.dimension == 3:
self.createGmshMesh()
else:
@@ -115,7 +118,7 @@ class Detector:
mesh = self.device_dict["mesh"]["2D_mesh"]
for mesh_line in mesh["mesh_line"]:
devsim.add_2d_mesh_line(mesh=mesh_name, **mesh_line)
- if self.control_dict["ac-weightfield"] == True:
+ if (self.control_dict["weightfield"] == True) :
mesh["region"]["material"] = "gas"
else:
pass
@@ -139,7 +142,7 @@ class Detector:
mesh_name = self.device
mesh = self.device_dict["mesh"]["gmsh_mesh"]
devsim.create_gmsh_mesh (mesh=mesh_name, file=mesh['file'])
- if self.control_dict["ac-weightfield"] == True:
+ if (self.control_dict["weightfield"] == True) :
mesh["region"][0]['material']="gas"
else:
pass
diff --git a/field/devsim_draw.py b/field/devsim_draw.py
index 4a859cb96f891e896a38f44b1a5dde3e04823fe1..b412fa329e34b81f07ce32d238144f471c80ca10 100644
--- a/field/devsim_draw.py
+++ b/field/devsim_draw.py
@@ -10,6 +10,7 @@ import os
from util.output import output
def draw_iv(device,V,I,path):
+
fig2=matplotlib.pyplot.figure()
matplotlib.pyplot.semilogy(V,I)
matplotlib.pyplot.xlabel('Voltage (V)')
@@ -58,6 +59,59 @@ def draw_iv(device,V,I,path):
canvas.SaveAs(os.path.join(path, "simIV{}to{}_picture.root".format(min(V),max(V))))
canvas.SaveAs(os.path.join(path, "simIV{}to{}_picture.pdf".format(min(V),max(V))))
+
+def draw_noise(device,V,noise,path):
+ fig2=matplotlib.pyplot.figure()
+ matplotlib.pyplot.semilogy(V,noise)
+ matplotlib.pyplot.xlabel('Voltage (V)')
+ matplotlib.pyplot.ylabel('Current (A)')
+ matplotlib.pyplot.yscale('log')
+ fig2.savefig(os.path.join(path, "{}_noise.png".format(device)))
+ fig2.clear()
+
+
+ file = ROOT.TFile(os.path.join(path, "simnoise{}to{}.root".format(min(V),max(V))), "RECREATE")
+ tree = ROOT.TTree("SicarTestnoise", "SicarTest with impactgen")
+ x = array('d', [0])
+ y = array('d', [0])
+
+ tree.Branch("voltage", x, "x/D")
+ tree.Branch("Current", y, "y/D")
+
+ for point in zip(V,noise):
+ x[0], y[0] = point
+ tree.Fill()
+
+ file.Write()
+ file.Close()
+
+ file = ROOT.TFile(os.path.join(path, "simnoise{}to{}.root".format(min(V),max(V))), "READ")
+ tree = file.Get("SicarTestnoise")
+
+ graph = ROOT.TGraph(tree.GetEntries())
+ for i, entry in enumerate(tree):
+ x = entry.x
+ y = entry.y
+ graph.SetPoint(i, x, y)
+
+ canvas = ROOT.TCanvas("canvas", "Graph", 800, 600)
+ graph.SetMarkerStyle(ROOT.kFullCircle)
+ graph.SetMarkerSize(0.5)
+ graph.SetMarkerColor(ROOT.kBlue)
+ graph.SetLineColor(ROOT.kWhite)
+ graph.Draw("AP")
+
+ graph.SetTitle("NoiseCurrent vs Voltage")
+ graph.GetXaxis().SetTitle("Voltage(V)")
+ graph.GetYaxis().SetTitle("NoiseCurrent(A)")
+
+ canvas.Update()
+ canvas.SaveAs(os.path.join(path, "simnoise{}to{}_picture.root".format(min(V),max(V))))
+ canvas.SaveAs(os.path.join(path, "simnoise{}to{}_picture.pdf".format(min(V),max(V))))
+
+
+
+
def draw_cv(device,V,C,path):
fig3=matplotlib.pyplot.figure(num=4,figsize=(4,4))
matplotlib.pyplot.plot(V, C)
diff --git a/field/devsim_solve.py b/field/devsim_solve.py
index 8b4ebb835c3a6c11a12ca73f9ae36ce4fe2d7527..9e7c5691f0a11c2a2206c8925b0199058b52c24d 100644
--- a/field/devsim_solve.py
+++ b/field/devsim_solve.py
@@ -2,7 +2,10 @@
# -*- encoding: utf-8 -*-
import devsim
+import subprocess
+import shlex
from .build_device import Detector
+from . import control_step
from . import model_create
from . import physics_drift_diffusion
from . import initial
@@ -38,21 +41,41 @@ paras = {
"voltage_step" : 1,
"acreal" : 1.0,
"acimag" : 0.0,
- "frequency" : 10.0,
+ "frequency" : 1000.0,
"Cylindrical_coordinate": False,
+
+
"ac-weightfield" : False,
+
+
+ "Voltage-step-model" : False,
+ "step":1,
+
+ "Mixed-model" : False,
}
def main(kwargs):
simname = kwargs['label']
is_cv = kwargs['cv']
is_wf = kwargs["wf"]
+ is_step = kwargs["step"]
+ is_Mixed = kwargs["mixed"]
if is_wf:
paras.update({"ac-weightfield": True})
else:
paras.update({"ac-weightfield": False})
+
+ if is_step:
+ paras.update({"Voltage-step-model": True})
+ else:
+ paras.update({"Voltage-step-model": False})
+ # Mix model is applied for weightingfield simulation of large work
+ if is_Mixed:
+ paras.update({"Mixed-model": True})
+ else:
+ paras.update({"Mixed-model": False})
devsim.open_db(filename="./output/field/SICARDB.db", permission="readonly")
@@ -79,6 +102,7 @@ def main(kwargs):
devsim.add_db_entry(material="Silicon", parameter="n_i", value=N_i, unit="/cm^3", description="Intrinsic Electron Concentration")
devsim.add_db_entry(material="SiliconCarbide", parameter="n_i", value=N_i, unit="/cm^3", description="Intrinsic Electron Concentration")
devsim.add_db_entry(material="gas", parameter="n_i", value="1e-9", unit="/cm^3", description="Intrinsic Electron Concentration")
+ devsim.add_db_entry(material="gas", parameter="Permittivity", value="1", unit="1", description="Permittivity")
devsim.add_db_entry(material="Silicon", parameter="n1", value=N_i, unit="/cm^3", description="n1")
devsim.add_db_entry(material="Silicon", parameter="p1", value=N_i, unit="/cm^3", description="p1")
@@ -145,7 +169,7 @@ def main(kwargs):
impact_label=MyDetector.device_dict['avalanche_model']
else:
impact_label=None
- if paras["ac-weightfield"] == True:
+ if (paras["ac-weightfield"] == True) or (paras["Mixed-model"] == True):
pass
else:
initial.DriftDiffusionInitialSolution(device, region, paras,irradiation_label=irradiation_label, irradiation_flux=irradiation_flux, impact_label=impact_label, set_contact_type=None,circuit_contacts=circuit_contacts)
@@ -185,7 +209,8 @@ def main(kwargs):
header_cv = ["Voltage","Capacitance"]
writer_cv = csv.writer(f_cv)
writer_cv.writerow(header_cv)
-
+ if is_step or is_Mixed == True:
+ v_goal = MyDetector.device_dict['bias']['voltage']
v_max = MyDetector.device_dict['bias']['voltage']
area_factor = MyDetector.device_dict['area_factor']
frequency = paras['frequency']
@@ -197,23 +222,56 @@ def main(kwargs):
voltage_step = -1 * paras['voltage_step']
while abs(v) <= abs(v_max):
voltage.append(v)
- if paras["ac-weightfield"]==True:
+ if (paras["ac-weightfield"]==True) or (paras["Mixed-model"]==True):
v=1
for contact in circuit_contacts:
- print(type(circuit_contacts))
print("+++++++++++++++++++++\n begin simulate Weight field\n +++++++++++++++++++++")
print(contact)
devsim.set_parameter(name="debug_level", value="info")
- devsim.set_parameter(device=device, name=physics_drift_diffusion.GetContactBiasName(contact), value=v)
- devsim.solve(type="dc", absolute_error=paras['absolute_error_VoltageSteps'], relative_error=paras['relative_error_VoltageSteps'], maximum_iterations=paras['maximum_iterations_VoltageSteps'])
- paras["milestone_step"] == 1
- paras.update({"milestone_step":paras["milestone_step"]})
- path = output(__file__, device,contact)
- milestone_save_wf_2D(device, region, v, path,contact)
- devsim.set_parameter(device=device, name=physics_drift_diffusion.GetContactBiasName(contact), value=0)
+ if (paras["ac-weightfield"]==True and paras["Mixed-model"]==False):
+ devsim.set_parameter(device=device, name=physics_drift_diffusion.GetContactBiasName(contact), value=v)
+ devsim.solve(type="dc", absolute_error=paras['absolute_error_VoltageSteps'], relative_error=paras['relative_error_VoltageSteps'], maximum_iterations=paras['maximum_iterations_VoltageSteps'])
+ paras["milestone_step"] == 1
+ paras.update({"milestone_step":paras["milestone_step"]})
+ path = output(__file__, device,contact)
+ milestone_save_wf_2D(device, region, v, path,contact)
+ devsim.set_parameter(device=device, name=physics_drift_diffusion.GetContactBiasName(contact), value=0)
+ if paras["Mixed-model"]==True:
+ print("This is design for large work in shell\n=============================\n test for refine mesh")
+ while v <=1:
+ # if v != 0:
+ # control_step.set_values(device, region)
+ # else:
+ # pass
+ command_list = ["devsim.set_parameter(device=device, name=physics_drift_diffusion.GetContactBiasName(contact), value=v)",
+ "devsim.solve(type='dc', absolute_error=paras['absolute_error_VoltageSteps'], relative_error=paras['relative_error_VoltageSteps'], maximum_iterations=paras['maximum_iterations_VoltageSteps'])",
+ "paras['milestone_step'] == 1",
+ "control_step.save_values(device,region)",
+ "paras.update({'milestone_step':paras['milestone_step']})",
+ "path = output(__file__, device,contact)"
+ "try_save(device, region, v, path,contact)"
+ ]
+ for command in command_list:
+ print(command)
+ process = subprocess.Popen(shlex.split(command),stdout=subprocess.PIPE)
+ while True:
+ output_info = process.stdout.readline().decode().strip() # type: ignore
+ if output_info:
+ print(output_info)
+ else:
+ break
+ process.wait()
+ # devsim.set_parameter(device=device, name=physics_drift_diffusion.GetContactBiasName(contact), value=v)
+ # devsim.solve(type="dc", absolute_error=paras['absolute_error_VoltageSteps'], relative_error=paras['relative_error_VoltageSteps'], maximum_iterations=paras['maximum_iterations_VoltageSteps'])
+ # paras["milestone_step"] == 1
+ # paras.update({"milestone_step":paras["milestone_step"]})
+ # path = output(__file__, device,contact)
+
+ print("+++++++++++++++++++++\n Simulation of Weight field is over !>.<!\n +++++++++++++++++++++")
exit()
elif paras["ac-weightfield"] ==False:
pass
+
devsim.set_parameter(device=device, name=physics_drift_diffusion.GetContactBiasName(circuit_contacts), value=v)
devsim.solve(type="dc", absolute_error=paras['absolute_error_VoltageSteps'], relative_error=paras['relative_error_VoltageSteps'], maximum_iterations=paras['maximum_iterations_VoltageSteps'])
physics_drift_diffusion.PrintCurrents(device, circuit_contacts)
@@ -241,7 +299,10 @@ def main(kwargs):
if MyDetector.dimension == 1:
milestone_save_1D(device, region, v, path)
elif MyDetector.dimension == 2:
- milestone_save_2D(device, region, v, path)
+ if MyDetector.device_dict.get("mesh", {}).get("gmsh_mesh", {}):
+ milestone_save_2D(device, region, v, path)
+ else:
+ milestone_save_2D(device, region, v, path)
elif MyDetector.dimension == 3:
milestone_save_3D(device, region, v, path)
else:
@@ -340,8 +401,8 @@ def milestone_save_2D(device, region, v, path):
metadata['dimension'] = 2
names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
- if v == 0:
- names.append('NetDoping')
+ # if v == 0:
+ # names.append('NetDoping')
for name in names: # scalar field on mesh point (instead of on edge)
with open(os.path.join(path, "{}_{}V.pkl".format(name,v)),'wb') as file:
@@ -410,6 +471,37 @@ def milestone_save_3D(device, region, v, path):
data['metadata'] = metadata
pickle.dump(data, file)
+
+
+def milestone_save_wf_3D(device, region, v, path,contact):
+ x=devsim.get_node_model_values(device=device,region=region,name="x")
+ y=devsim.get_node_model_values(device=device,region=region,name="y")
+ z=devsim.get_node_model_values(device=device,region=region,name="z")
+ Potential=devsim.get_node_model_values(device=device,region=region,name="Potential")
+
+ metadata = {}
+ metadata['voltage'] = v
+ metadata['dimension'] = 3
+
+ names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
+ if v == 0:
+ names.append('NetDoping')
+
+ for name in ['Potential']: # scalar field on mesh point (instead of on edge)
+ with open(os.path.join(path, "{}_{}_{}V.pkl".format(name,v,contact)),'wb') as file:
+ data = {}
+ data['values'] = eval(name) # refer to the object with given name
+ merged_list = [x, y]
+ transposed_list = list(map(list, zip(*merged_list)))
+ data['points'] = transposed_list
+ data['metadata'] = metadata
+ pickle.dump(data, file)
+def try_save(device, region, v, path,contact):
+ try:
+ milestone_save_wf_2D(device, region, v, path,contact)
+ except:
+ milestone_save_wf_3D(device, region, v, path,contact)
+
if __name__ == "__main__":
args = sys.argv[1:]
kwargs = {}
diff --git a/field/initial.py b/field/initial.py
index ce109d0e391bb7bcb0091d822b4bda285be6f7d0..330a71c7fd4d2cebb0e3ae386707164d7caab580 100644
--- a/field/initial.py
+++ b/field/initial.py
@@ -25,7 +25,7 @@ def switch_Cylindrical_coordinate(device,region):
devsim.set_parameter(name="element_node1_volume_model",value="ElementCylindricalNodeVolume@en1")
-def InitialSolution(device, region, circuit_contacts=None, set_contact_type=None, paras=None):
+def InitialSolution(device, region,paras, circuit_contacts, set_contact_type=None):
if paras["Cylindrical_coordinate"]==True:
switch_Cylindrical_coordinate(device,region)
else:
@@ -38,44 +38,43 @@ def InitialSolution(device, region, circuit_contacts=None, set_contact_type=None
devsim.edge_from_node_model(device=device,region=region,node_model="InitialElectron")
devsim.edge_from_node_model(device=device,region=region,node_model="InitialHole")
CreateSiliconPotentialOnly(device, region)
- if paras["ac-weightfield"]==True:
+ if paras["weightfield"]==True:
try:
CreateOxidePotentialOnly(device=device, region="SiO2", update_type="default")
for interface in devsim.get_interface_list(device=device):
CreateSiliconOxideInterface(device=device, interface=interface)
except:
- print("Waring info:++++++++++++++++++++++++++++++++++++++++++/nWarning: There is no SiO2 layer in your detector\n++++++++++++++++++++++++++++++++++++++")
+ print("================RASER info===============\nThere is no SiO2 layer in your detector\n===========warning============")
pass
+ else:
+ pass
# Set up the contacts applying a bias
for i in devsim.get_contact_list(device=device):
if set_contact_type and i in set_contact_type:
contact_type = set_contact_type[i]
else:
contact_type = {"type" : "Ohmic"}
-
- devsim.set_parameter(device=device, name=GetContactBiasName(i), value="0.0")
- #if circuit_contacts and i in circuit_contacts:
+ devsim.set_parameter(device=device, name=GetContactBiasName(i), value=0)
if str(circuit_contacts) in i :
CreateSiliconPotentialOnlyContact(device, region, i, contact_type,True)
- if paras["ac-weightfield"]==True:
+ if paras["weightfield"]==True:
try:
CreateOxideContact(device=device, region="SiO2", contact=i)
except:
- print("Waring info:++++++++++++++++++++++++++++++++++++++++++/nWarning: There is no SiO2 layer in your detector\n++++++++++++++++++++++++++++++++++++++")
+ print("===============RASER info===============\nThere is no SiO2 layer in your detector\n===========warning============")
pass
else:
###print "FIX THIS"
### it is more correct for the bias to be 0, and it looks like there is side effects
devsim.set_parameter(device=device, name=GetContactBiasName(i), value="0.0")
CreateSiliconPotentialOnlyContact(device, region, i, contact_type)
- if paras["ac-weightfield"]==True:
+ if paras["weightfield"]==True:
try:
CreateOxideContact(device=device, region="SiO2", contact=i)
except:
print("Waring info:++++++++++++++++++++++++++++++++++++++++++/nWarning: There is no SiO2 layer in your detector\n++++++++++++++++++++++++++++++++++++++")
pass
-
def DriftDiffusionInitialSolution(device, region, paras, irradiation_label=None, irradiation_flux=1e15, impact_label=None, circuit_contacts=None, set_contact_type=None):
if paras["Cylindrical_coordinate"]==True:
switch_Cylindrical_coordinate(device,region)
diff --git a/field/loop_section.py b/field/loop_section.py
new file mode 100644
index 0000000000000000000000000000000000000000..2a4ef00b6330b0dc5bc68bf3afb4ffa7bdca5437
--- /dev/null
+++ b/field/loop_section.py
@@ -0,0 +1,142 @@
+#!/usr/bin/env python3
+# -*- encoding: utf-8 -*-
+import devsim
+import pickle
+import os
+
+"""
+this part is for loop,
+control by paras
+
+"""
+from . import initial
+from . import physics_drift_diffusion
+
+import numpy as np
+
+
+class loop_section():
+ def __init__(self, paras,device,region,solve_model,irradiation):
+ self.paras = paras
+ self.step_model = False
+ self.solve_model = solve_model
+ self.device = device
+ self.region = region
+ self.irradiation = False
+ self.voltage = []
+ self.current = []
+ self.capacitance = []
+ self.noise = []
+ self.voltage_milestone = []
+ self.positions_mid = []
+ self.intensities = []
+
+ self.positions = []
+ self.electrons = []
+ self.holes = []
+
+ def initial_solver(self,contact,set_contact_type,irradiation_label,irradiation_flux,impact_label):
+ initial.InitialSolution(device=self.device, region=self.region, circuit_contacts=contact,paras=self.paras,set_contact_type=set_contact_type)
+ devsim.solve(type="dc", absolute_error=self.paras['absolute_error_Initial'], relative_error=self.paras['relative_error_Initial'], maximum_iterations=self.paras['maximum_iterations_Initial'])
+ print("======================\nFirst initialize successfully\n===============================")
+ if self.irradiation == False:
+ if self.solve_model == "wf":
+ pass
+ else:
+ print("======RASER info ===========\nNo radiation\n================info=================")
+
+ initial.DriftDiffusionInitialSolution(device=self.device, region=self.region, circuit_contacts=contact,paras=self.paras,set_contact_type=set_contact_type,
+ irradiation_label=None,irradiation_flux=0,impact_label=impact_label)
+ devsim.solve(type="dc", absolute_error=self.paras['absolute_error_Initial'], relative_error=self.paras['relative_error_Initial'], maximum_iterations=self.paras['maximum_iterations_Initial'])
+ elif self.irradiation == True:
+ print("======RASER info ===========\nradiation\n================info=================")
+ initial.DriftDiffusionInitialSolution(device=self.device, region=self.region, circuit_contacts=contact,paras=self.paras,set_contact_type=set_contact_type,
+ irradiation_label=irradiation_label,irradiation_flux=irradiation_flux,impact_label=impact_label)
+ devsim.solve(type="dc", absolute_error=self.paras['absolute_error_Initial'], relative_error=self.paras['relative_error_Initial'], maximum_iterations=self.paras['maximum_iterations_Initial'])
+ else:
+ print("======RASER info ===========\nirradiation should set as False or True\n================Error=================")
+ exit()
+ print("=====================\nDriftDiffusion initialize successfully\n======================")
+ print("=========RASER info =========\nAll initialization successfully\n=========info========== ")
+
+
+
+
+ def save_values(self,device, region):
+ Holes_values = devsim.get_node_model_values(device=device, region=region, name="Holes")
+ Electrons_values = devsim.get_node_model_values(device=device, region=region, name="Electrons")
+ Potential_values = devsim.get_node_model_values(device=device, region=region, name="Potential")
+ with open('./param_file/step-model/{}_Holes.pkl'.format(device), 'wb') as file:
+ file.truncate(0)
+ with open('./param_file/step-model/{}_Holes.pkl'.format(device), 'wb') as file:
+ pickle.dump(Holes_values, file)
+ with open('./param_file/step-model/{}_Electrons.pkl'.format(device), 'wb') as file:
+ file.truncate(0)
+ with open('./param_file/step-model/{}_Electrons.pkl'.format(device), 'wb') as file:
+ pickle.dump(Electrons_values, file)
+ with open('./param_file/step-model/{}_Potential.pkl'.format(device), 'wb') as file:
+ file.truncate(0)
+ with open('./param_file/step-model/{}_Potential.pkl'.format(device), 'wb') as file:
+ pickle.dump(Potential_values, file)
+
+ def load_values(self,values):
+ if values=="Holes":
+ with open('./param_file/step-model/{}_Holes.pkl'.format(self.device), 'rb') as file:
+ return pickle.load(file)
+ elif values=="Electrons":
+ with open('./param_file/step-model/{}_Electrons.pkl'.format(self.device), 'rb') as file:
+ return pickle.load(file)
+ elif values=="Potential":
+ with open('./param_file/step-model/{}_Potential.pkl'.format(self.device), 'rb') as file:
+ return pickle.load(file)
+
+ def set_values(self,device, region):
+ for i in ("Holes","Electrons","Potential"):
+ value = self.load_values(i)
+ devsim.set_node_values(device=device, region=region,name=i,values=value)
+
+
+ def loop_solver(self,circuit_contact,v_current,area_factor,path,device,region):
+ if self.solve_model =="step":
+ if v_current == 0:
+ pass
+ else:
+ print("=================RASER info==================\n Load last voltage successfully\n===============info===================")
+ self.set_values(device=device,region=region)
+
+ self.voltage.append(v_current)
+ devsim.set_parameter(device=self.device, name=physics_drift_diffusion.GetContactBiasName(circuit_contact), value=v_current)
+ devsim.solve(type="dc", absolute_error=self.paras['absolute_error_VoltageSteps'], relative_error=self.paras['relative_error_VoltageSteps'], maximum_iterations=self.paras['maximum_iterations_VoltageSteps'])
+ if self.solve_model !="wf":
+ physics_drift_diffusion.PrintCurrents(device=self.device, contact=circuit_contact)
+ electron_current= devsim.get_contact_current(device=self.device, contact=circuit_contact, equation="ElectronContinuityEquation")
+ hole_current = devsim.get_contact_current(device=self.device, contact=circuit_contact, equation="HoleContinuityEquation")
+ total_current = electron_current + hole_current
+ if(abs(total_current/area_factor)>105e-6):
+ print("==========RASER info===========\nCurrent is too large !\n==============Warning==========")
+ # break
+ self.current.append(total_current)
+ print(self.current)
+ if self.solve_model == "cv":
+ devsim.circuit_alter(name="V1", value=v_current)
+ devsim.solve(type="ac", frequency=self.paras["frequency"])
+ cap=1e12*devsim.get_circuit_node_value(node="V1.I", solution="ssac_imag")/ (-2*np.pi*self.paras["frequency"])
+ self.capacitance.append(cap*area_factor)
+ if self.solve_model == "noise":
+ devsim.solve(type="noise", frequency=self.paras["frequency"],output_node="V1.I")
+ noise=devsim.get_circuit_node_value(node="V1.I")
+ self.noise.append(noise)
+ self.save_values(device,region)
+
+
+
+
+ def get_voltage_values(self):
+ return self.voltage
+ def get_current_values(self):
+ return self.current
+ def get_cap_values(self):
+ return self.capacitance
+ def get_noise_values(self):
+ return self.noise
+
diff --git a/field/physics_drift_diffusion.py b/field/physics_drift_diffusion.py
index 70c75126cb7a41c2a84f6b71322e4e0fcf424423..42ac0c20042cb04567c5736a444c02a945c0225e 100644
--- a/field/physics_drift_diffusion.py
+++ b/field/physics_drift_diffusion.py
@@ -101,7 +101,7 @@ def CreateSiliconPotentialOnlyContact(device, region, contact, contact_type, is_
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
if is_circuit:
- # add_circuit_node(name=str(GetContactBiasName(contact)),value=0)
+ # add_circuit_node(name=str(GetContactBiasName(contact)),value=0)
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")
if is_circuit:
diff --git a/field/save_milestone.py b/field/save_milestone.py
new file mode 100644
index 0000000000000000000000000000000000000000..ed547d47d115e93417c701ff82db92e759d48911
--- /dev/null
+++ b/field/save_milestone.py
@@ -0,0 +1,201 @@
+
+
+import devsim
+import numpy as np
+import os
+import pickle
+from . import devsim_draw
+
+def milestone_save_1D(device, region, v, path):
+ x = np.array(devsim.get_node_model_values(device=device, region=region, name="x"))
+ Potential = np.array(devsim.get_node_model_values(device=device, region=region, name="Potential")) # get the potential dat
+ NetDoping= np.array(devsim.get_node_model_values(device=device, region=region, name="NetDoping"))
+ PotentialNodeCharge = np.array(devsim.get_node_model_values(device=device, region=region, name="PotentialNodeCharge"))
+ Electrons = np.array(devsim.get_node_model_values(device=device, region=region, name="Electrons"))
+ Holes = np.array(devsim.get_node_model_values(device=device, region=region, name="Holes"))
+ devsim.edge_average_model(device=device, region=region, node_model="x", edge_model="xmid")
+ x_mid = devsim.get_edge_model_values(device=device, region=region, name="xmid") # get x-node values
+ ElectricField = devsim.get_edge_model_values(device=device, region=region, name="ElectricField") # get y-node values
+ TrappingRate_n = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_n"))
+ TrappingRate_p = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_p"))
+
+ devsim_draw.draw1D(x,Potential,"Potential","Depth[cm]","Potential[V]", v, path)
+ devsim_draw.draw1D(x_mid,ElectricField,"ElectricField","Depth[cm]","ElectricField[V/cm]",v, path)
+ devsim_draw.draw1D(x,TrappingRate_n,"TrappingRate_n","Depth[cm]","TrappingRate_n[s]",v, path)
+ devsim_draw.draw1D(x,TrappingRate_p,"TrappingRate_p","Depth[cm]","TrappingRate_p[s]",v, path)
+
+ dd = os.path.join(path, str(v)+'V.dd')
+ devsim.write_devices(file=dd, type="tecplot")
+
+ metadata = {}
+ metadata['voltage'] = v
+ metadata['dimension'] = 1
+
+ names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
+ if v == 0:
+ names.append('NetDoping')
+
+ for name in names: # scalar field on mesh point (instead of on edge)
+ with open(os.path.join(path, "{}_{}V.pkl".format(name,v)),'wb') as file:
+ data = {}
+ data['values'] = eval(name) # refer to the object with given name
+ data['points'] = x
+ data['metadata'] = metadata
+ pickle.dump(data, file)
+
+def milestone_save_2D(device, region, v, path):
+ x = np.array(devsim.get_node_model_values(device=device, region=region, name="x")) # get x-node values
+ y = np.array(devsim.get_node_model_values(device=device, region=region, name="y")) # get y-node values
+ Potential = np.array(devsim.get_node_model_values(device=device, region=region, name="Potential")) # get the potential data
+ TrappingRate_n = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_n"))
+ TrappingRate_p = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_p"))
+
+ devsim.element_from_edge_model(edge_model="ElectricField", device=device, region=region)
+ devsim.edge_average_model(device=device, region=region, node_model="x", edge_model="xmid")
+ devsim.edge_average_model(device=device, region=region, node_model="y", edge_model="ymid")
+ ElectricField=np.array(devsim.get_edge_model_values(device=device, region=region, name="ElectricField"))
+ x_mid = np.array(devsim.get_edge_model_values(device=device, region=region, name="xmid"))
+ y_mid = np.array(devsim.get_edge_model_values(device=device, region=region, name="ymid"))
+
+ devsim_draw.draw2D(x,y,Potential,"Potential",v, path)
+ devsim_draw.draw2D(x_mid,y_mid,ElectricField,"ElectricField",v, path)
+ devsim_draw.draw2D(x,y,TrappingRate_n,"TrappingRate_n",v, path)
+ devsim_draw.draw2D(x,y,TrappingRate_p,"TrappingRate_p",v, path)
+
+ dd = os.path.join(path, str(v)+'V.dd')
+ devsim.write_devices(file=dd, type="tecplot")
+
+ metadata = {}
+ metadata['voltage'] = v
+ metadata['dimension'] = 2
+
+ names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
+ # if v == 0:
+ # names.append('NetDoping')
+
+ for name in names: # scalar field on mesh point (instead of on edge)
+ with open(os.path.join(path, "{}_{}V.pkl".format(name,v)),'wb') as file:
+ data = {}
+ data['values'] = eval(name) # refer to the object with given name
+ merged_list = [x, y]
+ transposed_list = list(map(list, zip(*merged_list)))
+ data['points'] = transposed_list
+ data['metadata'] = metadata
+ pickle.dump(data, file)
+
+
+def milestone_save_wf_2D(device, region, v, path,contact):
+ x = np.array(devsim.get_node_model_values(device=device, region=region, name="x")) # get x-node values
+ y = np.array(devsim.get_node_model_values(device=device, region=region, name="y")) # get y-node values
+ Potential = np.array(devsim.get_node_model_values(device=device, region=region, name="Potential")) # get the potential data
+
+ devsim.element_from_edge_model(edge_model="ElectricField", device=device, region=region)
+ devsim.edge_average_model(device=device, region=region, node_model="x", edge_model="xmid")
+ devsim.edge_average_model(device=device, region=region, node_model="y", edge_model="ymid")
+ ElectricField=np.array(devsim.get_edge_model_values(device=device, region=region, name="ElectricField"))
+ x_mid = np.array(devsim.get_edge_model_values(device=device, region=region, name="xmid"))
+ y_mid = np.array(devsim.get_edge_model_values(device=device, region=region, name="ymid"))
+
+ devsim_draw.draw2D(x,y,Potential,"Potential",v, path)
+ devsim_draw.draw2D(x_mid,y_mid,ElectricField,"ElectricField",v, path)
+
+ dd = os.path.join(path, str(v),str(contact)+'V.dd')
+ devsim.write_devices(file=dd, type="tecplot")
+
+ metadata = {}
+ metadata['voltage'] = v
+ metadata['dimension'] = 2
+
+ for name in ['Potential']: # scalar field on mesh point (instead of on edge)
+ with open(os.path.join(path, "{}_{}_{}V.pkl".format(name,v,contact)),'wb') as file:
+ data = {}
+ data['values'] = eval(name) # refer to the object with given name
+ merged_list = [x, y]
+ transposed_list = list(map(list, zip(*merged_list)))
+ data['points'] = transposed_list
+ data['metadata'] = metadata
+ pickle.dump(data, file)
+
+def milestone_save_3D(device, region, v, path):
+ x = np.array(devsim.get_node_model_values(device=device, region=region, name="x")) # get x-node values
+ y = np.array(devsim.get_node_model_values(device=device, region=region, name="y")) # get y-node values
+ Potential = np.array(devsim.get_node_model_values(device=device, region=region, name="Potential")) # get the potential data
+ TrappingRate_n = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_n"))
+ TrappingRate_p = np.array(devsim.get_node_model_values(device=device, region=region, name="TrappingRate_p"))
+
+ devsim.element_from_edge_model(edge_model="ElectricField", device=device, region=region)
+ devsim.edge_average_model(device=device, region=region, node_model="x", edge_model="xmid")
+ devsim.edge_average_model(device=device, region=region, node_model="y", edge_model="ymid")
+ ElectricField=np.array(devsim.get_edge_model_values(device=device, region=region, name="ElectricField"))
+ x_mid = np.array(devsim.get_edge_model_values(device=device, region=region, name="xmid"))
+ y_mid = np.array(devsim.get_edge_model_values(device=device, region=region, name="ymid"))
+
+ devsim_draw.draw2D(x,y,Potential,"Potential",v, path)
+ devsim_draw.draw2D(x_mid,y_mid,ElectricField,"ElectricField",v, path)
+ devsim_draw.draw2D(x,y,TrappingRate_n,"TrappingRate_n",v, path)
+ devsim_draw.draw2D(x,y,TrappingRate_p,"TrappingRate_p",v, path)
+
+ dd = os.path.join(path, str(v)+'V.dd')
+ devsim.write_devices(file=dd, type="tecplot")
+
+ metadata = {}
+ metadata['voltage'] = v
+ metadata['dimension'] = 2
+
+ names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
+ # if v == 0:
+ # names.append('NetDoping')
+
+ for name in names: # scalar field on mesh point (instead of on edge)
+ with open(os.path.join(path, "{}_{}V.pkl".format(name,v)),'wb') as file:
+ data = {}
+ data['values'] = eval(name) # refer to the object with given name
+ merged_list = [x, y]
+ transposed_list = list(map(list, zip(*merged_list)))
+ data['points'] = transposed_list
+ data['metadata'] = metadata
+ pickle.dump(data, file)
+
+def milestone_save_wf_3D(device, region, v, path,contact):
+ x=devsim.get_node_model_values(device=device,region=region,name="x")
+ y=devsim.get_node_model_values(device=device,region=region,name="y")
+ z=devsim.get_node_model_values(device=device,region=region,name="z")
+ Potential=devsim.get_node_model_values(device=device,region=region,name="Potential")
+
+ metadata = {}
+ metadata['voltage'] = v
+ metadata['dimension'] = 3
+
+ names = ['Potential', 'TrappingRate_p', 'TrappingRate_n']
+ if v == 0:
+ names.append('NetDoping')
+
+ for name in ['Potential']: # scalar field on mesh point (instead of on edge)
+ with open(os.path.join(path, "{}_{}_{}V.pkl".format(name,v,contact)),'wb') as file:
+ data = {}
+ data['values'] = eval(name) # refer to the object with given name
+ merged_list = [x, y]
+ transposed_list = list(map(list, zip(*merged_list)))
+ data['points'] = transposed_list
+ data['metadata'] = metadata
+ pickle.dump(data, file)
+
+
+
+def save_milestone(device, region, v, path,dimension,contact,is_wf):
+ if dimension ==1 :
+ milestone_save_1D(device, region, v, path)
+ if dimension == 2:
+ if is_wf == True:
+ milestone_save_wf_2D(device, region, v, path,contact)
+ elif is_wf == False:
+ milestone_save_2D(device, region, v, path)
+ else:
+ print("==========RASER info ==========\nis_wf only has 2 values, True or False\n==========Error=========")
+ if dimension == 3:
+ if is_wf == True:
+ milestone_save_wf_3D(device, region, v, path,contact)
+ elif is_wf == False:
+ milestone_save_3D(device, region, v, path)
+ else:
+ print("==========RASER info ==========\nis_wf only has 2 values, True or False\n==========Error=========")
\ No newline at end of file
diff --git a/field/solver_section.py b/field/solver_section.py
new file mode 100644
index 0000000000000000000000000000000000000000..3bea26c0c2bc03a3edda77fbc436097dfc7ca544
--- /dev/null
+++ b/field/solver_section.py
@@ -0,0 +1,269 @@
+#!/usr/bin/env python3
+# -*- encoding: utf-8 -*-
+import sys
+import os
+import subprocess
+import time
+import math
+
+import devsim
+import numpy as np
+
+from .build_device import Detector
+from . import model_create
+from . import save_milestone
+from . import loop_section
+from . import physics_drift_diffusion
+from util.output import output
+from .devsim_draw import *
+import multiprocessing
+
+
+v_current = 0
+V = []
+c = []
+Current = []
+noise = []
+
+paras = {
+ "absolute_error_Initial" : 1e10,
+ "relative_error_Initial" : 1e-3,
+ "maximum_iterations_Initial" : 1000,
+
+ "absolute_error_DriftDiffusion" : 1e20,
+ "relative_error_DriftDiffusion" : 1e-3,
+ "maximum_iterations_DriftDiffusion" : 1000,
+
+ "absolute_error_VoltageSteps" : 1e20,
+ "relative_error_VoltageSteps" : 1e-3,
+ "maximum_iterations_VoltageSteps" : 1000,
+
+ "milestone_mode" : True,
+ "milestone_step" : 50,
+
+ "voltage_step" : 0.1,
+ "acreal" : 1.0,
+ "acimag" : 0.0,
+ "frequency" : 1000.0,
+
+ "Cylindrical_coordinate": False,
+
+
+ "ac-weightfield" : False,
+
+
+ "Voltage-step-model" : False,
+ "step":1,
+
+}
+os.environ["OMP_NUM_THREADS"] = "1"
+def main (kwargs):
+ simname = kwargs['label']
+ is_cv = kwargs['cv']
+ is_wf = kwargs["wf"]
+ is_step = kwargs["step"]
+ is_noise = kwargs["noise"]
+
+ if is_wf:
+ paras.update({"weightfield": True})
+ else:
+ paras.update({"weightfield": False})
+
+ if is_step:
+ paras.update({"Voltage-step-model": True})
+ else:
+ paras.update({"Voltage-step-model": False})
+
+
+
+ device = simname
+ region = simname
+ MyDetector = Detector(device, devsim_solve_paras=paras)
+ MyDetector.mesh_define()
+
+ if "frequency" in MyDetector.device_dict:
+ paras.update({"frequency": MyDetector.device_dict['frequency']})
+ if "area_factor" in MyDetector.device_dict:
+ paras.update({"area_factor": MyDetector.device_dict['area_factor']})
+ if "default_dimension" in MyDetector.device_dict:
+ default_dimension =MyDetector.device_dict["default_dimension"]
+ if "irradiation" in MyDetector.device_dict:
+ irradiation = True
+ else:
+ irradiation = False
+
+ devsim.open_db(filename="./output/field/SICARDB.db", permission="readonly")
+
+ T = MyDetector.device_dict['temperature']
+ k = 1.3806503e-23 # J/K
+ q = 1.60217646e-19 # coul
+ devsim.add_db_entry(material="global", parameter="T", value=T, unit="K", description="T")
+ devsim.add_db_entry(material="global", parameter="k_T", value=k*T, unit="J", description="k*T")
+ devsim.add_db_entry(material="global", parameter="Volt_thermal", value=k*T/q, unit="J/coul", description="k*T/q")
+ N_c=2.82e19*pow(T/300,1.5)
+ N_v=1.83e19*pow(T/300,1.5)
+ devsim.add_db_entry(material="Silicon",parameter="N_c",value=N_c, unit="/cm^3", description="effective density of states in conduction band")
+ devsim.add_db_entry(material="Silicon",parameter="N_v",value=N_v, unit="/cm^3", description="effective density of states in valence band")
+ E_g=1.12*1.6e-19
+ N_i=pow(N_c*N_v,0.5)*math.exp(-E_g/(2*k*T))
+ devsim.add_db_entry(material="Silicon", parameter="n_i", value=N_i, unit="/cm^3", description="Intrinsic Electron Concentration")
+ devsim.add_db_entry(material="SiliconCarbide", parameter="n_i", value=N_i, unit="/cm^3", description="Intrinsic Electron Concentration")
+ devsim.add_db_entry(material="gas", parameter="n_i", value="1e-9", unit="/cm^3", description="Intrinsic Electron Concentration")
+ devsim.add_db_entry(material="gas", parameter="Permittivity", value="1", unit="1", description="Permittivity")
+ devsim.add_db_entry(material="Silicon", parameter="n1", value=N_i, unit="/cm^3", description="n1")
+ devsim.add_db_entry(material="Silicon", parameter="p1", value=N_i, unit="/cm^3", description="p1")
+
+ if "parameter_alter" in MyDetector.device_dict:
+ for material in MyDetector.device_dict["parameter_alter"]:
+ print (material)
+ for parameter in MyDetector.device_dict["parameter_alter"][material]:
+ print (parameter)
+ devsim.add_db_entry(material=material,
+ parameter=parameter['name'],
+ value=parameter['value'],
+ unit=parameter['unit'],
+ description=parameter['name'])
+ if "parameter" in MyDetector.device_dict:
+ devsim.add_db_entry(material=MyDetector.device_dict['parameter']['material'],parameter=MyDetector.device_dict['parameter']['name'],value=MyDetector.device_dict['parameter']['value'],unit=MyDetector.device_dict['parameter']['unit'],description=MyDetector.device_dict['parameter']['description'])
+ if "U_const" in MyDetector.device_dict:
+ U_const=MyDetector.device_dict["U_const"]
+ model_create.CreateNodeModel(device,region,"U_const",U_const)
+ else:
+ model_create.CreateNodeModel(device,region,"U_const",0)
+ if "irradiation" in MyDetector.device_dict:
+ irradiation_label=MyDetector.device_dict['irradiation']['irradiation_label']
+ irradiation_flux=MyDetector.device_dict['irradiation']['irradiation_flux']
+ else:
+ irradiation_label=None
+ irradiation_flux=None
+ if 'avalanche_model' in MyDetector.device_dict:
+ impact_label=MyDetector.device_dict['avalanche_model']
+ else:
+ impact_label=None
+
+
+ circuit_contacts=[]
+ if is_wf == True:
+ if MyDetector.device_dict.get("mesh", {}).get("2D_mesh", {}).get("ac_contact"):
+ print("=========RASER info===================\nACLGAD is simulating\n=============info====================")
+ for read_out_lsit in MyDetector.device_dict["mesh"]["2D_mesh"]["ac_contact"]:
+ circuit_contacts.append( read_out_lsit["name"])
+ for i,c in enumerate(circuit_contacts):
+ devsim.circuit_element(name="V{}".format(i), n1=physics_drift_diffusion.GetContactBiasName(c), n2=0,
+ value=0.0, acreal=paras['acreal'], acimag=paras['acimag'])
+ else:
+ print("===============RASER info===================\nNot AC detector\n===========info=============")
+ circuit_contacts.append( MyDetector.device_dict['bias']['electrode'])
+
+ else:
+ circuit_contacts=MyDetector.device_dict['bias']['electrode']
+ devsim.circuit_element(name="V1", n1=physics_drift_diffusion.GetContactBiasName(circuit_contacts), n2=0,
+ value=0.0, acreal=paras['acreal'], acimag=paras['acimag'])
+ T1 = time.time()
+ print("================RASER info============\nWelcome to RASER TCAD PART, mesh load successfully\n=============info===============")
+ devsim.set_parameter(name="debug_level", value="info")
+ devsim.set_parameter(name = "extended_solver", value=True)
+ devsim.set_parameter(name = "extended_model", value=True)
+ devsim.set_parameter(name = "extended_equation", value=True)
+
+ if is_cv ==True:
+ solve_model = "cv"
+ elif is_noise == True:
+ solve_model = "noise"
+ elif is_wf ==True:
+ solve_model = "wf"
+ elif is_step ==True:
+ solve_model = "step"
+ else :
+ solve_model = None
+ path = output(__file__, device)
+ loop=loop_section.loop_section(paras=paras,device=device,region=region,solve_model=solve_model,irradiation=irradiation)
+ def worker_function(queue, lock, circuit_contacts, v_current, area_factor, path, device, region, solve_model, irradiation, is_wf):
+ try:
+ print(f"è¿è¡Œ loop_solver,å‚æ•°:{circuit_contacts}, {v_current}, {area_factor}")
+ loop.loop_solver(circuit_contact=circuit_contacts, v_current=v_current, area_factor=area_factor, path=path, device=device, region=region)
+ result_message = "Execution completed successfully"
+
+ except Exception as e:
+ result_message = f"Error: {e}"
+ with lock:
+ queue.put(result_message)
+
+
+
+ if is_wf == True:
+ v_current=1
+ print("=======RASER info========\nBegin simulation WeightingField\n======================")
+ for contact in circuit_contacts:
+ folder_path = os.path.join(path,"/weghtingfield/{}_{}_device".format(contact,v_current))
+ if not os.path.exists(folder_path):
+ os.makedirs(folder_path)
+
+ paras["milestone_step"] == 1
+ paras.update({"milestone_step":paras["milestone_step"]})
+ devsim.circuit_element(name="V1", n1=physics_drift_diffusion.GetContactBiasName(contact), n2=0,
+ value=0.0, acreal=paras['acreal'], acimag=paras['acimag'])
+ loop.initial_solver(contact=contact,set_contact_type=None,irradiation_label=irradiation_label,irradiation_flux=irradiation_flux,impact_label=impact_label)
+
+ loop.loop_solver(circuit_contact=contact,v_current=v_current,area_factor=paras["area_factor"],path=path,device=device,region=region)
+ if(paras['milestone_mode']==True and v_current%paras['milestone_step']==0.0):
+ save_milestone.save_milestone(device=device, region=region, v=v_current, path=path,dimension=default_dimension,contact=contact,is_wf=is_wf)
+ devsim.write_devices(file=os.path.join(folder_path,"weightingfield.dat"), type="tecplot")
+ elif is_wf == False:
+ v_current = 0
+ loop.initial_solver(contact=circuit_contacts,set_contact_type=None,irradiation_label=irradiation_label,irradiation_flux=irradiation_flux,impact_label=impact_label)
+ v_current = 0
+ v_goal =MyDetector.device_dict['bias']['voltage']
+ if v_goal > 0:
+ voltage_step = paras['voltage_step']
+ else:
+ voltage_step = -1 * paras['voltage_step']
+ if is_step == False:
+ while abs(v_current) <= abs(v_goal):
+ loop.loop_solver(circuit_contact=circuit_contacts,v_current=v_current,area_factor=paras["area_factor"],path=path,device=device,region=region)
+ if(paras['milestone_mode']==True and v_current%paras['milestone_step']==0.0) or abs(v_current) == abs(v_goal) :
+ save_milestone.save_milestone(device=device, region=region, v=v_current, path=path,dimension=default_dimension,contact=circuit_contacts,is_wf=is_wf)
+ devsim.write_devices(file=os.path.join(path,"Ele_Characterization/{}".format(v_current)), type="tecplot")
+ v_current += voltage_step
+
+ if is_step:
+ lock = multiprocessing.Lock()
+ queue = multiprocessing.Queue()
+
+ while abs(v_current) <= abs(v_goal):
+ print("voltage_step:", voltage_step)
+ print("v_current:", v_current)
+ print("============================check======================")
+ print(queue)
+ print("============================queue======================")
+ time.sleep(5)
+ p = multiprocessing.Process(target=worker_function, args=(queue, lock, circuit_contacts, v_current, paras['area_factor'], path, device, region, solve_model, irradiation, is_wf))
+ p.start()
+ p.join()
+ while not queue.empty():
+ output_info = queue.get()
+ print("队列输出:", output_info) # 确认输出内容
+ if output_info is None:
+ print("è¦å‘Š: worker_function 返回了 None,å¯èƒ½å‘生了错误!")
+ # exit(1)
+ if (paras['milestone_mode'] and v_current % paras['milestone_step'] == 0.0) or abs(v_current) == abs(v_goal):
+ save_milestone.save_milestone(device=device, region=region, v=v_current, path=path, dimension=default_dimension, contact=circuit_contacts, is_wf=is_wf)
+ devsim.write_devices(file=os.path.join(path, "Ele_Characterization_{}.dd".format(v_current)), type="tecplot")
+
+ v_current += voltage_step
+
+ if is_wf != True:
+ if is_step ==False:
+ draw_iv(device, V=loop.get_voltage_values(), I=loop.get_current_values(),path=path)
+ if is_cv == True:
+ draw_cv(device, V=loop.get_voltage_values(), C=loop.get_cap_values(),path=path)
+ if is_noise == True:
+ draw_noise(device, V=loop.get_voltage_values(), noise=loop.get_noise_values(),path=path)
+ elif is_step == True:
+ draw_iv(device, V=V, I=Current,path=path)
+ if is_cv == True:
+ draw_cv(device, V=V, C=c,path=path)
+ if is_noise == True:
+ draw_noise(device, V=V, noise=noise,path=path)
+ T2 =time.time()
+ print("=========RASER info===========\nSimulation finish ,total used time: {}s !^w^!\n================".format(T2-T1))
\ No newline at end of file
diff --git a/field/test4hsic.py b/field/test4hsic.py
deleted file mode 100755
index 0f3399c0fb027de5d5afa64c40d7b6e52fe4bb08..0000000000000000000000000000000000000000
--- a/field/test4hsic.py
+++ /dev/null
@@ -1,158 +0,0 @@
-import ROOT
-import pickle
-import numpy as np
-from scipy.interpolate import Rbf
-from scipy.interpolate import griddata
-import math
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-
-# TODO: rewrite this after 3D field extraction verified
-
-def main(simname):
- if '3d' in simname:
- with open('./output/field/{}/x_500.0.pkl'.format(simname), 'rb') as file:
- x = pickle.load(file)
- with open('./output/field/{}/y_500.0.pkl'.format(simname), 'rb') as file:
- y = pickle.load(file)
- with open('./output/field/{}/z_500.0.pkl'.format(simname), 'rb') as file:
- z = pickle.load(file)
- with open('./output/field/{}/potential_500.0.pkl'.format(simname), 'rb') as file:
- potential = pickle.load(file)
-
- # 定义æ’å€¼ç½‘æ ¼
- nPoints = 100
- x_min, x_max = np.min(x), np.max(x)
- y_min, y_max = np.min(y), np.max(y)
- z_min, z_max = np.min(z), np.max(z)
-
- # 进行æ’值
- xi = np.linspace(x_min, x_max, nPoints)
- yi = np.linspace(y_min, y_max, nPoints)
- zi = np.linspace(z_min, z_max, nPoints)
- xi, yi, zi = np.meshgrid(xi, yi, zi)
- pi = griddata((x, y, z), potential, (xi, yi, zi), method='linear')
-
- # 计算梯度
- dx, dy, dz = np.gradient(pi, xi[0,0,0], yi[0,0,0], zi[0,0,0])
-
- # 计算总场强
- total_field = np.sqrt(dx**2 + dy**2 + dz**2)
- # 绘制电势的xzå¹³é¢æˆªé¢å›¾
- fig_potential_xz = plt.figure()
- ax_potential_xz = fig_potential_xz.add_subplot(111)
- ax_potential_xz.contourf(xi[:, nPoints//2, :], zi[:, nPoints//2, :], pi[:, nPoints//2, :], cmap='viridis')
- ax_potential_xz.set_xlabel('x')
- ax_potential_xz.set_ylabel('z')
- ax_potential_xz.set_title('Potential XZ Plane')
- plt.savefig('./output/field/{}/Potential_XZ_Plane.png'.format(simname))
-
- # 绘制电势的yzå¹³é¢æˆªé¢å›¾
- fig_potential_yz = plt.figure()
- ax_potential_yz = fig_potential_yz.add_subplot(111)
- ax_potential_yz.contourf(yi[:, nPoints//2, :], zi[:, nPoints//2, :], pi[:, nPoints//2, :], cmap='viridis')
- ax_potential_yz.set_xlabel('y')
- ax_potential_yz.set_ylabel('z')
- ax_potential_yz.set_title('Potential YZ Plane')
- plt.savefig('./output/field/{}/Potential_YZ_Plane.png'.format(simname))
- # 绘制电势的xyå¹³é¢æˆªé¢å›¾
- fig_xy_planes = plt.figure(figsize=(10, 10))
- for i in range(0, nPoints, 5):
- ax_xy_plane = fig_xy_planes.add_subplot(5, 5, i//5+1)
- ax_xy_plane.imshow(pi[:,:,i], extent=[x_min, x_max, y_min, y_max], cmap='viridis', origin='lower')
- ax_xy_plane.set_xlabel('x')
- ax_xy_plane.set_ylabel('y')
- ax_xy_plane.set_title('XY Plane')
- plt.tight_layout()
- plt.savefig('./output/{}/Potential_XY_Planes.png'.format(simname))
- #Plot total field XY plane
- fig_total_field_xy = plt.figure()
- ax_total_field_xy = fig_total_field_xy.add_subplot(111)
- ax_total_field_xy.contourf(xi[:, :, nPoints // 2], yi[:, :, nPoints // 2], total_field[:, :, nPoints // 2], cmap='viridis')
- ax_total_field_xy.set_xlabel('x')
- ax_total_field_xy.set_ylabel('y')
- ax_total_field_xy.set_title('Total Field XY Plane')
- plt.savefig('./output/field/{}/Total_Field_XY_Plane.png'.format(simname))
-
- #Plot total field YZ plane
- fig_total_field_yz = plt.figure()
- ax_total_field_yz = fig_total_field_yz.add_subplot(111)
- ax_total_field_yz.contourf(yi[:, nPoints // 2, :], zi[:, nPoints // 2, :], total_field[:, nPoints // 2, :], cmap='viridis')
- ax_total_field_yz.set_xlabel('y')
- ax_total_field_yz.set_ylabel('z')
- ax_total_field_yz.set_title('Total Field YZ Plane')
- plt.savefig('./output/field/{}/Total_Field_YZ_Plane.png'.format(simname))
-
- #Plot 20 dianchang XY planes
- fig_xy = plt.figure(figsize=(15, 15))
- nPlanes = 20
- for i in range(nPlanes):
- index = int(nPoints*i/(nPlanes-1)) % nPoints # Wrap index within the range of nPoints
- ax_xy = fig_xy.add_subplot(4, 5, i+1)
- ax_xy.imshow(pi[:, :, index], extent=[x_min, x_max, y_min, y_max], cmap='viridis', origin='lower')
- ax_xy.set_xlabel('x')
- ax_xy.set_ylabel('y')
- ax_xy.set_title('XY Plane {}'.format(i+1))
- plt.tight_layout()
- plt.savefig('./output/field/{}/Total_Field_XY_Planes.png'.format(simname))
-
-
-
- else:
- # 从pickle文件ä¸åŠ è½½xã€y和电势数æ®
- with open('./output/field/{}/x_500.0.pkl'.format(simname), 'rb') as file:
- x = pickle.load(file)
- with open('./output/field/{}/y_500.0.pkl'.format(simname), 'rb') as file:
- y = pickle.load(file)
- with open('./output/field/{}/potential_500.0.pkl'.format(simname), 'rb') as file:
- potential = pickle.load(file)
-
- # 定义æ’å€¼ç½‘æ ¼
- nPoints = 1000
- x_min, x_max = np.min(x), np.max(x)
- y_min, y_max = np.min(y), np.max(y)
-
- # 进行æ’值
- rbf = Rbf(x, y, potential, function='cubic')
-
- # 获å–æ’值结果
- xi = np.linspace(x_min, x_max, nPoints)
- yi = np.linspace(y_min, y_max, nPoints)
- xi, yi = np.meshgrid(xi, yi)
- zi = rbf(xi, yi)
-
- # 创建TH2F对象，用于å˜å‚¨æ’值结果
- hInterpolated = ROOT.TH2F('hInterpolated', 'Interpolated Potential', nPoints, x_min, x_max, nPoints, y_min, y_max)
-
- # å°†æ’值结果填充到TH2F对象ä¸
- for i in range(nPoints):
- for j in range(nPoints):
- hInterpolated.SetBinContent(i+1, j+1, zi[i, j])
-
- # ä¿å˜æ’值结果为ROOT对象
- output_file = ROOT.TFile('./output/{}/Interpolated_Potential_2D.root'.format(simname), 'recreate')
- hInterpolated.Write()
- output_file.Close()
-
- # 计算梯度
- dx = np.gradient(zi, xi[0,0], axis=0)
- dy = np.gradient(zi, yi[0,0], axis=1)
-
- # 计算总场强
- total_field = np.sqrt(dx**2 + dy**2)
-
- # 创建TH2F对象，用于å˜å‚¨æ€»åœºå¼º
- hTotalField = ROOT.TH2F('hTotalField', 'Total Field', nPoints, x_min, x_max, nPoints, y_min, y_max)
-
- # 将总场强填充到TH2F对象ä¸
- for i in range(nPoints):
- for j in range(nPoints):
- hTotalField.SetBinContent(i+1, j+1, total_field[i, j])
-
- # ä¿å˜æ€»åœºå¼ºä¸ºROOT对象
- output_file = ROOT.TFile('./output/field/{}/Total_Field_2D.root'.format(simname), 'recreate')
- hTotalField.Write()
- output_file.Close()
-
-if __name__ == "__main__":
- main("2dfield_4HSiC")
\ No newline at end of file