thetisproject / thetis Goto Github PK
View Code? Open in Web Editor NEWFinite element flow solver for simulating coastal and estuarine flows.
License: Other
Finite element flow solver for simulating coastal and estuarine flows.
License: Other
The new firedrake viz output format supports multiple fields in the same file (they are all interpolated/projected into the same space). I note that at the moment each field is dumped in its own file. This is partly because that's all that's supported for viz up to now. Is it better to put everything in one file?
Is it still the case that Zenodo requests are done manually for Thetis?
I'd like a Zenodo release for commit "424173012fbb6f60a85f544c560a18a061d5a5cb"
Thanks,
Simon
It seems that for any field that is defined on a continuous function space the corresponding vtk output produces broken vtu files. I think this is due to the self-proclaimed hack (line 123 of thetis/exporter.py) to avoid storing a separate dg coordinates function in every File
object. If I comment out that line, things seem to work again. This is presumably because vtu files with only continuous fields do not have a discontinuous mesh. Is the fix to just add a and not is_cg(function)
?
Reproducer:
from thetis import *
from thetis.exporter import ExportManager
mesh = UnitSquareMesh(2,2)
function_spaces = [
FunctionSpace(mesh, "DG", 0, name="P0"),
FunctionSpace(mesh, "CG", 1, name="P1"),
FunctionSpace(mesh, "CG", 2, name="P2"),
FunctionSpace(mesh, "DG", 1, name="P1DG"),
FunctionSpace(mesh, "DG", 2, name="P2DG")]
fields = {}
field_meta_data = {}
for fs in function_spaces:
name = fs.name + "_field"
field = Function(fs, name=name)
fields[name] = field
field_meta_data[name] = {'shortname': name, 'filename': name}
em = ExportManager('outputs', fields.keys(), fields, field_meta_data, export_type='vtk')
em.export()
There are two failures (see #183 for Jenkins errors):
test_adjoint/test_swe_adjoint.py::test_gradient_from_adjoint[setup_steady]
fails with a linear solver convergence error.test_adjoint/examples/test_examples.py::test_examples[tidalfarm.py]
also fails to converge. In this case we can see that it is because the functional is uniformly zero.Thetis implements a bunch of different time integrator schemes, but currently these are not tested properly.
We could at least solve an ODE with them (using a small mesh and a solution that is constant in space) and verify correct convergence rate.
Hi,
Odd problem that seems to occur on one of my simulations. I read in a netcdf file to create the bathymetry/topography. When I look at the output the bathymetry looks ok, except it seems to be clipped at 27m depth, despite the topography going to 40m+ (i.e. -40m). It otherwise looks OK (not skewed of "off" from my mesh).
So I then printed out the minimum in the bathymetry.data.dat (i.e. the highest land) variable at various stages of the simulation and get the following output:
Reading bathymetry/topography data from inputs/shetlands_palaeo_bathy_100m.nc
Interpolating bathymetry/topography data onto simulation mesh
Bathymetry/topography min:-196.91967167572423
Bathymetry/topography max:129.63601321824117
-196.91967167572423
before init
-196.91967167572423
before init cond
-196.91967167572423
dt = 4.0
Using time integrator: CrankNicolson
after init cond
26.610915552352306
before start
26.610915552352306
0 0 T= 0.00 eta norm: 0.1133 u norm: 0.1133 0.00
volume2d rel. error 0.0000e+00
Exporting: uv_2d elev_2d
Exporting:
time: 4.0
26.610915552352306
The lines of code are:
print("before init")
print(min(bathymetry2d.dat.data))
update_forcings(0.0)
# set initial condition for elevation, piecewise linear function
print("before init cond")
print(min(bathymetry2d.dat.data))
solverObj.assign_initial_conditions(uv=Constant((1.0e-6,0.0)), elev=Constant(1.0e-6))
print("after init cond")
print(min(bathymetry2d.dat.data))
#uv, elev = solverObj.timestepper.solution.split()
# Run model
print("before start")
print(min(bathymetry2d.dat.data))
You can see the bathymetry gets altered after assign_initial_conditions. This has been driving me mad for a couple of days now. I've dug around the code and can't see anything obvious. Any ideas from someone who knows the code better?
Currently releasing a version of Thetis through zenodo is a manual process. It would be nice if this could be done automatically, the same way as already occurs using firedrake-zenodo with firedrake.
We need a portable binary output format for sharing model results. Presumably hdf5 or netcdf are the the preferred formats. As Thetis may use some esoteric function spaces, I guess we'll need to convert the outputs to P0, P1, or P1DG before exporting.
Memory consumption seems to grow boundlessly with large simulations. Here's an example of the idealized estuary example (turbulence branch):
$pidstat -h -r -p 18846 30
Linux 4.2.0-35-generic (hyakutake) 04/09/2016 _x86_64_ (4 CPU)
# Time UID PID minflt/s majflt/s VSZ RSS %MEM Command
1460219513 1000 18846 514.13 0.00 2312504 1374720 8.42 python
1460219543 1000 18846 368.17 0.07 2534136 1592508 9.75 python
1460219573 1000 18846 399.67 0.00 2787808 1846216 11.30 python
1460219603 1000 18846 310.50 0.00 2970008 2028532 12.42 python
1460219633 1000 18846 183.33 0.00 3077760 2136284 13.08 python
1460219663 1000 18846 0.00 0.00 3077760 2136284 13.08 python
1460219693 1000 18846 1.30 0.00 3077760 2136388 13.08 python
1460219723 1000 18846 0.00 0.00 3077760 2136388 13.08 python
1460219753 1000 18846 8.73 0.00 3107268 2165872 13.26 python
1460219783 1000 18846 98.73 0.00 3262180 2320784 14.21 python
1460219813 1000 18846 306.13 0.00 3469808 2528008 15.47 python
1460219843 1000 18846 363.87 0.00 3715588 2773976 16.98 python
1460219873 1000 18846 338.30 0.00 3903300 2961496 18.13 python
1460219903 1000 18846 380.43 0.00 4152840 3210848 19.65 python
1460219933 1000 18846 298.73 0.00 4333568 3391808 20.76 python
1460219963 1000 18846 318.77 0.00 4555468 3613784 22.12 python
1460219993 1000 18846 354.37 0.00 4772844 3831080 23.45 python
1460220023 1000 18846 276.20 0.00 4951260 4009200 24.54 python
1460220053 1000 18846 353.80 0.00 5171488 4229532 25.89 python
1460220083 1000 18846 365.97 0.00 5356540 4414732 27.02 python
1460220113 1000 18846 258.50 0.00 5567644 4625752 28.32 python
1460220143 1000 18846 386.77 0.00 5787396 4845760 29.66 python
1460220173 1000 18846 383.00 0.00 6004996 5063300 30.99 python
1460220203 1000 18846 345.73 0.00 6223184 5280888 32.33 python
1460220233 1000 18846 280.57 0.00 6403664 5461848 33.43 python
If I manually invoke garbage collector (e.g. at the end of each export) there's no increase:
$pidstat -h -r -p 21984 30
Linux 4.2.0-35-generic (hyakutake) 04/09/2016 _x86_64_ (4 CPU)
# Time UID PID minflt/s majflt/s VSZ RSS %MEM Command
1460221257 1000 21984 3064.67 0.00 2137176 1200828 7.35 python
1460221287 1000 21984 404.30 0.00 2451024 1514800 9.27 python
1460221317 1000 21984 110.47 0.00 2498376 1558064 9.54 python
1460221347 1000 21984 0.00 0.00 2498376 1558064 9.54 python
1460221377 1000 21984 46.37 0.00 2498220 1558212 9.54 python
1460221407 1000 21984 0.17 0.00 2498220 1558212 9.54 python
1460221437 1000 21984 23.37 0.00 2497416 1557428 9.53 python
1460221467 1000 21984 0.17 0.00 2497416 1557428 9.53 python
1460221497 1000 21984 30.23 0.00 2497416 1557468 9.53 python
1460221527 1000 21984 0.00 0.00 2497416 1557468 9.53 python
1460221557 1000 21984 35.63 0.00 2497416 1557468 9.53 python
1460221587 1000 21984 0.00 0.00 2497416 1557468 9.53 python
1460221617 1000 21984 56.87 0.00 2537072 1596768 9.77 python
1460221647 1000 21984 35.97 0.00 2497416 1557468 9.53 python
1460221677 1000 21984 11.33 0.00 2527436 1587444 9.72 python
1460221707 1000 21984 23.83 0.00 2497416 1557468 9.53 python
1460221737 1000 21984 0.00 0.00 2497416 1557468 9.53 python
1460221767 1000 21984 35.13 0.00 2497416 1557468 9.53 python
1460221797 1000 21984 0.00 0.00 2497416 1557468 9.53 python
1460221827 1000 21984 35.17 0.00 2497416 1557468 9.53 python
1460221857 1000 21984 0.00 0.00 2497416 1557468 9.53 python
1460221887 1000 21984 0.00 0.00 2497416 1557468 9.53 python
If system memory starts to run out, garbage collector seems to be called automatically, but you'll have to wait until then.
Why is it necessary to call the garbage collector manually?
And even so, why is there so much memory being allocated at run time? I thought that once you've created the necessary functions, there should be no need to allocate for more.
I'd need a Zenodo release for the current Firedrake stack for the 3D baroclinic paper. See the JSON file
firedrake.json
In many cases users want to compute custom diagnostics of simulation state at run time. It would be useful to have a generic method for supporting this.
I think all the existing callbacks (volume and mass conservation etc) could then be ported to the same interface.
It appears that the default model options are in global scope:
from thetis import *
import thetis.options
a = thetis.options.ModelOptions2d()
a.timestepper_options.solver_parameters = {'for': 'bar'}
print('a: {:}'.format(a.timestepper_options.solver_parameters))
b = thetis.options.ModelOptions2d()
print('b: {:}'.format(b.timestepper_options.solver_parameters))
Prints
a: {'foo': 'bar'}
b: {'foo': 'bar'}
Probably not what we want. At least this affects the test suite.
I think this is because we define the default options as trait objects in attach_paired_options
decorator, rather than a class that should be initiated separately for each instance.
In order to do dual weighted residual (DWR) error estimation in Thetis, we need strong residuals. In the (closed) PR https://github.com/thetisproject/thetis/pull/146, I started working on this by defining strong residuals for each term in the shallow water and 2d tracer equations. However, the code structure needed some major changes. For error estimation, the idea was to use callbacks to periodically output strong residual fields, which are then weighted by adjoint solutions to give the DWR estimator.
The strong residual is also needed for SUPG stabilisation. Suppose we have the strong residual R
for the steady tracer equation. Then SUPG amounts to adding a stabilisation term
tau * dot(u, grad(psi)) * R * dx
to the weak form, where tau
is the stabilisation parameter, u
is the fluid velocity and psi
is the test function used. Instead of having a separate term for SUPG, perhaps it would make sense to simply add
tau * dot(u, grad(psi)) * R_T * dx
to the weak residual in each term T
, where R_T
is the corresponding term in the strong residual.
To do list:
tracer_eq_2d
.tracer_eq
.shallowwater_eq
.momentum_eq
.Also for DWR we need edge residuals to take into account fluxes:
tracer_eq_2d
.tracer_eq
.shallowwater_eq
.momentum_eq
.Hi
I was wondering if it would be possible to do a release of Thetis at the following commit id "bab2234b4f84de6c3e5244ca0e466ec15a62d991".
Thanks for your help!
Use python logging module to handle all print statements etc.
In addition to the current support for sigma coordinates, we would like to be able to use z-layers. This is important where we have shallow bathymetry and would like to exploit this to reduce the number of degrees of freedom we're using. Also probably useful for wetting and drying.
The implementation will have cross-cutting implications across both Firedrake and PyOP2, but the aim here is to iron out the user-facing interface (and hence the necessary abstractions) which will guide the implementation.
Currently, Firedrake only supports a single global number of extruded layers across the whole domain. The first step is providing the ability to have a variable number of columns in the domain. Proposed API:
base = Mesh(...)
start = MeshFunction(base, dtype=int)
stop = MeshFunction(base, dtype=int)
# The convenience version that just has a globally constant layer number remains.
mesh = ExtrudedMesh(base, start_layer=start, stop_layer=stop)
Supporting this should be reasonably straightforward. Note that we do not expect to lose performance by making the number of layers per-column. The performance advantage of extruded meshes comes from the data layout, not that we have hard-coded one of the loop bounds.
We would like not to have internal vertical facets that are "hanging".
+---+---+---+
| | | |
+---+---+---+
| | | |
+---+---+---+
| | <- This causes a problem
+---+
The reason for this is that we don't want to have to treat these specially in the code (all the logic becomes much more complicated). The proposed solution is to create degenerate cells where there is a step in the number of layers, squashing the volume of the hanging facet to zero.
+---+---+---+
| | | |
+---+---+---+
| | | |
+---+---+---+
| /
+/
We will need to be a little careful here, because now all the terms in the Jacobian will be zero, and if we're don't take care, computations of J/detJ on the facet will create a NaN.
I think the model is that "top" layer number is the same in all columns, and so the data model encodes the "start" layer in each column. Since hanging vertical facets are squashed to zero, and therefore do not contribute non-zero values at any point, when we integrate over the vertical facets we have to inspect the number of layers in both of the columns and choose the starting layer corresponding to the shortest column.
In updating to the new Firedrake api (finat, not fiat elements). I fixed element_continuity
. In doing so, I realised that there was a bug in its definition.
def element_continuity(fiat_element):
"""Return an :class:`ElementContinuity` instance with the
continuity of a given element.
:arg fiat_element: The fiat element to determine the continuity
of.
:returns: A new :class:`ElementContinuity` instance.
"""
import FIAT
cell = fiat_element.get_reference_element()
...
edofs = fiat_element.entity_dofs()
dim = cell.get_spatial_dimension()
dg = True
for i in range(dim - 1):
if any(len(k) for k in edofs[i].values()):
dg = False
break
return ElementContinuity(dg, dg, dg)
This incorrectly reports elements with edge dofs, but not vertex dofs as dg (oops!).
Fixing this, leads to a bunch of tests failing. So I wonder what's gone wrong. (this fix is to loop:
for i in range(dim)
, not for i in range(dim-1)
.
The api-changes
branch does this using the new finat elements. But fails tests because it's not bug-for-bug compatible.
I think we could simplify the time integrator API a little:
Omit solution function from TimeIntegrator.advance
call
Currently TimeIntegrator.advance
takes the solution function as an argument. This is redundant as in practice TimeIntegrator
needs to store the solution function as a property in any case. It's also difficult to ensure meaningful operation in case the user happens to pass a different function to the advance
method.
Omit dt
from TimeIntegrator.advance
+ add set_dt
method
We normally use a constant time step so there's no need to pass the time step to every call. In case the user wants to change the time step, I would add set_dt
method instead. Then one could make sure that the time step gets updated correctly. In cases where changing time step is not allowed (e.g. linear multistep methods) we'd raise an exception.
Here's a draft of the base classes
class TimeIntegratorBase(object):
"""Base class for all time integrator objects."""
__metaclass__ = ABCMeta
@abstractproperty
def advance(self, t, update_forcings=None):
"""Advances equations for one time step."""
pass
@abstractproperty
def initialize(self, *args, **kwargs):
"""Initializes the time integrator"""
pass
class TimeIntegrator(TimeIntegratorBase):
"""Base class for all time integrator objects that solve a single equation"""
def __init__(self, equation, solution, fields, dt, solver_parameters={}):
super(TimeIntegrator, self).__init__()
self.equation = equation
self.solution = solution
self.fields = fields
self.dt_const = Constant(dt)
# unique identifier for solver
self.name = '-'.join([self.__class__.__name__,
self.equation.__class__.__name__])
self.solver_parameters = {}
self.solver_parameters.update(solver_parameters)
def set_dt(self, dt):
"""Updates time step"""
self.dt_const.assign(dt)
Here I'm separating TimeIntegratorBase
and TimeIntegrator
classes: All time integrators are derived from the former (e.g. 3D coupled time integrators) while all "normal" time integrators (that iterate a single equation) inherit the latter.
Here I haven't changed the signature for update_forcings
function which is still passed to the advance
method. There might be a better way of handling that too.
Currently function spaces are chosen with the boolean options.mimetic
: False
corresponds to P(n)DG-P(n)DG space and True
to RT(n+1)-P(n)DG pair.
#11 adds an additional options.continuous_pressure
that provides P(n)DG-P(n+1) pair.
All cases use options.order
as the base degree (n).
This needs to be generalized to support more function spaces. Easy way is to define names for each combination, for example
options.element_family = 'dg-dg' # p(n)dg-p(n)dg
options.element_family = 'rt-dg' # rt(n+1)dg-p(n)dg
options.element_family = 'dg-cg' # p(n)dg-p(n+1)
The above names do not contain the degrees as sane combinations are presumably dictated by the stability constraints of shallow water equations. (Are there exceptions?)
Alternatively we could use human readable names (if such exist)?
Apparently a few people are looking for a steady state passive tracer example in Thetis. So I thought I would code up the TELEMAC-2D 'point discharge with diffusion' test case which I considered using a standalone SUPG P1 solver in this paper: https://zenodo.org/record/3653373#.Xlkdaumnw5k
I made the new branch passive-tracer-example
, where there is a coupled 'SteadyState'
solver. If I run examples/passiveTracer/passiveTracer.py
, it gives a very encouraging plot:
Visually, it agrees nicely with the asymptotic solution and SUPG approximation as plotted in the paper (Figure 1, p.92). However, the values are significantly negative! Does anyone have an idea what is going wrong?
#29 fixes the top/bottom conditions for p1dg slope limiters but lateral boundaries are still broken. That is, the elements at boundaries only see the interior neighbors which often leads to too strong limiting.
When running some of the tests I get the following error:
py.test test/operations/
============================================================================================== test session starts ===============================================================================================
platform linux2 -- Python 2.7.9, pytest-2.8.2, py-1.4.30, pluggy-0.3.1
rootdir: /data/lmitche1/models/cofs, inifile:
plugins: pep8-1.0.6, ipdb-0.1.dev2, xdist-1.11, benchmark-2.4.1, cov-1.8.1
collected 32 items
test/operations/test_operations_2d-3d.py ................................
================================================================================== 32 passed, 4 pytest-warnings in 2.81 seconds ==================================================================================
Error in atexit._run_exitfuncs:
Traceback (most recent call last):
File "/usr/lib/python2.7/atexit.py", line 24, in _run_exitfuncs
func(*targs, **kargs)
File "/homes/lmitche1/.local/lib/python2.7/site-packages/colorama/initialise.py", line 19, in reset_all
AnsiToWin32(orig_stdout).reset_all()
File "/homes/lmitche1/.local/lib/python2.7/site-packages/colorama/ansitowin32.py", line 67, in __init__
strip = conversion_supported or not is_a_tty(wrapped)
File "/homes/lmitche1/.local/lib/python2.7/site-packages/colorama/ansitowin32.py", line 17, in is_a_tty
return hasattr(stream, 'isatty') and stream.isatty()
ValueError: I/O operation on closed file
Error in sys.exitfunc:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/usr/lib/python2.7/atexit.pyc in _run_exitfuncs()
22 func, targs, kargs = _exithandlers.pop()
23 try:
---> 24 func(*targs, **kargs)
25 except SystemExit:
26 exc_info = sys.exc_info()
/homes/lmitche1/.local/lib/python2.7/site-packages/colorama/initialise.pyc in reset_all()
17
18 def reset_all():
---> 19 AnsiToWin32(orig_stdout).reset_all()
20
21
/homes/lmitche1/.local/lib/python2.7/site-packages/colorama/ansitowin32.pyc in __init__(self, wrapped, convert, strip, autoreset)
65 # should we strip ANSI sequences from our output?
66 if strip is None:
---> 67 strip = conversion_supported or not is_a_tty(wrapped)
68 self.strip = strip
69
/homes/lmitche1/.local/lib/python2.7/site-packages/colorama/ansitowin32.pyc in is_a_tty(stream)
15
16 def is_a_tty(stream):
---> 17 return hasattr(stream, 'isatty') and stream.isatty()
18
19
ValueError: I/O operation on closed file
This appears to be because in this mode pytest has already closed stdout before colorama gets to it in atexit. colorama.init is claimed to be a no-op on all platforms but windows. But this is not the case, because it registers an atexit handler in all cases. I think you're only using colorama for symbolic ansi escape sequences. So I think one can get around this problem by just not calling colorama.init. We're unlikely to run on windows any time soon.
Hear is the Code
##########################################
from thetis import *
mesh2d = RectangleMesh(20, 20, 10e3, 10e3)
fs_p1 = FunctionSpace(mesh2d, 'CG', 1)
bathymetry_2d = Function(fs_p1, name='Bathymetry').assign(10.0)
solver_obj = solver.FlowSolver(mesh2d, bathymetry_2d, n_layers=6)
options = solver_obj.options
options.element_family = 'dg-dg'
options.polynomial_degree = 1
options.timestepper_type = 'SSPRK22'
options.timestepper_options.use_automatic_timestep = False
options.solve_salinity = False
options.solve_temperature = False
options.simulation_export_time = 50.0
options.simulation_end_time = 3600.
options.timestep = 25.0
solver_obj.create_function_spaces()
init_elev = Function(solver_obj.function_spaces.H_2d)
coords = SpatialCoordinate(mesh2d)
init_elev.project(2.0*exp(-((coords[0] - 4e3)**2 + (coords[1] - 4.5e3)2)/2.2e32))
solver_obj.assign_initial_conditions(elev=init_elev)
solver_obj.iterate()
###########################################
the output
What am I doing Wrong
##########################################
/home/qismath/firedrake/lib/python3.6/site-packages/pyproj/crs.py:77: FutureWarning: '+init=:' syntax is deprecated. ':
' is the preferred initialization method.
return _prepare_from_string(" ".join(pjargs))
UFL:WARNING Discontinuous Lagrange element requested on triangle * interval, creating DQ element.
Coupled time integrator: CoupledTwoStageRK
2D time integrator: ESDIRKTrapezoid
3D time integrator: SSPRK22ALE
3D implicit time integrator: BackwardEuler
Traceback (most recent call last):
File "test_flow.py", line 23, in
solver_obj.assign_initial_conditions(elev=init_elev)
File "/home/qismath/firedrake/src/thetis/thetis/solver.py", line 877, in assign_initial_conditions
self.create_equations()
File "/home/qismath/firedrake/src/thetis/thetis/solver.py", line 800, in create_equations
elem_height=self.fields.v_elem_size_2d)
File "/home/qismath/firedrake/src/thetis/thetis/utility.py", line 950, in init
assert (len(nodes) == out_nodes)
AssertionError
The Equation
class as it is implemented now is too rigid for developing custom time integrators. We need a more flexible API where the user can choose which terms are treated implicitly/explicitly/differently.
Equation class defines all the terms in the equation. It takes a bunch of options and the function space where the solution belongs:
function_space = solver.function_spaces.P1DG
options = {'use_grad_div_viscosity': True}
eq = MyEquation(function_space, **options)
class MyEquation(Equation):
def __init__(self, function_space, ...)
self.function_space = function_space
self.test = TestFunction(self.function_space)
...
Options are typically booleans like above.
The terms are implemented as methods that return the corresponding form, including interface terms and boundary conditions
class MyEquation(Equation):
...
def viscosity(self, fields, bnd_conditions):
"""viscosity operator with interior penalty terms."""
grad_div = self.options.use_grad_div_viscosity
f = 0
f += ... self.test * dx
return f
Sign convention: all terms are assumed to be on the left hand side of the equation A + term = 0.
The terms depend on a fields
dictionary
fields={'uv_3d': solver.fields.uv_3d,
'elev_3d': solver.fields.elev_3d,
'viscosity': None}
Fields should be Coefficient
s, e.g. Function
or Constant
instances, or UFL expressions(?). Some fields are optional, e.g. if viscosity
is omitted (or it is None
) the corresponding terms will be ignored. (Alternatively we could have option.use_viscosity
.) Some fields are required and should be asserted.
The equation object knows all the terms that it implements
class MyEquation(Equation):
def __init__(...):
...
# add some python magic to do this
self.terms = {}
self.terms['viscosity'] = self.viscosity
self.terms['elevation_gradient'] = self.elevation_gradient
To differentiate how the terms should be treated, the user can add labels:
eq.set_label('implicit', ('viscosity', 'elevation_gradient'))
eq.set_label('explicit', 'advection')
eq.set_label('source', 'wind_stress')
The equation class provides a method for evaluating the forms
f_implicit = eq.get_form('implicit', fields)
# returns eq.viscosity(fields) + eq.elevation_gradient(fields)
f_explicit = eq.get_form('explicit', fields)
# returns eq.advection(fields)
Time integrator then takes the equation, fields and labels as an argument
labels = {'implicit': ('viscosity', 'elevation_gradient'),
'explicit': 'advection'}
ti = MyTimeIntegrator(eq, fields, labels, bnd_conditions)
The labels depend on the time integrator, e.g. backward Euler only has implicit terms.
Using the Equation.get_form
method the time integrator can then construct the solvers it needs:
class MyTimeIntegrator(Equation):
def __init__(self, equation, fields, labels, bnd_conditions):
self.equation = equation
...
mass_term = ...
f_k1 = mass_term + self.dt*self.equation.get_form('implicit', fields_k1)
self.problem_k1 = NonlinearVariationalProblem(f_k1, self.solution_k1)
self.solver_k1 = NonlinearVariationalSolver(self.problem_k1)
def set_timestep(self, dt):
self.dt.assign(dt)
def advance(self, time, update_forcings=None):
...
NOTE: Currently the equation object knows the fields it uses and also the solution function. These however may depend on the time integrator so it's better to define them here.
Boundary conditions are defined in a dictionary that define which fields should be replaced at the boundary (for weak boundary conditions). There can be additional boundary specific options.
ocean_bnd_id = 1
river_bnd_id = 2
ocean_uv = Constant((1, 0, 0))
river_uv = Constant((0.05, 0, 0))
ocean_funcs = {'uv_3d': ocean_uv}
river_funcs = {'uv_3d': river_uv, 'clamped': True}
bnd_advection = {ocean_bnd_id: ocean_funcs,
river_bnd_id: river_funcs}
# apply to term
f_adv = eq.advection(fields, bnd_advection)
# or time integrator
bcs = {'advection': bnd_advection}
ti = MyTimeIntegrator(eq, fields, labels, bnd_conditions)
It should be possible to support Dirichlet boundary conditions in a similar fashion,
ocean_funcs = DirichletBC(...)
this just will be passed to the solver instead of the term.
For the sake of clarity the field names should follow field_metadata
(defined in field_defs.py ) as closely as possible. Likewise, options should be similar to options.py where applicable.
With recent firedrake updates, the 3D test cases fail if there's 2D-3D coupling with non-zero coupling term.
For example, the lockExchange
example manifests this error:
$ python lockExchange.py -r coarse
UFL:WARNING Discontinuous Lagrange element requested on triangle * interval, creating DQ element.
Coupled time integrator: CoupledTwoStageRK
2D time integrator: ESDIRKTrapezoid
3D time integrator: SSPRK22ALE
3D implicit time integrator: DIRK22
{'ksp_type': 'gmres', 'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'multiplicative'}
2D mesh: 99 nodes, 128 triangles
3D mesh: 10 layers, 1280 prisms
Horizontal element size: 1414.21 ... 1414.21 m
Vertical element size: 2.000 ... 2.000 m
Element family: dg-dg, degree: 1
Number of tracer DOFs: 7680
Number of cores: 1
Tracer DOFs per core: ~7680.0
- dt 2d swe: 4.711801835659298
- dt h. advection: 70.71067639983731
- dt v. advection: 20.83333333333321
- dt viscosity: 19.999999027690038
- CFL adjusted dt: 2D: inf 3D: 19.999999027690038
- chosen dt: 2D: 19.999999027690038 3D: 19.999999027690038
- adjusted dt: 2D: 19.565217391304348 3D: 19.565217391304348
dt = 19.565217
max h diff 851.8518104386496 - 851.8518896211225
Running lock exchange problem with options:
Resolution: coarse
Reynolds number: 1.0
Use slope limiters: True
Horizontal viscosity: 1000.0
Lax-Friedrichs factor vel: 0.0
Lax-Friedrichs factor trc: 0.0
Exporting to outputs_coarse_dg-dg_tri_p1_visc-const_Re1.0_lf-vel0.0_lf-trc0.0_lim
Elem size: 1414.2135279967463 1414.2135937246817
0 0 T= 0.00 eta norm: 0.0000 u norm: 0.0000 0.00
temp_3d mass rel. error 0.0000e+00
temp_3d overshoot 0 0
volume2d rel. error 0.0000e+00
volume3d rel. error 0.0000e+00
Traceback (most recent call last):
File "PETSc/petscdmshell.pxi", line 322, in petsc4py.PETSc.DMSHELL_CreateFieldDecomposition
File "/home/karna/temp/firedrake/src/firedrake/firedrake/dmhooks.py", line 367, in create_field_decomposition
call_setup=True)
File "/home/karna/temp/firedrake/src/firedrake/firedrake/dmhooks.py", line 172, in add_hook
raise ValueError("Expecting non-empty stack")
ValueError: Expecting non-empty stack
The above exception was the direct cause of the following exception:
SystemError: <built-in method insert of list object at 0x7fd499423c48> returned a result with an error set
Exception ignored in: 'petsc4py.PETSc.traceback'
SystemError: <built-in method insert of list object at 0x7fd499423c48> returned a result with an error set
Traceback (most recent call last):
File "lockExchange.py", line 261, in <module>
parse_options()
File "lockExchange.py", line 257, in parse_options
run_lockexchange(**args_dict)
File "lockExchange.py", line 215, in run_lockexchange
solver_obj.iterate()
File "/home/karna/sources/thetis/thetis/solver.py", line 1164, in iterate
update_forcings, update_forcings3d)
File "/home/karna/sources/thetis/thetis/coupled_timeintegrator.py", line 692, in advance
self.timesteppers.swe2d.solve_stage(i_stage, t, update_forcings)
File "/home/karna/sources/thetis/thetis/rungekutta.py", line 559, in solve_stage
self.solve_tendency(i_stage, t, update_forcings)
File "/home/karna/sources/thetis/thetis/rungekutta.py", line 551, in solve_tendency
self.solver[i_stage].solve()
File "/home/karna/temp/firedrake/src/firedrake/firedrake/variational_solver.py", line 264, in solve
self.snes.solve(None, work)
File "PETSc/SNES.pyx", line 557, in petsc4py.PETSc.SNES.solve
petsc4py.PETSc.Error: error code 77
[0] SNESSolve() line 4437 in /tmp/pip-req-build-zdul6tge/src/snes/interface/snes.c
[0] SNESSolve_NEWTONLS() line 225 in /tmp/pip-req-build-zdul6tge/src/snes/impls/ls/ls.c
[0] KSPSolve() line 707 in /tmp/pip-req-build-zdul6tge/src/ksp/ksp/interface/itfunc.c
[0] KSPSetUp() line 377 in /tmp/pip-req-build-zdul6tge/src/ksp/ksp/interface/itfunc.c
[0] PCSetUp() line 894 in /tmp/pip-req-build-zdul6tge/src/ksp/pc/interface/precon.c
[0] PCSetUp_FieldSplit() line 580 in /tmp/pip-req-build-zdul6tge/src/ksp/pc/impls/fieldsplit/fieldsplit.c
[0] PCFieldSplitSetDefaults() line 537 in /tmp/pip-req-build-zdul6tge/src/ksp/pc/impls/fieldsplit/fieldsplit.c
[0] Petsc has generated inconsistent data
[0] Unhandled case, must have at least two fields, not 0
The error occurs in self.timesteppers.swe2d.solve_stage
on the second iteration (iteration 1, RK stage 0) because the coupling term is zero during the first iteration. The coupling term is a source term in the 2D momentum equation (function solver.fields.split_residual_2d
).
The error only occurs if fieldsplit is used in the 2D solver. The default solver_parameters are {'ksp_type': 'gmres', 'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'multiplicative', 'snes_type': 'newtonls'}
.
We need better documentation of the source code, at least
We should also update the website to contain
Term
, Equation
, TimeIntegrator
abstractionsHi all! I am having a very odd issue with firedrake/thetis. I am trying to run the same code which runs normally in my other machine. The same script used to run in this machine until yesterday, so in order to be sure I reinstalled the firedrake environment and thetis since there was no change in the script I am trying to run.
Could you please guide me as to what might have gone wrong in case you came across something similar?
I attach the error message produced. :
In order to simplify this, I went back to the examples and run channel2d. I tried this both in the master branch and the wetting and drying branch I normally use. I added the output in a text file:
Also, since this is my first bug report, can you please advise if that was clear enough for future reference?
Kind regards,
Thanasis
Hi,
Ran a tidal model using 'dg-dg' and the 2D SWE. I saved the Velocity2D as h5 files, planning to do some post-processing afterwards.
The post processing code is:
...
V = VectorFunctionSpace(mesh, "DG", 1)
uv = Function(V, name='vel_2d')
print uv.dat.data.shape
for i in range(start_file,int(t_end/t_export)+1):
print('Reading h5 files. Time ',i,i*t_export)
checkpoint_file = DumbCheckpoint(file_location + '/Velocity2d_{:05}'.format(i), mode=FILE_READ)
checkpoint_file.load(uv)
checkpoint_file.close()
The model was run on 10 cores, so running my post-processing code on 10 cores also, I get:
Reading Velocity h5 files. Time 10 9000
(163497, 2)
Reading Velocity h5 files. Time 10 9000
Traceback (most recent call last):
File "post_processing_vel.py", line 56, in <module>
checkpoint_file.load(uv)
File "/home/firedrake/src/firedrake/firedrake/checkpointing.py", line 238, in load
self.vwr.popGroup()
File "/usr/lib/python3.5/contextlib.py", line 77, in __exit__
self.gen.throw(type, value, traceback)
File "/home/firedrake/src/PyOP2/pyop2/petsc_base.py", line 347, in vec_context
yield self._vec
File "/home/firedrake/src/firedrake/firedrake/checkpointing.py", line 236, in load
v.load(self.vwr)
File "PETSc/Vec.pyx", line 445, in petsc4py.PETSc.Vec.load
petsc4py.PETSc.Error: error code 76
[3] VecLoad() line 931 in /tmp/pip-req-build-udnn6xhw/src/vec/vec/interface/vector.c
[3] VecLoad_Default() line 404 in /tmp/pip-req-build-udnn6xhw/src/vec/vec/utils/vecio.c
[3] VecLoad_HDF5() line 264 in /tmp/pip-req-build-udnn6xhw/src/vec/vec/utils/vecio.c
[3] Error in external library
[3] Error in HDF5 call H5Dopen2() Status -1
I'm assuming I've done something dumb and this is not a bug, but any help appreciated.
Currently test/swe2d/test_thacker.py
only tests the CrankNicolson
time integrator.
Indeed, all rungekutta.DIRKGeneric
based time integrators in the 2D solver (BackwardEuler
, DIRK22
, DIRK33
) are not wetting-drying aware as the do not assemble the modified bathymetry term (ShallowWaterTerm.wd_bathymetry_displacement
).
I've renamed the repo to thetis and changed the files accordingly in pr #17. If travis shows green light I think we're good.
utility.get_facet_mask
appears to be broken, breaking the 3D model. Below is a minimal example:
from firedrake import *
mesh2d = UnitSquareMesh(10, 10)
mesh = ExtrudedMesh(mesh2d, 10, 0.1)
V = FunctionSpace(mesh, 'DG', 1)
from thetis import get_facet_mask
get_facet_mask(V, 'geometric', 'top')
# returns [0, 1, 2, 3, 4, 5]
get_facet_mask(V, 'geometric', 'bottom')
# also returns [0, 1, 2, 3, 4, 5]
# expected: the top/bottom nodes of a prism, e.g. [0, 1, 2] and [3, 4, 5]
Please could I have Zenodo release for Thetis master? (commit 3161124)
Cheers
It would be useful to have a generic way to export diagnostic fields at run time. For example, as part of some earlier work it became useful to export vorticity at regular intervals (I include the vorticity callback script below).
As it stands, there does not appear to be a standard way to export diagnostic fields like that. One suggestion raised (Tuomas) has been that the user should be able to declare a diagnostic field ( e.g. vorticity_2d
) and then a CallBack function that is used to evaluate it. Once registering the field should be accessible as any other field (solver.fields.vorticity_2d
) and exported with options.fields_to_export = ['vorticity_2d']
class Vorticity2DCallback(DiagnosticCallback):
"""Callback that calculates the vertical vorticity from a DG velocity field"""
def __init__(self, solver_obj, time=0, name="vorticity2d", **kwargs):
"""
:arg solver_obj: Thetis solver object
:arg time: Starting time
:arg **kwargs: any additional keyword arguments, see DiagnosticCallback
"""
super(Vorticity2DCallback,self).__init__(solver_obj)
kwargs.setdefault('append_to_log', False)
kwargs.setdefault('export_to_hdf5', False)
self.outputdir = create_directory(os.path.join(self.solver_obj.options.output_directory,name))
self.solver_obj = solver_obj
self.vtk_file = File(os.path.join(self.outputdir,'Vorticity2d.pvd'))
self.time = time
V = FunctionSpace(solver_obj.mesh2d, 'CG', 1)
# Define trial and test fuction
self.zeta = TrialFunction(V)
self.psi = TestFunction(V)
self.vorticity_cg = Function(V, name=name) # this is where the solution will be stored
self.a = self.psi * self.zeta * dx
self.n = FacetNormal(solver_obj.mesh2d)
self._name = name
self._variable_names = ["min_max"]
@property
def name(self):
return self._name
@property
def variable_names(self):
return self._variable_names
def __call__(self):
time = self.solver_obj.simulation_time + self.time
uv, elev = self.solver_obj.timestepper.solution.split()
# omega_z = du/dy - dv/dx
L = (- Dx(self.psi, 1) * uv[0] + Dx(self.psi, 0) * uv[1]) * dx + \
(self.psi * uv[0] * self.n[1] - self.psi * uv[1] * self.n[0]) * ds
solve(self.a == L, self.vorticity_cg, solver_parameters={'ksp_type': 'cg'})
# Exporting vtk File
self.vtk_file.write(self.vorticity_cg)
# Exporting hdf5 File
chk = DumbCheckpoint(self.outputdir + "/Vorticity_%d" % int(time), mode=FILE_CREATE)
chk.store(self.vorticity_cg, name="Vorticity_2d")
omega_min, omega_max = np.min(self.vorticity_cg.dat.data), np.max(self.vorticity_cg.dat.data)
return [omega_min,omega_max]
def message_str(self, *args):
line = 'Exporting vorticity, minimum: {:f} , maximum: {:f}'.format(args[0],args[1])
return line
In current master test_adjoint
doesn't pass travis tests. I tried re-running a couple of earlier commits that used to be green but they fail now:
https://travis-ci.org/thetisproject/thetis/builds/143139296
Any idea what's going on? Out of memory? Doesn't seem to run out of time.
I think I did blow the caches on travis recently so that could be a part of the reason.
Not sure this is fixable - but it would be good to discuss.
If thetis is imported before thetis_adjoint the exported firedrake object such as solve, Function, etc. are the standard firedrake ones, and not the firedrake_adjoint wrappers, which obviously breaks the adjoint capability. For instance:
$ python -c "from thetis_adjoint import *; print solve.__module__"
dolfin_adjoint.solving
$ python -c "from thetis import *; from thetis_adjoint import *; print solve.__module__"
firedrake.solving
The reason is that if thetis is already imported, the from thetis import *
at the last line of thetis_adjoint/init.py just reuses the already loaded thetis module, and thus never executes the from firedrake_adjoint import *
in thetis/firedrake.py. Just adding an extra from firedrake_adjoint import *
after the from thetis_adjoint import *
(or equivalently adding this as the last line of thetis_adjoint/init.py) does not help as the symbols used by thetis itself will still be the original firedrake ones.
One of the things that this breaks is the py.test tests. If I run the adjoint tests (which have from thetis_adjoint import *
) individually, they run fine; if I run them together with the other tests they break, as they will have been preceded by tests that have imported thetis.
With recent updates, 3D problems run ~6x slower. For example, examples tracerBox
and lockExchange
are affected. 2D example stommel2d
seems to be unaffected.
Bottom roughness length z_0 is currently stored in physical_constants
dict as a scalar value. As z_0 may be spatially varying it should be stored in options as a Constant/Function/Coefficient.
If the user provides maximum velocity scale and maximum viscosity, we can compute maximum stable time step based on the discretization and the time integration scheme. This kind of approach is already implemented for the 3D model.
After changes to timeintegrators and swe equations+terms the channel2d test no longer works. Fixes to update the calls to new API are in fix-swe-sspimex branch but the example still blows up during the run.
Vertical velocity solver appears to be broken. E.g. test continuity3d/test_continuity_mes.py::test_setup1
fails.
There have been no related changes in Thetis recently so it's likely that this is a Firedrake bug.
In #6 travis dies with
==================================== ERRORS ====================================
__________ ERROR collecting test/continuity3d/test_continuity_mes.py ___________
test/continuity3d/test_continuity_mes.py:8: in <module>
from scipy import stats
E ImportError: No module named scipy
__________ ERROR collecting test/swe2d/test_steady_state_basin_mms.py __________
test/swe2d/test_steady_state_basin_mms.py:15: in <module>
from scipy import stats
E ImportError: No module named scipy
__________ ERROR collecting test/tracerEq/test_steady_adv-diff_mms.py __________
test/tracerEq/test_steady_adv-diff_mms.py:8: in <module>
from scipy import stats
E ImportError: No module named scipy
=============== 53 passed, 2 skipped, 3 error in 988.77 seconds ================
The command "py.test --travis -v test" exited with 1.
although python-scipy
is in .travis.yml
. Any ideas?
In setup9:
grad_h_x = '''h0/(2.0*lx*sqrt(0.3*x[0]*x[0] + 0.2*x[1]*x[1]))*0.6*x[0] - 3.0*pi/lx*sin(pi*(3.0*x[0] + 1.0*x[1])/lx)'''
grad_h_y = '''h0/(2.0*lx*sqrt(0.3*x[0]*x[0] + 0.2*x[1]*x[1]))*0.4*x[1] - 1.0*pi/lx*sin(pi*(3.0*x[0] + 1.0*x[1])/lx)'''
If x
is the origin (entirely plausible). Then we end up with:
grad_h_x = h0/(2*lx*sqrt(0))*0.6*0 - ...
=> BOOM.
The cofs FFC branch is currently identical to fd_bendy which is somewhat outdated. For example it predates this fix for parallel compilation:
https://bitbucket.org/mapdes/ffc/pull-requests/79/sort-entry-ir-to-ensure-deterministic/diff
We could cherry-pick the fix or rebase cofs branch on a newer version of master. I'd go for the latter but to me it's unclear what is the latest bendy compatible commit in FFC.
These look like OS failures:
https://jenkins.ese.ic.ac.uk:8000/thetisproject/thetis/58
https://jenkins.ese.ic.ac.uk:8000/thetisproject/thetis/59
In addition, login button doesn't seem to work - I cannot stop/restart builds.
(firedrake) jh1889@rheic ~/work/software/firedrake/src/thetis/examples/tidal_barrage $ python lagoon2d.py
dt = 25
Using time integrator: CrankNicolson
Traceback (most recent call last):
File "lagoon2d.py", line 88, in
number_timesteps=5, name="lagoon_1", export_to_hdf5=True)
TypeError: Can't instantiate abstract class LagoonCallback with abstract methods message_str
Up-to-date firedrake and thetis (updated this morning). Let me know if you want any other info.
Thetis is not happy with the namespaces
In [1]: from thetis import *
In [2]: limiter.__name__
Out[2]: 'firedrake.slope_limiter.limiter'
In [3]: from thetis import limiter
In [4]: limiter.__name__
Out[4]: 'firedrake.slope_limiter.limiter'
How should I import thetis.limiter
for this to work?
I think I understand why these are here, however, they can cause some slightly strange side effects.
I think the function cache is easier to address. Effectively you're using this for "work vectors". When you want a temporary function, you ask for one given a particular function space.
I think we could add something like this to Firedrake. Here's how I would do it.
The FunctionSpace object gains a get_function(name=None)
, restore_function(function)
pair. You can then get a temporary function that you can use:
f = V.get_function()
... do stuff with f
V.restore_function(f)
We can wrap this up in a context manger so one could do:
with V.work_space() as f:
... do stuff with f
This is neater, because if the only reference to the work space is inside V, once V gets collected, so do the temporary functions.
The linear problem cache is a little trickier, however, I think it is still possible to fix. The entries in the problem cache become invalid when the mesh they were built on top of go out of scope. So we could make these collectable by adding a level on top using a WeakKeyDictionary
. Sketch:
import weakref
class ProblemCache(object):
def __init__(self):
self._cache = weakref.WeakKeyDictionary()
def __getitem__(self, key):
mesh, key = key
return self._cache[mesh][key]
def __setitem__(self, key, val):
mesh, key = key
d = self._cache.get(mesh, {})
d[key] = val
self._cache[mesh] = d
def __contains__(self, key):
mesh, key = key
return mesh in self._cache and self._cache[mesh].get(key) is not None
def clear(self):
self._cache.clear()
So the idea is to use a weakref to the Mesh as a first-level key, the second level acts as before. Now if the Mesh goes out of scope, the cached problems are no longer strongly referenced by the cache and can be collected.
Running the Columbia plume example on our Linux cluster runs out of memory.
Here's a small test case that reproduces the issue. Mesh file is located in examples/columbia_plume/
example
from thetis import *
meshfile = 'mesh_cre-plume_03_normal.msh'
mesh2d = Mesh(meshfile)
print_output('Loaded mesh ' + mesh2d.name)
# mesh2d = RectangleMesh(90, 90, 100e3, 100e3)
nlayers = 20
dt = 15.
t_export = 150.
t_end = 50*t_export
P1_2d = FunctionSpace(mesh2d, 'CG', 1)
bathymetry_2d = Function(P1_2d, name='bathymetry')
bathymetry_2d.assign(200.)
solver_obj = solver.FlowSolver(mesh2d, bathymetry_2d, nlayers)
options = solver_obj.options
options.element_family = 'dg-dg'
options.timestepper_type = 'SSPRK22'
options.solve_salinity = True
options.solve_temperature = False
options.use_implicit_vertical_diffusion = False
options.use_bottom_friction = False
options.use_turbulence = False
options.use_baroclinic_formulation = False
options.simulation_export_time = t_export
options.simulation_end_time = t_end
options.timestepper_options.use_automatic_timestep = False
options.timestep = dt
options.no_exports = True
solver_obj.assign_initial_conditions(salt=Constant(30.0))
solver_obj.iterate()
Running this with 96 cores, show ~linear increase in memory consumption (of rank 0):
pidstat -h -r -p 19096 15
Linux 3.10.0-514.26.2.el7.x86_64 (compute-0-0) 05/03/2019 _x86_64_ (16 CPU)
# Time UID PID minflt/s majflt/s VSZ RSS %MEM Command
1556868669 51118 19096 11164.60 0.13 2113196 378416 0.58 python
1556868684 51118 19096 4139.53 0.00 2113732 379028 0.58 python
1556868699 51118 19096 2289.73 0.00 2115304 380472 0.58 python
1556868714 51118 19096 416.73 0.00 2116508 381656 0.58 python
1556868729 51118 19096 1909.07 0.00 2117632 382872 0.58 python
1556868744 51118 19096 395.40 0.00 2118832 384096 0.58 python
1556868759 51118 19096 78.07 0.00 2119928 385260 0.59 python
1556868774 51118 19096 1885.40 0.00 2121128 386408 0.59 python
1556868789 51118 19096 43.33 0.00 2122364 387516 0.59 python
1556868804 51118 19096 401.07 0.00 2123428 388728 0.59 python
1556868819 51118 19096 57.40 0.00 2124628 389808 0.59 python
1556868834 51118 19096 1868.27 0.00 2125828 391136 0.59 python
1556868849 51118 19096 264.40 0.00 2126956 392296 0.60 python
1556868864 51118 19096 789.93 0.00 2128052 393464 0.60 python
1556868879 51118 19096 747.53 0.00 2129508 394624 0.60 python
1556868894 51118 19096 2065.40 0.00 2130708 395788 0.60 python
1556868909 51118 19096 33.87 0.00 2131832 396780 0.60 python
1556868924 51118 19096 447.13 0.00 2133068 398096 0.61 python
1556868939 51118 19096 60.67 0.00 2134168 399272 0.61 python
1556868954 51118 19096 3031.33 0.00 2135300 400472 0.61 python
1556868969 51118 19096 42.60 0.00 2136532 401472 0.61 python
1556868984 51118 19096 472.20 0.00 2137596 402816 0.61 python
1556868999 51118 19096 63.07 0.00 2138792 403960 0.61 python
1556869014 51118 19096 1782.80 0.00 2139892 405208 0.62 python
1556869029 51118 19096 47.60 0.00 2141120 406320 0.62 python
1556869044 51118 19096 465.93 0.00 2142320 407600 0.62 python
1556869059 51118 19096 76.80 0.00 2143520 408640 0.62 python
1556869074 51118 19096 1779.67 0.00 2144620 409912 0.62 python
On a static mesh a forward Euler update of field T
can be expressed as
M*T_{n+1} - M*T_{n} = dt*L
where M
is the mass matrix and L
is the linear form of the equation.
If the mesh is time dependent, the conservative arbitrary Lagrangian-Eulerian (ALE) formulation becomes
M_{n+1}*T_{n+1} - M_{n}*T_{n} = dt*L
where the mass matrices now depend on time. Here L
is also modified to take into account the mesh velocity.
What is the easiest way to implement this in firedrake? My guess is that it would make sense to compute and store M_{n} T_{n}
in the previous time step and then use it to compute T_{n+1}
the next iteration. Thus there would be no need to store the mass matrices.
Alternatively you could cast the old solution to the new mesh and then do a conventional solve on a static mesh:
M_{n+1}*T_{n+1} - M_{n+1}*P*T_{n} = dt*L
with P = inv(M_{n+1})*M_{n}
. But this is probably costly.
This should work on arbitrary function spaces: The formulation applies to both tracers (which are most likely DG) and velocity (which can be partially continuous).
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.