From b349147ea2ad0c294085dedd0f5f86013c40c502 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jose=20Luis=20Cerc=C3=B3s=20Pita?= Date: Fri, 10 Aug 2012 16:57:58 +0200 Subject: [PATCH] Improved free surface evolving method in order to support long simulations. --- src/Mod/Ship/simRun/Sim/fsEvolution.py | 70 +++++++++++++----------- src/Mod/Ship/simRun/Sim/fsEvolution.pyc | Bin 4941 -> 0 bytes src/Mod/Ship/simRun/Simulation.py | 8 --- src/Mod/Ship/simRun/Simulation.pyc | Bin 4343 -> 0 bytes 4 files changed, 38 insertions(+), 40 deletions(-) delete mode 100644 src/Mod/Ship/simRun/Sim/fsEvolution.pyc delete mode 100644 src/Mod/Ship/simRun/Simulation.pyc diff --git a/src/Mod/Ship/simRun/Sim/fsEvolution.py b/src/Mod/Ship/simRun/Sim/fsEvolution.py index a4d2e0161..14fdf4ba3 100644 --- a/src/Mod/Ship/simRun/Sim/fsEvolution.py +++ b/src/Mod/Ship/simRun/Sim/fsEvolution.py @@ -45,33 +45,19 @@ class simFSEvolution: @param t Actual time (without adding dt). """ self.fs = fs - # Allocate memory nx = self.fs['Nx'] ny = self.fs['Ny'] nF = nx*ny - # Evaluate potential gradients grad = self.evaluateGradient() - # Integrate variables + # In order to improve results in really long simulations free surface + # will performed considering external waves and second order effects + # in two different ways. First external waves at time t will be + # substracted, then second order waves will be computed, and finally + # external waves at t+dt will be added. for i in range(0,nx): for j in range(0,ny): - # Get value at pos using characteristics method - gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j]) - gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal) - self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + dt*gradVal - # Velocity potential - self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + \ - dt*self.fs['accPot'][i,j] + \ - 0.5*dt*dt*grav*self.fs['pos'][i,j][2] - # Acceleration potential. This is really hard to simulate - # accurately due to numerical diffusion of the function, so - # external waves, and diffracted waves will be computed - # in two different ways: - # * External waves will be considered analitically, - # substracting waves at t, and adding waves at t+dt - # * Second order waves will be computed substracting external - # waves to free surface height, and then imposing boundary - # condition. pos = np.copy(self.fs['pos'][i,j]) + # Substract external waves at time t. for w in waves['data']: A = w[0] T = w[1] @@ -81,16 +67,39 @@ class simFSEvolution: k = 2.0*np.pi/wl frec = 2.0*np.pi/T l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading) - # Substract external waves height in order to know second - # order waves free surface amplitude. - amp = A*np.sin(k*l - frec*(t+dt) + phase) + amp = A*np.sin(k*l - frec*t + phase) pos[2] = pos[2] - amp - # Compute analitic external waves acceleration potential - amp0 = grav*A*np.cos(k*l - frec*t + phase) - amp1 = grav*A*np.cos(k*l - frec*(t+dt) + phase) - self.fs['accPot'][i,j] = self.fs['accPot'][i,j] - amp0 + amp1 - # Now impose free surface boundary condition - # self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + grav*pos[2] + amp = - grav/frec*A*np.sin(k*l - frec*t + phase) + self.fs['velPot'][i,j] = self.fs['velPot'][i,j] - amp + amp = grav*A*np.cos(k*l - frec*t + phase) + self.fs['accPot'][i,j] = self.fs['accPot'][i,j] - amp + # Now compute second order waves using position copy, + # where external waves are excluded, in order impose + # free surface boundary condition relative to second + # order phenomena. + self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + \ + dt*self.fs['accPot'][i,j] + # self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + \ + # grav*pos[2] + # Restore external waves to velocity and acceleration + # potentials. + for w in waves['data']: + A = w[0] + T = w[1] + phase = w[2] + heading = np.pi*w[3]/180.0 + wl = 0.5 * grav / np.pi * T*T + k = 2.0*np.pi/wl + frec = 2.0*np.pi/T + l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading) + amp = - grav/frec*A*np.sin(k*l - frec*(t+dt) + phase) + self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + amp + amp = grav*A*np.cos(k*l - frec*(t+dt) + phase) + self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + amp + # Update free surface point position + gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j]) + gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal) + self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + dt*gradVal # Impose values at beach (far free surface) for i in range(0,nx): for j in [0,ny-1]: @@ -107,15 +116,12 @@ class simFSEvolution: ny = self.fs['Ny'] nF = nx*ny grad = np.ndarray((nF,3), dtype=np.float32) - FF = open('gradient', 'w') for i in range(0,nx): for j in range(0,ny): pos = self.fs['pos'][i,j] grad[i*ny+j] = self.gradientphi(pos) gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j]) gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal) - FF.write('%g\t%g\n' % (pos[1], gradVal)) - FF.close() return grad def gradientphi(self, pos): diff --git a/src/Mod/Ship/simRun/Sim/fsEvolution.pyc b/src/Mod/Ship/simRun/Sim/fsEvolution.pyc deleted file mode 100644 index e8cf2ab57f85012e37793c30e682c17326e5a245..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 4941 zcmcgwUvJ~a5#Oa`$&&5+Q^=Oxv^fTs;+)huw@n{{)@Ynua!JtCccp{6HEC`_DiB^!A~-MrMNE{o}XK<8WYtI9l-XCf)uqgv#v0EdGV~1diA$ zT%4B4FEwzh#Er_-`)O4;nO~1$rpv;050~D87)vo0qJ@UV4Lkt_N@MgF5ca(|(k2~v zCQk3x?dj(*aZ|TvdvRpc(AW&x&2OPU0NSghR`Z( zI}#0tmSPuICiMGg?RVpz+D}vhk;Ws1#nG*Hs+8yY`}F+^HFWCj|19?RkGny#j~#Xf z(f)DJ+i&aHPTx&NW)%iotsn|atMxqEn%s+OO?nj8IxZT)8}@s|GT$&E&KP(*lNH)J z>R4iDQW0bHE21%}vRV=XC#z!k7ik!oxc$TuZ$0{4PHJNKhg@69w4cgJU6_(ME3<^D zRWYuy!Es$Qcb)(+OBzc|)_8D*Wi_9;CNP-BqrqdQQI^FHS*(jT;su&_#wDJzD$bx@ zGA>1GWgk{hHYB&=%hk*h~s$jrS03oN^2)izh2W z_V=ul{-#y1HpKv})pKiLl!7(v&!yHrhBePR@$2A^H9%JX6$1RxO9n>SbS!K=NZYQ5 z4)`D<<}g7ITBDMuPJUb4_lr{k(2yH>PgBH&d#<#7!t&BVzOg^J)RnQ{@XWvsvuJiN z0@IBL#&&%_h&r}!UYZAoRTtuDn0<``bO^aZB(V{9cnhw*5(G<>f8+nu2Ju#t{nW zPcx5U5V>K8by1x5xY0@7b8{0@AIF2pcheC_-47_$a3~I0b?Ol8m=}~C`^HefCW|6_ z>S&bp?&p?$9uu%UV_Jnf{}*m}m841Rt%k}2vdE3zg8*k*8*)R|t#xZlZdtF&7v;8m z&3d==3SmDxs8N~G1O0X!W#hg-cXwbJnsW`7vobe&TQ3HP= z$CIW6o|BWZ7B)+wA#WNXTEjV2J@|0OW^rXZTf zU^qHpG&14;8f!Bn{Ctx$1o}z`sqz6)&BQYY*tC|!4FEBSkP~7%C&a?2OR+j1vjkP= z1%xw4Df5+^DtCdXk5jLMVf@LkNIaO~n2`UmN@MNZCLP0m zGl%gaf;k_cCVU8c-@!oUa`ds-*u*O&OXo#6c!W#S;SwtZ?QH_KWdqvJ$);78FXQs) zkkez$AbyAH%?y6v#Eb)(rBGr3?~o;NK&@;H)pg917*kMFfj3hR|0Yp@_#nNv9{ouY zh=Eoh7zOGSwRbh0k|^(j3!#_JflQB>5R zxRyt{{513^2!br)0C3$(jGG@LoUTmQoLr|9MdP5`5aa7O7_f>?Su`8lI4?~cx5nhg zGd6MAvqie`)YgZ?R~U<2uh zUt`|^mS+U8Q1Q8V$IFV3KdJI+C#HAA*$ZJ1d4Z3FWK$txdu~K;hsOO# zMKuSZBzpGG&L>;aeFG3#j5yy@_%7`yypKdR? zN&?(On3^tjQftaOh_5Mc;+y_9`k()J;J}BV8c9=R6?wg0yHmSX+sgb>RJ$^vOs!)wnDB+uRfI<(wm#g^U|B2$_vt4P`NJIyyQ*k zH6?FJuO-1`#>g&52o-11vi)KZv{(`KGrbg9R~x z3~y*hbJlT_R(3lr*q;6ezfJ(jvHEDYQ={BK^+jxrn}_g~NnBZV=(C$cwHLc1>)Md^ ztGKXPj9wk_3>}WMcGKc8Gp;OV0psLfpxC!j8JPqKmiB549)-F@U0prl>MUGr8^wNr zY9MkZ(htPdWLT$s7RXsm&g$~Ap@HhEAG;k6ReN^J!njqu^=Ah4`)OfaTqH&{-08@; zCN-HjGlajPcYCvld+d#)5+V>8!_t3!)ui#&hv zhWuT|zYoKuym)D4?#w0uMi9o^WSa zxE0$d_HU>hHXK04X$Q6g32+R##SO&(ggHGDI_>kCR6CrX3ZNA)G^O}NfXyLjM32XC zP7VX`MBvrfbxkK(-~w!HD!L#n5@T#xk~6Fb;kumZiv49do7Wa_4Jn02+X!(8_m3(Q z55lAb+^R7NWV0Pkeh)@*73bl=9H)s1ANhwX!?872luR+$PjeGKi%)=*JOpeecnKZa zujrR!HyXR}+bDKsV;}JTF8tCrRcr~_^!X=YdJyjHbvMFKDLviWxiSydc{*W=&~LSZ zG@R1Ml>s~5+}NX#CnU7v>L7+|+Cs@*QGXF18ykMDwT8nuQ`xZdaASXCbd+u!A$nlZ zpU;Bkcqn|x>~7}<4~61a;|xuuY9Dn)h;r6WoJ~kPPF#BI`P@&U6rhu0kDNJiI^g@U z85_lXl4%lF^{!ba0-jfmZ3R?vw z0^J0HhZ3%^MIr9)6I=j7){y-|k;a!f**Eht1Qu~s1#%vUY06Nsh%1b6HXeEa8gLGj zp$5cIFY|6YuZdatabr!+>tb$PX>5qOd8H9h#g+M%oCB&5R^@y_g(W#xj1O7_>7k%}fH(%j7m}Pi}3xP*o z)({D=-g=dK<0sYIuTr0zs>9Ki+O?+H0xg7{Ex8S;zjm=91NgcoYuuFYZvevZW^{GI02et|1YLjrTHk_X zzWa6f6j?|bXX&57HpSuQ=4ShlI{YarcTfD|)bcu-jRJSf23TIGt!dX2Q%w58EG*Bq1#=u;T3Q(?{5-sW7~1x3=RKOu zTW!A@o9HeC4|?SR-W70K=^@D=Mbz}9HpwfEQqbhmozjslky;Xvdli;{W~_}5jUyo< z*q7s178e{6t)oA>#|1bPnn+gRKf*H|3UkFkuWM&qK}?9)t-rbPFt#bD(2P=m!xRQVHI zMEgiQaQfI(Djrx`6#bIJn4Il&l_jHJO&bb-MvBunL!4@BH7<5()($GEw`$A5Qmqa4 zTnXL}Zt35YKU+Z?f9zS;dgW4m7p<$U2Do;`XS&Rw=+iceC7h7>cmRS}{nvrzWf5*8 zUk#BmSXYihA9SYBv)j^t$M|jFxzkSr4fPue3ePT$`7@t%LFx&;nHKek=N|4LclQ~& z{RRzTWJdUa0$pgvgs&54W_~=Vj2l-)NB}WanJA?;#10g|=|-q`d{}~Qsu$_y49ffb zfAOBUPc2xFCjA93J_LAZ1qlEv0eGkQiSdc^86s*7z>>9yYc ztC*jkQI;T_50Krjw?1plb!eg@5WqX~=%xLm({f~rB#X#@isxdSk4_bQ#VGoit(jK))!Ik1Z$}8)z;eo1^<1eK>z>%