Dear all,
I’m trying to run a simulation of a GEM to quantify gain and transparency. I used the biconical gem example in Garfield/Examples/neBEM, added the gas volume to the script and run the for loop to just 5 events and it takes too long to run. I dont know if im using some element size or configuration differently in nebem. Also, i would like to ask if is there any way to save the resolved field on a separate file and then just import it, like when using ansys or elmer, to run faster. I pasted my code below:
#include <TApplication.h>
#include <iostream>
#include "Garfield/ComponentNeBem3d.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/MediumConductor.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/MediumPlastic.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/SolidHole.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/Sensor.hh"
using namespace Garfield;
int main(int argc, char* argv[]) {
TApplication app("app", &argc, argv);
MediumMagboltz gas("Ar", 70.0, "co2", 30);
gas.SetTemperature(300);
gas.SetPressure(760);
gas.EnableDrift();
MediumConductor Cu;
MediumPlastic Kp;
Kp.SetDielectricConstant(4.0);
// Geometry.
GeometrySimple geo;
double kptx = 0.0;
double kpty = 0.0;
double kptz = 0.0;
double lenLX = 0.0120;
double lenLY = 0.0120;
double kptLZ = 0.0050;
double maxdia = 0.0070;
double mindia = 0.0050;
double lwcprx = 0.0;
double lwcpry = 0.0;
double lwcprLZ = 0.0005;
double lwcprz = kptz - (kptLZ / 2.0) - (lwcprLZ / 2.0);
double indgap = 0.05;
double anodex = 0.0;
double anodey = 0.0;
double anodeLZ = 0.0;
double anodez = lwcprz - (lwcprLZ / 2.0) - indgap - (anodeLZ / 2.0);
double upcprx = 0.0;
double upcpry = 0.0;
double upcprLZ = 0.0005;
double upcprz = kptz + (kptLZ / 2.0) + (upcprLZ / 2.0);
double drftgap = 0.05;
double drftx = 0.0;
double drfty = 0.0;
double drftLZ = 0.0;
double drftz = upcprz + (upcprLZ / 2.0) + drftgap + (drftLZ / 2.0);
double drftV = -750;
double upcprV = -250;
double lwcprV = 250;
double anodeV = 750;
std::cout << "kptz: " << kptz << std::endl;
std::cout << "lwcprz: " << lwcprz << std::endl;
std::cout << "anodez: " << anodez << std::endl;
std::cout << "upcprz: " << upcprz << std::endl;
std::cout << "drftz: " << drftz << std::endl;
SolidBox cathode(drftx, drfty, drftz, lenLX / 2.0, lenLY / 2.0, drftLZ / 2.0);
cathode.SetBoundaryPotential(drftV);
geo.AddSolid(&cathode, &Cu);
SolidHole topCu(upcprx, upcpry, upcprz, maxdia / 2.0, maxdia / 2.0,
lenLX / 2.0, lenLY / 2.0, upcprLZ / 2.0);
topCu.SetSectors(3);
topCu.SetBoundaryPotential(upcprV);
geo.AddSolid(&topCu, &Cu);
SolidHole topKapton(kptx, kpty, (kptz + (kptLZ / 4.0)), maxdia / 2.0,
mindia / 2.0, lenLX / 2.0, lenLY / 2.0, kptLZ / 4.0);
topKapton.SetSectors(3);
topKapton.SetBoundaryDielectric();
geo.AddSolid(&topKapton, &Kp);
SolidHole btmKapton(kptx, kpty, (kptz - (kptLZ / 4.0)), mindia / 2.0,
maxdia / 2.0, lenLX / 2.0, lenLY / 2.0, kptLZ / 4.0);
btmKapton.SetSectors(3);
btmKapton.SetBoundaryDielectric();
geo.AddSolid(&btmKapton, &Kp);
SolidHole btmCu(lwcprx, lwcpry, lwcprz, maxdia / 2.0, maxdia / 2.0,
lenLX / 2.0, lenLY / 2.0, lwcprLZ / 2.0);
btmCu.SetSectors(3);
btmCu.SetBoundaryPotential(lwcprV);
geo.AddSolid(&btmCu, &Cu);
SolidBox anode(anodex, anodey, anodez, lenLX / 2.0, lenLY / 2.0,
anodeLZ / 2.0);
anode.SetBoundaryPotential(anodeV);
geo.AddSolid(&anode, &Cu);
geo.SetMedium(&gas);
double tgtElSize = 10.e-4;
int minEl = 3, maxEl = 5;
int xcopy = 10, ycopy = 10, zcopy = 0;
ComponentNeBem3d nebem;
nebem.SetGeometry(&geo);
nebem.SetTargetElementSize(tgtElSize);
nebem.SetMinMaxNumberOfElements(minEl, maxEl);
nebem.SetPeriodicityX(lenLX);
nebem.SetPeriodicityY(lenLY);
nebem.SetPeriodicCopies(xcopy, ycopy, zcopy);
nebem.UseLUInversion();
// nebem.EnableDebugging();
nebem.Initialise();
/*
ViewGeometry geomView2dClose(&geo);
geomView2dClose.SetArea(-lenLX, -lenLY, upcprz, lenLX, lenLY, lwcprz);
geomView2dClose.SetPlane(0, 1, 0, 0, 0, 0.0);
geomView2dClose.Plot2d();
ViewField EfieldView(&nebem);
EfieldView.SetArea(-lenLX / 2.0, -lenLY / 2.0, lwcprz - 10.0e-4, lenLX / 2.0,
lenLY / 2.0, upcprz + 10.0e-4);
EfieldView.SetPlane(0, 1, 0, 0, 0, 0);
EfieldView.PlotContour("e");
std::cout << "plot feito?" << std::endl;
*/
Sensor* sensor = new Sensor();
sensor->AddComponent(&nebem);
sensor->SetArea(-3*lenLX, -3*lenLY, -0.02, 3*lenLX, 3*lenLY, 0.02);
AvalancheMicroscopic* aval = new AvalancheMicroscopic();
aval->SetSensor(sensor);
double acumulador = 0;
double x0,y0,z0,t0,e0;
double x1,y1,z1,t1,e1;
int status;
const int nAval = 5;
for (int i=0; i < nAval; i++){
double x = 0;//rng.Uniform(-pitch, pitch);
double y = 0;//rng.Uniform(-pitch, pitch);
double z = 180e-4;//copperthickness + 0.5*kaptonthickness + 1*pitch; // 140 um acima do furo do gem
double t = 0;
double e = 0.0;
std::cout << i+1 << "/" << nAval << " Elétrons soltos" << std::endl;
aval->AvalancheElectron(x,y,z,t,e,0,0,0);
int ne, ni;
aval->GetAvalancheSize(ne,ni);
std::cout << "Avalanche Size :" << ne << std::endl;
aval->GetElectronEndpoint(0,x0,y0,z0,t0,e0,x1,y1,z1,t1,e1,status);
std::cout<< "posição inicial: x = " << x0 << ", y = " << y0 << ", z = " << z0 << " at " << t0 << std::endl;
std::cout << "posição final: x = " << x1 << ", y = " << y1 << ", z = " << z1 << " at " << t1 << std::endl;
std::cout << "\n" << std::endl;
if (status == -1 && z1 < 1e-4){
std::cout << "electron passed through" << std::endl;
acumulador++;
}
if( status == -5 || status == -7 || status == -9){
std::cout << "electron got lost "<< std::endl;
std::cout << "\n" << std::endl;
}
}
std::cout << " transparency " << acumulador/nAval<< std::endl;
return 0;
}