From c88661d770e459b3c544b4d7f962e678adab9c52 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Wed, 9 Aug 2017 18:12:55 +0200
Subject: [PATCH] Bugfixes to 2d shallowwater
After a lot of experimenting and stuff...
---
test/hyperbolic/shallowwater.mini | 11 +++++------
test/hyperbolic/shallowwater.ufl | 17 +++++++++--------
2 files changed, 14 insertions(+), 14 deletions(-)
diff --git a/test/hyperbolic/shallowwater.mini b/test/hyperbolic/shallowwater.mini
index 374b2e89..24f4942b 100644
--- a/test/hyperbolic/shallowwater.mini
+++ b/test/hyperbolic/shallowwater.mini
@@ -3,14 +3,13 @@ __exec_suffix = {diff_suffix}
diff_suffix = numdiff, symdiff | expand diff
-extension = 10.0 10.0
-upperright = 10.0 10.0
-cells = 20 20
-elementType = simplical
-elements = 20 20
+extension = 50.0 50.0
+upperright = 50.0 50.0
+cells = 200 200
+elementType = hexahedral
[instat]
-T = 0.5
+T = 10.0
dt = 0.005
[wrapper.vtkcompare]
diff --git a/test/hyperbolic/shallowwater.ufl b/test/hyperbolic/shallowwater.ufl
index 076ab49f..49edd94b 100644
--- a/test/hyperbolic/shallowwater.ufl
+++ b/test/hyperbolic/shallowwater.ufl
@@ -1,7 +1,7 @@
cell = quadrilateral
x = SpatialCoordinate(cell)
-f = 0.2 * exp(-(pow((x[0] - 5), 2) / 2 + pow(x[1] - 5, 2) / 2)) + 1.
+f = conditional(sqrt(pow(x[0] - 25., 2) + pow(x[1] - 25., 2)) <= 11, 2., 1.)
V = FiniteElement("DG", cell, 1)
MV = MixedElement(V, V, V)
@@ -18,18 +18,19 @@ mass = inner(u, v)*dx
h, q0, q1 = split(u)
flux = as_matrix([[q0, q1],
- [h*q0*q0 + 0.5*g*h*h, h*q0*q1],
- [h*q0*q1, h*q1*q1 + 0.5*g*h*h]])
+ [q0*q0/h + 0.5*g*h*h, q0*q1/h],
+ [q0*q1/h, q1*q1/h + 0.5*g*h*h]])
-bflux = as_matrix([[-q0, -q1],
- [h*q0*q0 + 0.5*g*h*h, h*q0*q1],
- [h*q0*q1, h*q1*q1 + 0.5*g*h*h]])
+bflux = as_matrix([[-q0, -q1],
+ [q0*q0/h + 0.5*g*h*h, q0*q1/h],
+ [q0*q1/h, q1*q1/h + 0.5*g*h*h]])
# Define numerical fluxes to choose from
-alpha = Max(abs(n[0]*q0('+') + n[1]*q1('+')) + sqrt(g*h('+')), abs(n[0]*q0('-') + n[1]*q1('-')) + sqrt(g*h('-')))
+alpha = Max(abs(n[0]*q0('+') + n[1]*q1('+')) / h('+') + sqrt(g*h('+')), abs(n[0]*q0('-') + n[1]*q1('-')) / h('-') + sqrt(g*h('-')))
llf_flux = dot(avg(flux), n) - 0.5*alpha*jump(u)
-boundary_flux = 0.5*dot(flux + bflux, n) + as_vector([0., q0, q1])
+alpha_b = abs(n[0]*q0 + n[1]*q1) / h + sqrt(g*h)
+boundary_flux = 0.5*dot(flux + bflux, n) + alpha_b * as_vector([0., q0, q1])
numerical_flux = llf_flux
r = -1. * inner(flux, grad(v))*dx \
--
GitLab