Merge pull request #360 from thetisproject/359_fix_output_api

Firedrake outputs API has changed

examples/balzano/balzano.py

@@ -26,6 +26,7 @@
     doi: 10.1016/j.advwatres.2009.09.005.
 """
 from thetis import *
+from firedrake.output.vtk_output import VTKFile
 
 outputdir = 'outputs'
 mesh2d = RectangleMesh(12, 6, 13800, 7200)
@@ -81,9 +82,9 @@
 }
 
 # User-defined output: moving bathymetry and eta_tilde
-wd_bathfile = File(os.path.join(outputdir, 'moving_bath.pvd'))
+wd_bathfile = VTKFile(os.path.join(outputdir, 'moving_bath.pvd'))
 moving_bath = Function(P1_2d, name="moving_bath")
-eta_tildefile = File(os.path.join(outputdir, 'eta_tilde.pvd'))
+eta_tildefile = VTKFile(os.path.join(outputdir, 'eta_tilde.pvd'))
 eta_tilde = Function(P1_2d, name="eta_tilde")
 
 

examples/baroclinic_channel/baroclinic_channel.py

@@ -151,7 +151,7 @@
 mask_temp_relax_3d.dat.data[:] = numpy.maximum(mask_numpy_y0, mask_numpy_y1)
 ix = mask_temp_relax_3d.dat.data < 0
 mask_temp_relax_3d.dat.data[ix] = 0.0
-# File('mask.pvd').write(mask_temp_relax_3d)
+# VTKFile('mask.pvd').write(mask_temp_relax_3d)
 # NOTE must accept ufl expressions as forcing term!
 options.temperature_source_3d = mask_temp_relax_3d/t_temp_relax*(temp_relax - solver_obj.fields.temp_3d)
 

examples/baroclinic_eddies/baroclinic_eddies.py

@@ -29,6 +29,7 @@
 """
 
 from thetis import *
+from firedrake.output.vtk_output import VTKFile
 from diagnostics import *
 
 
@@ -242,7 +243,7 @@ def compute_2d_baroc_head(solver_obj, output):
     # custom export of surface temperature field
     surf_temp_2d = Function(solver_obj.function_spaces.H_2d, name='Temperature')
     extract_surf_temp = SubFunctionExtractor(solver_obj.fields.temp_3d, surf_temp_2d)
-    surf_temp_file = File(options.output_directory + '/SurfTemperature2d.pvd')
+    surf_temp_file = VTKFile(options.output_directory + '/SurfTemperature2d.pvd')
 
     def export_func():
         extract_surf_temp.solve()

examples/bottomFriction/ekman_bottom.py

@@ -100,7 +100,7 @@ def bottom_ekman_test(layers=50, verify=True, iterate=True,
             uv_ana = Function(solver_obj.function_spaces.P1DGv, name='solution')
             uv_ana.interpolate(uv_ana_expr)
 
-            out = File(options.output_directory + '/uv_analytical/uv_analytical.pvd')
+            out = VTKFile(options.output_directory + '/uv_analytical/uv_analytical.pvd')
             out.write(uv_ana)
 
             l2_tol = 0.05

examples/bottomFriction/ekman_surface.py

@@ -7,6 +7,7 @@
 solution.
 """
 from thetis import *
+from firedrake.output.vtk_output import VTKFile
 
 depth = 20.0
 
@@ -89,7 +90,7 @@ def surface_ekman_test(layers=50, verify=True, iterate=True,
     uv_ana = Function(solver_obj.function_spaces.P1DGv, name='solution')
     uv_ana.interpolate(uv_ana_expr)
 
-    out = File(options.output_directory + '/uv_analytical/uv_analytical.pvd')
+    out = VTKFile(options.output_directory + '/uv_analytical/uv_analytical.pvd')
     out.write(uv_ana)
 
     # initialize with a linear v profile to speed-up convergence

examples/bottomFriction/steadyChannel.py

@@ -126,7 +126,7 @@ def bottom_friction_test(layers=25, gls_closure='k-omega',
             z_0 = float(options.bottom_roughness)
             log_uv = Function(solver_obj.function_spaces.P1DGv, name='log velocity')
             log_uv.project(as_vector((u_b / kappa * ln((xyz[2] + depth + z_0)/z_0), 0, 0)))
-            out = File(options.output_directory + '/log_uv/log_uv.pvd')
+            out = VTKFile(options.output_directory + '/log_uv/log_uv.pvd')
             out.write(log_uv)
 
             uv_p1_dg = Function(solver_obj.function_spaces.P1DGv, name='velocity p1dg')

examples/channel_inversion/inverse_problem.py

@@ -133,4 +133,4 @@
     oc.rename(name)
     print_function_value_range(oc, prefix='Optimal')
     if not no_exports:
-        File(f'{options.output_directory}/{name}_optimised.pvd').write(oc)
+        VTKFile(f'{options.output_directory}/{name}_optimised.pvd').write(oc)

examples/columbia_plume/atm_forcing.py

@@ -89,8 +89,8 @@ def test():
     pres_fn = 'tmp/AtmPressure2d.pvd'
     wind_fn = 'tmp/WindStress2d.pvd'
     print('Saving output to {:} {:}'.format(pres_fn, wind_fn))
-    out_pres = File(pres_fn)
-    out_wind = File(wind_fn)
+    out_pres = VTKFile(pres_fn)
+    out_wind = VTKFile(wind_fn)
     hours = 24*1.5
     granule = 4
     simtime = numpy.arange(granule*hours)*3600./granule

examples/columbia_plume/cre-plume.py

@@ -175,7 +175,7 @@
                               visc_bnd_dist, scalar=visc_bnd_value)
 viscosity_bnd_2d.assign(viscosity_bnd_2d + options.horizontal_viscosity)
 ExpandFunctionTo3d(viscosity_bnd_2d, viscosity_bnd_3d).solve()
-# File('bnd_visc.pvd').write(viscosity_bnd_2d)
+# VTKFile('bnd_visc.pvd').write(viscosity_bnd_2d)
 solver_obj.function_spaces.P1_2d.restore_work_function(viscosity_bnd_2d)
 options.horizontal_viscosity = viscosity_bnd_3d
 
@@ -313,7 +313,7 @@ def interp_ocean_bnd(time):
                               lx_relax, scalar=1.0/t_bnd_relax)
 ExpandFunctionTo3d(mask_tmp_2d, mask_uv_relax_3d).solve()
 solver_obj.function_spaces.P1_2d.restore_work_function(mask_tmp_2d)
-# File('mask.pvd').write(mask_tracer_relax_3d)
+# VTKFile('mask.pvd').write(mask_tracer_relax_3d)
 # options.temperature_source_3d = mask_tracer_relax_3d*(temp_bnd_3d - solver_obj.fields.temp_3d)
 # options.salinity_source_3d = mask_tracer_relax_3d*(salt_bnd_3d - solver_obj.fields.salt_3d)
 # options.momentum_source_3d = mask_uv_relax_3d*(uv_bnd_3d - solver_obj.fields.uv_3d)

examples/columbia_plume/ncom_forcing.py

@@ -122,11 +122,11 @@ def test_interpolator():
     vvel_fn = 'tmp/vvel.pvd'
     elev_fn = 'tmp/elev.pvd'
     print('Saving output to {:} {:} {:} {:} {:}'.format(salt_fn, temp_fn, uvel_fn, vvel_fn, elev_fn))
-    out_salt = File(salt_fn)
-    out_temp = File(temp_fn)
-    out_uvel = File(uvel_fn)
-    out_vvel = File(vvel_fn)
-    out_elev = File(elev_fn)
+    out_salt = VTKFile(salt_fn)
+    out_temp = VTKFile(temp_fn)
+    out_uvel = VTKFile(uvel_fn)
+    out_vvel = VTKFile(vvel_fn)
+    out_elev = VTKFile(elev_fn)
 
     out_salt.write(salt)
     out_temp.write(temp)

examples/columbia_plume/roms_forcing.py

@@ -113,8 +113,8 @@ def test_interpolator():
     salt_fn = 'tmp/salt_roms.pvd'
     temp_fn = 'tmp/temp_roms.pvd'
     print('Saving output to {:} {:}'.format(salt_fn, temp_fn))
-    out_salt = File(salt_fn)
-    out_temp = File(temp_fn)
+    out_salt = VTKFile(salt_fn)
+    out_temp = VTKFile(temp_fn)
 
     out_salt.write(salt)
     out_temp.write(temp)

examples/columbia_plume/test_bathy_smoothing.py

@@ -26,7 +26,7 @@
     alpha=1e2, exponent=2.5,
     minimum_depth=3.5, niter=30)
 
-out = File('bath.pvd')
+out = VTKFile('bath.pvd')
 out.write(bathymetry_2d)
 out.write(new_bathymetry_2d)
 
@@ -58,7 +58,7 @@ def compute_hcc(bathymetry_2d, nlayers):
     return f_hcc
 
 
-out = File('hcc.pvd')
+out = VTKFile('hcc.pvd')
 hcc = compute_hcc(bathymetry_2d, nlayers)
 out.write(hcc)
 hcc = compute_hcc(new_bathymetry_2d, nlayers)

examples/columbia_plume/tidal_forcing.py

@@ -31,8 +31,8 @@ def test():
     elev_outfn = 'tmp/tidal_elev.pvd'
     uv_outfn = 'tmp/tidal_uv.pvd'
     print('Saving to {:} {:}'.format(elev_outfn, uv_outfn))
-    elev_out = File(elev_outfn)
-    uv_out = File(uv_outfn)
+    elev_out = VTKFile(elev_outfn)
+    uv_out = VTKFile(uv_outfn)
     for t in numpy.linspace(0, 12*3600., 49):
         tbnd.set_tidal_field(t)
         if elev_field.function_space().mesh().comm.rank == 0:

examples/discrete_turbines/tidal_array.py

@@ -6,6 +6,7 @@
 """
 
 from thetis import *
+from firedrake.output.vtk_output import VTKFile
 
 # Set output directory, load mesh, set simulation export and end times
 outputdir = 'outputs'
@@ -28,7 +29,7 @@
 bathymetry_2d.assign(Constant(50.0))
 x = SpatialCoordinate(mesh2d)
 h_viscosity = Function(P1_2d).interpolate(conditional(le(x[0], 50), 51-x[0], 1.0))
-File(outputdir + '/viscosity/viscosity.pvd').write(h_viscosity)
+VTKFile(outputdir + '/viscosity/viscosity.pvd').write(h_viscosity)
 
 
 # Turbine options

examples/dome/dome.py

@@ -23,6 +23,7 @@
     z-coordinate, isopycnal and non-hydrostatic models. Ocean Modelling, 11(1-2):69-97.
 """
 from thetis import *
+from firedrake.output.vtk_output import VTKFile
 import dome_setup as setup
 import diagnostics
 
@@ -150,7 +151,7 @@
 mask_temp_relax_3d.dat.data[:] = numpy.maximum(mask_numpy_x0, mask_numpy_x1)
 ix = mask_temp_relax_3d.dat.data < 0
 mask_temp_relax_3d.dat.data[ix] = 0.0
-# File('mask.pvd').write(mask_temp_relax_3d)
+# VTKFile('mask.pvd').write(mask_temp_relax_3d)
 options.temperature_source_3d = mask_temp_relax_3d/t_temp_relax*(temp_relax - solver_obj.fields.temp_3d)
 
 # use salinity field as a passive tracer for tracking inflowing waters
@@ -256,7 +257,7 @@ def compute_depth_av_inflow(uv_inflow_3d, uv_inflow_2d):
 # Export bottom salinity
 bot_salt_2d = Function(solver_obj.function_spaces.H_2d, name='Salinity')
 extract_bot_salt = SubFunctionExtractor(solver_obj.fields.salt_3d, bot_salt_2d, boundary='bottom')
-bot_salt_file = File(options.output_directory + '/BotSalinity2d.pvd')
+bot_salt_file = VTKFile(options.output_directory + '/BotSalinity2d.pvd')
 
 
 def export_func():

examples/tidalfarm/tidalfarm.py

@@ -210,4 +210,4 @@ def update_forcings(t):
     # options, such as maxiter and pgtol can be passed on.
     td_opt = minimize(rf, bounds=[0, max_density],
                       options={'maxiter': 100, 'pgtol': 1e-3})
-    File('optimal_density.pvd').write(td_opt)
+    VTKFile('optimal_density.pvd').write(td_opt)

examples/tohoku_inversion/inverse_problem.py

@@ -120,6 +120,6 @@
 source = get_source(mesh2d, source_model, initial_guess=inv_manager.control_coeff_list)
 if source_model == "okada":
     source.subfault_variables = args.okada_parameters
-outfile = File(f"{options.output_directory}/elevation_optimised.pvd")
+outfile = VTKFile(f"{options.output_directory}/elevation_optimised.pvd")
 print_function_value_range(source.elev_init, prefix="Optimal")
 outfile.write(source.elev_init)

examples/tohoku_inversion/plot_optimized_source.py

@@ -45,7 +45,7 @@
 
 # Plot in vtu format
 print_output("Writing optimized elevation to vtu")
-outfile = File(f"{output_dir}/elevation_optimized.pvd")
+outfile = VTKFile(f"{output_dir}/elevation_optimized.pvd")
 print_function_value_range(elev, prefix="Optimal")
 outfile.write(elev)
 

test/bottomFriction/test_bottom_friction.py

@@ -106,7 +106,7 @@ def run_bottom_friction(do_assert=True, do_export=False, **model_options):
         log_uv = Function(solver_obj.function_spaces.P1DGv, name='log velocity')
         log_uv.project(as_vector((u_b / kappa * ln((xyz[2] + depth + z_0)/z_0), 0, 0)))
         if do_export:
-            out = File(outputdir + '/log_uv/log_uv.pvd')
+            out = VTKFile(outputdir + '/log_uv/log_uv.pvd')
             out.write(log_uv)
 
     solver_obj.iterate()

test/continuity3d/test_continuity_mes.py

@@ -127,8 +127,8 @@ def run(setup, refinement, order, do_export=True):
     w_analytical = setup_obj.w(x, y, z, lx, ly)
 
     if do_export:
-        out_w = File(os.path.join(options.output_directory, 'w.pvd'))
-        out_w_ana = File(os.path.join(options.output_directory, 'w_ana.pvd'))
+        out_w = VTKFile(os.path.join(options.output_directory, 'w.pvd'))
+        out_w_ana = VTKFile(os.path.join(options.output_directory, 'w_ana.pvd'))
 
     solver_obj.w_solver.solve()
     uvw = solver_obj.fields.uv_3d + solver_obj.fields.w_3d

test/firedrake/test_divergence_2d.py

@@ -43,8 +43,8 @@ def compute(refinement=1, order=1, do_export=False):
     if do_export:
         print('analytical div {:} {:}'.format(div_uv.dat.data.min(), div_uv.dat.data.max()))
         # export analytical solutions
-        out_uv = File('uv.pvd')
-        out_div = File('div.pvd')
+        out_uv = VTKFile('uv.pvd')
+        out_div = VTKFile('div.pvd')
         out_uv.write(uv)
         out_div.write(div_uv)
 

test/firedrake/test_implicit_diffusion.py

@@ -36,8 +36,8 @@ def test_implicit_diffusion(do_export=False, do_assert=True):
     mesh = ExtrudedMesh(mesh2d, layers=layers, layer_height=-depth/layers)
 
     if do_export:
-        sol_file = File('sol.pvd')
-        ana_file = File('ana_sol.pvd')
+        sol_file = VTKFile('sol.pvd')
+        ana_file = VTKFile('ana_sol.pvd')
 
     # define function spaces
     fam = 'DG'

test/firedrake/test_implicit_friction.py

@@ -32,7 +32,7 @@ def test_implicit_friction(do_export=False, do_assert=True):
     mesh = ExtrudedMesh(mesh2d, layers=50, layer_height=-depth/layers)
 
     if do_export:
-        out_file = File('implicit_bf_sol.pvd')
+        out_file = VTKFile('implicit_bf_sol.pvd')
 
     # ----- define function spaces
     deg = 1

test/momentumEq/test_h-viscosity_mes.py

@@ -87,7 +87,7 @@ def run(refinement, **model_options):
 
     # export analytical solution
     if not options.no_exports:
-        out_uv_ana = File(os.path.join(options.output_directory, 'uv_ana.pvd'))
+        out_uv_ana = VTKFile(os.path.join(options.output_directory, 'uv_ana.pvd'))
 
     def export_func():
         if not options.no_exports:

test/momentumEq/test_v-viscosity_mes.py

@@ -92,7 +92,7 @@ def run(refinement, **model_options):
     solverobj.fields.uv_3d.project(ana_uv_expr)
     # export analytical solution
     if not options.no_exports:
-        out_uv_ana = File(os.path.join(options.output_directory, 'uv_ana.pvd'))
+        out_uv_ana = VTKFile(os.path.join(options.output_directory, 'uv_ana.pvd'))
 
     def export_func():
         if not options.no_exports:

test/pressure_grad/test_baroc_head_mes.py

@@ -72,8 +72,8 @@ def compute_l2_error(refinement=1, fs_type='P1DGxP1DG', no_exports=True):
     mesh.coordinates.dat.data[:, 2] = new_z
 
     if not no_exports:
-        out_density = File('density.pvd')
-        out_baroc_head = File('baroc_head.pvd')
+        out_density = VTKFile('density.pvd')
+        out_baroc_head = VTKFile('baroc_head.pvd')
 
     # project initial density
     beta = -1.5/depth

test/pressure_grad/test_int_pg_mes.py

@@ -82,9 +82,9 @@ def compute_l2_error(refinement=1, quadratic=False, no_exports=True):
     mesh.coordinates.dat.data[:, 2] = new_z
 
     if not no_exports:
-        out_density = File('density.pvd')
-        out_bhead = File('baroc_head.pvd')
-        out_pg = File('int_pg.pvd')
+        out_density = VTKFile('density.pvd')
+        out_bhead = VTKFile('baroc_head.pvd')
+        out_pg = VTKFile('int_pg.pvd')
 
     # project initial scalar
     density_3d.project(density_expr)

test/pressure_grad/test_pg-stack_mes.py

@@ -177,10 +177,10 @@ def compute_l2_error(refinement=1, quadratic_pressure=False, quadratic_density=F
     print_output('Int.PG  L2 error: {:}'.format(l2_err_pg))
 
     if not no_exports:
-        out_temp = File('temperature.pvd')
-        out_density = File('density.pvd')
-        out_bhead = File('baroc_head.pvd')
-        out_pg = File('int_pg.pvd')
+        out_temp = VTKFile('temperature.pvd')
+        out_density = VTKFile('density.pvd')
+        out_bhead = VTKFile('baroc_head.pvd')
+        out_pg = VTKFile('int_pg.pvd')
 
         # export numerical solution
         out_temp.write(temp_3d)

test/slopelimiter/test_slopelimiter.py

@@ -39,7 +39,7 @@ def vertex_limiter_test(dim=3, type='linear', direction='x', export=False):
     tracer.project(tracer_original)
 
     if export:
-        tracer_file = File('tracer.pvd')
+        tracer_file = VTKFile('tracer.pvd')
         tracer_file.write(tracer)
 
     limiter = VertexBasedP1DGLimiter(p1dg)

test/solver3d/test_baroclinic_mms.py

@@ -279,27 +279,27 @@ def run(setup, refinement, polynomial_degree, do_export=True, **options):
     ana_w_3d = sdict['w_3d']
     f = Function(solver_obj.function_spaces.U, name='ana w').project(ana_w_3d)
     if do_export:
-        File(options.output_directory + '/ana_w.pvd').write(f).write(f)
+        VTKFile(options.output_directory + '/ana_w.pvd').write(f).write(f)
 
     ana_temp_3d = sdict['temp_3d']
     f = Function(solver_obj.function_spaces.H, name='ana temp').project(ana_temp_3d)
     if do_export:
-        File(options.output_directory + '/ana_temp.pvd').write(f).write(f)
+        VTKFile(options.output_directory + '/ana_temp.pvd').write(f).write(f)
 
     ana_rho_3d = sdict['density_3d']
     f = Function(solver_obj.function_spaces.H, name='ana rho').project(ana_rho_3d)
     if do_export:
-        File(options.output_directory + '/ana_rho.pvd').write(f).write(f)
+        VTKFile(options.output_directory + '/ana_rho.pvd').write(f).write(f)
 
     ana_bhead_3d = sdict['baroc_head_3d']
     f = Function(solver_obj.function_spaces.H, name='ana bhead').project(ana_bhead_3d)
     if do_export:
-        File(options.output_directory + '/ana_bhead.pvd').write(f).write(f)
+        VTKFile(options.output_directory + '/ana_bhead.pvd').write(f).write(f)
 
     ana_intpg_3d = sdict['int_pg_3d']
     f = Function(solver_obj.function_spaces.U, name='ana int pg').project(ana_intpg_3d)
     if do_export:
-        File(options.output_directory + '/ana_intpg.pvd').write(f).write(f)
+        VTKFile(options.output_directory + '/ana_intpg.pvd').write(f).write(f)
 
     options.volume_source_2d = sdict['vol_source_2d']
     options.momentum_source_2d = sdict['mom_source_2d']

test/sphere/test_williamson.py

@@ -191,7 +191,7 @@ def run(refinement, cell='triangle', setup=setup_williamson2, **model_options):
     solver_obj.create_function_spaces()
     if not options.no_exports:
         # Store analytical elevation to disk
-        out = File(outputdir + '/Elevation2d_ana/Elevation2d_ana.pvd')
+        out = VTKFile(outputdir + '/Elevation2d_ana/Elevation2d_ana.pvd')
         ana_elev = Function(solver_obj.function_spaces.H_2d, name='Elevation')
 
     def export():