1. No category

Modelado Numérico de Tsunamis
un barniz…y algunos desafíos
Patricio Catalan
Profesor Adjunto
Departamento de Obras Civiles
UTFSM
2
¿Por Dónde Comenzar?
Un ejercicio de modelación numérica en realidad es un ejercicio de simulación, donde el
objetivo final es lograr una reproducción lo más precisa de la realidad o fenómeno de
interés
Es un problema complejo, que cuenta al menos con tres elementos claves
Algoritmos y
Software
Matemática
Modelado
- Física
Objetivos:
1. Desarrollar modelos matemáticos apropiados que describan la física de los tsunami
2. Desarrollar modelos numéricos precisos y eficientes para resolver estos modelos matemáticos
Ridgway Scott, PASI 2013
3
Partamos por la física
Partamos por la física
Después de todo, el objetivo/problema es reproducir la realidad
Previo… que representa para Uds. esto:
@p
@x
5
Definición de un tsunami
Perturbación repentina de la superficie de un cuerpo de agua, de gran escala, que genera
ondas de períodos muy largos
La definición formal
2004 Indian Ocean Tsunamis
Comparisons between model results and Jason-1 measurements
Serie de ondas
oceánicas
período
entre 5 y 60 minutos generadas por una
(Wang
and Liu, de
JHR
44, 2006)
perturbación
(left) and TOPEX measurements
(right) a gran escala del océano.
En ninguna parte se menciona la amplitud del tsunami
6
Ok…pero…
Esa definición no nos ayuda mucho, pero si se analizan los
datos se tiene lo siguiente
El tsunami es una onda de período muy largo…luego
su longitud es grande también: L
L > +100 km
Se propagan en el oceáno (principalmente), cuya
profundidad media es h = 4 km
Luego, la relación h/L es ≪ .
¿Por qué es relevante esto?
Por que en teoría de oleaje, si esa razón es pequeña, se
dice que es una onda de aguas someras
En ese caso, la velocidad de propagación de la onda es
C=√ (gh)
Resultados de campo sugieren que al menos el frente
principal se comporta así
1/17/13%
7
Algunas Ventajas de la Escala de los Tsunamis
Dos medidas de linealidad de los procesos
Profundidad/Longitud: h/L
Muy pequeño en el oceáno abierto : lineal
Permite aproximación de aguas someras (velocidad uniforme en la
vertical)
Amplitud/Profundidad: A/h
Relevante en aguas someras y zona de inundación : No lineal
8
Luego
Luego, es aparente que los tsunamis pueden ser modelados por teorías de onda larga…
¿Que significa esto?
Vamos a revisar un poco la mecánica de fluidos fundamental.
Requerimos poder representar el comportamiento cinemático y dinámico de un fluido.
Cinemático: Cómo se mueve, sin importar las fuerzas
Dinámico: Incorporando las fuerzas.
9
Cinemática: Conservación de la masa
d
dt
I
⇢d8 +
8
D
+
Dt
I
A
⇢~v · n
ˆ dA = 0
· ⇥v = 0
10
Dinámica
Esta es fácil:
F = m*a
11
Dinámica: Ecuación de Navier-Stokes
Supuestos: Fluido newtoniano, flujo incompresible
@u
+ (~u · r)~u = ~g
@t
1
2
rp + ⌫r ~u
⇢
Aceleración
Local
Aceleración
Advectiva
Fuerza de gravedad
(por unidad de
masa)
Fuerza de superficie
debida a la presión
Fuerza de superficie
debida a la
esfuerzos de corte
de origen viscoso
12
El sistema así no nos es útil…
Nos aprovechamos de la condición de aguas someras
Hipótesis:
Velocidad vertical despreciable: w=0
Velocidad horizontal prácticamente uniforme en la vertical: u(z)=u , v(z)=0
Distribución de presión hidrostática (depende linealmente de la profundidad
solamente)
Esto permite integrar las ecuaciones entre el fondo del mar y la superficie, provistas las
condiciones de contorno siguientes
Velocidad paralela al fondo
Velocidad paralela a la superficie
Esfuerzos en la superficie despreciables
Esfuerzos de fondo representados a través de la fricción
La integración requiere del uso de la Regla de Leibniz para evaluar las cantidades en la
frontera
13
El sistema final
1 @H
+ r · (H~u) = O(µ2 )
✏ @t
A
h
✓ ◆2
h
2
µ =
L
✏=
H = h + ✏⌘
@~u
+ ✏~u · r~u + r⌘ = O(µ2 )
@t
⌘ : Desplazamiento de la superficie libre
~u : Vector velocidad
MOMENTUM
NON LINEAR SHALLOW WATER EQUATIONS
CONTINUIDAD
14
NLSWE
he Shallow Water Equations
A Hyperbolic
differential
equation
EDP de tipopartial
hiperbólico
✔😄
Permite
uso schemes
de expresiones explícitas
• Enables
explicit
Puede presentar discontinuidades/shocks
Solutions
form discontinuities / shocks
Requiere de gran precisión sin oscilaciones cerca de las discontinuidades ✗
• Require
high accuracy in smooth parts
Requiere de cuidado en zonas secas/mojadas ✗
without oscillations near discontinuities
Solutions
include
dry
areasad-hoc
Requiere
de otros
términos
Corioliswater depths ruin simulations
• Negative
Rotura
Often high requirements to accuracy
• Order of spatial/temporal discretization
• Floating point rounding errors
COMMUN. MATH. SCI.
c 2007 International Press
⃝
Vol. 5, No. 1, pp. 133–160
A SECOND-ORDER WELL-BALANCED POSITIVITY PRESERVING
CENTRAL-UPWIND SCHEME FOR THE SAINT-VENANT SYSTEM∗
Can be difficult to capture "lake at rest"
ALEXANDER KURGANOV† AND GUERGANA PETROVA‡
A standing wave or shock
15
Modelo ampliado
Aceleración
de Coriolis
Presión no
hidrostática
Fricción de
Fondo
Difusión
Horizontal
Con esto, en principio el problema de encontrar un modelo matemático está solucionado
Algoritmos y
Software
Matemática
Modelado
- Física
16
¿Cómo se comportan estos modelos?
Probamos 7 modelos de uso frecuente en Chile que tienen como ecuaciones gobernantes
las presentadas anteriormente.
Para minimizar los efectos propios del tsunami, tratamos de mantener las mismas
configuraciones y dominio
Mismas condiciones iniciales
Mismo dominio de cálculo
Sin embargo, fuimos variando la resolución de la discretización
17
coefficient((thin(lines).(The(dark(lines(are(the(average(time(series(and(the(dotted(line(is(the(standard(deviation(among(runs.(
4
Comparación IntraModelos
2
η [m]
IntraCmodel(results(could(be(grouped(in(three(typical(model(response(to(changes(in(grid(resolution:(
0
-2
-4
6
Case A
-6
-8
-10
0
2
η [m]
• (Case* A)( High( sensitivity( in( free( surface( displacement( but( negligible(
dependency( on( phase( (arrival( time).( Variability,( estimated( by( the(
standard(deviation,(was(comparable(to(the(average(signal.(
4
• 2(models((fell(into(this(group.
0
-2
-4
-6
0
50
100
150
Elapsed Time [Min]
200
6
• 2(models(fell(into(this(group.(
300
Case B
•
•
4
2
η [m]
• (Case* B)( Large( sensitivity( in( free( surface( displacement( and( a(
notorious( dependency( on( phase( (arrival( time).( Variability( was(
comparable(to(the(average(signal.(
250
•
0
-2
•
-4
-6
0
50
100
150
Elapsed Time [Min]
200
250
300
•
6
• 3models(fell(into(this(group.(
Case C
4
2
η [m]
• (Case* C)( Reduced( sensitivity( in( free( surface( displacement( but( ( it(
increases( significantly( at( crest( and( troughs.( Variability( was( greatly(
reduced( when( compared( to( other( cases,( showing( a( remarkable(
consistency(in(the(phase.(
0
-2
Né
Ag
-4
-6
0
50
100
150
Elapsed Time [Min]
200
250
300
Sy
Ap
18
Sample*InterGModel*Results
Comparación InterModelos
These(figures(show(the(time(series(of(all(models(for(a(fixed(Manning(roughness(factor(and(comparable(grid(
settings.(
10
8
η [m]
6
4
2
η [m]
• (Coarse* Grid)( A( high( degree( of( variability( is(
observed( both( in( amplitude( and( arrival( times.(
Peak( amplitudes( vary( between( 1C6,( with( large(
variations(in(temporal(response.((
• Overall(structure(of(the(time(series(is(controlled(by(
resonant(conditions(at(bay(level.
0
-2
-4
-6
-8
-10
0
50
100
150
Elapsed Time [Min]
200
250
300
10
• (Finer* Grid)( Most( models( tend( to( yield( a( similar(
temporal( response,( with( comparable( peak((
amplitudes( and( phases,( although( the( location( of(
the(maxima(can(vary(significantly(in(time.((
8
6
4
2
0
-2
-4
• A( couple( of( models( deviate( notoriously,( which( is(
under(further(investigation.
-6
-8
-10
0
50
100
150
200
250
300
Elapsed Time [Min]
19
¿Cómo se Explica Esto?
Algoritmos y
Software
Matemática
Modelado
- Física
20
Implementación Numérica
El desafío es pasar de un esquema continuo (como las derivadas parciales), a un esquema
discreto como lo es una malla de cálculo numérica
Sin perder precisión
Con bajo costo
Sin introducir errores
La ecuación es hiperbólica y resuelta como un problema de condición inicial o de contorno.
Universidad Técnica Federico Santa María
El método típico es de marcha, en el espacio y el tiempo.
FONDEF D11I1119
de Aguas Someras: Métodos de resolución
La forma de resolver la ecuaciónEcuación
parcial en
diferencias finitas puede tener efectos en la
Se resumirá en la siguiente tabla los métodos empleados para ambos la CPU y GPU.
respuesta del problema
Difusión Numérica
CPU
Convección Numérica
Comcot
Tunami
Neowave
Tusunawi
GeoC
Dispersión Numérica
Conservative Conservative
NonNonConse
Formulation
Conservative
Conservative
Method
Finite
Difference
Finite
Difference
Finite
Difference
Finite
Elements
Fi
Diffe
Scheme
Explicit
Explicit
Implicit
Explicit
Exp
Tab. 2: Análisis de los métodos empleados en cada Software.
21
Esquemas numéricos
Tab. 2: Análisis de los métodos empleados en cada Software.
iscusión:
Diferencia Central
Ventajas:
s importante
destacar que el método explícitos de diferencias finitas son preferibles porque sólamente evaluaciones
da paso de tiempo son requeridas. También es fácil de ver que los esquemas de diferencia finitas son mas populares
Es más sencillo de programar y requieren menos tiempo en la computadora por paso de
la evolución en el dominio del tiempo y el espacio. La razón principal se puede atribuir al hecho de que los sistemas
tiempo.
⋄
ferencia finitas
son altamente
paralelizables y pueden ser utilizada con openMP para la CPU y openCL or CUDA
la GPU, (e.g.
Brodtkorb et
(2010,
2012)).que
A demás,
como la disminución
el tiempo
requerido para completar
El esquema
esal.más
preciso
un Central–Upwind
deenprimer
orden.
mulaciónes es uno de nuestro objetivo principales, un esquema de diferencia finitas explícito debería ser el enfoque
Un parámetro de disipación es necesario para acercarse a un estado de equilibrio.
ebe adoptarse.
Permite el uso de múltiples grillas (grillas anidadas) para regiones extensas.
aciónDesventajas:
de Aguas Someras: Esquemas
Es unsemétodo
altamente
inestable.
siguiente tabla
resumen los
esquemas de
solución para la CPU y GPU empleadas tanto para la continuidad y la
ión de momentum. El concepto de esquema guarda relación con la forma en que se realiza la discretización de las
Es algo más disipativo de manera numérica,
lo queelsignifica
queenlaelonda
decae
debido
a
bles y operadores para pasar de una variable continua, comoFig.
por
ejemplo,
nivel
del
mar
caso
de
un
tsunami,
a
19 Two-dimensional structured mesh for finite difference approximations.
la precisión
delel método
y no
sólo debido
procesos
ariable discreta,
como lo sería
nivel del mar
calculado
ahora ena los
puntos defísicos.
la malla de cálculo. La definición de
s fundamental,
pues
distintos
esquemas de para
cálculoobtener
pueden favorecer
distintos aspectos,
como pororden,
ejemplo,
la estabiliTiene
que
ser modificado
las propiedades
de segundo
para
mejorar
umérica y/o la propagación de errores de truncación. Adicionalmente, en la tabla también se incluye algunas
formas
!
la
precisión.
Ui+1, j −Ui−1, j
∂ U !!y ciudades.:
adas paralela aumentar la la precisión conforme se acerca la onda a las zonas de interés como bahías
≡ δ2x
! ≈
∂
x
2
∆
x
Depende en gran medida del ∆t y la discretización espacial ∆x.
i, j
Comcot
CPU
Tunami
Neowave
TsunAWI
GeoClaw
GPU
!
Ui, j+1 −Ui, j−1
∂ U !!
≡ δ2y
≈
!
∂ y i, j
2∆ y
SWE
Leap-frog
CentralDifference
First-order In general,
Leap-frog
First-order
Centralfinite difference
approximation
will involve a stencil of point
Upwind
Upwind
Upwind
Momentum
First-order
Upwind
CentralDifference
First-order
Upwind
Accuracy
NestedGrids
NestedGrids
Continuity
Leap-frog
First-order
Upwind
CentralUpwind
Mesh
Refinement
Adaptative Mesh
Refinement
Variable
Time-step
22
e adoptarse.
Esquemas numéricos
Ux i ≈
ión de Aguas Someras: Esquemas
U(xi + ∆ x) −U(xi − ∆ x) Ui+1 −Ui−1
=
≡δ
2∆ x
2∆ x
The finite difference operator δ2x is called a central difference operator. Finite dif
one-sided.
For example,
a backward
difference
approximation
guiente tabla se resumen los esquemas de solución para la CPU
y GPU
empleadas
tanto para
la continuidad
y lais,
n de momentum. El concepto de esquema guarda relación con la forma en que se realiza la discretización de las
1
(Ui −Ui−1a) ≡ δx−Ui ,
s y operadores para pasar de una variable continua, como por ejemplo, el nivel del mar en el casoUde
un tsunami,
xi ≈
∆x
able discreta, como lo sería el nivel del mar calculado ahora en los puntos de la malla de cálculo. La definición de
and a forward
difference
approximation
undamental, pues distintos esquemas de cálculo pueden favorecer
distintos
aspectos,
como poris,
ejemplo, la estabilimérica y/o la propagación de errores de truncación. Adicionalmente, en la tabla también se incluye algunas
formas
1
(Ui+1 −Ui ) ≡ δx+Ui .
Ux i ≈
as paralela aumentar la la precisión conforme se acerca la onda a las zonas de interés como bahías
y ciudades.:
∆x
Comcot
CPU
Tunami
Neowave
Leap-frog
CentralDifference
First-order
Upwind
Momentum
First-order
Upwind
CentralDifference
First-order
Upwind
Accuracy
NestedGrids
NestedGrids
Continuity
GPU
Exercise 1. Write
a M ATLAB function
which computes the central differen
TsunAWI
GeoClaw
SWE
x1 , x2 , . . . , xNx on the domain [0, 1] with periodic boundary conditions (i.e U0
Leap-frog
First-order
Centralhave the form function
dU = centraldiff(U)
where U = (U1 , U2 ,
(δ2xU1 , δ2xU2 , . . . , δ2xUpwind
UNx −1 , δ2xUNx )T . Upwind
Leap-frog
First-order
Upwind
CentralUpwind
Mesh
Adaptative Mesh
Variable
Refinement
Refinement
Time-step
Exercise 2. Write the
function forwarddiff
which uses a forward difference appro
Tab. 3: Análisis de los esquemas empleados por cada Software.
cusión:
anilo S. Kusanovic, Patricio Catalán
Fig. 19 Two-dimensional structured mesh for finite difference approximations.
Exercise 3. Write the function backwarddiff
which uses a backward difference
!
Ui+1, j −Ui−1, j
∂ U !!
input.
≡ δ2xUi, j .
! ≈
∂x
Tsumami Modelling
2∆ x
i, j
!
Ui, j+1 −Ui, j−1
∂ U !!
≡ δ2yUi, j .
≈
∂ y !i, j
2∆ y
7/11
In general, finite difference approximation will involve a stencil of points surrounding Ui, j .
47.1 Local Truncation Error for a Derivative Approxim
23
¿Por qué la sensibilidad?
FONDEF D11I1119
A comparison between the time series fot the hydrostatic solution in the presented formulation and
by the non-hydrostatic solution is shown in figure(11). It can be seen that the hydrostatic model does
Diferentes
esquemas
tienen
distintos
errores
detravels.
truncado
shape and
amplitude
of the
solitary wave
as it
However, the propagation profile seems to travel
only changing
its cada
shape.paso, un cierto nivel de error se genera que de acumularse, puede
Esto significa
que en
generar deficiencias en la modelación
Ejemplo:
Disipación numérica
(a)
(b)
24
Desafío
Identificar el origen de las discrepancias entre los modelos.
25