diff --git a/src/SimNDT/engine/efit2d.py b/src/SimNDT/engine/efit2d.py new file mode 100644 index 0000000..54b3db3 --- /dev/null +++ b/src/SimNDT/engine/efit2d.py @@ -0,0 +1,575 @@ +from SimNDT.engine.engineBase import EngineBase + +try: + import SimNDT.engine.efit2dcython as EFIT2DCython +except: + print("error in importation for Serial EFIT2D") + +from SimNDT.core.material import Material +from SimNDT.core.constants import BC + +import copy +import numpy as np + +global ErrorImportCL + +try: + import pyopencl as cl + + ErrorImportCL = False +except ImportError: + ErrorImportCL = True + +c11 = 4350 * 4350 * 7750 +c12 = 4350 * 4350 * 7750 - 2260 * 2260 * 7750 +c22 = c11 +c44 = 2260 * 2260 * 7750 +pzt = Material("pzt", 7750.0, c11, c12, c22, c44) + + +class EFIT2D(EngineBase): + def __init__(self, simPack, Platform): + EngineBase.__init__(self, simPack, Platform) + + self.max_value = 0.0 + self._Receiver = False + + def setup_CL(self): + if self.Platform == "OpenCL": + self.initFieldsCL() + + def materialSetup(self): + + Materials = copy.deepcopy(self.simPack.Materials) + MRI, NRI = np.shape(self.simPack.Simulation.Im) + + # Condicion de vacio si uno de los materiales es aire + for n in range(len(Materials)): + if Materials[n].Rho < 2.0: + Materials[n].Rho = 10e23 + Materials[n].C11 = 1e-20 + Materials[n].C12 = 1e-20 + Materials[n].C22 = 1e-20 + Materials[n].C44 = 1e-20 + Materials[n].Eta_v = 1e-20 + Materials[n].Eta_s = 1e-20 + # verificar que no sea cero exactamente + if Materials[n].C44 == 0: + Materials[n].C44 = 1e-30 + if Materials[n].Eta_s == 0: + Materials[n].Eta_s = 1e-30 + + self.Rho = np.zeros((MRI, NRI), dtype=np.float32) + self.C11 = np.zeros((MRI, NRI), dtype=np.float32) + self.C12 = np.zeros((MRI, NRI), dtype=np.float32) + self.C22 = np.zeros((MRI, NRI), dtype=np.float32) + self.C44 = np.ones((MRI, NRI), dtype=np.float32) * 1e-30 + self.Eta_v = np.zeros((MRI, NRI), dtype=np.float32) + self.Eta_s = np.zeros((MRI, NRI), dtype=np.float32) + + for n in range(len(Materials)): + ind = np.nonzero(self.simPack.Simulation.Im == Materials[n].Label) + self.Rho[ind] = Materials[n].Rho + self.C11[ind] = Materials[n].C11 + self.C12[ind] = Materials[n].C12 + self.C22[ind] = Materials[n].C22 + self.C44[ind] = Materials[n].C44 + self.Eta_v[ind] = Materials[n].Eta_v + self.Eta_s[ind] = Materials[n].Eta_s + + # Find the base materials + for i, item in enumerate(Materials): + if item.Label == 0: + base_material = i + + # Set material in boundaries + ind = np.nonzero(self.simPack.Simulation.Im == 255.0) + self.Rho[ind] = Materials[base_material].Rho + self.C11[ind] = Materials[base_material].C11 + self.C12[ind] = Materials[base_material].C12 + self.C22[ind] = Materials[base_material].C22 + self.C44[ind] = Materials[base_material].C44 + self.Eta_v[ind] = Materials[base_material].Eta_v + self.Eta_s[ind] = Materials[base_material].Eta_s + + def initFields(self): + + MRI, NRI = np.shape(self.simPack.Simulation.Im) + Signal = self.simPack.Signal + t = self.simPack.Simulation.t + self.source = Signal.generate(t) + + self.Vx = np.zeros((MRI, NRI), dtype=np.float32) + self.Vy = np.zeros((MRI, NRI), dtype=np.float32) + self.Txx = np.zeros((MRI, NRI), dtype=np.float32) + self.Txy = np.zeros((MRI, NRI), dtype=np.float32) + self.Tyy = np.zeros((MRI, NRI), dtype=np.float32) + self.DVx = np.zeros((MRI, NRI), dtype=np.float32) + self.DVy = np.zeros((MRI, NRI), dtype=np.float32) + self.ABS = np.ones((MRI, NRI), dtype=np.float32) + self.SV = np.zeros((MRI, NRI), dtype=np.float32) + + def staggeredProp(self): + + MRI, NRI = np.shape(self.simPack.Simulation.Im) + + BXtemp = np.zeros((MRI, NRI), dtype=np.float32) + BYtemp = np.zeros((MRI, NRI), dtype=np.float32) + self.BX = np.zeros((MRI, NRI), dtype=np.float32) + self.BY = np.zeros((MRI, NRI), dtype=np.float32) + self.C44_Eff = np.zeros((MRI, NRI), dtype=np.float32) + + BXtemp[:, :] = 1.0 / self.Rho[:, :] + BYtemp[:, :] = 1.0 / self.Rho[:, :] + + self.BX[:-2, :] = 0.5 * (BXtemp[1:-1, :] + BXtemp[:-2, :]) + self.BX[-2, :] = np.copy(BXtemp[-2, :]) + + self.BY[:, :-2] = 0.5 * (BYtemp[:, 1:-1] + BYtemp[:, :-2]) + self.BY[:, -2] = np.copy(BYtemp[:, -2]) + + self.C44_Eff[:-2, :-2] = 4. / ( + (1. / self.C44[:-2, :-2]) + (1. / self.C44[1:-1, :-2]) + (1. / self.C44[:-2, 1:-1]) + ( + 1. / self.C44[1:-1, 1:-1])) + + self.Eta_vs = np.zeros((MRI, NRI), dtype=np.float32) + self.Eta_ss = np.zeros((MRI, NRI), dtype=np.float32) + self.Eta_vs = self.Eta_v + 2 * self.Eta_s + self.Eta_ss[:-2, :-2] = 4. / ( + (1. / self.Eta_s[:-2, :-2]) + (1. / self.Eta_s[1:-1, :-2]) + (1. / self.Eta_s[:-2, 1:-1]) + ( + 1. / self.Eta_s[1:-1, 1:-1])) + + def applyBoundaries(self): + Materials = copy.deepcopy(self.simPack.Materials) + MRI, NRI = np.shape(self.simPack.Simulation.Im) + + APARA = 0.015 + + Top = int(self.simPack.Simulation.TapGrid[0]) + Bottom = int(self.simPack.Simulation.TapGrid[1]) + Left = int(self.simPack.Simulation.TapGrid[2]) + Right = int(self.simPack.Simulation.TapGrid[3]) + + _Top = Top - np.arange(0, Top) + _Bottom = np.arange(MRI - Bottom + 1, MRI) - MRI + Bottom - 1 + _Left = Left - np.arange(0, Left) + _Right = np.arange(NRI - Right + 1, NRI) - NRI + Right - 1 + + Boundaries = self.simPack.Boundary + for boundary in Boundaries: + if boundary.Name == "Top" and boundary.BC == BC.AirLayer: + self.BX[0:Top, :] = 0.0 + self.BY[0:Top, :] = 0.0 + self.C11[0:Top, :] = 0.0 + self.C12[0:Top, :] = 0.0 + self.C22[0:Top, :] = 0.0 + self.C44_Eff[0:Top, :] = 0.0 + self.Eta_vs[0:Top, :] = 0.0 + self.Eta_ss[0:Top, :] = 0.0 + + if boundary.Name == "Bottom" and boundary.BC == BC.AirLayer: + self.BX[MRI - Bottom:, :] = 0.0 + self.BY[MRI - Bottom:, :] = 0.0 + self.C11[MRI - Bottom:, :] = 0.0 + self.C12[MRI - Bottom:, :] = 0.0 + self.C22[MRI - Bottom:, :] = 0.0 + self.C44_Eff[MRI - Bottom:, :] = 0.0 + self.Eta_vs[MRI - Bottom:, :] = 0.0 + self.Eta_ss[MRI - Bottom:, :] = 0.0 + + if boundary.Name == "Left" and boundary.BC == BC.AirLayer: + self.BX[:, 0:Left] = 0.0 + self.BY[:, 0:Left] = 0.0 + self.C11[:, 0:Left] = 0.0 + self.C12[:, 0:Left] = 0.0 + self.C22[:, 0:Left] = 0.0 + self.C44_Eff[:, 0:Left] = 0.0 + self.Eta_vs[:, 0:Left] = 0.0 + self.Eta_ss[:, 0:Left] = 0.0 + if boundary.Name == "Right" and boundary.BC == BC.AirLayer: + self.BX[:, NRI - Right:] = 0.0 + self.BY[:, NRI - Right:] = 0.0 + self.C11[:, NRI - Right:] = 0.0 + self.C12[:, NRI - Right:] = 0.0 + self.C22[:, NRI - Right:] = 0.0 + self.C44_Eff[:, NRI - Right:] = 0.0 + self.Eta_vs[:, NRI - Right:] = 0.0 + self.Eta_ss[:, NRI - Right:] = 0.0 + + for j in range(0, NRI): + self.ABS[0:Top, j] = np.exp(- ((APARA * _Top) ** 2)) + self.ABS[MRI - Bottom + 1:, j] = np.exp(- ((APARA * _Bottom) ** 2)) + for i in range(0, MRI): + self.ABS[i, 0:Left] = np.exp(- ((APARA * _Left) ** 2)) + self.ABS[i, NRI - Right + 1:] = np.exp(- ((APARA * _Right) ** 2)) + + def sourceSetup(self): + + longitudinal = self.simPack.Source.Longitudinal + shear = self.simPack.Source.Shear + + Pressure = self.simPack.Source.Pressure + Displacement = self.simPack.Source.Displacement + + if Pressure and not Displacement: + self.typeSource = 0 + elif not Pressure and Displacement: + self.typeSource = 1 + elif Pressure and Displacement: + self.typeSource = 2 + else: + self.typeSource = 0 + + if longitudinal and not shear: + self.typeWave = np.int32(0) + elif not longitudinal and shear: + self.typeWave = np.int32(1) + elif longitudinal and shear: + self.typeWave = np.int32(2) + else: + self.typeWave = np.int32(0) + + self.Amplitude = self.simPack.Signal.Amplitude + XL = self.simPack.Inspection.XL + YL = self.simPack.Inspection.YL + NX = np.size(XL, 0) + + IsPZT = self.simPack.Transducers[0].PZT + + size1 = int((self.simPack.Simulation.TapGrid[0])) + size2 = int((self.simPack.Simulation.TapGrid[1])) + + for IT in range(-size1, 0): + for m in range(0, NX): + xl = int(XL[m, 0]) + yl = int(YL[m, 0]) + self.BX[xl + IT, yl] = 1.0 / pzt.Rho if IsPZT else 0 + self.BY[xl + IT, yl] = 1.0 / pzt.Rho if IsPZT else 0 + self.C11[xl + IT, yl] = pzt.C11 if IsPZT else 0 + self.C12[xl + IT, yl] = pzt.C12 if IsPZT else 0 + self.C22[xl + IT, yl] = pzt.C22 if IsPZT else 0 + self.C44[xl + IT, yl] = pzt.C44 if IsPZT else 0 + + if self.simPack.Inspection.Name == "Transmission": + for IT in range(1, size2): + for m in range(0, NX): + xl = int(XL[m, 1]) + yl = int(YL[m, 1]) + self.BX[xl + IT, yl] = 1.0 / pzt.Rho if IsPZT else 0 + self.BY[xl + IT, yl] = 1.0 / pzt.Rho if IsPZT else 0 + self.C11[xl + IT, yl] = pzt.C11 if IsPZT else 0 + self.C12[xl + IT, yl] = pzt.C12 if IsPZT else 0 + self.C22[xl + IT, yl] = pzt.C22 if IsPZT else 0 + self.C44[xl + IT, yl] = pzt.C44 if IsPZT else 0 + + def simSetup(self): + self.MRI, self.NRI = np.shape(self.simPack.Simulation.Im) + XL = self.simPack.Inspection.XL + YL = self.simPack.Inspection.YL + + self.NX = np.size(XL, 0) + self.XL = np.copy(np.int32(XL[:, 0])) + self.YL = np.copy(np.int32(YL[:, 0])) + + if self.simPack.Transducers[0].Window and not self.simPack.Transducers[0].PointSource: + self.win = np.float32(1.0 - np.cos(2 * np.pi * np.arange(0, self.NX) / (self.NX + 1))) + else: + self.win = np.ones((self.NX,), dtype=np.float32) + + def initFieldsCL(self): + + self.Txx_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Txx) + self.Tyy_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Tyy) + self.Txy_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Txy) + self.Vx_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Vx) + self.Vy_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.Vy) + + self.DVx_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.DVx) + self.DVy_buf = cl.Buffer(self.ctx, self.mf.READ_WRITE | self.mf.COPY_HOST_PTR, hostbuf=self.DVy) + + self.ABS_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.ABS) + self.BX_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.BX) + self.BY_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.BY) + self.C11_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.C11) + self.C12_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.C12) + self.C44_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.C44_Eff) + + self.ETA_vs_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.Eta_vs) + self.ETA_s_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.Eta_s) + self.ETA_ss_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.Eta_ss) + + self.XL_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.XL) + self.YL_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.YL) + + self.WIN_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.win) + + if not self.simPack.Inspection.Name == "Tomography": + self._Receiver = True + XL = self.simPack.Inspection.XL + YL = self.simPack.Inspection.YL + self.receiver_buf = cl.Buffer(self.ctx, self.mf.WRITE_ONLY | self.mf.COPY_HOST_PTR, + hostbuf=self.receiver_signals) + self.XXL = np.copy(np.int32(XL[:, 1])) + self.YYL = np.copy(np.int32(YL[:, 1])) + self.XXL_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.XXL) + self.YYL_buf = cl.Buffer(self.ctx, self.mf.READ_ONLY | self.mf.COPY_HOST_PTR, hostbuf=self.YYL) + + self.program = cl.Program(self.ctx, self.kernel()).build() + + def receiverSetup(self): + TimeSteps = int(self.simPack.Simulation.TimeSteps) + N_IR = np.size(self.simPack.Inspection.IR, 1) + self.receiver_signals = np.zeros((TimeSteps, N_IR - 1), dtype=np.float32) + + # @blab+ + def clEnqueue(self, Var, Var_buf): + # external access to Txx .. for snapshoting needs enqueueing of CL buffers + cl.enqueue_copy(self.queue, Var, Var_buf) + + def receivers(self, Var, Var_buf): + if self._Receiver: + self.program.Receiver_EFIT2D(self.queue, (self.NX,), None, self.Txx_buf, self.receiver_buf, + np.int32(self.n), + self.XXL_buf, self.YYL_buf).wait() + else: + cl.enqueue_copy(self.queue, Var, Var_buf) + self.receiver_signals[self.n, :] = self.simPack.Inspection.getReceivers(Var) + + def saveOutput(self): + if self._Receiver: + cl.enqueue_copy(self.queue, self.receiver_signals, self.receiver_buf).wait() + + def receiversSerial(self, Var): + self.receiver_signals[self.n, :] = self.simPack.Inspection.getReceivers(Var) + + def run(self): + + y = np.float32(self.source[self.n]) + self.program.Velocity_EFIT2D_Voigt(self.queue, (self.NRI, self.MRI,), None, + self.Txx_buf, self.Txy_buf, self.Tyy_buf, + self.Vx_buf, self.Vy_buf, self.DVx_buf, self.DVy_buf, + self.BX_buf, self.BY_buf, self.ABS_buf, self.ddx).wait() + + if self.typeSource == 1 or self.typeSource == 2: + self.program.Source_Vel_EFIT2D(self.queue, (self.NX,), None, + self.Vx_buf, + self.Vy_buf, + self.XL_buf, + self.YL_buf, + y, self.typeWave, self.WIN_buf).wait() + + self.program.Stress_EFIT2D_Voigt(self.queue, (self.NRI, self.MRI,), None, + self.Txx_buf, self.Txy_buf, self.Tyy_buf, + self.Vx_buf, self.Vy_buf, self.DVx_buf, self.DVy_buf, + self.C11_buf, self.C12_buf, self.C44_buf, self.ETA_vs_buf, + self.ETA_s_buf, self.ETA_ss_buf, self.ABS_buf).wait() + + if self.typeSource == 0 or self.typeSource == 2: + self.program.Source_EFIT2D(self.queue, (self.NX,), None, + self.Txx_buf, + self.Tyy_buf, + self.Txy_buf, + self.XL_buf, + self.YL_buf, + y, self.typeWave, self.WIN_buf).wait() + + self.receivers(self.Txx, self.Txx_buf) + + def runGL(self): + cl.enqueue_copy(self.queue, self.Vx, self.Vx_buf) + cl.enqueue_copy(self.queue, self.Vy, self.Vy_buf) + self.SV = np.sqrt(self.Vx ** 2 + self.Vy ** 2) + value = np.max(np.abs(self.SV)) + self.max_value = self.max_value if self.max_value >= value else value + self.SV = 20. * np.log10(np.abs(self.SV) / self.max_value + 1e-40) + + def runGLSerial(self, step=50): + + vx = np.copy(self.Vx) + vy = np.copy(self.Vy) + + self.SV = np.sqrt(vx ** 2 + vy ** 2) + value = np.max(np.abs(self.SV)) + self.max_value = self.max_value if self.max_value >= value else value + self.SV = 20. * np.log10(np.abs(self.SV) / self.max_value + 1e-40) + + def runSerial(self): + + y = np.float32(self.source[self.n]) + + self.Vx, self.Vy, self.DVx, self.DVy = EFIT2DCython.velocityVoigt(self.Txx, self.Txy, self.Tyy, + self.Vx, self.Vy, self.DVx, self.DVy, + self.BX, self.BY, self.ABS, self.ddx, + self.dt) + + if self.typeSource == 1 or self.typeSource == 2: + self.Vx, self.Vy = EFIT2DCython.sourceVel(self.Vx, self.Vy, self.XL, self.YL, y, self.typeWave, self.win, + self.dtdxx) + + self.Txx, self.Tyy, self.Txy = EFIT2DCython.stressVoigt(self.Txx, self.Txy, self.Tyy, + self.Vx, self.Vy, self.DVx, self.DVy, + self.C11, self.C12, self.C44_Eff, + self.Eta_vs, self.Eta_s, self.Eta_ss, self.ABS, + self.dtx) + + if self.typeSource == 0 or self.typeSource == 2: + self.Txx, self.Tyy = EFIT2DCython.sourceStress(self.Txx, self.Tyy, self.XL, self.YL, y, self.typeWave, + self.win, self.dtdxx) + + self.receiversSerial(self.Txx) + + def kernel(self): + + macro = """ + #define MRI %s + #define NRI %s + #define ind(i, j) ( ( (i)*NRI) + (j) ) + #define dtx %gf + #define dtdxx %gf + #define Stencil 2 + #define NX %s + #define dt %gf + + + + """ % (str(self.MRI), str(self.NRI), self.dtx, self.dtdxx, + str(self.NX), self.dt) + + return macro + """ + + + __kernel void Receiver_EFIT2D(__global float *Buffer, __global float *receiver, const int t, + __global const int *XXL, __global const int *YYL){ + + + __private float _tmp = 0.0f; + + + for (int i=0; i0){ + dvx[ind(i,j)] = (BX[ind(i,j)]*ddx)*( Txx[ind(i+1,j)] - Txx[ind(i,j)] + Txy[ind(i,j)] - Txy[ind(i,j-1)] ); + vx[ind(i,j)] += dt*dvx[ind(i,j)]; + } + + if (i > 0 && j< NRI-1){ + dvy[ind(i,j)] = (BY[ind(i,j)]*ddx)*( Txy[ind(i,j)] - Txy[ind(i-1,j)] + Tyy[ind(i,j+1)] - Tyy[ind(i,j)] ); + vy[ind(i,j)] += dt*dvy[ind(i,j)]; + + } + + barrier(CLK_GLOBAL_MEM_FENCE); + // Apply absorbing boundary conditions + vx[ind(i,j)] *= ABS[ind(i,j)]; + vy[ind(i,j)] *= ABS[ind(i,j)]; + + + } + + + __kernel void Stress_EFIT2D_Voigt(__global float *Txx, __global float *Txy, __global float *Tyy, + __global float *vx, __global float *vy, + __global float *dvx, __global float *dvy, + __global const float *C11, __global const float *C12, __global const float *C44, + __global const float *ETA_VS, __global const float *ETA_S, __global const float *ETA_SS, + __global const float *ABS) { + + + int j = get_global_id(0); + int i = get_global_id(1); + + + if (i>0 && j>0 ){ + Txx[ind(i,j)] += ( ( C11[ind(i,j)]* dtx )*(vx[ind(i,j)] - vx[ind(i-1,j)]) + + ( C12[ind(i,j)]* dtx )*(vy[ind(i,j)] - vy[ind(i,j-1)]) + + ( (ETA_VS[ind(i,j)])*dtx) * (dvx[ind(i,j)] - dvx[ind(i-1,j)]) + + ( (ETA_S [ind(i,j)])*dtx) * (dvy[ind(i,j)] - dvy[ind(i,j-1)]) ); + + Tyy[ind(i,j)] += ( ( C12[ind(i,j)]* dtx )*(vx[ind(i,j)] - vx[ind(i-1,j)]) + + ( C11[ind(i,j)]* dtx )*(vy[ind(i,j)] - vy[ind(i,j-1)]) + + ( (ETA_VS[ind(i,j)])*dtx) *(dvy[ind(i,j)] - dvy[ind(i,j-1)]) + + ( (ETA_S [ind(i,j)])*dtx) * (dvx[ind(i,j)] - dvx[ind(i-1,j)]) ); + } + + if (i