调研概述：基于激光加速器的小型化自由电子激光 本调研以中国科学院上海光学精密机械研究所发表于Nature的论文为中心，介绍了我国在基于激光加速器的小型化自由电子激光领域取得的成就。从自由电子激光器（FEL）、激光尾场加速两个核心模块，原理、系统构成、应用三个方面简要介绍了该激光技术的特色点与先进处。自由电子激光器为我们展现了产生激光的另一种原理，激光尾场加速则向我们揭示了非线性系统独特的特性。

换热器排管助手V1.0: 正式版使用说明及优化说明 增速7400倍!换热器排管助手V1.0正式版更新！Github仓库地址为：https://github.com/DMCXE/Pipe_arrangement

换热器排管助手V0.9：极低运行效率下的换热器排管及绘制 项目源起大三下学期，涉及了两个关于换热器的课程设计&核动力设备大作业，不出意外，过段时间的课设三也是直立式蒸汽发生器。无论在课设还是大作业，排管始终是一个令人头疼的问题。课设中1.1倍根号下总管数确定中心管数属实令人费解（考虑优化的时候知道是怎么回事了）。核动力设备大作业中也因为没有合理的排布（类圆形）而被老师批评。对于数百个乃至数千根管，如何快速、简便、合理的绘制成了萦绕在我心头解不开的问题。于是从暑假拖到了现在。程序目标程序将要实现：对于给定的管间距参数，在任意排管方式（正方形、正三角形）下，对于任意的管数N，都能够给出一种排布方式及其排布图。且能够被其它程序方便引用。设计思路由于编者能力水平低下，此篇文章中，将介绍一种目前效率极其低下但容易实现的思路。依据排布要求，对于总数N，生成第一象限中（含坐标轴）范围为(N,N)的圆心位置。假定一个半径为R的圆，遍历所有圆心位置，取出所有在圆内且满足边界条件的圆心位置。计算R下符合条件的圆心位置数量是否近似满足总管数在第一象限分布的数量。若不满足，则令R递增，直至满足为止。坐标变换，通过第一象限的点生产全象限的完整分布绘制排布图。思路非常暴力，非常简单，非常低效。但是正在优化，有了初步思路。程序设计总览为了便于调用，整个排布图将设计为类形式，它包含：class Pip_arrangement: def __init__(self,s1,s2,s3,e,r,N) #传入参数 def point_squar(self)#正方形排布 def pos_stard(self)#第一象限范围筛选 def arrangement(self)#全排布 def visualize(possiton,r)#可视化 def test_visualize()排布类型基本思路为，对于总数N，生成第一象限中（含坐标轴）范围为(N,N)的圆心位置。实际上，可以通过小圆与大圆面积近似相等进一步缩短生成范围以大幅度减少复杂度。粗略计算，集合内元素含量可以由目前的N^2简化到N。这里给出正方形排布的生产程序。def point_squar(self): pos = np.zeros([1,2]) for i in range(self.N): for j in range(self.N): pos = np.append(pos,[[self.s1*j,self.s3+self.s2*(i+1)]],axis=0) pos = pos[1:] return pos范围判断这里判断圆心以及圆是否处于半径范围内的方法仍然是遍历。更改生成范围和采用矩阵化操作都能够大幅度的加快计算速度。遍历的思路是：取出坐标，判断条件满足，满足即放入新数组，不满足扩展大圆半径。def pos_stard(self): pos = copy.deepcopy(self.point_squar()) R = self.s1 + self.e count_tot = 0 count_Y = 0 while count_tot < (self.N-2*count_Y)/4+count_Y: count_Y = 0 pos_stard = np.zeros([1, 2]) for p in pos: if p[0]**2 + p[1]**2 <= (R-self.e-self.r)**2: pos_stard = np.append(pos_stard,[p],axis=0) for pp in pos_stard: if pp[0] == 0.: count_Y = count_Y+1 pos_stard = copy.deepcopy(pos_stard[1:]) count_tot = len(pos_stard) R = R + 1 #mm return R,pos_stard,count_tot坐标扩展仍然是遍历，将第一象限的坐标通过对称操作变换到全部象限空间中。需要注意的是，python对数组的操作会影响原数组，因此需要使用deepcopy的方法。def arrangement(self): pos = copy.deepcopy(self.pos_stard()[1]) pos_first_with_axis_y = copy.deepcopy(self.pos_stard()[1]) R = self.pos_stard()[0] pos_second_without_axis_y = np.zeros([1,2]) for ppp in pos_first_with_axis_y: if ppp[0] != 0: ppp[0] = -ppp[0] pos_second_without_axis_y = np.append(pos_second_without_axis_y,[ppp],axis=0) pos_second_without_axis_y = pos_second_without_axis_y[1:] pos_up = np.append(pos,pos_second_without_axis_y,axis=0) ppos = np.append(pos,pos_second_without_axis_y,axis=0) pos_down = pos_up for q in pos_down: q[1] = -q[1] pos = np.append(ppos,pos_down,axis=0) L =len(pos) return R,pos,L可视化可视化需要先使用matplotlib.artist中的Artist函数，可以帮助我们在确定圆心坐标和半径的情况下批量画出圆。def visualize(possiton,r): dr = [] figure, axes = plt.subplots() for xx in possiton: dr.append(plt.Circle(xx, r, fill=False, lw=0.25)) axes.set_aspect(1) for pic in dr: axes.add_artist(pic) Artist.set(axes, xlabel='X-Axis', ylabel='Y-Axis', xlim=(-300, 300), ylim=(-300, 300), title='Arrangement of heat transfer tubes') plt.savefig('test111.png', dpi=1200, bbox_inches='tight') plt.title('Arrangement of heat transfer tubes') plt.show()用法传入需要的变量，通过Pip_arrangement(s1,s2,s3,e,r,N).arrangement()可以得到半径、位置坐标等信息。def test_visualize(): s1 = 3.2 s2 = s1 * np.sqrt(3)/2 s3 = 2 * s1 e = 0.1 r = 1.1 N = 1000 pipe = Pip_arrangement(s1,s2,s3,e,r,N) pos = pipe.arrangement()[1] visualize(pos,r)全部代码import numpy as np import matplotlib.pyplot as plt from matplotlib.artist import Artist import copy import time class Pip_arrangement: def __init__(self,s1,s2,s3,e,r,N): self.s1 = s1 self.s2 = s2 self.s3 = s3 self.e = e self.r = r self.N = N def point_squar(self): pos = np.zeros([1,2]) for i in range(self.N): for j in range(self.N): pos = np.append(pos,[[self.s1*j,self.s3+self.s2*(i+1)]],axis=0) pos = pos[1:] return pos def pos_stard(self): pos = copy.deepcopy(self.point_squar()) R = self.s1 + self.e count_tot = 0 count_Y = 0 while count_tot < (self.N-2*count_Y)/4+count_Y: count_Y = 0 pos_stard = np.zeros([1, 2]) for p in pos: if p[0]**2 + p[1]**2 <= (R-self.e-self.r)**2: pos_stard = np.append(pos_stard,[p],axis=0) for pp in pos_stard: if pp[0] == 0.: count_Y = count_Y+1 pos_stard = copy.deepcopy(pos_stard[1:]) count_tot = len(pos_stard) R = R + 1 #mm return R,pos_stard,count_tot def arrangement(self): pos = copy.deepcopy(self.pos_stard()[1]) pos_first_with_axis_y = copy.deepcopy(self.pos_stard()[1]) R = self.pos_stard()[0] pos_second_without_axis_y = np.zeros([1,2]) for ppp in pos_first_with_axis_y: if ppp[0] != 0: ppp[0] = -ppp[0] pos_second_without_axis_y = np.append(pos_second_without_axis_y,[ppp],axis=0) pos_second_without_axis_y = pos_second_without_axis_y[1:] pos_up = np.append(pos,pos_second_without_axis_y,axis=0) ppos = np.append(pos,pos_second_without_axis_y,axis=0) pos_down = pos_up for q in pos_down: q[1] = -q[1] pos = np.append(ppos,pos_down,axis=0) L =len(pos) return R,pos,L def visualize(possiton,r): dr = [] figure, axes = plt.subplots() for xx in possiton: dr.append(plt.Circle(xx, r, fill=False, lw=0.25)) axes.set_aspect(1) for pic in dr: axes.add_artist(pic) Artist.set(axes, xlabel='X-Axis', ylabel='Y-Axis', xlim=(-3000, 3000), ylim=(-3000, 3000), title='Arrangement of heat transfer tubes') plt.savefig('test111.png', dpi=1200, bbox_inches='tight') plt.title('Arrangement of heat transfer tubes') plt.show() def test_visualize(): s1 = 3.2 s2 = s1 * np.sqrt(3)/2 s3 = 2 * s1 e = 0.1 r = 1.1 N = 1000 pipe = Pip_arrangement(s1,s2,s3,e,r,N) pos = pipe.arrangement()[1] visualize(pos,r) if __name__ == '__main__': st = time.time() test_visualize() stt = time.time() print(stt-st)运行结果在正方形排布下（六边形懒得做了，优化后会做的），数据来源来自于核动力设备课程大作业——直立式蒸汽发生器中相关参数，1000根管子下的排布，仅仅是1000根，就用时1620s，27分钟！！！！不过得出的结果还算令人满意优劣分析及展望更多的数学准备以减少数据量。对于总数N，生成第一象限中（含坐标轴）范围为(N,N)的圆心位置。实际上，可以通过小圆与大圆面积近似相等进一步缩短生成范围以大幅度减少复杂度。这一步已经算出来了。重复利用矩阵运算的优势。 以Numpy为代表的一系列优秀计算库优化了矩阵运算，大幅度的提高运算速度。相比于遍历的土方法，显然矩阵运算更为优秀。更优的算法。算法灵感来源于蒙特卡洛（特指扔硬币算圆周率那个）（可能因为都有圆），肯定存在更优秀的算法。灵活度。实际存在更严苛的限制条件，可以进一步的拓展可拓展性。算是本篇唯一一个优点了罢，正六边形排布改一下生成就可以实现。总结假期拖延症严重。明天又要线上了。

单粒子在场中的运动（一）：一种粗糙的计算 项目源起在亚布力的时候，偶然的一次滑水，在GitHub上找到了某个描述在某个感兴趣装置中的粒子运动模拟。欣喜之余仔细阅读了源码与作者的项目介绍，不免大失所望，不仅缺乏很多必要的头文件，且作者直言并没有足够的知识去解决这个问题。因此处于好奇，不免想要尝试一下。（虽然经过短暂的构思，发现确实存在很多问题）。初期尝试由于先前想要计算的问题十分复杂，而且发生的过程极快，因此难以在现有知识水平下进行程序的搭建。突然想到了某些圣遗物和某个被封了很多年知识科普团队，于是定下了第一步计划：在不考虑粒子间相互作用、不考虑相对论效应、不考虑Maxwell方程的前提下，单粒子或多粒子在常电磁场中的三维运动可视化。且求解过程基于物理方程描述。（即忘记粒子是如何在场中运动的结论）目录{timeline}{timeline-item color="#19be6b"}物理过程的数学表达{/timeline-item}{timeline-item color="#19be6b"}程序化准备{/timeline-item}{timeline-item color="#19be6b"}Codeing{/timeline-item}{timeline-item color="#19be6b"}运行结果{/timeline-item}{timeline-item color="#ed4014"}缺点分析及展望{/timeline-item}{timeline-item color="#ed4014"}结束语{/timeline-item}{/timeline}物理过程的数学表达“当听者以为自己已有所得时，真正的内行会觉得连皮毛也不是。” ――克莱因，《数学在19 世纪的发展》 我想，单粒子轨道这个看似简单的问题可能也深邃到我们尚只理解其皮毛 。 ——《计算等离子体物理导论》单粒子在场中运动，从入门角度而言，任何上过高中的理科生都可以回答这个问题：磁场改变运动方向，电场改变粒子速度。通过高中的大量训练，我们知道如何联立洛伦兹力与向心力求解平面中的粒子运动，复杂一点的情况下，小球受轨道电场力平抛等等。再最开始，我也想直接进行一番物理计算，把运动模式给出来，即给出闭式方程。但这显然是不科学的。场在空间的复杂分布会超出我们的计算能力，且并不是所有物理过程的数学描述都有通解。在电磁场中，带电粒子所受的洛伦兹力为：运动轨迹由Newton‘s Law给出：通过以上方程组，理论上可以求解出空间中粒子在常电磁场作用下的任意运动。程序化准备接下来我们简要的分析一下，为程序化做准备（离散化方程）。在笛卡尔坐标下，任何物理量都可以由三维分量表示，比如：电场E = [E_x,E_y,E_z]，磁场B =[B_x,B_y,B_z]，场力F = [F_x,F_y,F_z]粒子速度V= [v_x,v_y,v_z]，粒子位置Position= [X,Y,Z]对于任意时刻t，都可以计算出时刻t下的F，由：采用向后差分（方程的解可由初始条件确定），可离散化为：整理为：得到当前时刻下的速度，进一步可以得到当前时刻的位置，由整理为：最后，我们能够得到代求解的方程组为：当然，用前朝的剑斩本朝的官，肯定存在大问题。后文会总结，未来会改进。在程序（Python）中，这段表示为：F = p.q * (E + np.cross(p.v,B)) p.v = p.v + F * (timestep/p.m) p.pos = p.pos + p.v * timestep如果采用比如数学功能较强的软件，只需要给出最开始的数学描述微分方程即可。Codeing剩下的工作就是Coding，让它能够在符合语法规范的情况下运行。总的来看，主要包含以下几个模块：class Particle #粒子 class Field #场 class ParticleSimulator #粒子运动 def visualize(simulator) #可视化 def test_visualize() #运行粒子粒子是粒子仿真的重点，具有质量与电荷，需要定义初始位置与初始速度。定义粒子类：class Particle: def __init__(self,start_pos,v,mass,Quantity): #start_pos , v is 3 dim self.pos = start_pos self.m = mass self.q = Quantity self.v = v self.sum = np.zeros([1,3]) #储存粒子的历史位置场场是粒子在场中运动的重点，目前虽然是单一指向的常场，但长远考虑，需要建立一个场类：class Field: def __init__(self): ''' self.B = np.array([20, 20, 30]) self.E = np.array([10., 20.,-15]) ''' self.B = np.array([0, 0, 0.5]) self.E = np.array([0, 0, 0.2])粒子仿真粒子仿真按时步进行，对于每一个时步求解上述方程组并更新粒子位置与记录路径总和。在一个ParticleSimulator类中，对于时间dt的更新方式为：def evolve(self,dt): #参数传入更新时间dt timestep = 1e-3 #确定更新的最小时步 nsteps = int(dt / timestep) #分割时间内计算空间 #传入场参数 f = self.f B = f.B E = f.E for i in range(nsteps): for p in self.particles: #控制方程求解 F = p.q * (E + np.cross(p.v,B)) p.v = p.v + F * (timestep/p.m) p.pos = p.pos + p.v * timestep p.sum = np.append(p.sum,[p.pos],axis=0)可视化可视化利用matplotlib库中的animation.FuncAnimation()函数，对于指定的ParticleSimulator.evolve粒子运动更新方式，按照函数传入的设定参数进行绘制更新。这里更多的是函数的使用方式。比如绘制时.set_data()方法只允许传入两个坐标参数，第三维需要额外使用.set_3d_properties()传入。更多的细节不再赘述。line.set_data(X[i],Y[i]) line.set_3d_properties(Z[i])使用方式在一个定义的test_visualize()中，首先设定粒子信息：#单粒子时 particles = [Particle(position, v, mass, q)] #多粒子时 particles = [Particle(position1, v1, mass1, q1), Particle(position2, v2, mass2, q2), Particle(position3, v3, mass3, q3)]接着设定场field = Field()传入仿真参数simulator = ParticleSimulator(particles,field)可视化visualize(simulator)运行时调用test_visualize()即可。代码总览import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation class Particle: def __init__(self,start_pos,v,mass,Quantity): #start_pos , v is 3 dim self.pos = start_pos self.m = mass self.q = Quantity self.v = v self.sum = np.zeros([1,3]) class Field: def __init__(self): self.B = np.array([0, 0, 0.5]) self.E = np.array([0, 0, 0.]) class ParticleSimulator: def __init__(self,particles,Field): self.particles = particles self.f = Field #最小时间dt内更新轨迹 def evolve(self,dt): timestep = 1e-3 nsteps = int(dt / timestep) f = self.f B = f.B E = f.E for i in range(nsteps): for p in self.particles: F = p.q * (E + np.cross(p.v,B)) p.v = p.v + F * (timestep/p.m) p.pos = p.pos + p.v * timestep p.sum = np.append(p.sum,[p.pos],axis=0) def visualize(simulator): X = [p.pos[0] for p in simulator.particles] Y = [p.pos[1] for p in simulator.particles] Z = [p.pos[2] for p in simulator.particles] fig = plt.figure() ax = fig.add_subplot(projection = "3d") line, = ax.plot(X, Y, Z) ax.set(xlim3d=(-20,10), xlabel='X') ax.set(ylim3d=(-20,10), ylabel='Y') ax.set(zlim3d=(-10,20), zlabel='Z') def init(): line.set_data([], []) line.set_3d_properties([]) return line, def init2(): line.set_data([], []) line.set_3d_properties([]) return line def animate(i): simulator.evolve(0.1) X = [p.sum[:,0] for p in simulator.particles] Y = [p.sum[:,1] for p in simulator.particles] Z = [p.sum[:,2] for p in simulator.particles] line.set_data(X[i],Y[i]) line.set_3d_properties(Z[i]) return line, anim = animation.FuncAnimation(fig, animate, frames=len(X), init_func=init, blit= True, interval=1) plt.show() def test_visualize(): p1 = np.array([0, 0, 0]) v = np.array([1,0,0]) particles = [Particle(p1, v, 3, 1)] field = Field() simulator = ParticleSimulator(particles,field) visualize(simulator) if __name__ == '__main__': test_visualize()运行结果运行环境处理器：基于Apple Silicon的Arm64构架的M1内存：16GB DDR4操作系统：macOS Monterey 12.4语言：miniforge虚拟环境下的Python3.9具有初速度的粒子在恒定磁场下运动带电粒子在具有一定沿X方向初速度，在沿Z方向恒定磁场下运动，所得的图像如下图所示：根据高中知识可以知道，是圆形轨道。具有初速度的粒子在恒定电磁场下运动带电粒子在具有一定沿X方向初速度，在沿Z方向恒定电磁场下运动，所得的图像如下图所示：初速度为0的粒子在任意电磁场下的运动带电粒子在某一任意指向电磁场下的运动如下图所示：单粒子在受控电磁场下被约束理想情况下，单粒子在受控的电磁场作用下被约束在特定范围内。多粒子程序在二维平面内支持多粒子运动，但当拓展至三维空间内，由于程序逻辑，粒子是交替产生的。缺点分析及展望目前，该项目存在非常多的缺陷。计算逻辑低效。本次演示的计算逻辑为绘图驱动计算。即计算过程是伴随着动画的渲染而不断更新的。因此运行低效、处理困难、规模繁琐。——改进方向：先进行粒子运动计算，再由数据驱动路径绘制。数值解法能量耗散。向后差分的方法会造成体系能量的耗散，会造成系统随着时间演化精度不断下降。——改进方向：采用更先进的差分方法，比如通过哈密尔顿量计算，确保体系完备。运行效率低。随着运行时间增长，Particle中存储了大量的解，且每次绘图都需要重新遍历，时间开销逐步增大。通过改进缺点1可以克服物理场复杂化计算可靠性未验证。当物理场随着空间分布后，计算是否任然可靠未验证。拓宽粒子相互作用与Maxwell方程组结束语又是不务正业的几天时间，别人该写实习报告的写实习报告，该学习预习的预习，该搞项目的搞项目；能力不足还想学人造轮子，只能是自讨苦吃！

DeepLearning学习笔记之一：重新认识 第一次接触DL还是之前提过的DDPG。但实际上就如同学习最优策略的数学理论一样，对那些泛函内容只能照葫芦画瓢般的模仿。所以还是决定从头认识一下他，看看他到底离我多远，或者离我多近。这一切的前提是，学过高数、线性代数与概率论，以及一些计算机语言和思维。这些大一就已经经历过了。要理解DeepLearning，我们不妨从一个线性回归说起——这是初中的知识，在一些点之间划一条直线。呃呃，我又发现我现在还是打不了公式代码，我再也不用这个数据库了！那我先介绍几个基本的语言实现。由于第一次接触DL使用的是Tensorflow 1.X——一个开源的DL框架，自然的我换到MacOs平台后第一个想要验证的就是它。巧合的是Tensorflos以及Apple Inc.已经为M1提供了原生ARM般程序与基于METAL的GPU加速，在经历一番艰难的Tensorflow 1.X->2.X代码修改后，老程序还是跑起来了。首先介绍第一个重要的代码实现：它能够算出任意可导函数的导数。（还在研究争取造轮子）with tf.GradientTape() as t: y = 2*x*x x_grad = t.gradient(y,x)之后是随机梯度下降算法SGD暂时看不懂没关系，等我睡一觉起来把公式模块弄出来就好了，困了，晚安！手写实现为：def sgd(params,grads,lr,batch_size): for param, grads in zip(params,grads):#zip使得w和b拆分成列 param.assign_sub(lr*grads/batch_size)#这里就是SGD算法，.assign_sub表示将原来存在的param元素减去（）内 API实现一行即可：trainer = tf.keras.optimizers.SGD(learning_rate=0.03)随便查个图