text:
# The two envelopes - Revisited

Palynka
Posers and Puzzles 19 Mar '09 11:18
1. 21 Mar '09 23:22
THE DIFFERENCE
The difference in picking the second envelop to begin with, and switching after we open the first is this.

With the latter we have information we did not have before, specifically a dollar amount which limits our possibilities to two equally likely cases.

PARADOX #1 - Common Sense =VS= The Math
It is counter to our common sense that this would make a difference, but the math works out differently than we naturally expect from a common sense perspective. This is the paradox.

It shouldn't matter but it does, similar in vein to (albeit in a distinctly different manner from) the Monty Hall problem.

PARADOX #2 - Is the amount really important?
Or perhaps the paradox is in the fact that regardless of what we actually find, it is better to switch. In this respect, the fact we know the amount is important, but the actual amount is not.

The difficulty of testing the actual returns via a series of random trials is that the original terms of the trial is unbounded, making it impossible to select the random numbers required to test it.
2. Palynka
Upward Spiral
24 Mar '09 11:43
I'm having trouble replicating my results. 😕

Ok, here's what I'm getting now and you can tell me how to proceed or how to correct what I made wrong.

We know that the amounts x cannot come from a uniform for all possible values of x (no such distribution exists) so let's set assume that the draws come a uniform distribution with an unknown bound of X.

The question is then if we can infer anything about X from x. Using a similar reasoning as in the doomsday argument, we set a prior for X: P(X) = k/X (with k a normalizing parameter). The unconditional distribution of x is then such that P(x) = k/X.

By Bayes' theorem, we get the posterior distribution (note that the k's cancel out as they should):
P(X|x) = P(x|X)*P(X)/P(x) = x/X^2

We know that for any x > X/2 then we must hold the larger amount an we can approximate the probability by a continuous one (X large) and integrate P(X|x) between x and 2x. This gives us a value of exactly one half. P(X<2x) = 1/2. Conversely, P(X>2x) = 1/2

But in the first case (X<2x), we necessarily hold the larger envelope and swapping leads to a loss of -1/2x. In the second case (X>2x), then we might hold the larger or the smaller envelope. A naive interpretation of the last sentence would assume that we could be holding any of the two envelopes with equal probability and therefore the expected value of switching would be:
1/2 * -x/2 + 1/2* (1/2*-x/2 + 1/2*x) = -1/8 x

The expected value of switching would then be negative! The problem is that in the second case the probabilities may depend on X and x. I was able to find the correct expression to the expected value for switching when x<N/2 and it gave me a zero expected value of switching. I can't seem to be able to find it now, so it's probably wrong.

The key intuition here is that if your amount in the original envelope is in the upper half of the distribution, then you cannot double. But if the distribution is uniform over a given interval, then you have 50% chance of being in that upper half!
3. 05 Apr '09 03:57
I've now been thinking about this one for a while, and I think I have a little insight...

In terms of running a test with millions of trials and an infinite ceiling, the reason why one would not be able to see a measured positive difference in the outcome from switching envelopes is the very lack of a "highest" value. since the envelopes themselves could be infinitely large in value, any data from a particular envelope choice is outweighed by an infinite number of infinitely larger cases. and this is true no matter how large the envelope is... there are an infinite number of cases that outweigh the current choice to the point where by magnitude, an expectation of 5/4x vs x is meaningless and immeasurable. and furthermore, since there is no known "highest value" envelope, there is no single exception case where the player would be making a bad play by switching envelopes... as you say it SEEMS paradoxical.

it has also been said that the apparent paradox is that if one took the strategy of "picking" an envelope, and switching before it is opened, and THEN opening the envelope, that the person should get the same results as the person switching after seeing the envelope. yet the finite answer says that no matter what value you see you should switch.

i contend that it is the ACT of observation of the envelope that defines it at a finite value. akin to the box containing "schroedinger's cat", the envelope you open has undefined value until you actually open it to see how much money is in it. until opening it has equal probability of being any value in the list of possible envelope values. but then AFTER opening the envelope, it is defined as being valued "x" and we are given a choice of choosing between the possibility of the other envelope being valued "2x" or "1/2x." and it turns out that for any FINITE value x, and with no promise of there being a "highest" valued envelope, there is more to gain from switching then there is to lose [with an expected value of (5/4)x].

thoughts?
4. 05 Apr '09 04:12
Originally posted by Palynka
I'm having trouble replicating my results. 😕

Ok, here's what I'm getting now and you can tell me how to proceed or how to correct what I made wrong.

We know that the amounts x cannot come from a uniform for all possible values of x (no such distribution exists) so let's set assume that the draws come a uniform distribution with an unknown bound of X. ...[text shortened]... is uniform over a given interval, then you have 50% chance of being in that upper half!
i haven't really gone through your calculations for accuracy, but i think i understand your basic reasoning... however, i think most of it is based on an assumption i don't necessarily agree we can make.

essentially, you are assuming there is a highest possible value M that the envelope could be, but we do not know how large that is. in fact, this assumption is exactly against the problem as proposed, and if it were true, of course there would be an error in the (6/5)x expected value interpretation i've been discussing: that if the envelope in question is the "highest possible" value, the player would be a fool to switch envelopes.

then you propose (and please correct me if i'm wrong) that based on the assumption there is a largest possible value, but we don't know how large it is, we should make a guess as to whether our value is in the top half of a normal distribution that has a ceiling of M. this would be an educated guess of course, a probabilistic calculation that is analogous to the question "how likely is it that the host of this game show has so many envelopes that the value i see is one of the 'smaller' envelopes rather than one of the 'larger' ones."

i don't know the "doomsday argument" but it seems to me that these assumptions are largely flawed considering that the posed problem gives us no indication of the magnitude of the highest possible value in an envelope, and in fact explicitly states that the values are drawn from a list that tends towards infinity. even if there WERE a ceiling value, our guess at its magnitude would be like throwing darts in the wind, in my opinion.
5. Palynka
Upward Spiral
05 Apr '09 10:20
Thanks for the comments Aethereal, I'll try to go by them when I have more time , I just managed to skim them for the moment.

Just so I know how to frame this, are you familiar (very, loosely, not at all) with the use of priors in probability?
6. 05 Apr '09 20:24
Originally posted by Palynka
I'm having trouble replicating my results. 😕

Ok, here's what I'm getting now and you can tell me how to proceed or how to correct what I made wrong.

We know that the amounts x cannot come from a uniform for all possible values of x (no such distribution exists) so let's set assume that the draws come a uniform distribution with an unknown bound of X. ...[text shortened]... is uniform over a given interval, then you have 50% chance of being in that upper half!
i've been examining this post and the info in the "doomsday argument" wikipedia page to help me better understand the math behind what you're doing here...

but i think there's a logical error that is somewhat important: you discuss two cases, one in which X (the upper bound on x) is > 2x and one in which X < 2x and say that they have equal probability.

however, the envelopes must always be double one another. so if X < 2x, the only possible case is that we are holding the biggest envelope in the set... but every envelope that is NOT the biggest envelope in the set contains an x such that X > 2x. so how can these cases have equal probability? in fact the larger the set gets, the less likely it is that we are holding the largest envelope in the set (the case with a (-1/2)x expected value) and the more likely we are holding an envelope that is not the largest in the set (the case with a [1/2*(-1/2*x) + 1/2*(x)] = +1/4*(x) expected value).

and i think that as the set of possible envelopes tends towards infinity, the probability that we are holding the very largest envelope in the set tends towards zero, making the probability that we are NOT holding the very largest envelope tend toward 1.
7. 06 Apr '09 00:42
Originally posted by Palynka
Thanks for the comments Aethereal, I'll try to go by them when I have more time , I just managed to skim them for the moment.

Just so I know how to frame this, are you familiar (very, loosely, not at all) with the use of priors in probability?
i would say i'm between very and loosely familiar with the use of priors... what i am very good well versed in is extrapolation and inductive reasoning. i.e. i pick things up usually the first time through an explanation, if it's clearly explained, and am very fluent in the mechanics of algebra and calculus.

so if it seems like i don't get the way in which you're framing the problem, a brief rundown of intent or a brief mathematical step-by-step walkthrough is usually enough to get me on the same page (for example, the maths you used in a previous post defining a prior and analyzing posterior was fine, and i agreed with the method except for the key point i outlined in my previous post)

cheers!
8. 06 Apr '09 10:27
I don't really see the problem

I read somewhere in this thread you have 50% chance of gaining 100% and 50% chance of loosing 50%.
( 200% + 50% ) / 2 = 125% per switch

So after 1.000.000 switches, you have X * 1,25 ^ 1.000.000 dollar?
I think it doesn't matter what you choose.
9. Palynka
Upward Spiral
06 Apr '09 11:12
Originally posted by Aetherael
i contend that it is the ACT of observation of the envelope that defines it at a finite value. akin to the box containing "schroedinger's cat", the envelope you open has undefined value until you actually open it to see how much money is in it. until opening it has equal probability of being any value in the list of possible envelope values. but then AF ...[text shortened]... more to gain from switching then there is to lose [with an expected value of (5/4)x].
This is exactly why I think the use of priors must be the correct way to formulate this. The act of opening the envelope updates the prior and so such a way would be consistent with indifference before opening but not after opening. The opening of the envelope reveals information.

The problem I have with your final remarks in this post is it uses information from an impossible posterior (unbounded uniform). Although it is possible for that ex-post the probabilities are 1/2, this cannot be because the values came from an unbounded uniform!

Another beauty of using priors is that we can use an unbounded uniform as a prior (reflecting our complete lack of information beyond the fact that the amounts are positive), which will imply a well-defined posterior after the envelope is opened.
10. Palynka
Upward Spiral
06 Apr '09 11:22
Originally posted by Aetherael
i haven't really gone through your calculations for accuracy, but i think i understand your basic reasoning... however, i think most of it is based on an assumption i don't necessarily agree we can make.

essentially, you are assuming there is a highest possible value [b]M
that the envelope could be, but we do not know how large that is. in fact, th ...[text shortened]... guess at its magnitude would be like throwing darts in the wind, in my opinion.[/b]
Well, you can call it an assumption, but for me it is a logical necessity. The fact that unbounded uniform distributions do not exist leaves us with only two choices, if we don't want simply to answer: "the problem is incorrectly posed". These are:

1) The distribution is not uniform.
2) The distribution is not unbounded.

Personally I think that weakening 2) is more true to the spirit of the question. Of course, arbitrarily fixing the limit would also not be true to the spirit of the question. So what I try to do is to impose as little as possible by not excluding ANY possible value ex-ante by using an unbounded prior.

Like this, all values are possible ex-ante and the "true" distribution is still uniform.
11. Palynka
Upward Spiral
06 Apr '09 11:322 edits
Originally posted by Aetherael
i've been examining this post and the info in the "doomsday argument" wikipedia page to help me better understand the math behind what you're doing here...

but i think there's a logical error that is somewhat important: you discuss two cases, one in which X (the upper bound on x) is > 2x and one in which X < 2x and say that they have equal probabilit ...[text shortened]... ro, making the probability that we are NOT holding the very largest envelope tend toward 1.
however, the envelopes must always be double one another. so if X < 2x, the only possible case is that we are holding the biggest envelope in the set... but every envelope that is NOT the biggest envelope in the set contains an x such that X > 2x. so how can these cases have equal probability? in fact the larger the set gets, the less likely it is that we are holding the largest envelope in the set (the case with a (-1/2)x expected value) and the more likely we are holding an envelope that is not the largest in the set (the case with a [1/2*(-1/2*x) + 1/2*(x)] = +1/4*(x) expected value).

I'm not sure I understand... I got the value of 1/2 through integrating: [;\int_n^{2n} P(N|n)dN ;].

where n is the value in the envelope and N is the upper bound. I integrate over the upper bound.

In fact, getting a "low" value of n is more likely when N is "low", so I don't think one can make an unconditionally claim. This is why I integrate over P(N|n).
12. Palynka
Upward Spiral
06 Apr '09 11:54
But now you're making me think that there might be a way to apply the prior directly to n...

I think my mistake is that when I'm in the P(N>2n), I assume equal likelihood of holding each envelope. Conditionally, this might not be true (and my intuition is that it isn't).
13. 06 Apr '09 14:081 edit
Originally posted by Palynka
But now you're making me think that there might be a way to apply the prior directly to n...

I think my mistake is that when I'm in the P(N>2n), I assume equal likelihood of holding each envelope. Conditionally, this might not be true (and my intuition is that it isn't).
yeah this is the part where i was straying from your argument - the integration was fine, and the splitting into two cases to examine P(N>2n) and P(N<2n)

but my point is that as the set of possible values in the envelope grows, P(N<2n) shrinks. Think of what the terms N and n mean in the setup of the problem we've been dealing with... we are trying to determine what the probability that the upper bound is less than double the value in our envelope. but this is ONLY true if we have chosen the largest possible envelope of all allowed in the set; that is, that we have chosen the envelope that sets the upper bound.

just to use a concrete example... imagine the bounded set contains envelopes {1,2,4,8,16,32,64}. the player doesn't know of course that 64 is the ceiling value, but you can see that when N=64, the only time N<2n is when n=N, or rather when we are looking at the largest envelope. in this concrete example, we see that the probability we chose n=N is 1/7... as the number of envelopes in the set grows this probability will tend toward 0, and the probability that we did NOT choose the biggest envelope will tend toward 1.

all other envelopes contribute to the P(N>2n), which has a positive expectation of +(1/4)n from the math in your earlier post.

i guess what i mean is, though we are using a uniform distribution for the possible values of n, that these possible values are discrete numbers, and are ONLY the powers of 2... so to think that half of them would be larger than N/2 and half of them smaller than N/2 is folly. in fact only one value in the distribution will be ever be larger than N/2 and that envelope will be the upper bound N.
14. Palynka
Upward Spiral
06 Apr '09 14:521 edit
Originally posted by Aetherael
yeah this is the part where i was straying from your argument - the integration was fine, and the splitting into two cases to examine P(N>2n) and P(N<2n)

but my point is that as the set of possible values in the envelope grows, P(N<2n) shrinks. Think of what the terms N and n mean in the setup of the problem we've been dealing with... we are trying to the distribution will be ever be larger than N/2 and that envelope will be the upper bound N.
but my point is that as the set of possible values in the envelope grows, P(N<2n) shrinks.

The set of possible values is not bounded, so I don't know what you mean by that. My prior already includes all possible values of N up to infinity. And I'm integrating over that support (up to infinity). Think about it. If the integral of n up to 2n is 1/2, then the integral of 2n to 3n is smaller than 1/2. So your idea that P(N<2n) "shrinks" when N increases* is already included.

*I'm interpreting your statement like that, but note that "possible values" are up to infinity. There is no value for which I impose a zero-probability of happening.

{1,2,4,8,16,32,64}
Ha! Ok, I may have misled you with my initial posts. In my replicating calculations I was thinking of the set of possible values being the reals, not just a sequence of the powers of 2. This was where we were miscommunicating (my bad, sorry...)

i guess what i mean is, though we are using a uniform distribution for the possible values of n, that these possible values are discrete numbers, and are ONLY the powers of 2... so to think that half of them would be larger than N/2 and half of them smaller than N/2 is folly.
Yes, absolutely correct. See above.

Somewhere along the way I realized exactly what you're saying now and dropped the powers of 2 assumption. It is too restricting and I suspect solves the problem exactly in the way you also suspect. Before opening the envelope, switching is infinity-infinity (indeterminate) and after opening the envelope, switching is optimal. The paradox disappears because the indeterminacy solves the conundrum.

It doesn't seem very satisfying, right? Well, to me it isn't. I think it imposes too much which isn't in the question (the powers of 2), along with the unknown interval size. So I'm now thinking how to work it out in a more general setting...
15. 06 Apr '09 19:311 edit
Originally posted by Palynka
[b]but my point is that as the set of possible values in the envelope grows, P(N<2n) shrinks.

The set of possible values is not bounded, so I don't know what you mean by that. My prior already includes all possible values of N up to infinity. And I'm integrating over that support (up to infinity). Think about it. If the integral of n up to 2n is 1/2, n interval size. So I'm now thinking how to work it out in a more general setting...[/b]
ok now i see where we separated in terms of my understanding of your argument. now it seems ok, except i will ask you to consider this:

1) wouldn't setting the problem over the reals contain an infinitude of irrational numbers that, in the spirit of the problem, wouldn't really work as a "dollar amount" in the envelope? since they are uniformly distributed throughout the real line, this wouldn't hurt our ability to uniformly distribute the RATIONALS, and so i can see an argument for using this as a setting for the problem. however,

2) using the rationals requires an infinitesimally accurate system of "value," that is to say the number you see in the envelope could be as tiny as .0000....0001, which i think again is not within the spirit of the problem. also the rationals include things like 5/7 dollars and 41/3 dollars, which it seems to me goes against the spirit of the problem as a dollar amount. if it were posed as "you open an envelope and there is a number in it... and are told it is either exactly half or exactly double the number in the other envelope" etc. then i could see this as plausible. but i think these restrictions are not flippant, but really within the spirit of the problem, and as such

3) i believe the most "natural" setting (har har) for the problem is over the natural numbers. but then since we get into the issue that any non-power of 2 would give you too much information about the other envelope, so we restrict our scope to {1, 2, 4, 8, 16, ...}

perhaps these restrictions are too strong, but i personally think it is more in line with the spirit of the problem. that being said, infinitesimally accurate decimals and a placement over the rationals is entirely arguable, and i think that all the maths you've done so far would work over the rationals...i just think it doesn't seem right to tell the player "in your envelope is 147 and 93756/142019274 dollars... do you think there is double that or half that in the other envelope?" 🙂 just being cheeky.


text:
# Are there full cycle cryptographic/one-way hash primitives?

I'm looking for behavior similiar to that of LCGs, (i.e. input and output sizes are same). Full cycle of $2^{32}$ different inputs generates full cycle of $2^{32}$ different outputs, distribution of which depends on large key. Obviously combined LCGs (with generated parameters based on the perturb seed key) seem trivially invertible through modular inverse, but you get the idea.

The hash has to accept small input and produce never-colliding equal sized output, and is keyed with large 'salt'/'nonce' key, to make precomputing generic table of all $2^{32}$ hashes for trivial inversion infeasible. The hash has to remain one way when 'seed' key and output is known (ie input should be difficult to compute out of those, short of bruteforcing 2^31 inputs for each particular seed+output and meeting in the middle for birthday).

This exotic requirement sets it apart from classical bijective trapdoor functions (RSA, MQ et al), as the output has to be much smaller than some input parameters (ie the "salt" seed).

Is this even feasible? It does not appear to be with traditional ARX/feistel, however I suspect there are still useable trapdoors in number theory - for example FSB and SWIFFTX. However I'm not sure if those can be adapted to the requirement above.

I'm OK with security parameter well below nbits/2, as long it is still reasonably difficult to invert algebraically.

Bottom line: Albeit related, perfect hash functions are not acceptable answer as there does not seem to be examples of cryptographically secure ones, much less any working as a trapdoor.

Block ciphers such as RC6 or AES rounds show some promise (by using the 'decrypt' part as our hash), but naive implementation seems to be reversible by running the rounds backwards (remember, the cipher key, our "seed" sbox which introduces non-linearity in those is publicly known, whereas the tiny "plaintext" is what is to be protected).

I presume block ciphers which break linearity through heavy input data dependency might work, however I have no idea what candidates I should be looking for.

-
$2^{32}$ brute force can be done in seconds (on a GPU) unless the 'hash' is very slow, so avoiding attacks faster than brute force does not seem to gain you anything much. You could just use format-preserving encryption and accept the fact that knowing the key allows decryption. – otus Oct 10 at 4:14
I'm well aware. To avoid further discussions on the example security parameter, imagine 512 everywhere in the text instead of 32. – kat Oct 10 at 4:29
If you change 32->512 the answer changes to: just use a collision resistant 512-bit hash. – otus Oct 10 at 4:35
@otus Which hash function do you recommend - it has to be still full cycle (bijective). I do realize universe would end before enumerating all states with 512 bits, however my application is not to avert bruteforce through overkill like this - that's what the seed is for - it's an interactive protocol with constrained bandwidth and short lifetimes of this weird MAC I am trying to do. – kat Oct 10 at 4:41
You should probably leave this question as is and ask another question about the specific problem you are trying to solve. It is not at all clear to me that you need a permutation. – otus Oct 10 at 5:08

I'm sorry, but I don't believe that there are any known 'one way permutations' for small inputs (that stay one way if you know the 'key').

When we create a cryptographic permutation (or bijection; that's what we call a 'full cycle' function), we have two basic ways to do it:

• We can take a series of easy to invert permutations ("round functions"), and glue them together to form a function that we hope is hard to invert.

This is the general design principle behind most block ciphers; we believe that if the key is secret, this can work very well. However, if you give someone the key, it becomes easy to invert (by just inverting each individual round function).

• We can take advantage of some mathematical relation which is one-to-one (that is, it avoids collisions), but it doesn't do it for any naively obvious reasons.

We often find these sorts of things within public key cryptography; one such is $f(x) = g^x \bmod p$ (for appropriate $g, p$); this function can map the inputs $(1, p-1)$ into the range $(1, p-1)$ such that every output is possible.

One problem with this is that these types of functions are large; in the specific case I mentioned, we need $p$ to be perhaps 2048 bits to be really secure, and that's hard to describe as a small range.

-
RSA is probably what I'm going to end up with (though not happy about it). On a related note, in the event g happens to be a weak composite (can't afford primality tests) and p is always hardcoded value, wouldn't there be a shortcut to derive x for an adversary? – kat Oct 10 at 4:57
Actually, the function $f(x) = g^x$ that I mentioned is not RSA; instead, it's modular exponentiation. As for $p$ and $g$, you could select $p$ to be the group14 prime from datatracker.ietf.org/doc/rfc3526, and $g = p - key^2$; as long as $\sqrt{p} > key > 1$, that'll always be a strong value for which $f(x)$ is a bijection.. – poncho Oct 10 at 11:14
My apologies for incorrect lingo (i refer to it as RSA as RSA simply exploits similiar congruence). My configuration at the moment is (in C-ish notation): G = primitive_roots_of_p[seed1]^seed2; p = 2^31-1; x = text to be hashed; Unfortunately I can't afford to square key, as that would result having "seed" even smaller than it already is. – kat Oct 10 at 11:17
If $p$ is the 31 bit prime you appear to mention, well, that's quite weak (that is, easy to invert). I suspect that otus is right; you need to rethink the problem. Why do you think you need a public permutation anyways? What problem are you actually trying to solve? – poncho Oct 10 at 11:32
It is the only one which appears to have bijective mapping with "hashing" property. Ie forward step is one powmul, the inverse is 2^30 powmuls. I'd rather avoid it if I could... As for the exact protocol, it's byzantine fault tolerance/quorum consensus, where client has to prove ownership of data it previously said it owns. The idea is to prevent clients from announcing us data they can can easily generate at runtime. This trapdoor forces client to not lie to us, as it would turn into fairly difficult PoW problem and would have to solve 2^30 steps for each 4 bytes it didnt store. – kat Oct 10 at 11:38

format: plain


text:
# Multigrid preconditioners for the hybridized Discontinuous Galerkin discretisation of the shallow water equations

Numerical climate- and weather-prediction models require the fast solution of the equations of fluid dynamics. Discontinuous Galerkin (DG) discretisations have several advantageous properties. They can be used for arbitrary domains and support a structured data layout, which is particularly important on modern chip architectures. For smooth solutions, higher order approximations can be particularly efficient since errors decrease exponentially in the polynomial degree. Due to the wide separation of timescales in atmospheric dynamics, semi-implicit time integrators are highly efficient, since the implicit treatment of fast waves avoids tight constraints on the time step size, and can therefore improve overall efficiency. However, if implicit-explicit (IMEX) integrators are used, a large linear system of equations has to be solved in every time step. A particular problem for DG discretisations of velocity-pressure systems is that the normal Schur-complement reduction to an elliptic system for the pressure is not possible since the numerical fluxes introduce artificial diffusion terms. For the shallow water equations, which form an important model system, hybridized DG methods have been shown to overcome this issue. However, no attention has been paid to the efficient solution of the resulting linear system of equations. In this paper we show that the elliptic system for the flux unknowns can be solved efficiently by using a non-nested multigrid algorithm. The method is implemented in the Firedrake library and we demonstrate the superior performance of the algorithm for a representative model problem.

## Authors

• 2 publications
• 5 publications
• 9 publications
• 2 publications
• ### IMEX HDG-DG: a coupled implicit hybridized discontinuous Galerkin (HDG) and explicit discontinuous Galerkin (DG) approach for shallow water systems

We propose IMEX HDG-DG schemes for planar and spherical shallow water sy...
11/07/2017 ∙ by Shinhoo Kang, et al. ∙ 0

• ### A Semi-Implicit Meshless Method for Incompressible Flows in Complex Geometries

We present an exponentially convergent semi-implicit meshless algorithm ...
05/29/2021 ∙ by Shantanu Shahane, et al. ∙ 0

• ### High-order numerical solutions to the shallow-water equations on the rotated cubed-sphere grid

High-order numerical methods are applied to the shallow-water equations ...
01/14/2021 ∙ by Stéphane Gaudreault, et al. ∙ 0

• ### Exponential Integrators with Parallel-in-Time Rational Approximations for the Shallow-Water Equations on the Rotating Sphere

High-performance computing trends towards many-core systems are expected...
05/16/2018 ∙ by Martin Schreiber, et al. ∙ 0

• ### Exponential Integrators with Parallel-in-Time Rational Approximations for Climate and Weather Simulations

High-performance computing trends towards many-core systems are expected...
05/16/2018 ∙ by Martin Schreiber, et al. ∙ 0

• ### Local time stepping for the shallow water equations in MPAS-Ocean

We assess the performance of a set of local time-stepping schemes for th...
06/14/2021 ∙ by Giacomo Capodaglio, et al. ∙ 0

• ### Mixed-precision for Linear Solvers in Global Geophysical Flows

Semi-implicit time-stepping schemes for atmosphere and ocean models requ...
03/30/2021 ∙ by Jan Ackmann, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The dynamics of geophysical fluids is characterised by phenomena which occur at very different time scales. While large-scale flow in the atmosphere is described by the slow evolution of quasi-geostrophic modes, fast acoustic waves play an important rôle in incorporating non-linear dynamics and rapid adjustment to geostrophic balance. However, acoustic waves carry very little energy and hence do not have to be represented to high accuracy in computer simulations. Efficient numerical models for weather- and climate-prediction often exploit this scale-separation by using semi-implicit time integrators. Treating the fast waves implicitly (and thereby removing any instabilities due to acoustic modes) avoids tight restrictions on the time step size. It is sufficient to chose the timestep such that it allows the stable and accurate representation of the slow modes, which are treated explicitly. Examples of popular time-integrators are semi-implicit semi-Lagrangian methods used by the UK Met Office [1] and implicit-explicit (IMEX) methods [2, 3, 4, 5] (see also [6] for applications in atmospheric modelling).

In a computer model the equations of fluid-dynamics have to be represented on a grid with finite resolution. While simple finite-difference- or finite-volumes schemes on structured meshes (such as latitude-longitude grids) work well for local area models, they suffer from the convergence of grid lines at the poles on global grids [7]. This leads to artificially tight constraints on the timestep size due to the CFL condition because of the small grid cells at the poles. Furthermore, all-to-all communications limit scalability in massively parallel implementations. Finite element methods on quasi-uniform unstructured global grids overcome this problem, but care has to be taken when solving the resulting linear systems which arise in semi-implicit timestepping. The standard approach uses a Schur-complement reduction to a pressure system. Naively this results in a dense Schur-complement operator in pressure space. For conforming finite element methods this can be traced back to the fact that the mass-matrices are non-diagonal. In this case the problem is easily overcome by under-integrating and co-locating the quadrature points with nodal basis points, or forming an approximate Schur-complement [8].

Discontinuous Galerkin (DG) discretisations have recently received significant attention in atmospheric modelling [9]

. DG methods have several advantages, especially on modern chip architectures: since the unknowns are associated with cells of the mesh, data access can be coalesced, which is particularly important for chips with wide vector units. For sufficiently smooth solutions, higher order methods are numerically efficient since errors decrease exponentially as the polynomial degree of the DG approximation increases. Matrix-free implementations are particularly efficient on chips with limited memory bandwidth since they have a very high computational intensity; sum factorisation techniques can further reduce the computational complexity

[10]. Unfortunately, for DG methods the numerical flux introduces off-diagonal matrix elements which cause issues for semi-implicit discretisations. The flux generates artificial diffusion terms, which are numerical artifacts and can not be removed in an obvious way. This rules out the naive Schur-complement approach for solving the implicit linear system.

In [11] a solution to this problem is presented for the shallow water equations (SWEs), a popular two-dimensional model system for atmospheric dynamics. The authors use the hybridized Discontinuous Galerkin (HDG) method (see e.g. [12, 13, 14, 15]), which is adapted for the SWEs in [16, 17]. Introducing a set of flux-unknowns on the facets of the mesh allows the reduction to a sparse Schur-complement system in the flux space. In [11] this problem is assembled and solved with a direct method via LU factorisation. While this works for small problems, the cost of this direct solver grows rapidly as the problem size increases, even if regularities in the sparsity pattern are exploited (see [18] for a general overview of direct methods). This approach is therefore clearly unfeasible for large scale applications, especially if they are eventually to be used in time critical operational forecasting systems.

In this paper we present an alternative solution which overcomes those issues and we demonstrate the computational performance of our approach. We introduce a new non-nested multigrid method for solving the linear systems arising in the IMEX approach for the SWEs, when discretised using a hybridized (high order) DG method. The method is strongly influenced by the work of Gopalakrishnan and Tan [19] and Cockburn, Dubois, Gopalakrishnan and Tan [20]. In particular [19] introduces a non-nested multigrid V-cycle for the hybridised mixed method for standard second order elliptic PDEs discretised using Raviart Thomas elements for the velocity approximation. This has been extended in [20] to the multigrid solution of the HDG discretisation for the same elliptic PDE. In this paper we adapt this approach to solve the linear and non-linear shallow water equations. The non-linear SWEs present additional challenges since they require the construction of a Lax-Friedrich flux, which results in a vector-valued flux space [17] for which the method in [20] is not directly applicable.

More specifically, the new contributions of this paper are:

1. The non-nested multigrid method of [20] is extended to the HDG solution of the linearised SWEs with IMEX methods.

2. To solve the non-linear SWEs by this approach we construct novel intergrid- and coarse-level operators for vector-valued flux spaces.

3. Our HDG solvers and multigrid algorithms are implemented in the Slate language [21]

, which forms part of the open source Firedrake library

[22] for the solution of finite element problems based on composable abstractions.

4. Based on this, the superior performance of the non-nested multigrid approach for both the linear- and non-linear SWEs is demonstrated for a parallel implementation.

Our approach differs from the more algebraic method presented in [23, 24], which constructs a multigrid hierarchy directly on the trace-space defined on the skeleton of the mesh; this is achieved by using -refinement and statically condensing out unknowns on fine edges which are not represented on the coarser mesh.

##### Structure.

This paper is organised as follows: In Section 2 we briefly review the physics of the shallow water equations. The hybridized DG discretisation is described in Section 3 and the non-nested multigrid preconditioner for the resulting linear systems in IMEX timestepping is presented in Section 4. The specific integrators considered in this paper are introduced in Section 5 with focus on the linear solvers used for implicit methods. The implementation in the Firedrake library is described in Section 6 and results are presented in Section 7. We finally conclude and outline directions for future work in Section 8. Some more technical issues such as the explicit expressions for key bilinear forms are relegated to the appendices.

## 2 Shallow water equations

The shallow water equations (SWEs) describe the evolution of a layer of fluid under gravitational and inertial forces. In a rotating frame a Coriolis term has to be included as well. The SWEs can be written in conservative form for the geopotential height and the column-integrated momentum , which are defined for all points in a bounded domain . The time evolution of and is described by the following two equations:

 ∂ϕ∂t+cg∇⋅u =0(mass conservation) (1) ∂u∂t+cg∇⋅(u⊗uϕ+ϕB+12(ϕB+ϕ)2) =R:=cg(ϕB+ϕ)∇ϕB−fu⊥(momentum conservation)

Suitable boundary conditions have to be included at the boundary of . The geopotential height is defined as where denotes the perturbation of the fluid relative to its state at rest and is the gravitational acceleration. Furthermore where is the bathymetry of the sea-floor, measured from the surface at rest (see Fig. 1). The momentum is related to the column-averaged velocity of the fluid at each point by . is the latitude-dependent Coriolis parameter, and the vector is obtained by rotating by 90 degrees in the counter-clockwise direction.

Eq. (1) has been non-dimensionalised as follows: let be a reference height, and define . Further assume that any horizontal distances are measured in units of some reference length scale and times in units of a reference time scale . , are given in units of and the momentum in units of . The dimensionless reference gravity wave speed which multiplies the flux terms in Eq. (1) is then given by

 cg=¯¯¯¯T¯¯¯¯L¯¯¯ϕ1/2B.

For convenience, we introduce the variable , which allows us to write Eq. (1) in compact form as

 ∂q∂t+∇⋅F(q)=s(q),or % equivalently∂qi∂t+∂jFij(q)=si(q). (2)

Here and in the following summation over repeated indices is implicitly assumed. The exact form of and can be deduced from Eq. (1) and is given by (see [11])

 F(q) =⎛⎜ ⎜ ⎜ ⎜⎝uvu2ϕ+ϕB+ϕBϕ+12ϕ2uvϕ+ϕBuvϕ+ϕBv2ϕ+ϕB+ϕBϕ+12ϕ2⎞⎟ ⎟ ⎟ ⎟⎠, s(q) =⎛⎜⎝0cgϕ∂1ϕB+fvcgϕ∂2ϕB−fu⎞⎟⎠.

The SWEs can be linearised about (this is a good approximation if at every point in the domain) to obtain

 ∂ϕ∂t+cg∇⋅u =0, ∂u∂t+cgϕB∇ϕ =−fu⊥. (3)

It is important to note that in a non-rotating frame () fast gravity waves are supported by Eq. (3). This can be seen by taking the divergence of the second equation and inserting it into the time derivative of the first, to obtain the wave equation

 ∂2ϕ∂t2−c2g∇⋅(ϕB∇ϕ)=0.

The (local) speed of the gravity waves is given by . More generally, in a rotating frame the fundamental solutions of the linear SWEs in Eq. (3) are inertia-gravity (Poincaré) waves. In this case the frequency is related to the wave vector by the dispersion relation . In analogy to Eq. (2), the linearised SWEs in Eq. (3) can be written as

 ∂q∂t+∇⋅FL(q)=s(q)⇔∂qi∂t+∂j(FL)ij(q)=si(q) (4)

with and the linear operator given by

 FL(q)=cg⎛⎜⎝uvϕBϕ00ϕBϕ⎞⎟⎠.

## 3 IMEX HDG methods for the shallow water equations

### 3.1 Spatial discretisation

To discretise the SWEs with the discontinuous Galerkin (DG) method, we first construct a mesh of non-overlapping cells that partition the domain such that . The mesh spacing is defined as . Let be the set of facets of the mesh; is often called the “skeleton of ”. We further assume that the mesh consists of triangular or quadrilateral cells and has no hanging nodes. More specifically, each facet is a straight line which forms the side of exactly two triangles or quadrilaterals. For simplicity we use periodic boundary conditions in this paper, which implies that all facets are internal.

Since functions in the DG space which we introduce below are not continuous across cell-boundaries, care has to be taken when defining quantities on facets. For two neighbouring cells and that share a facet , let be the outward unit normal on the boundary of element , and the corresponding outward normal of cell (see Fig. 2), such that on the shared facet . On each facet , the following average , difference and jump operators can be defined for scalar () and vector valued () quantities:

 {{s}} =12(s++s−), {{a}} =12(a++a−), ⟨⟨s⟩⟩ =s+−s−, ⟨⟨a⟩⟩ =a+−a−, [[s]] =s+n++s−n−, [[a]] =a+⋅n++a−⋅n−.

At any point on the facet and for any scalar- or vector-valued piecewise polynomial functions , we define and . Note that the jump of a scalar quantity is vector-valued and vice versa.

Write for the space of polynomials of degree at most on an arbitrary domain . The (scalar-valued) DG space of order on the mesh is defined as

 Vh:={ψ∈L2(Ωh):ψ|K∈Pp(K),∀K∈Ωh}. (5)
##### Case 1: Linear SWEs.

It is instructive to first discuss the DG discretisation of the linear SWEs in Eq. (4). This can be found by assuming that111In a slight abuse of notation, we use the same variable names for the continuous fields , , , etc. in Section 2 and for the DG fields discussed from this section onwards. , then multiplying by a test function and integrating by parts in each cell to obtain

 (∂qi∂t,vi)K−((FL)ij(q),∂jvi)K+⟨(F∗L)ij(q+,q−)nj,vi⟩∂K=(si(q),vi)Kfor % all v∈Wh. (6)

Here is the outward normal on each facet of cell , with . Round and angular brackets denote scalar products in function spaces defined on the cell and its boundary respectively:

 (f,g)K =(f,g)L2(K)=∫Kf(x)g(x)dx, ⟨f,g⟩∂K =(f,g)L2(∂K)=∫∂Kf(x)g(x)ds.

The numerical flux appearing in Eq. (6) (still to be defined) will depend on the values , of the field on both sides of the facets comprising and will couple the unknowns between neighbouring cells. For the linear problem both the local Lax-Friedrichs (Rusanov) flux [25] and the upwind flux (see e.g. [26]) are candidates for the numerical flux . As shown in [16, 17], those are defined as

 (F∗(LF)L)ij(q+,q−)nj ={{(FL)ij(q)}}nj+cg2√ϕB⟨⟨qi⟩⟩ (7a) (F∗(up)L)ij(q+,q−)nj ={{(FL)ij(q)}}nj+cg2√ϕBBik⟨⟨qk⟩⟩ (7b)

where the matrix defined on is given by

 B=⎛⎜⎝1000n2xnxny0nxnyn2y⎞⎟⎠.

If Eq. (6) is integrated in time with a semi-implicit method, a coupled sparse linear system for the fields has to be solved in each timestep. The standard approach in many atmospheric models is to eliminate the momentum unknowns from the discrete equations and solve the resulting elliptic Schur-complement system for the height perturbation , followed by reconstruction of the momentum field (see e.g. [1] for the equivalent procedure applied to the full Navier Stokes equations). Unfortunately, for the DG discretisation considered here, the numerical flux in Eq. (6) introduces artificial diffusion terms in the momentum equation. Although those terms disappear in the limit , they result in a dense Schur-complement operator which renders this naive Schur-complement approach unfeasible here. It is possible to address this issue by forming an approximate Schur-complement and using the resulting approximate solver as a preconditioner for a Krylov subspace method. A similar approach is pursued to deal with the non-diagonal velocity mass matrix in conforming mixed finite element discretisations in [8, 27]. However, this requires additional iterations over the large -system. As our numerical results in Section 7 demonstrate, this approximate Schur-complement method (which essentially drops the artificial diffusion terms in the preconditioner) is not competitive with the hybridized DG method pursued in this paper since the latter allows the formulation of an exact Schur-complement in a suitably defined finite element space on the facets.

### 3.2 Hybridised DG discretisation

The key idea of the hybridized DG method described in [16, 17] is to introduce a new variable , which is associated with the skeleton of the mesh, see Fig. 2 (right). As we will show below, it is possible to form an exact Schur-complement in . For this, recall that the facets are straight lines bounding the cells. On these we define the piecewise polynomial space

 Λh ={λ∈L2(Eh):λ|e∈Pp(e),∀e∈Eh} and setˆΛh =Λh×Λh×Λh. (8)

This allows the introduction of a flux which only depends on the value inside the cell and the field on its boundary. To see this, recall that for hyperbolic conservation laws the numerical upwind flux can be constructed by solving a Riemann problem on the interface between two cells (see e.g. section 2.4 of [28]). The solution contains both the field inside the cell and an intermediate state , which is usually eliminated with a jump condition. Since this jump condition involves the value of fields on both sides of the interface, the resulting flux in Eq. (7b) depends on and . As argued in [16, 17], the HDG method represents the intermediate state as in the space instead of eliminating ; the jump condition is enforced weakly.

replaces the numerical flux in Eq. (6) and is given by either (Lax-Friedrichs) or (upwind) defined on the boundary of each cell as

 (ˆF(LF)L)ij(q,ˆq)nj =(FL)ij(q)nj+cg√ϕB(qi−ˆqi) (9a) (ˆF(up)L)ij(q,ˆq)nj =(FL)ij(q)nj+cg√ϕBBik(qk−ˆqk). (9b)

To close the system of equations, continuity of the flux is enforced weakly by requiring that on each facet

 ⟨{{(ˆFL)ij(q,ˆq)nj}},ˆvi⟩e=0for all test functions ˆv∈ˆΛh. (10)

To simplify notation, introduce the following abbreviations:

 L(q,v) =∑K∈Ωh{((FL)ij(q),∂jvi)K−⟨(F∗L)ij(q+,q−)nj,vi⟩∂K} (11) ˆL(q,ˆq,v) =∑K∈Ωh{((FL)ij(q),∂jvi)K−⟨(ˆFL)ij(q,ˆq)nj,vi⟩∂K} N0(q,v) =∑K∈Ωh(si(q),vi)K(see Eq. % (???)) M(q,v) =∑K∈Ωh(qi,vi)K Ξ(q,ˆq,ˆv) =∑e∈Eh⟨{{(ˆFL)ij(q,ˆq)nj}},ˆvi⟩e

With the definitions in Eq. (11), after summing over all cells the weak form of the linear SWEs discretised with DG in Eq. (6) can be written in condensed form as

 ∂M(q,v)∂t=N0(q,v)+L(q,v)for all v∈Wh. (12)

The equivalent HDG discretisation of Eq. (6) is obtained by replacing by in Eq. (12) and enforcing the continuity condition in Eq. (10). This results in

 ∂M(q,v)∂t =N0(q,v)+ˆL(q,ˆq,v), subject toΞ(q,ˆq,ˆv) =0for all v∈Wh, ˆv∈ˆΛh. (13)

As has been shown in [17], the solution of the original, non-hybridized system in Eq. (12) is identical to the solution of the HDG problem given in Eq. (13). In general, the additional field introduced in Eq. (13) on the facets has three components. As will be discussed in Section 4.1 below, can be eliminated for the upwind flux, whereas does not enter the equations if the Lax-Friedrichs method is used.

##### Case 2: non-linear SWEs.

Similarly, the weak form of the non-linear SWEs in Eq. (2) can be written for each cell as

 (∂qi∂t,vi)K−(Fij(q),∂jvi)K+⟨F∗ij(q+,q−)nj,vi⟩∂K=(si(q),vi)Kfor all v∈Wh. (14)

where is the non-linear Lax-Friedrichs flux

 F(LF)ij(q+,q−)nj={{Fij(q)}}nj+cg2τ∗⟨⟨qi⟩⟩ (15)

with . Further define

 N(q,v)=∑K∈Ωh{((F−FL)ij(q),∂jvi)K−⟨(F∗−F∗L)ij(q+,q−)nj,vi⟩∂K+(si(q),vi)K} (16)

With the definitions in Eqs. (11) and (16), after summing over all cells the weak form of the DG discretisation in Eq. (14) can be written in condensed form as

 ∂M(q,v)∂t=N(q,v)+L(q,v)for all v∈Wh. (17)

As above, the equivalent HDG discretisation of Eq. (14) is obtained by replacing by in the linear term in Eq. (17) and enforcing the continuity condition in Eq. (10). This results in

 ∂M(q,v)∂t =N(q,v)+ˆL(q,ˆq,v), subject toΞ(q,ˆq,ˆv) =0for all v∈Wh, ˆv∈ˆΛh. (18)

By splitting the right hand side into two terms, we have isolated the fast dynamics due to gravity waves in , whereas any slower modes are described by . This observation is crucial for the construction of suitable semi-implicit time integrators, which avoid numerical instabilities due to the fast gravity waves.

### 3.3 Time discretisation

To integrate Eq. (18) (Eq. (13) can be dealt with in exactly the same way) in time, the term is integrated explicitly, whereas is treated implicitly. Denoting the solution at time by , a simple scheme would for example be the following “-method”:

 M(q(n+1),v)−M(q(n),v)Δt=N(q(n),v)+θˆL(q(n+1),ˆq(n+1),v)+(1−θ)L(q(n),v)for all v∈Wh (19)

with and subject to the flux condition Eq. (10). The equivalent expression for the linear SWEs in Eq. (13) can be obtained by replacing in Eq. (19). At each time step, calculating the field requires the solution of the linear system

 M(q(n+1),v)−θΔtˆL(q(n+1),ˆq(n+1),v) =R(q(n),ˆq(n),v) (20) :=M(q(n),v)+Δt(N(q(n),v)+(1−θ)L(q(n),v)) subject toΞ(q,ˆq,ˆv) =0for all v∈Wh, ˆv∈ˆΛh.

Note that for a purely explicit method (), the system matrix arising from the left hand side of the system in Eq. (20) is simply the mass matrix (obviously, in this case Eq. (19) reduces to the explicit Euler integrator). On the other hand, -stage IMEX methods [2, 3, 4, 5, 6] are a generalisation of the method in Eq. (19). They require the construction of intermediate states which are obtained by solving the system of linear equations

 M(Q(i),v) =M(q(n),v)+Δti−1∑j=1aijN(Q(j),v)+Δti∑j=1˜aijˆL(Q(j),ˆQ(j),v), (21) withΞ(Q(i),ˆQ(i),ˆv) =0for each stage i=1,…,s and all v∈Wh, % ˆv∈ˆΛh.

The state at the next timestep satisfies

 M(q(n+1),v)=M(q(n),v)+Δts∑i=1biN(Q(i),v)+Δts∑i=1˜biˆL(Q(i),ˆQ(i),v). (22)

and is obtained by an additional mass-solve. Each particular IMEX method is determined by a choice of coefficients , , , which are commonly known as the Butcher tableau coefficients for the explicit and implicit terms. The equations in Eq. (21) can be solved stage-by-stage. At every stage a linear system of the abstract general form

 ˆA(q,ˆq,v,ˆv)=M(q,v)−αΔtˆL(q,ˆq,v)+Ξ(q,ˆq,ˆv)=R(v)for all v∈Wh, ˆv∈ˆΛh. (23)

with some positive constant has to be solved222 Note that the -method” in Eq. (19) could be written as a 2-stage IMEX method with

but this would be inefficient since it requires an additional mass-solve in the first stage.. Typically the value of is different in each stage of the IMEX method and the right hand side depends on previously calculated fields. To see that in Eq. (23) is indeed linear in the combined field , use the definitions of , and in Eqs. (9a), (9b), (11) and observe that

 ˆA(c1q1+c2q2,c1ˆq1+c2ˆq2,v,ˆv)=c1ˆA(q1,ˆq1,v,ˆv)+c2ˆA(q2,ˆq2,v,ˆv)for all v∈Wh, ˆv∈ˆΛh

and for all , .

## 4 Schur-complement multigrid solver

Constructing efficient solvers for Eq. (23) is the topic of this paper. In the following section we give explicit expressions for for the Lax-Friedrichs and upwind flux. We then describe a Schur complement approach which reduces the problem to solving an elliptic system in the flux variable on the facets. This flux system is solved with a new non-nested multigrid preconditioner based on ideas in [20].

### 4.1 Linear equation

Recall that in the domain, on the facets, and write the corresponding test-functions as , . Using the definitions in Eqs. (7a) and (7b), the bilinear form on the left hand side of Eq. (23) can be written down explicitly for the upwind- and Lax-Friedrichs flux. As shown in [17], the resulting bilinear form for the upwind flux does not depend on . For the Lax-Friedrichs flux the form does not depend on . Explicit expressions for those forms are given as follows:

Upwind flux:

 ˆA(up)(q,ˆq,v,ˆv)=ˆA(up)(ϕ,u,ˆϕ,ψ,w,ˆψ) :=(ϕψ+u⋅w)Ωh−cgαΔt[(u⋅∇ψ+ϕBϕ∇⋅w)Ωh (24) −([[uψ]]+2√ϕB({{ϕψ}}−ˆϕ{{ψ}})+ϕBˆϕ[[w]])Eh] +(ˆψ[[[u]]+2√ϕB({{ϕ}}−ˆϕ)])Eh

Lax-Friedrichs flux:

 ˆA(LF)(q,ˆq,v,ˆv)=ˆA(LF)(ϕ,u,ˆu,ψ,w,ˆw) :=(ϕψ+u⋅w)Ωh−cgαΔt[(u⋅∇ψ+ϕBϕ∇⋅w)Ωh (25) −(ˆu[[ψ]]+2√ϕB({{u⋅w}}−ˆu⋅{{w}})+ϕB[[ϕw]])Eh] +(ˆw⋅[ϕB[[ϕ]]+2√ϕB({{u}}−ˆu)])Eh

The system in Eq. (23) can be solved with a Schur-complement approach by eliminating the field defined on the elements and leaving a problem for defined on the facets. For this write

 (26)

where is the Riesz representer of . The scalar products in Eq. (26) are defined as

 (q,v)Ωh =(ϕψ+u⋅w)Ωh:=∫Ω(ϕ(x)ψ(x)+u(x)⋅w(x))dx, (ˆq,ˆv)Eh

and the linear operators

 A :Wh→Wh, B :ˆΛh→Wh, C :Wh→ˆΛh, M :ˆΛh→ˆΛh. (27)

Explicit expressions for the operators in Eq. (27) can be deduced from Eqs. (25) and (24) and are given both for the Lax-Friedrichs- and for the upwind flux in B.

### 4.2 Schur complement solve

Let the inverse of be defined in the obvious way as

 (FAq,v)Ωh=(q,v)Ωhfor all q,v∈Wh. (28)

It can now be shown that Eq. (23) can be solved in three steps.

Step 1.

Given in Eq. (26), calculate the modified right hand side

 ˆR(ˆv)=(CFf,ˆv)Ωhfor all ˆv∈ˆΛh.
Step 2.

Solve the Schur-complement system

 ˆS(ˆq,ˆv)=(Mˆq,ˆv)Eh+(cgαΔt)2(CFBˆq,ˆv)Eh=ˆR(ˆv)for all ˆv∈ˆΛh (29)

for the flux unknown .

Step 3.

Reconstruct the solution as

 q=F(f−(cgαΔt)2Bˆq). (30)

The crucial observation is that both and only couple unknowns within one grid cell and therefore can both be represented by a set of cell-local matrices. Since further and couple the unknowns in one grid cell to the flux unknowns on its boundary, the matrix representation of the operator in Eq. (29) is sparse. Recall further that for the Lax-Friedrichs flux only and appear in Eq. (23) whereas for the upwind flux, the equation only contains and . Therefore Eq. (29) reduces to

 ˆS(up)(ˆϕ,ˆψ) =ˆR(up)(ˆψ) for all ˆψ∈Λh (upwind) (31a) ˆS(LF)(ˆu,ˆw) =ˆR(LF)(ˆw) for all ˆw∈Λh×Λh (Lax-Friedrichs) (31b)

in those cases. Before discussing the solution of Eq. (29) in Section 4.4, we briefly discuss the bilinear form in Eq. (24) in the context of previous work.

### 4.3 Relationship to previous work in [20]

For the HDG upwind flux, the Schur-complement system in Eq. (29) is very similar to the elliptic problem written down in Theorem 2.1 of [20]. To see this, observe that for constant , Eq. (2.5) in [20] can be written in our notation as

 ˆA(CGDT)(ϕ,u,ˆϕ,ψ,w,ˆψ)=b(w,ψ)

for some right hand side . The bilinear form is given by

 ˆA(CGDT)(ϕ,u,ˆϕ,ψ,w,ˆψ) :=(u⋅w)Ωh−c−1[(κu⋅∇ψ+ϕ∇⋅w)Ωh (32) −(κ[[uψ]]+2κτK({{ϕψ}}−ˆϕ{{ψ}})+ˆϕ[[w]])Eh] +(ˆψ[[[u]]+2τK({{ϕ}}−ˆϕ)])Eh

where is a constant. If we assume spatially non-varying bathymetry () for the shallow water problem and set , , in Eq. (32), the resulting bilinear form differs from the one given in our Eq. (24) only by a zero order term, namely:

 ˆA(up)(ϕ,u


text:
# Solid of Revolution – Comparing Methods #2

To get my solid, I rotated the area between the curves $$y=x$$ and $$y=x^2$$ around the line $$y=1.25$$.  I then approximated it by cylindrical shells.  The width of each of my shells was 1/16 and the height was $$x^2-x$$, where the $$x$$-value was chosen from the right endpoint of each interval.

One of the issues with the print comes from the fact aht  the top shell does not connect smoothly to the shell below it.  Initially I thought that this resulted from an error in the printer.  When I looked at the Cinema 4D file, however, I discovered that the top shell does not make contact with the shell below it because its height is less than 1/16 (the bowl shape is very flat there).  Additionally, the shell second from the top does not make contact with the shell below it, but the height difference is not noticable.  Another issue with the print was that the leftover supports look bad.  The inside of the object was filled with them, and even after clearing out most of the supports, it still looks terrible.  We are thinking about printing the object on a different printer, such as the Form Labs printer (which uses stereolithography) or the ProJet 260 (which uses powder), in order to get the inside part of the object to look clean.

# Solid of Revolution – Comparing Methods #1

We decided to find a solid of revolution for which both the washer method and cylindrical shell method worked and to model it with both methods. Laura Taalman has done similar designs for shell approximations which can be found here and here

We decided to rotated the area between $$y=x$$ and $$y=x^2$$ about the line $$y=1.25$$. This created a bowl-like object which is curved on the outside, straight on the inside and with a hole at the bottom of the bowl. We split this project up and I created a model of the object using the washer method, while Ryan used the cylindrical shell method. We planned to each use 16 slices/shells so we could compare our two models later.

Creating this object was a very similar process to that of creating the sphere with washers. Since I was using 16 washers and the object has a height of 1, each washer was 1/16 thick. In Cinema4D I used the preprogramed tube object and made each of them to have a height of 1/16. In order to calculate the inner and outer radius of each washer I used the value of middle of each washer, 1/32, 3/32, 5/32, 7/32, etc. I plugged these values into the equation $$f(x)=1.25-x$$ for the inner radius value and $$g(x)=1.24-x^2$$ for the outer radius value. The only other thing I had to change for each was was its height coordinate. My washers in Cinema4D were parallel to the $$xz$$-plane so I adjusted each washer’s $$y$$-coordinate so they would be spread out to the correct height. I increased the $$y$$-coordinate of each was by 1/16 (the height of each washer) from the previous.

Since I used the preprogramed tube object I was unable to use Magic Merge easily to connect these object and instead just used the ‘Connect + Delete’ command under the ‘Mesh’ tab. This means that there are inner walls in my object, but since each washer is relatively small this is not hugely problematic and adds extra support.

I ran into one issue when I was finishing up with my object. The smallest washer (and the last one I created) did not touch the previous washer. All the other washers had been able to lie on top of the previous one since their outer radius was larger than the inner radius of the previous washer. The 15th washer had an inner radius of 0.34375 cm and the 16th washer had an outer radius of 0.31152cm and thus there was a gap. This is due to the fact the bowl is very flat at the bottom.

We decided to remedy this we would just delete the 16th washer and make a note of it here as well as in the description of the object on Thingiverse.

We printed Ryan’s object with cylindrical shells using the Makerbot and the supports on the inside were difficult to remove and didn’t look good. We are currently discussing printing on another printer for these objects and I will update later on how that goes!

# Volume by Cylindrical Shells

On Monday, June 15, I modeled a volume by cylindrical shells from Calculus II.  I used Example 1 in 7.3 of Stewart’s Essential Calculus, which is a volume of revolution of the curve $$y=2x^2-x^3$$ about the y-axis. This is shaped a bit like a stadium. The plan is to approximate this volume using 16 cylindrical shells. I first sketched out the curve in 2-dimensions to get a feel for the profile of the shape. Next, I wrote out each point on the curve from $$x=0$$ to $$x=2$$ in intervals of length 1/8, namely $$(0,0), … , (1.875,0.439), (2,0)$$.

Each cylindrical shell is determined by its height, the thickness of the shell and either the inner or outer radius. By construction, each shell has thickness 1/8 and for each $$k=0,1,2, \dots, 15$$, the inner radius of a shell is $$k/8$$, while the outer radius is $$(k+1)/8$$. I chose the height of each shell to be $$f(k/8)$$, which is the the $$x$$-value closer to the origin.

To summarize, if $$r=$$inner radius of shell, $$R=$$outer radius of shell, and $$h=$$height of tube, then for $$k = 0, 1, 2, \dots , 15$$,

$$r = k/8, \quad R = (k+1)/8, \quad h=2r^2 – r^3.$$

Since the height function is increasing between $$x=0$$ and $$x=4/3$$, the cylindrical shells lie inside the volume of revolution. Between $$x=4/3$$ and $$x=2$$, the function is decreasing and the shells lie outside the volume.

I then made a model of this Riemann approximation for the solid in Cinema 4D by inserting tubes (Cinema 4D’s name for “cylindrical shells”) of inner radius $$r$$, outer radius $$R$$, and height $$h$$. When Cinema 4D inserts a tube, it places half of the tube above the $$xy$$-plane (the $$xz$$-plane in Cinema 4D) and half of it below. Therefore, I needed to add half of each tube’s height to put each shell onto the $$xy$$-plane.

The object was first printed with the supports setting on just in case, although I thought that they did not need supports. The result was that a couple superfluous strands of plastic running along the shells, which needed to be removed. Given the geometry of the shape, I would advise against using supports in building this object.

I would advise that you use a raft to build this object, however, because it was very difficult to remove the object from the build plate. (It took 5 minutes of very careful tugging by David Pfaff.)

You can find this objects on Thingiverse here.


text:
# Second order partial derivatives

Suppose that $$z=g(x,y),\, x=s+t,$$ and $$y=st$$, where all first and second order partial derivatives of $$g$$ exist and are continuous.

Show that $$\frac{\partial ^2 z}{\partial s\partial t}=\frac{\partial ^2 g}{\partial x^2}+x\frac{\partial^2g}{\partial x\partial y}+y\frac{\partial^2 g}{\partial y^2}+\frac{\partial g}{\partial y}.$$

Someone told me I have to use the chain rule twice, but I still don't quite understand what I'm meant to do.

• Well, I expanded all the second order terms (e.g. the term on the LHS would be d/ds(dz/dt), and I assumed that dg/dy = dz/dy since z = g(x,y). And I found the partial derivatives dx/ds, dx/dt, dy/ds, dy/dt etc, but I'm having trouble understanding how to put it all together – user190322 Nov 6 '14 at 11:51
• Could you edit the question to show how you expanded $\frac{\partial^2g}{\partial x^2}$ for example? It will give us a better look on the troubles you're having and it will likely produce more helpfull answers. – gebruiker Nov 6 '14 at 11:57
• If you are having trouble formatting the maths, there is a tutorial here on how to format on this site. – gebruiker Nov 6 '14 at 12:00

EDIT: @Kevin is correct, I needed to replace full derivaties w.r.t $s, t$ with partials. The remaining parts of the answer is correct. I would reply to Kevin as a comment about this but the "add comment" function on my PC seems to be disabled.


text:
# All Seminars

Show:
Title: Congruences from quaternion algebras
Seminar: Algebra
Speaker: Kimball Martin of University of Oklahoma
Contact: John Duncan, john.duncan@emory.edu
Date: 2017-11-07 at 4:00PM
Venue: W306
Abstract:
Studying congruences between modular forms is a prosperous avenue in number theory. One approach to obtaining congruences involves computations on the Jacobian (Mazur, Ribet, ...). For instance, Mazur uses the Jacobian to determine when there is a weight 2 cusp form of prime level congruent to an Eisenstein series, which has various applications. We will explore another approach to obtaining congruences of modular forms using the arithmetic of quaternion algebras and the Jacquet-Langlands correspondence. This will lead to (1) generalizations work of Mazur and Ribet on weight 2 Eisenstein congruences, and (2) a phenomenon of many mod 2 congruences between weight k cusp forms.
Title: Extremal number of configurations in a grid
Seminar: Combinatorics
Speaker: Marcelo Sales of University of Sao Paulo
Contact: Dwight Duffus, dwight@mathcs.emory.edu
Date: 2017-11-06 at 4:00PM
Venue: W302
Abstract:
A configuration is a finite set of points with no three collinear. Two configurations have the same order type if there exists a bijection between these two configurations that preserves the orientation of every ordered triple. A configuration A contains a copy of a configuration B some subset of A has the same order type of B and we denote by B \subset A. For a configuration B and an integer m, the extremal number ex(m,B)= max {|A| : B is not a subset of A, A \subset [m]^2} is the maximum size of a subset of the grid $[m]^2$ without a copy of $B$. We discuss some bounds on this function for general B.
Title: Insights from computational fluid dynamic modelling for aortic arch pathologies
Seminar: Numerical Analysis and Scientific Computing
Speaker: Massimiliano Marrocco-Trischitta of San Donato Hospital in Milan, Italy
Date: 2017-11-03 at 2:00PM
Venue: W301
Abstract:
Objectives. To assess whether the geometrical and hemodynamic reappraisal of the Ishimarus Aortic Arch Map according to Aortic Arch Classification in Type I, II, and III, may provide valuable information regarding the suitability for thoracic endovascular aortic repair (TEVAR), and the risk of aortic dissection. Methods. Anonymized thoracic computed tomography scans of healthy aortas were reviewed, and stratified according to the Aortic Arch Classification. Twenty patients of each Type of Arch were selected. Further processing allowed calculation of angulation and tortuosity of each proximal landing zones. Data were described indicating both proximal landing zone and Type of Arch (e.g. 0/I). Also, among these 60 CT angiography scans, 15 were selected, 5 per Type of Arch, for further analysis. Computational fluid dynamics were performed to compute displacement forces, exerted by pulsatile blood flow on the aortic wall in the defined landing areas. Equivalent surface tractions were computed dividing the displacement forces magnitude of each proximal landing zone by the corresponding area. The three-dimensional orientation (x,y,z) of displacement forces was described as an upward (z direction), and a sideways component (x-y plane).
Title: Improving Question Answering by Bridging Linguistic Structures and Statistical Learning
Defense: Dissertation
Speaker: Tomasz Jurczyk of Emory University
Contact: TBA
Date: 2017-11-02 at 4:00PM
Venue: W301
Abstract:
Question answering (QA) has lately gained lots of interest from both academic and industrial research. No matter the question, search engine users expect the machines to provide answers instantaneously, even without searching through relevant websites.\\ \\While a significant portion of these questions ask for concise and well known facts, more complex questions do exist and they often require dedicated approaches to provide robust and accurate systems.\\ \\This thesis explores linguistically-oriented approaches for both factoid and non-factoid question answering and applications to cross-genre tasks. The contributions include new annotation schemes for the question answering oriented corpora, extracting linguistic structures and performing matching, and early exploration of applications to conversation dialog tasks.
Title: 576 Fermions
Seminar: Algebra
Speaker: Theo Johnson-Freyd of Perimeter Institute
Contact: John Duncan, john.duncan@emory.edu
Date: 2017-10-24 at 4:00PM
Venue: W306
Abstract:
The Stolz--Teichner conjectures predict that the generalized cohomology theory called Topological Modular Forms has a geometric model in terms of the space of 2-dimensional supersymmetric quantum field theories, and that holomorphic vertex operator superalgebras provide the geometric model for nontrivial degrees of TMF. Since TMF is periodic with period 576, these conjectures in particular predict an equivalence between holomorphic VOSAs of different central charge that had not been discovered by physicists. I will report on progress constructing this "periodicity" equivalence geometrically. Specifically, I will explain the solution to the warm-up problem of constructing geometrically the 8-fold periodicity of real K-theory: my solution realizes this periodicity as an example of super symplectic reduction. I will then explain why I believe the Conway group Co0 will play a role in the 576-fold periodicity problem, and why my recent computation of H^4(Co0) provides evidence for this belief.
Title: Survey on recent results on maximal tori of algebraic groups
Seminar: Algebra
Speaker: Philippe Gille of CNRS, Lyon
Contact: John Duncan, john.duncan@emory.edu
Date: 2017-10-24 at 5:00PM
Venue: W306
Abstract:
Prasad and Rapinchuk investigated the isopectrality problem for certain Riemannian varieties by analysing in what extent a semisimple algebraic group defined over a number field is determined by its maximal tori. We shall report advances on this topic by Chernousov/Rapinchuk/Rapinchuk, Bayer-Fluckiger/Lee/Parimala and others by discussing the case of non-archimedean fields and local-global principles.
Title: Packing nearly optimal Ramsey $R(3,t)$ graphs
Seminar: Combinatorics
Speaker: Lutz Warnke of Georgia Institute of Technology
Contact: Dwight Duffus, dwight@mathcs.emory.edu
Date: 2017-10-23 at 4:00PM
Venue: W302
Abstract:
TBA
Title: Backtracking-Based Accelerated Descent Methods for Large-Scale Linear Inverse Problems
Seminar: Numerical Analysis and Scientific Computing
Speaker: Xianqi Li of University of Florida
Contact: Lars Ruthotto, lruthotto@emory.edu
Date: 2017-10-20 at 2:00PM
Venue: W301
Abstract:
Large-scale linear inverse problems arise in a wide range of applications such as image processing and statistical inference. However, the high dimensional (possibly dense and ill-conditioned) matrix in data fidelity term often brings significantly computational challenges when solving the formulated optimization problem and hence hindered the applicability of the sophisticated interior point method and second-order optimization methods. To tackle those challenges, first-order gradient descent method turns into a good choice. In this talk, we first review some classical first-order accelerated descent methods, then introduce our proposed backtracking based accelerated descent methods, which are capable of hunting for more aggressive stepsize via conducting fewer number of line searches. A brief convergence analysis will be presented. The numerical results on structured (low rank and/or sparsity and/or group sparsity) network learning and total-variation based image reconstruction problems indicate the efficiency and effectiveness of the proposed algorithms.
Title: You Are Already Living Inside a Computer
Seminar: Computer Science
Speaker: Ian Bogost of Georgia Tech
Contact: TBA
Date: 2017-10-20 at 3:00PM
Venue: W201
Abstract:
Futurists and philosophers have made dramatic predictions about the future of computers. Artificial intelligence might end the need for human work, or it might enslave humanity, or it might facilitate a kind of rapture into machines, where people's consciousnesses could upload and, through simulation, become immortal.\\ \\Realized or not, those future visions offer a stark contrast with the reality of computing today: One where digital machinery is embedded in the ordinary, human world rather than leading away from it. The computational aspects of ordinary things from smartphone apps to internet-connected toasters have become goals unto themselves, rather than just a means to an end. As it spreads from desktops and back-offices to pockets, cameras, cars, and door locks, the affection people have with computers transfers onto other, even more ordinary objects. And the more people love using computers for everything, the more life feels incomplete unless it takes place inside them. Reality might have beaten the futurists to the punch, by turning computing into a way of life.
Title: The density of squarefree values taken by a polynomial
Colloquium: N/A
Speaker: Manjul Bhargava of Princeton
Contact: John Duncan, john.duncan@emory.edu
Date: 2017-10-20 at 4:00PM
Venue: W201
Abstract:
It is well known that the density of integers that are squarefree is $6/\pi^2$, giving one of the more intriguing occurrences of $\pi$ where one might not a priori expect it! A natural next problem that has played an important role in number theory is that of understanding the density of squarefree values taken by an integer polynomial. We survey a number of recent results on this problem for various types of polynomials - some of which use the ABC Conjecture'' and some of which do not.


text:
# Managing subplot layouts in MATLAB2

Posted by Doug Hull,

Many people use subplot to layout axes, but do not know there is a second way of using it so that all the axes can be sized differently. Also, a method for having a common x-axis is shown.

h(4) = subplot(6,1,4:6);
h(4) = subplot(2,1,2);
linkaxes(h,'x')


text:
Volume 397 - The Ninth Annual Conference on Large Hadron Collider Physics (LHCP2021) - Joint sessions
Searches for additional Higgs bosons at the LHC
J. Tao*  on behalf of the ATLAS and CMS Collaborations
Full text: pdf
Pre-published on: October 20, 2021
Published on: November 17, 2021
Abstract
Some latest results of the direct searches for additional Higgs bosons at ATLAS and CMS are presented. The searches for heavy neutral Higgs bosons, light neutral Higgs bosons, singly and doubly charged Higgs bosons are performed using proton-proton collision data at $\sqrt{s}$ = 13 TeV collected by the ATLAS and CMS detectors. No significant sign of any beyond standard model Higgs boson is observed. For each search the upper limit on the Higgs boson production cross section times branching ratio is set.
DOI: https://doi.org/10.22323/1.397.0077
How to cite

Metadata are provided both in "article" format (very similar to INSPIRE) as this helps creating very compact bibliographies which can be beneficial to authors and readers, and in "proceeding" format which is more detailed and complete.

Open Access


text:
Free Version
Moderate

# Second Order Homogeneous - General Solution

DIFFEQ-GQG27J

Find the general solution of the differential equation:

$$x'''+2x''-x'-2x=0$$

A

$x(t)=c_1e^t+c_2te^{t}+c_3e^{-2t}$

B

$x(t)=c_1e^{-t}+c_2te^{-t}+c_3e^{-2t}$

C

$x(t)=c_1e^t+c_2e^{-t}+c_3te^{-t}$

D

$x(t)=c_1e^t+c_2e^{-t}+c_3e^{-2t}$