--- trunk/doc/user/wave.tex 2006/03/08 05:45:51 580 +++ trunk/doc/user/wave.tex 2006/03/08 05:48:51 581 @@ -20,25 +20,18 @@ \begin{eqnarray} \label{WAVE natural} \sigma\hackscore{ij}n\hackscore{j}=0 \end{eqnarray} -for all time $t>0$. At initial time $t=0$ the displacement -$u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are give: +for all time $t>0$. + +At initial time $t=0$ the displacement +$u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given: \begin{eqnarray} \label{WAVE initial} -u\hackscore{i}(0,x)=0 & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \; +u\hackscore{i}(0,x)= \left\{\begin{array}{cc} U\hackscore{0} & \mbox{for $x$ at point charge, $x\hackscore{C}$} \\ 0 & \mbox{elsewhere}\end{array} \right. & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \; \end{eqnarray} for all $x$ in the domain. -The points on the left face of the domain, ie. points $x=(x\hackscore{i})$ with $x\hackscore{0}=0$, -are moved into the $x\hackscore{2}$ direction by $u^{max}$. This is expressed by the constraint -\begin{eqnarray} \label{WAVE impact} -u\hackscore{2}(x,t)=s(t) \mbox{ for } x=(x\hackscore{i}) with x\hackscore{0}=0 -\end{eqnarray} -where -\begin{eqnarray} \label{WAVE impact s} -s(t)=u^{max} (1-e^{-(\tau^{-1}t)^3}) -\end{eqnarray} -$\tau$ measures the time needed to reach the maximum displacement $u^{max}$ ( -actually $\tau (ln 2)^{\frac{1}{3}}) \approx 0.9 \tau$ is the time until $\frac{u^{max}}{2}$ is reached). -The smaller $\tau$ the faster $u^{max}$ is reached. +Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we + set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point + $x\hackscore C$. We use an explicit time-integration scheme to calculate the displacement field $u$ at certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n-1)}+h$ with time step size $h>0$. In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. We use the Verlet scheme \index{Verlet scheme} with constant time step size $h$ @@ -46,7 +39,7 @@ \begin{eqnarray} \label{WAVE dyn 2} u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ \end{eqnarray} -for all $n=1,2,\ldots$. It is designed to solve a system of equations of the form +for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form \begin{eqnarray} \label{WAVE dyn 2b} u\hackscore{,tt}=G(u) \end{eqnarray} @@ -62,21 +55,17 @@ \end{eqnarray} derived from \eqn{WAVE natural} where \begin{eqnarray} \label{WAVE dyn 3a} -\sigma\hackscore{ij}^{(n-1)} & = & \lambda u^{(n-1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n-1)}\hackscore{i,j} + u^{(n-1)}\hackscore{j,i}) -\end{eqnarray}. -Additionally $a^{(n)}$ has to fullfill the constraint -\begin{eqnarray}\label{WAVE dyn 4} -a^{(n)}\hackscore{2}(x)= s\hackscore{,tt}(t^{(n)})\mbox{ where } s\hackscore{,tt}(t)=\frac{u^{max}}{\tau^2} -(6\tau^{-1}t- 9(\tau^{-1}t)^4) e^{-(\tau^{-1}t)^3} +\sigma\hackscore{ij}^{(n-1)} & = & \lambda u^{(n-1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n-1)}\hackscore{i,j} + u^{(n-1)}\hackscore{j,i}). \end{eqnarray} -derived from \eqn{WAVE impact}. - +Now we have converted our problem for displacement, $u^{(n)}$, into a problem for +acceleration, $a^(n)$, which now depends +on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$. In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through -the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which are relevant here are +the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The forms which are relevant here are $$\label{WAVE dyn 100} -D\hackscore{ij} u\hackscore{j} = - X\hackscore{ij,j}\; . +D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .$$ Implicitly, the natural boundary condition \label{WAVE dyn 101} @@ -86,48 +75,70 @@ \label{WAVE dyn 101} u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0 -where $r$ and $q$ are given functions. +where $r$ and $q$ are given functions. However in this problem we don't need to define +any constraints for our problem. -With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$, $r$ and $q$: +With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$: $$\label{WAVE dyn 6} \begin{array}{ll} D\hackscore{ij}=\rho \delta\hackscore{ij}& X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\ -r\hackscore{i}=\delta\hackscore{i2} \cdot s\hackscore{,tt}(t^{(n)}) & q\hackscore{i}=\delta\hackscore{i2} \cdot x\hackscore{0}.\texttt{whereZero()} \end{array}$$ -Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero -(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero. + + +%Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero +%(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero. -In the following script defines a the function \function{wavePropagation} which +The following script defines a the function \function{wavePropagation} which implements the Verlet scheme to solve our wave propagation problem. The argument \var{domain} which is a \Domain class object -defines the domain of the problem. \var{h} and \var{tend} is the step size -and the time mark after which the simulation will be terminated. \var{lam}, \var{mu} and -\var{rho} are material properties. \var{s_tt} is a function which returns $s\hackscore{,tt}$ at a given time: +defines the domain of the problem. \var{h} and \var{tend} are the step size +and the time mark after which the simulation will be terminated respectively. \var{lam}, \var{mu} and +\var{rho} are material properties. \begin{python} from esys.linearPDEs import LinearPDE -from numarray import identity +from numarray import identity,zeros,ones from esys.escript import * -def wavePropagation(domain,h,tend,lam,mu,rho,s_tt): +from esys.escript.pdetools import Locator +def wavePropagation(domain,h,tend,lam,mu,rho,U0): x=domain.getX() # ... open new PDE ... mypde=LinearPDE(domain) - mypde.setLumpingOn() + mypde.setSolverMethod(LinearPDE.LUMPING) kronecker=identity(mypde.getDim()) - mypde.setValue(D=kronecker*rho, \ - q=x[0].whereZero()*kronecker[1,:]) + # spherical source at middle of bottom face + xc=[width/2.,width/2.,0.] + # define small radius around point xc + # Lsup(x) returns the maximum value of the argument x + src_radius = 0.1*Lsup(domain.getSize()) + print "src_radius = ",src_radius + dunit=numarray.array([1.,0.,0.]) # defines direction of point source + mypde.setValue(D=kronecker*rho) # ... set initial values .... n=0 - u=Vector(0,ContinuousFunction(domain)) - u_last=Vector(0,ContinuousFunction(domain)) + # initial value of displacement at point source is constant (U0=0.01) + # for first two time steps + u=U0*whereNegative(length(x-xc)-src_radius)*dunit + u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit t=0 + # define the location of the point source + L=Locator(domain,xc) + # find potential at point source + u_pc=L.getValue(u) + print "u at point charge=",u_pc + u_pc_x = u_pc[0] + u_pc_y = u_pc[1] + u_pc_z = u_pc[2] + # open file to save displacement at point source + u_pc_data=open('./data/U_pc.out','w') + u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) while t