diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..d0574b5 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,15 @@ +services: + - docker + +branches: + only: + - master + +before_install: +- docker pull firedrakeproject/firedrake +- docker run -dti -v `pwd`:/home/firedrake/sapphire -w /home/firedrake --name fd firedrakeproject/firedrake +- docker exec fd ls +- docker exec fd pip3 install pytest + +script: +- docker exec fd bash -c "source firedrake/bin/activate && python3 -m pytest -v sapphire/tests" diff --git a/README.md b/README.md index a77867a..7d304b7 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ is an engine for constructing PDE-based simulations discretized in space with mixed finite elements and in time with finite differences. - +![](https://travis-ci.org/geo-fluid-dynamics/sapphire.svg?branch=master) [![DOI](https://zenodo.org/badge/157389237.svg)](https://zenodo.org/badge/latestdoi/157389237) diff --git a/sapphire/benchmarks/freeze_water_in_cavity.py b/sapphire/benchmarks/freeze_water_in_cavity.py index 99abd69..a45bd07 100644 --- a/sapphire/benchmarks/freeze_water_in_cavity.py +++ b/sapphire/benchmarks/freeze_water_in_cavity.py @@ -70,10 +70,8 @@ def dirichlet_boundary_conditions(sim): W = sim.function_space - dim = sim.mesh.geometric_dimension() - return [fe.DirichletBC( - W.sub(1), (0.,)*dim, "on_boundary"), + W.sub(1), (0., 0.), "on_boundary"), fe.DirichletBC(W.sub(2), sim.hot_wall_temperature, 1), fe.DirichletBC(W.sub(2), sim.cold_wall_temperature, 2)] @@ -90,22 +88,22 @@ def initial_values(sim): sim.prandtl_number = sim.prandtl_number.assign(Pr) - dim = sim.mesh.geometric_dimension() + w = fe.Function(sim.function_space) + + p, u, T = w.split() + + p.assign(0.) + + ihat, jhat = sim.unit_vectors() - T_c = sim.cold_wall_temperature.__float__() + u.assign(0.*ihat + 0.*jhat) - w = fe.interpolate( - fe.Expression( - (0.,) + (0.,)*dim + (T_c,), - element = sim.element), - sim.function_space) + T.assign(sim.cold_wall_temperature) F = heat_driven_cavity_variational_form_residual( sim = sim, solution = w)*fe.dx(degree = sim.quadrature_degree) - T_h = sim.hot_wall_temperature.__float__() - problem = fe.NonlinearVariationalProblem( F = F, u = w, @@ -158,7 +156,7 @@ def variational_form_residual(sim, solution): class Simulation(sapphire.simulations.convection_coupled_phasechange.Simulation): - def __init__(self, *args, spatial_dimensions, meshsize, **kwargs): + def __init__(self, *args, meshsize, **kwargs): self.reference_temperature_range__degC = fe.Constant(10.) @@ -166,17 +164,9 @@ def __init__(self, *args, spatial_dimensions, meshsize, **kwargs): self.cold_wall_temperature = fe.Constant(0.) - if spatial_dimensions == 2: - - Mesh = fe.UnitSquareMesh - - elif spatial_dimensions == 3: - - Mesh = fe.UnitCubeMesh - super().__init__( *args, - mesh = Mesh(*(meshsize,)*spatial_dimensions), + mesh = fe.UnitSquareMesh(meshsize, meshsize), variational_form_residual = variational_form_residual, initial_values = initial_values, dirichlet_boundary_conditions = dirichlet_boundary_conditions, diff --git a/sapphire/benchmarks/melt_octadecane_in_cavity.py b/sapphire/benchmarks/melt_octadecane_in_cavity.py index 2e3a804..88d1278 100644 --- a/sapphire/benchmarks/melt_octadecane_in_cavity.py +++ b/sapphire/benchmarks/melt_octadecane_in_cavity.py @@ -4,13 +4,21 @@ def initial_values(sim): - return fe.interpolate( - fe.Expression( - (0., 0., 0., sim.cold_wall_temperature.__float__()), - element = sim.element), - sim.function_space) - - + w = fe.Function(sim.function_space) + + p, u, T = w.split() + + p.assign(0.) + + ihat, jhat = sim.unit_vectors() + + u.assign(0.*ihat + 0.*jhat) + + T.assign(sim.cold_wall_temperature) + + return w + + def dirichlet_boundary_conditions(sim): W = sim.function_space diff --git a/sapphire/postprocessing/__init__.py b/sapphire/postprocessing/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/sapphire/postprocessing/vtk.py b/sapphire/postprocessing/vtk.py deleted file mode 100644 index 53e7996..0000000 --- a/sapphire/postprocessing/vtk.py +++ /dev/null @@ -1,186 +0,0 @@ -import numpy -import scipy.interpolate -from matplotlib import pyplot as plt -import vtk -from vtk.util.numpy_support import vtk_to_numpy - - -def read_vtk_data(vtk_filepath): - - reader = vtk.vtkXMLUnstructuredGridReader() - - reader.SetFileName(vtk_filepath) - - reader.Update() - - data = reader.GetOutput() - - return data - - -def get_connectivity(vtk_data): - - cell_data = vtk_to_numpy(vtk_data.GetCells().GetData()) - - # See https://perso.univ-rennes1.fr/pierre.navaro/read-vtk-with-python.html - connectivity = numpy.take( - cell_data, - [i for i in range(cell_data.size) if i%4 != 0])\ - .reshape(cell_data.size//4, 3) - - return connectivity - - -def plot_mesh(vtk_data, axes = None): - """ Plot the mesh """ - connectivity = get_connectivity(vtk_data) - - point_data = vtk_to_numpy(vtk_data.GetPoints().GetData()) - - if axes is None: - - _, axes = plt.subplots() - - plt.triplot(point_data[:,0], point_data[:,1], connectivity, axes = axes) - - axes.set_aspect("equal") - - return axes - - -def plot_scalar_field_contours( - vtk_data, - scalar_solution_component = 0, - filled = False, - colorbar = True, - axes = None, - **kwargs): - """ Plot contours of a scalar field """ - point_data = vtk_to_numpy(vtk_data.GetPoints().GetData()) - - x = point_data[:,0] - - y = point_data[:,1] - - u = vtk_to_numpy(vtk_data.GetPointData().GetArray( - scalar_solution_component)) - - if axes is None: - - _, axes = plt.subplots() - - args = (x, y, get_connectivity(vtk_data), u) - - if filled: - - mappable = axes.tricontourf(*args, **kwargs) - - else: - - mappable = axes.tricontour(*args, **kwargs) - - axes.set_aspect("equal") - - axes.set_xlabel("$x$") - - axes.set_ylabel("$y$") - - if colorbar: - - return axes, plt.colorbar(mappable = mappable, ax = axes) - - else: - - return axes - - -def plot_scalar_field(vtk_data, scalar_solution_component = 0, axes = None): - """ Plot a scalar field """ - return plot_scalar_field_contours( - vtk_data = vtk_data, - scalar_solution_component = scalar_solution_component, - filled = True, - axes = axes, - levels = 128) - - -def plot_vector_field(vtk_data, vector_solution_component = 0, axes = None, - **kwargs): - """ Plot a vector field """ - point_data = vtk_to_numpy(vtk_data.GetPoints().GetData()) - - x = point_data[:,0] - - y = point_data[:,1] - - u = vtk_to_numpy(vtk_data.GetPointData().GetArray( - vector_solution_component)) - - if axes is None: - - _, axes = plt.subplots() - - plt.quiver(x, y, u[:,0], u[:,1], axes = axes, **kwargs) - - axes.set_aspect("equal") - - axes.set_xlabel("$x$") - - axes.set_ylabel("$y$") - - return axes - - -def format_data_for_streamplot(x, y, u, v): - - x_unique = numpy.unique(x) - - y_unique = numpy.unique(y) - - x_grid, y_grid = numpy.meshgrid(x_unique, y_unique) - - u_grid = scipy.interpolate.griddata((x, y), u, (x_grid, y_grid)) - - v_grid = scipy.interpolate.griddata((x, y), v, (x_grid, y_grid)) - - return x_unique, y_unique, u_grid, v_grid - - -def plot_streamlines(vtk_data, vector_solution_component = 0, axes = None, - max_linewidth = 3., - **kwargs): - """ Plot streamlines of a vector field """ - points = vtk_to_numpy(vtk_data.GetPoints().GetData())[:,:2] - - velocities = vtk_to_numpy(vtk_data.GetPointData().GetArray( - vector_solution_component)) - - if axes is None: - - _, axes = plt.subplots() - - u = velocities[:,0] - - v = velocities[:,1] - - x_unique, y_unique, u_grid, v_grid = format_data_for_streamplot( - x = points[:,0], - y = points[:,1], - u = u, - v = v) - - speed_grid = numpy.sqrt(u_grid**2 + v_grid**2) - - stream = axes.streamplot(x_unique, y_unique, u_grid, v_grid, - linewidth = max_linewidth*speed_grid/speed_grid.max(), - arrowsize = 1.e-8, - **kwargs) - - axes.set_aspect("equal") - - axes.set_xlabel("$x$") - - axes.set_ylabel("$y$") - - return axes, stream - \ No newline at end of file diff --git a/sapphire/simulation.py b/sapphire/simulation.py index 72f1b06..779e147 100644 --- a/sapphire/simulation.py +++ b/sapphire/simulation.py @@ -177,6 +177,10 @@ def assign_parameters(self, parameters): return self + def unit_vectors(self): + + return unit_vectors(self.mesh) + def unit_vectors(mesh): diff --git a/setup.py b/setup.py index 4190f06..e8f7287 100644 --- a/setup.py +++ b/setup.py @@ -3,11 +3,6 @@ setuptools.setup( name = 'sapphire', - version = '0.4.3a0', + version = '0.4.4a0', packages = setuptools.find_packages(), - install_requires = [ - 'matplotlib', - 'numpy', - 'scipy', - 'vtk'] ) diff --git a/tests/model_validation/test__convection_coupled_phasechange.py b/tests/model_validation/test__convection_coupled_phasechange.py index 0170bab..3c180ca 100644 --- a/tests/model_validation/test__convection_coupled_phasechange.py +++ b/tests/model_validation/test__convection_coupled_phasechange.py @@ -73,7 +73,7 @@ def test__validate__melt_octadecane__regression(): assert(abs(sim.liquid_area - expected_liquid_area) < tolerance) -def freeze_water(endtime, s, tau, rx, nx, rt, nt, q, dim = 2, outdir = ""): +def freeze_water(endtime, s, tau, rx, nx, rt, nt, q, outdir = ""): mu_l__SI = 8.90e-4 # [Pa s] @@ -97,11 +97,9 @@ def freeze_water(endtime, s, tau, rx, nx, rt, nt, q, dim = 2, outdir = ""): quadrature_degree = q, element_degree = rx - 1, time_stencil_size = rt + 1, - spatial_dimensions = dim, meshsize = nx, output_directory_path = str(outdir.join( "freeze_water/" - + "dim{0}/".format(dim) + "s{0}_tau{1}/".format(s, tau) + "rx{0}_nx{1}_rt{2}_nt{3}/".format(rx, nx, rt, nt) + "q{0}/".format(q)))) diff --git a/tests/postprocessing/test__vtk.py b/tests/postprocessing/test__vtk.py deleted file mode 100644 index f2d0ab7..0000000 --- a/tests/postprocessing/test__vtk.py +++ /dev/null @@ -1,168 +0,0 @@ -import xml.etree.ElementTree -import sapphire.postprocessing.vtk -import sapphire.test - - -datadir = sapphire.test.datadir - -vtk_filename = "freeze_water/solution_4.vtu" - -def test__plot_mesh(datadir): - - data = sapphire.postprocessing.vtk.read_vtk_data( - vtk_filepath = str(datadir.join(vtk_filename))) - - axes = sapphire.postprocessing.vtk.plot_mesh(vtk_data = data) - - outpath = datadir.join("mesh.png") - - print("Saving {0}".format(outpath)) - - axes.get_figure().savefig(str(outpath)) - - -def test__plot_scalar_field_contours(datadir): - - data = sapphire.postprocessing.vtk.read_vtk_data( - vtk_filepath = str(datadir.join(vtk_filename))) - - - for filled in (False, True): - - axes, colorbar = sapphire.postprocessing.vtk.plot_scalar_field_contours( - vtk_data = data, - scalar_solution_component = 2, - filled = filled, - levels = 8) - - colorbar.ax.set_title("$T$") - - name = "temperature_contours" - - if filled: - - name += "_filled" - - outpath = datadir.join(name + ".png") - - print("Saving {0}".format(outpath)) - - axes.get_figure().savefig(str(outpath)) - - -def test__plot_scalar_field(datadir): - - data = sapphire.postprocessing.vtk.read_vtk_data( - vtk_filepath = str(datadir.join(vtk_filename))) - - axes, colorbar = sapphire.postprocessing.vtk.plot_scalar_field( - vtk_data = data, - scalar_solution_component = 2) - - colorbar.ax.set_title("$T$") - - outpath = datadir.join("temperature.png") - - print("Saving {0}".format(outpath)) - - axes.get_figure().savefig(str(outpath)) - - -def test__plot_vector_field(datadir): - - data = sapphire.postprocessing.vtk.read_vtk_data( - vtk_filepath = str(datadir.join(vtk_filename))) - - axes = sapphire.postprocessing.vtk.plot_vector_field( - vtk_data = data, - vector_solution_component = 1, - headwidth = 5) - - outpath = datadir.join("velocity.png") - - print("Saving {0}".format(outpath)) - - axes.get_figure().savefig(str(outpath)) - - -def test__plot_streamlines(datadir): - - data = sapphire.postprocessing.vtk.read_vtk_data( - vtk_filepath = str(datadir.join(vtk_filename))) - - axes, _ = sapphire.postprocessing.vtk.plot_streamlines( - vtk_data = data, - vector_solution_component = 1) - - axes.set_xlim((0., 1.)) - - axes.set_ylim((0., 1.)) - - outpath = datadir.join("velocity_streamlines.png") - - print("Saving {0}".format(outpath)) - - axes.get_figure().savefig(str(outpath)) - - -def test__plot_superposed_scalar_and_vector_fields(datadir): - - data = sapphire.postprocessing.vtk.read_vtk_data( - vtk_filepath = str(datadir.join(vtk_filename))) - - axes, colorbar = sapphire.postprocessing.vtk.plot_scalar_field( - vtk_data = data, - scalar_solution_component = 2) - - colorbar.ax.set_title("$T$") - - axes = sapphire.postprocessing.vtk.plot_vector_field( - vtk_data = data, - vector_solution_component = 1, - axes = axes, - headwidth = 5) - - outpath = datadir.join("temperature_and_velocity.png") - - print("Saving {0}".format(outpath)) - - axes.get_figure().savefig(str(outpath)) - - -vtk_dir = "freeze_water/" - -pvd_filepath = vtk_dir + "solution.pvd" - -def test__plot_unsteady_superposed_scalar_and_vector_fields(datadir): - - etree = xml.etree.ElementTree.parse( - str(datadir.join(pvd_filepath))).getroot() - - for time, vtu_filename in [ - (element.attrib["timestep"], element.attrib["file"]) - for element in etree[0]]: - - data = sapphire.postprocessing.vtk.read_vtk_data( - vtk_filepath = str(datadir.join(vtk_dir + vtu_filename))) - - axes, colorbar = sapphire.postprocessing.vtk.plot_scalar_field( - vtk_data = data, - scalar_solution_component = 2) - - colorbar.ax.set_title("$T$") - - axes = sapphire.postprocessing.vtk.plot_vector_field( - vtk_data = data, - vector_solution_component = 1, - axes = axes, - headwidth = 5) - - axes.set_title("$t = {0}$".format(time)) - - outpath = datadir.join( - "temperature_and_velocity__t{0}.png".format(time)) - - print("Saving {0}".format(outpath)) - - axes.get_figure().savefig(str(outpath)) - \ No newline at end of file diff --git a/tests/postprocessing/test__vtk/freeze_water/solution.pvd b/tests/postprocessing/test__vtk/freeze_water/solution.pvd deleted file mode 100644 index ef8d768..0000000 --- a/tests/postprocessing/test__vtk/freeze_water/solution.pvd +++ /dev/null @@ -1,8 +0,0 @@ - - - - - - - - diff --git a/tests/postprocessing/test__vtk/freeze_water/solution_0.vtu b/tests/postprocessing/test__vtk/freeze_water/solution_0.vtu deleted file mode 100644 index 174aceb..0000000 Binary files a/tests/postprocessing/test__vtk/freeze_water/solution_0.vtu and /dev/null differ diff --git a/tests/postprocessing/test__vtk/freeze_water/solution_2.vtu b/tests/postprocessing/test__vtk/freeze_water/solution_2.vtu deleted file mode 100644 index 51298f4..0000000 Binary files a/tests/postprocessing/test__vtk/freeze_water/solution_2.vtu and /dev/null differ diff --git a/tests/postprocessing/test__vtk/freeze_water/solution_4.vtu b/tests/postprocessing/test__vtk/freeze_water/solution_4.vtu deleted file mode 100644 index d44df98..0000000 Binary files a/tests/postprocessing/test__vtk/freeze_water/solution_4.vtu and /dev/null differ