From 8087f95720108640e4b78a9b6ef8ebf0f47c6fb7 Mon Sep 17 00:00:00 2001 From: Emilien Bauer Date: Tue, 1 Aug 2023 13:32:32 +0100 Subject: [PATCH 01/25] Add tiling. --- devito/operator/xdsl_operator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/devito/operator/xdsl_operator.py b/devito/operator/xdsl_operator.py index 2ffe000863..ff2283166a 100644 --- a/devito/operator/xdsl_operator.py +++ b/devito/operator/xdsl_operator.py @@ -53,9 +53,9 @@ # gpu-launch-sink-index-computations seemed to have no impact MLIR_GPU_PIPELINE = '"builtin.module(test-math-algebraic-simplification,scf-parallel-loop-tiling{parallel-loop-tile-sizes=128,1,1},func.func(gpu-map-parallel-loops),convert-parallel-loops-to-gpu,fold-memref-alias-ops,expand-strided-metadata,lower-affine,gpu-kernel-outlining,canonicalize,cse,convert-arith-to-llvm{index-bitwidth=64},finalize-memref-to-llvm{index-bitwidth=64},convert-scf-to-cf,convert-cf-to-llvm{index-bitwidth=64},canonicalize,cse,gpu.module(convert-gpu-to-nvvm,reconcile-unrealized-casts,canonicalize,gpu-to-cubin),gpu-to-llvm,canonicalize,cse)"' -XDSL_CPU_PIPELINE = "stencil-shape-inference,convert-stencil-to-ll-mlir,printf-to-llvm" +XDSL_CPU_PIPELINE = "stencil-shape-inference,convert-stencil-to-ll-mlir{tile-sizes=64,64},printf-to-llvm" XDSL_GPU_PIPELINE = "stencil-shape-inference,convert-stencil-to-ll-mlir{target=gpu},printf-to-llvm" -XDSL_MPI_PIPELINE = lambda decomp: f'"dmp-decompose-2d{decomp},canonicalize-dmp,convert-stencil-to-ll-mlir,dmp-to-mpi{{mpi_init=false}},lower-mpi,printf-to-llvm"' +XDSL_MPI_PIPELINE = lambda decomp: f'"dmp-decompose-2d{decomp},canonicalize-dmp,convert-stencil-to-ll-mlir{{tile-sizes=64,64}},dmp-to-mpi{{mpi_init=false}},lower-mpi,printf-to-llvm"' class XDSLOperator(Operator): From 63ab369f3f0ab9923b745f9f40b3646f9a644a2d Mon Sep 17 00:00:00 2001 From: Emilien Bauer Date: Tue, 1 Aug 2023 13:34:06 +0100 Subject: [PATCH 02/25] Add proper quoting. --- devito/operator/xdsl_operator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/devito/operator/xdsl_operator.py b/devito/operator/xdsl_operator.py index ff2283166a..292efa3b62 100644 --- a/devito/operator/xdsl_operator.py +++ b/devito/operator/xdsl_operator.py @@ -53,8 +53,8 @@ # gpu-launch-sink-index-computations seemed to have no impact MLIR_GPU_PIPELINE = '"builtin.module(test-math-algebraic-simplification,scf-parallel-loop-tiling{parallel-loop-tile-sizes=128,1,1},func.func(gpu-map-parallel-loops),convert-parallel-loops-to-gpu,fold-memref-alias-ops,expand-strided-metadata,lower-affine,gpu-kernel-outlining,canonicalize,cse,convert-arith-to-llvm{index-bitwidth=64},finalize-memref-to-llvm{index-bitwidth=64},convert-scf-to-cf,convert-cf-to-llvm{index-bitwidth=64},canonicalize,cse,gpu.module(convert-gpu-to-nvvm,reconcile-unrealized-casts,canonicalize,gpu-to-cubin),gpu-to-llvm,canonicalize,cse)"' -XDSL_CPU_PIPELINE = "stencil-shape-inference,convert-stencil-to-ll-mlir{tile-sizes=64,64},printf-to-llvm" -XDSL_GPU_PIPELINE = "stencil-shape-inference,convert-stencil-to-ll-mlir{target=gpu},printf-to-llvm" +XDSL_CPU_PIPELINE = '"stencil-shape-inference,convert-stencil-to-ll-mlir{tile-sizes=64,64},printf-to-llvm"' +XDSL_GPU_PIPELINE = '"stencil-shape-inference,convert-stencil-to-ll-mlir{target=gpu},printf-to-llvm"' XDSL_MPI_PIPELINE = lambda decomp: f'"dmp-decompose-2d{decomp},canonicalize-dmp,convert-stencil-to-ll-mlir{{tile-sizes=64,64}},dmp-to-mpi{{mpi_init=false}},lower-mpi,printf-to-llvm"' From e15ce9709edc65bd156ad067b15482ba39a68124 Mon Sep 17 00:00:00 2001 From: Emilien Bauer Date: Wed, 2 Aug 2023 13:03:59 +0100 Subject: [PATCH 03/25] Add dimensionality-1 tiling dimensions logic. --- devito/operator/xdsl_operator.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/devito/operator/xdsl_operator.py b/devito/operator/xdsl_operator.py index 292efa3b62..1450c95917 100644 --- a/devito/operator/xdsl_operator.py +++ b/devito/operator/xdsl_operator.py @@ -53,9 +53,9 @@ # gpu-launch-sink-index-computations seemed to have no impact MLIR_GPU_PIPELINE = '"builtin.module(test-math-algebraic-simplification,scf-parallel-loop-tiling{parallel-loop-tile-sizes=128,1,1},func.func(gpu-map-parallel-loops),convert-parallel-loops-to-gpu,fold-memref-alias-ops,expand-strided-metadata,lower-affine,gpu-kernel-outlining,canonicalize,cse,convert-arith-to-llvm{index-bitwidth=64},finalize-memref-to-llvm{index-bitwidth=64},convert-scf-to-cf,convert-cf-to-llvm{index-bitwidth=64},canonicalize,cse,gpu.module(convert-gpu-to-nvvm,reconcile-unrealized-casts,canonicalize,gpu-to-cubin),gpu-to-llvm,canonicalize,cse)"' -XDSL_CPU_PIPELINE = '"stencil-shape-inference,convert-stencil-to-ll-mlir{tile-sizes=64,64},printf-to-llvm"' +XDSL_CPU_PIPELINE = lambda nb_tiled_dims: f'"stencil-shape-inference,convert-stencil-to-ll-mlir{{tile-sizes={",".join(["64"]*nb_tiled_dims)}}},printf-to-llvm"' XDSL_GPU_PIPELINE = '"stencil-shape-inference,convert-stencil-to-ll-mlir{target=gpu},printf-to-llvm"' -XDSL_MPI_PIPELINE = lambda decomp: f'"dmp-decompose-2d{decomp},canonicalize-dmp,convert-stencil-to-ll-mlir{{tile-sizes=64,64}},dmp-to-mpi{{mpi_init=false}},lower-mpi,printf-to-llvm"' +XDSL_MPI_PIPELINE = lambda decomp, nb_tiled_dims: f'"dmp-decompose-2d{decomp},canonicalize-dmp,convert-stencil-to-ll-mlir{{tile-sizes={",".join(["64"]*nb_tiled_dims)}}},dmp-to-mpi{{mpi_init=false}},lower-mpi,printf-to-llvm"' class XDSLOperator(Operator): @@ -114,7 +114,9 @@ def _jit_compile(self): Printer(stream=module_str).print(self._module) module_str = module_str.getvalue() - xdsl_pipeline = XDSL_CPU_PIPELINE + to_tile = len(list(filter(lambda s : str(s) in ["x", "y", "z"], self.dimensions)))-1 + + xdsl_pipeline = XDSL_CPU_PIPELINE(to_tile) mlir_pipeline = MLIR_CPU_PIPELINE if is_omp: @@ -126,7 +128,7 @@ def _jit_compile(self): # reduce the domain of the computation (as devito has already done that for us) slices = ','.join(str(x) for x in shape) decomp = f"{{strategy=2d-grid slices={slices} restrict_domain=false}}" - xdsl_pipeline = XDSL_MPI_PIPELINE(decomp) + xdsl_pipeline = XDSL_MPI_PIPELINE(decomp, to_tile) elif is_gpu: xdsl_pipeline = XDSL_GPU_PIPELINE mlir_pipeline = MLIR_GPU_PIPELINE From bacf1afa49934adef7745280ef5adb52026faf46 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Wed, 2 Aug 2023 18:55:51 +0300 Subject: [PATCH 04/25] mpi: Init effort for serial modelling on wave operator --- devito/ir/ietxdsl/cluster_to_ssa.py | 8 ++++---- devito/ir/ietxdsl/ietxdsl_functions.py | 13 ++++++------- devito/ir/ietxdsl/lowering.py | 2 +- devito/mpi/distributed.py | 4 +++- fast/wave2d.py | 20 ++++++++++++++------ fast/wave3d.py | 19 +++++++++---------- 6 files changed, 37 insertions(+), 29 deletions(-) diff --git a/devito/ir/ietxdsl/cluster_to_ssa.py b/devito/ir/ietxdsl/cluster_to_ssa.py index 2ffb0d1887..1cb9c347d3 100644 --- a/devito/ir/ietxdsl/cluster_to_ssa.py +++ b/devito/ir/ietxdsl/cluster_to_ssa.py @@ -272,7 +272,7 @@ def _ensure_same_type(self, *vals: SSAValue): new_vals.append(val) continue # insert an integer to float cast op - conv = arith.SIToFPOp.get(val, builtin.f32) + conv = arith.SIToFPOp(val, builtin.f32) self.block.add_op(conv) new_vals.append(conv.result) return new_vals @@ -340,7 +340,7 @@ def match_and_rewrite(self, op: func.FuncOp, rewriter: PatternRewriter): self.seen_ops.add(op) rewriter.insert_op_at_start([ - t0 := func.Call.get('timer_start', [], [builtin.f64]) + t0 := func.Call('timer_start', [], [builtin.f64]) ], op.body.block) ret = op.get_return_op() @@ -348,8 +348,8 @@ def match_and_rewrite(self, op: func.FuncOp, rewriter: PatternRewriter): rewriter.insert_op_before([ timers := iet_ssa.LoadSymbolic.get('timers', llvm.LLVMPointerType.typed(builtin.f64)), - t1 := func.Call.get('timer_end', [t0], [builtin.f64]), - llvm.StoreOp.get(t1, timers), + t1 := func.Call('timer_end', [t0], [builtin.f64]), + llvm.StoreOp(t1, timers), ], ret) rewriter.insert_op_after_matched_op([ diff --git a/devito/ir/ietxdsl/ietxdsl_functions.py b/devito/ir/ietxdsl/ietxdsl_functions.py index 20a4d105d8..b4a9498fb1 100644 --- a/devito/ir/ietxdsl/ietxdsl_functions.py +++ b/devito/ir/ietxdsl/ietxdsl_functions.py @@ -20,12 +20,11 @@ # XDSL specific imports from xdsl.irdl import AnyOf, Operation, SSAValue from xdsl.dialects.builtin import (ContainerOf, Float16Type, Float32Type, - Float64Type, Builtin, i32, f32) + Float64Type, i32, f32) -from xdsl.dialects.arith import Muli, Addi from devito.ir.ietxdsl import iet_ssa -from xdsl.dialects import memref, arith, builtin, llvm +from xdsl.dialects import memref, arith, builtin from xdsl.dialects.experimental import math import devito.types @@ -74,7 +73,7 @@ def print_calls(cgen, calldefs): print("Call not translated in calldefs") return - call = Call.get(call_name, C_names, C_typenames, C_typeqs, prefix, retval) + call = Call(call_name, C_names, C_typenames, C_typeqs, prefix, retval) cgen.printCall(call, True) @@ -180,10 +179,10 @@ def add_to_block(expr, arg_by_expr: dict[Any, Operation], result): # reconcile differences if isinstance(rhs.typ, builtin.IntegerType): - rhs = arith.SIToFPOp.get(rhs, lhs.typ) + rhs = arith.SIToFPOp(rhs, lhs.typ) result.append(rhs) else: - lhs = arith.SIToFPOp.get(lhs, rhs.typ) + lhs = arith.SIToFPOp(lhs, rhs.typ) result.append(lhs) @@ -426,7 +425,7 @@ def myVisit(node, block: Block, ssa_vals={}): print(f"Call {node.name} instance translated as comment") return - call = Call.get(call_name, C_names, C_typenames, C_typeqs, prefix, retval) + call = Call(call_name, C_names, C_typenames, C_typeqs, prefix, retval) block.add_ops([call]) print(f"Call {node.name} translated") diff --git a/devito/ir/ietxdsl/lowering.py b/devito/ir/ietxdsl/lowering.py index 515fe2f7a1..8ccf6b6d25 100644 --- a/devito/ir/ietxdsl/lowering.py +++ b/devito/ir/ietxdsl/lowering.py @@ -371,7 +371,7 @@ def match_and_rewrite(self, op: memref.Store, rewriter: PatternRewriter, ssa_indices=[idx], result_type=llvm.LLVMPointerType.typed(op.memref.memref.element_type) ), - store := llvm.StoreOp.get(op.value, gep), + store := llvm.StoreOp(op.value, gep), ], [], ) diff --git a/devito/mpi/distributed.py b/devito/mpi/distributed.py index 5b0c890ebe..50577659cb 100644 --- a/devito/mpi/distributed.py +++ b/devito/mpi/distributed.py @@ -196,7 +196,9 @@ def __init__(self, shape, dimensions, input_comm=None, topology=None): # mpi4py takes care of that when the object gets out of scope self._input_comm = (input_comm or MPI.COMM_WORLD).Clone() - topology = ('*', '*', 1) + if len(shape) == 3: + topology = ('*', '*', 1) + if topology is None: # `MPI.Compute_dims` sets the dimension sizes to be as close to each other # as possible, using an appropriate divisibility algorithm. Thus, in 3D: diff --git a/fast/wave2d.py b/fast/wave2d.py index 409bbfb33b..4d5a7eed5b 100644 --- a/fast/wave2d.py +++ b/fast/wave2d.py @@ -1,8 +1,11 @@ # Based on the implementation of the Devito acoustic example implementation # Not using Devito's source injection abstraction import sys + import numpy as np -from devito import TimeFunction, Eq, Operator, solve, norm, XDSLOperator + +from devito import (TimeFunction, Eq, Operator, solve, norm, + XDSLOperator, configuration) from examples.seismic import RickerSource from examples.seismic import Model, TimeAxis @@ -92,7 +95,7 @@ def plot_2dfunc(u): # Define the wavefield with the size of the model and the time dimension u = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) - +# Another one to clone data u2 = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) # We can now write the PDE @@ -101,18 +104,19 @@ def plot_2dfunc(u): pde = u.dt2 - u.laplace # The PDE representation is as on paper -pde +# pde stencil = Eq(u.forward, solve(pde, u.forward)) -stencil +# stencil # Finally we define the source injection and receiver read function to generate # the corresponding code print(time_range) -print("Init norm:", norm(u)) +print("Init norm:", np.linalg.norm(u.data[:])) src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) op0 = Operator([stencil] + src_term, subs=model.spacing_map, name='SourceDevitoOperator') + # Run with source and plot op0.apply(time=time_range.num-1, dt=model.critical_dt) @@ -125,8 +129,10 @@ def plot_2dfunc(u): print("Init Devito linalg norm 2 :", np.linalg.norm(u.data[2])) print("Norm of initial data:", norm(u)) -# import pdb;pdb.set_trace() + +configuration['mpi'] = 0 u2.data[:] = u.data[:] +configuration['mpi'] = 'basic' # Run more with no sources now (Not supported in xdsl) op1 = Operator([stencil], name='DevitoOperator') @@ -146,7 +152,9 @@ def plot_2dfunc(u): # Reset initial data +configuration['mpi'] = 0 u.data[:] = u2.data[:] +configuration['mpi'] = 'basic' #v[:, ..., :] = 1 diff --git a/fast/wave3d.py b/fast/wave3d.py index abf70d3cdc..1497d3bea2 100644 --- a/fast/wave3d.py +++ b/fast/wave3d.py @@ -2,7 +2,8 @@ # Not using Devito's source injection abstraction import sys import numpy as np -from devito import TimeFunction, Eq, Operator, solve, norm, XDSLOperator +from devito import (TimeFunction, Eq, Operator, solve, norm, + XDSLOperator, configuration) from examples.seismic import RickerSource from examples.seismic import Model, TimeAxis @@ -124,9 +125,10 @@ def plot_3dfunc(u): print("Init linalg norm 1 :", np.linalg.norm(u.data[1])) print("Init linalg norm 2 :", np.linalg.norm(u.data[2])) -print("Norm of initial data:", norm(u)) -import pdb;pdb.set_trace() +print("Norm of initial data:", np.linalg.norm(u.data[:])) +configuration['mpi'] = 0 u2.data[:] = u.data[:] +configuration['mpi'] = 'basic' # Run more with no sources now (Not supported in xdsl) op1 = Operator([stencil], name='DevitoOperator') @@ -142,15 +144,13 @@ def plot_3dfunc(u): print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) -import pdb;pdb.set_trace() - - # Reset initial data +configuration['mpi'] = 0 u.data[:] = u2.data[:] +configuration['mpi'] = 'basic' #v[:, ..., :] = 1 - -print("Reinitialise data: Devito norm:", norm(u)) +print("Reinitialise data: Devito norm:", np.linalg.norm(u.data[:])) print("Init XDSL linalg norm:", np.linalg.norm(u.data[0])) print("Init XDSL linalg norm:", np.linalg.norm(u.data[1])) print("Init XDSL linalg norm:", np.linalg.norm(u.data[2])) @@ -160,10 +160,9 @@ def plot_3dfunc(u): xdslop.apply(time=time_range.num-1, dt=model.critical_dt) xdsl_output = u.copy() + print("XDSL norm:", norm(u)) print(f"xdsl output norm: {norm(xdsl_output)}") -import pdb;pdb.set_trace() - print("XDSL output linalg norm:", np.linalg.norm(u.data[0])) print("XDSL output linalg norm:", np.linalg.norm(u.data[1])) print("XDSL output linalg norm:", np.linalg.norm(u.data[2])) From 7f3b37ef46ed4bce1d9b04a1a4eae92b0251880e Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Wed, 2 Aug 2023 19:34:25 +0300 Subject: [PATCH 05/25] mpi: wip --- fast/wave2d.py | 1 - fast/wave3d.py | 10 +++++----- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/fast/wave2d.py b/fast/wave2d.py index 4d5a7eed5b..9fb0758f9b 100644 --- a/fast/wave2d.py +++ b/fast/wave2d.py @@ -142,7 +142,6 @@ def plot_2dfunc(u): if args.plot: plot_3dfunc(u) -#devito_output = u.data[:] print("After Operator 1: Devito norm:", norm(u)) print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) diff --git a/fast/wave3d.py b/fast/wave3d.py index 1497d3bea2..3a485fb3c7 100644 --- a/fast/wave3d.py +++ b/fast/wave3d.py @@ -93,7 +93,7 @@ def plot_3dfunc(u): # Define the wavefield with the size of the model and the time dimension u = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) - +# Another one to clone data u2 = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) # We can now write the PDE @@ -111,7 +111,7 @@ def plot_3dfunc(u): # the corresponding code print(time_range) -print("Init norm:", norm(u)) +print("Init norm:", np.linalg.norm(u.data[:])) src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) op0 = Operator([stencil] + src_term, subs=model.spacing_map, name='SourceDevitoOperator') # Run with source and plot @@ -124,8 +124,8 @@ def plot_3dfunc(u): print("Init linalg norm 0 :", np.linalg.norm(u.data[0])) print("Init linalg norm 1 :", np.linalg.norm(u.data[1])) print("Init linalg norm 2 :", np.linalg.norm(u.data[2])) - print("Norm of initial data:", np.linalg.norm(u.data[:])) + configuration['mpi'] = 0 u2.data[:] = u.data[:] configuration['mpi'] = 'basic' @@ -139,7 +139,7 @@ def plot_3dfunc(u): plot_3dfunc(u) #devito_output = u.data[:] -print("After Operator 1: Devito norm:", norm(u)) +print("After Operator 1: Devito norm:", np.linalg.norm(u.data[:])) print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) @@ -150,7 +150,7 @@ def plot_3dfunc(u): configuration['mpi'] = 'basic' #v[:, ..., :] = 1 -print("Reinitialise data: Devito norm:", np.linalg.norm(u.data[:])) +print("Reinitialise data for XDSL:", np.linalg.norm(u.data[:])) print("Init XDSL linalg norm:", np.linalg.norm(u.data[0])) print("Init XDSL linalg norm:", np.linalg.norm(u.data[1])) print("Init XDSL linalg norm:", np.linalg.norm(u.data[2])) From e54cd7b5a411030598048177c742398f04aa4103 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Wed, 2 Aug 2023 21:28:58 +0300 Subject: [PATCH 06/25] mpi: wip --- fast/wave3d.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/fast/wave3d.py b/fast/wave3d.py index 3a485fb3c7..4dd3647cd6 100644 --- a/fast/wave3d.py +++ b/fast/wave3d.py @@ -36,7 +36,7 @@ def plot_3dfunc(u): import pyvista as pv cmap = plt.colormaps["viridis"] values = u.data[0, :, :, :] - vistagrid = pv.UniformGrid() + vistagrid = pv.ImageData() vistagrid.dimensions = np.array(values.shape) + 1 vistagrid.spacing = (1, 1, 1) vistagrid.origin = (0, 0, 0) # The bottom left corner of the data set @@ -102,14 +102,14 @@ def plot_3dfunc(u): pde = u.dt2 - u.laplace # The PDE representation is as on paper -pde +# pde stencil = Eq(u.forward, solve(pde, u.forward)) stencil # Finally we define the source injection and receiver read function to generate # the corresponding code -print(time_range) +# print(time_range) print("Init norm:", np.linalg.norm(u.data[:])) src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) @@ -121,6 +121,8 @@ def plot_3dfunc(u): if args.plot: plot_3dfunc(u) +devito_norm = norm(u) +print("Init linalg norm 0 (inlined) :", norm(u)) print("Init linalg norm 0 :", np.linalg.norm(u.data[0])) print("Init linalg norm 1 :", np.linalg.norm(u.data[1])) print("Init linalg norm 2 :", np.linalg.norm(u.data[2])) @@ -151,18 +153,19 @@ def plot_3dfunc(u): #v[:, ..., :] = 1 print("Reinitialise data for XDSL:", np.linalg.norm(u.data[:])) -print("Init XDSL linalg norm:", np.linalg.norm(u.data[0])) -print("Init XDSL linalg norm:", np.linalg.norm(u.data[1])) -print("Init XDSL linalg norm:", np.linalg.norm(u.data[2])) +print("Init XDSL linalg norm 0:", np.linalg.norm(u.data[0])) +print("Init XDSL linalg norm 1:", np.linalg.norm(u.data[1])) +print("Init XDSL linalg norm 2:", np.linalg.norm(u.data[2])) # Run more with no sources now (Not supported in xdsl) -xdslop = XDSLOperator([stencil], name='xDSLOperator') +xdslop = Operator([stencil], name='Operator') xdslop.apply(time=time_range.num-1, dt=model.critical_dt) xdsl_output = u.copy() print("XDSL norm:", norm(u)) print(f"xdsl output norm: {norm(xdsl_output)}") -print("XDSL output linalg norm:", np.linalg.norm(u.data[0])) -print("XDSL output linalg norm:", np.linalg.norm(u.data[1])) -print("XDSL output linalg norm:", np.linalg.norm(u.data[2])) +print("XDSL output linalg norm 0:", np.linalg.norm(u.data[0])) +print("XDSL output linalg norm 1:", np.linalg.norm(u.data[1])) +print("XDSL output linalg norm 2:", np.linalg.norm(u.data[2])) +print("devito-norm:", devito_norm) \ No newline at end of file From 625e9764238784f398df015f72af926f02b36c78 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Thu, 3 Aug 2023 11:55:13 +0300 Subject: [PATCH 07/25] wave3d: cleanup --- fast/wave3d.py | 59 ++++++++++++++++++++++---------------------------- 1 file changed, 26 insertions(+), 33 deletions(-) diff --git a/fast/wave3d.py b/fast/wave3d.py index 4dd3647cd6..d094a0982b 100644 --- a/fast/wave3d.py +++ b/fast/wave3d.py @@ -88,45 +88,41 @@ def plot_3dfunc(u): # First, position source centrally in all dimensions, then set depth src.coordinates.data[0, :] = np.array(model.domain_size) * .5 -# We can plot the time signature to see the wavelet -# src.show() - # Define the wavefield with the size of the model and the time dimension u = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) # Another one to clone data u2 = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) +ub = TimeFunction(name="ub", grid=model.grid, time_order=to, space_order=so) + # We can now write the PDE # pde = model.m * u.dt2 - u.laplace + model.damp * u.dt # import pdb;pdb.set_trace() pde = u.dt2 - u.laplace -# The PDE representation is as on paper -# pde - stencil = Eq(u.forward, solve(pde, u.forward)) -stencil # Finally we define the source injection and receiver read function to generate # the corresponding code # print(time_range) -print("Init norm:", np.linalg.norm(u.data[:])) +# print("Init norm:", np.linalg.norm(u.data[:])) src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) op0 = Operator([stencil] + src_term, subs=model.spacing_map, name='SourceDevitoOperator') # Run with source and plot op0.apply(time=time_range.num-1, dt=model.critical_dt) + if len(shape) == 3: if args.plot: plot_3dfunc(u) -devito_norm = norm(u) -print("Init linalg norm 0 (inlined) :", norm(u)) -print("Init linalg norm 0 :", np.linalg.norm(u.data[0])) -print("Init linalg norm 1 :", np.linalg.norm(u.data[1])) -print("Init linalg norm 2 :", np.linalg.norm(u.data[2])) -print("Norm of initial data:", np.linalg.norm(u.data[:])) +# devito_norm = norm(u) +# print("Init linalg norm 0 (inlined) :", norm(u)) +# print("Init linalg norm 0 :", np.linalg.norm(u.data[0])) +# print("Init linalg norm 1 :", np.linalg.norm(u.data[1])) +# print("Init linalg norm 2 :", np.linalg.norm(u.data[2])) +# print("Norm of initial data:", np.linalg.norm(u.data[:])) configuration['mpi'] = 0 u2.data[:] = u.data[:] @@ -136,36 +132,33 @@ def plot_3dfunc(u): op1 = Operator([stencil], name='DevitoOperator') op1.apply(time=time_range.num-1, dt=model.critical_dt) +configuration['mpi'] = 0 +ub.data[:] = u.data[:] +configuration['mpi'] = 'basic' + if len(shape) == 3: if args.plot: plot_3dfunc(u) -#devito_output = u.data[:] -print("After Operator 1: Devito norm:", np.linalg.norm(u.data[:])) -print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) -print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) -print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) +# print("After Operator 1: Devito norm:", np.linalg.norm(u.data[:])) +#print("Devito norm 0:", np.linalg.norm(u.data[0])) +#print("Devito norm 1:", np.linalg.norm(u.data[1])) +#print("Devito norm 2:", np.linalg.norm(u.data[2])) # Reset initial data configuration['mpi'] = 0 u.data[:] = u2.data[:] configuration['mpi'] = 'basic' -#v[:, ..., :] = 1 -print("Reinitialise data for XDSL:", np.linalg.norm(u.data[:])) -print("Init XDSL linalg norm 0:", np.linalg.norm(u.data[0])) -print("Init XDSL linalg norm 1:", np.linalg.norm(u.data[1])) -print("Init XDSL linalg norm 2:", np.linalg.norm(u.data[2])) +# print("Reinitialise data for XDSL:", np.linalg.norm(u.data[:])) +# print("Init XDSL linalg norm 0:", np.linalg.norm(u.data[0])) +# print("Init XDSL linalg norm 1:", np.linalg.norm(u.data[1])) +# print("Init XDSL linalg norm 2:", np.linalg.norm(u.data[2])) # Run more with no sources now (Not supported in xdsl) -xdslop = Operator([stencil], name='Operator') +xdslop = XDSLOperator([stencil], name='XDSLOperator') xdslop.apply(time=time_range.num-1, dt=model.critical_dt) -xdsl_output = u.copy() - -print("XDSL norm:", norm(u)) -print(f"xdsl output norm: {norm(xdsl_output)}") -print("XDSL output linalg norm 0:", np.linalg.norm(u.data[0])) -print("XDSL output linalg norm 1:", np.linalg.norm(u.data[1])) -print("XDSL output linalg norm 2:", np.linalg.norm(u.data[2])) -print("devito-norm:", devito_norm) \ No newline at end of file +print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) +print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) +print("XDSL output norm 2:", np.linalg.norm(u.data[2]), "vs:", np.linalg.norm(ub.data[2])) From b73489a48caed6275e085b87b952a22558463b21 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Thu, 3 Aug 2023 12:03:14 +0300 Subject: [PATCH 08/25] wave2d: cleanup --- fast/wave2d.py | 52 ++++++++++++++++++++------------------------------ fast/wave3d.py | 5 ----- 2 files changed, 21 insertions(+), 36 deletions(-) diff --git a/fast/wave2d.py b/fast/wave2d.py index 9fb0758f9b..40bfb41020 100644 --- a/fast/wave2d.py +++ b/fast/wave2d.py @@ -1,7 +1,6 @@ # Based on the implementation of the Devito acoustic example implementation # Not using Devito's source injection abstraction import sys - import numpy as np from devito import (TimeFunction, Eq, Operator, solve, norm, @@ -97,21 +96,19 @@ def plot_2dfunc(u): u = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) # Another one to clone data u2 = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) +ub = TimeFunction(name="ub", grid=model.grid, time_order=to, space_order=so) # We can now write the PDE # pde = model.m * u.dt2 - u.laplace + model.damp * u.dt # import pdb;pdb.set_trace() pde = u.dt2 - u.laplace -# The PDE representation is as on paper -# pde - stencil = Eq(u.forward, solve(pde, u.forward)) # stencil # Finally we define the source injection and receiver read function to generate # the corresponding code -print(time_range) +# print(time_range) print("Init norm:", np.linalg.norm(u.data[:])) src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) @@ -124,11 +121,10 @@ def plot_2dfunc(u): if args.plot: plot_2dfunc(u) -print("Init Devito linalg norm 0 :", np.linalg.norm(u.data[0])) -print("Init Devito linalg norm 1 :", np.linalg.norm(u.data[1])) -print("Init Devito linalg norm 2 :", np.linalg.norm(u.data[2])) - -print("Norm of initial data:", norm(u)) +#print("Init Devito linalg norm 0 :", np.linalg.norm(u.data[0])) +#print("Init Devito linalg norm 1 :", np.linalg.norm(u.data[1])) +#print("Init Devito linalg norm 2 :", np.linalg.norm(u.data[2])) +# print("Norm of initial data:", norm(u)) configuration['mpi'] = 0 u2.data[:] = u.data[:] @@ -138,14 +134,14 @@ def plot_2dfunc(u): op1 = Operator([stencil], name='DevitoOperator') op1.apply(time=time_range.num-1, dt=model.critical_dt) -if len(shape) == 2: - if args.plot: - plot_3dfunc(u) +configuration['mpi'] = 0 +ub.data[:] = u.data[:] +configuration['mpi'] = 'basic' -print("After Operator 1: Devito norm:", norm(u)) -print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) -print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) -print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) +#print("After Operator 1: Devito norm:", norm(u)) +#print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) +#print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) +#print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) # import pdb;pdb.set_trace() @@ -157,21 +153,15 @@ def plot_2dfunc(u): #v[:, ..., :] = 1 -print("Reinitialise data: Devito norm:", norm(u)) -print("Init XDSL linalg norm:", np.linalg.norm(u.data[0])) -print("Init XDSL linalg norm:", np.linalg.norm(u.data[1])) -print("Init XDSL linalg norm:", np.linalg.norm(u.data[2])) +#print("Reinitialise data: Devito norm:", norm(u)) +#print("Init XDSL linalg norm:", np.linalg.norm(u.data[0])) +#print("Init XDSL linalg norm:", np.linalg.norm(u.data[1])) +#print("Init XDSL linalg norm:", np.linalg.norm(u.data[2])) # Run more with no sources now (Not supported in xdsl) -xdslop = XDSLOperator([stencil], name='xDSLOperator') +xdslop = Operator([stencil], name='xDSLOperator') xdslop.apply(time=time_range.num-1, dt=model.critical_dt) -xdsl_output = u.copy() -print("XDSL norm:", norm(u)) -print(f"xdsl output norm: {norm(xdsl_output)}") - -print("XDSL output linalg norm:", np.linalg.norm(u.data[0])) -print("XDSL output linalg norm:", np.linalg.norm(u.data[1])) -print("XDSL output linalg norm:", np.linalg.norm(u.data[2])) - - +print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) +print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) +print("XDSL output norm 2:", np.linalg.norm(u.data[2]), "vs:", np.linalg.norm(ub.data[2])) diff --git a/fast/wave3d.py b/fast/wave3d.py index d094a0982b..d7db2ebeba 100644 --- a/fast/wave3d.py +++ b/fast/wave3d.py @@ -102,11 +102,6 @@ def plot_3dfunc(u): stencil = Eq(u.forward, solve(pde, u.forward)) -# Finally we define the source injection and receiver read function to generate -# the corresponding code -# print(time_range) - -# print("Init norm:", np.linalg.norm(u.data[:])) src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) op0 = Operator([stencil] + src_term, subs=model.spacing_map, name='SourceDevitoOperator') # Run with source and plot From e4db7f9049db31aced3ef7bb3ccf3c5c775a8acc Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Thu, 3 Aug 2023 12:32:36 +0300 Subject: [PATCH 09/25] mpi-mfe: Add --- fast/mfe_2D.py | 73 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) create mode 100644 fast/mfe_2D.py diff --git a/fast/mfe_2D.py b/fast/mfe_2D.py new file mode 100644 index 0000000000..42e0ec213b --- /dev/null +++ b/fast/mfe_2D.py @@ -0,0 +1,73 @@ +# A 2D heat diffusion using Devito +# BC modelling included +# PyVista plotting included + +import argparse +import numpy as np + +from devito import Grid, TimeFunction, Eq, solve, Operator, Constant, norm, XDSLOperator +from examples.seismic import plot_image +from examples.cfd import init_hat + +parser = argparse.ArgumentParser(description='Process arguments.') + +parser.add_argument("-d", "--shape", default=(11, 11), type=int, nargs="+", + help="Number of grid points along each axis") +parser.add_argument("-so", "--space_order", default=2, + type=int, help="Space order of the simulation") +parser.add_argument("-to", "--time_order", default=1, + type=int, help="Time order of the simulation") +parser.add_argument("-nt", "--nt", default=40, + type=int, help="Simulation time in millisecond") +parser.add_argument("-bls", "--blevels", default=2, type=int, nargs="+", + help="Block levels") +parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot3D") +args = parser.parse_args() + +# Some variable declarations +nx, ny = args.shape +nt = args.nt +nu = .5 +dx = 1. / (nx - 1) +dy = 1. / (ny - 1) +sigma = .25 + +dt = sigma * dx * dy / nu +so = args.space_order +to = args.time_order + +print("dx %s, dy %s" % (dx, dy)) + +grid = Grid(shape=(nx, ny), extent=(2., 2.)) +u = TimeFunction(name='u', grid=grid, space_order=so) + +# Reset our data field and ICs +#init_hat(field=u.data[0], dx=dx, dy=dy, value=1.) +u.data[:, 2:3, 2:3] = 1 + +a = Constant(name='a') +# Create an equation with second-order derivatives +# eq = Eq(u.dt, a * u.laplace, subdomain=grid.interior) +eq = Eq(u.dt, a * u.laplace) +stencil = solve(eq, u.forward) +eq_stencil = Eq(u.forward, stencil) + +# Create boundary condition expressions +x, y = grid.dimensions +t = grid.stepping_dim + +initdata = u.data[:] +op = Operator([eq_stencil], name='DevitoOperator') +op.apply(time=nt, dt=dt, a=nu) +print(u.data[0, :]) +print("Devito Field norm is:", norm(u)) + +u.data[:, : , :] = 0 +u.data[:, 2:3 , 2:3] = 1 +# Reset data and run XDSLOperator +#init_hat(field=u.data[0], dx=dx, dy=dy, value=1.) +xdslop = Operator([eq_stencil], name='XDSLOperator') +xdslop.apply(time=nt, dt=dt, a=nu) +print(u.data[0, :]) + +print("XDSL Field norm is:", norm(u)) From ea7fe19f283aa47b1515899896b2519800981176 Mon Sep 17 00:00:00 2001 From: Anton Lydike Date: Thu, 3 Aug 2023 11:03:52 +0100 Subject: [PATCH 10/25] hacky fix for row major dmp.grid --- devito/operator/xdsl_operator.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/devito/operator/xdsl_operator.py b/devito/operator/xdsl_operator.py index 2ffe000863..1bf6ebda67 100644 --- a/devito/operator/xdsl_operator.py +++ b/devito/operator/xdsl_operator.py @@ -85,7 +85,10 @@ def _make_interop_o(self): @property def mpi_shape(self) -> tuple: dist = self.functions[0].grid.distributor - return dist.topology, dist.myrank + # temporary fix: + # swap dim 0 and 1 in topology because dmp.grid is row major and not column major + + return (dist.topology[1], dist.topology[0], *dist.topology[2:]), dist.myrank def _jit_compile(self): """ From d45c113f17d74cd42e14b5e1a8530bfde7fbc8ba Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Fri, 4 Aug 2023 17:18:05 +0300 Subject: [PATCH 11/25] bench: Conditional execution heat2d --- fast/diffusion_2D_wBCs.py | 26 +++++++++++++++----------- fast/wave2d.py | 1 + 2 files changed, 16 insertions(+), 11 deletions(-) diff --git a/fast/diffusion_2D_wBCs.py b/fast/diffusion_2D_wBCs.py index 0de0340254..5afe187445 100644 --- a/fast/diffusion_2D_wBCs.py +++ b/fast/diffusion_2D_wBCs.py @@ -21,7 +21,9 @@ type=int, help="Simulation time in millisecond") parser.add_argument("-bls", "--blevels", default=2, type=int, nargs="+", help="Block levels") -parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot3D") +parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot2D") +parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run") +parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run") args = parser.parse_args() # Some variable declarations @@ -55,15 +57,17 @@ x, y = grid.dimensions t = grid.stepping_dim -initdata = u.data[:] -op = Operator([eq_stencil], name='DevitoOperator') -op.apply(time=nt, dt=dt, a=nu) +import pdb;pdb.set_trace() -print("Devito Field norm is:", norm(u)) +# initdata = u.data[:] +if args.devito: + op = Operator([eq_stencil], name='DevitoOperator') + op.apply(time=nt, dt=dt, a=nu) + print("Devito Field norm is:", norm(u)) -# Reset data and run XDSLOperator -init_hat(field=u.data[0], dx=dx, dy=dy, value=1.) -xdslop = XDSLOperator([eq_stencil], name='XDSLOperator') -xdslop.apply(time=nt, dt=dt, a=nu) - -print("XDSL Field norm is:", norm(u)) +if args.xdsl: + # Reset data and run XDSLOperator + #init_hat(field=u.data[0], dx=dx, dy=dy, value=1.) + xdslop = XDSLOperator([eq_stencil], name='XDSLOperator') + xdslop.apply(time=nt, dt=dt, a=nu) + print("XDSL Field norm is:", norm(u)) diff --git a/fast/wave2d.py b/fast/wave2d.py index 40bfb41020..47d08a5d52 100644 --- a/fast/wave2d.py +++ b/fast/wave2d.py @@ -27,6 +27,7 @@ parser.add_argument("-bls", "--blevels", default=1, type=int, nargs="+", help="Block levels") parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot2D") +parser.add_argument("-mode", "--mode", default='devito', type=str, help="Operator mode") args = parser.parse_args() From ae3d586dbbb74fe6436f447643085632e372aa2a Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Fri, 4 Aug 2023 17:22:26 +0300 Subject: [PATCH 12/25] bench: Conditional execution heat3d --- fast/diffusion_2D_wBCs.py | 12 +++--------- fast/diffusion_3D_wBCs.py | 38 ++++++++++++++++++++------------------ 2 files changed, 23 insertions(+), 27 deletions(-) diff --git a/fast/diffusion_2D_wBCs.py b/fast/diffusion_2D_wBCs.py index 5afe187445..b5f850db85 100644 --- a/fast/diffusion_2D_wBCs.py +++ b/fast/diffusion_2D_wBCs.py @@ -53,21 +53,15 @@ stencil = solve(eq, u.forward) eq_stencil = Eq(u.forward, stencil) -# Create boundary condition expressions -x, y = grid.dimensions -t = grid.stepping_dim - -import pdb;pdb.set_trace() - -# initdata = u.data[:] if args.devito: op = Operator([eq_stencil], name='DevitoOperator') op.apply(time=nt, dt=dt, a=nu) print("Devito Field norm is:", norm(u)) +# Reset data +init_hat(field=u.data[0], dx=dx, dy=dy, value=1.) + if args.xdsl: - # Reset data and run XDSLOperator - #init_hat(field=u.data[0], dx=dx, dy=dy, value=1.) xdslop = XDSLOperator([eq_stencil], name='XDSLOperator') xdslop.apply(time=nt, dt=dt, a=nu) print("XDSL Field norm is:", norm(u)) diff --git a/fast/diffusion_3D_wBCs.py b/fast/diffusion_3D_wBCs.py index eebd13a7c5..9cb39d5179 100644 --- a/fast/diffusion_3D_wBCs.py +++ b/fast/diffusion_3D_wBCs.py @@ -20,6 +20,8 @@ parser.add_argument("-bls", "--blevels", default=2, type=int, nargs="+", help="Block levels") parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot3D") +parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run") +parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run") args = parser.parse_args() @@ -77,21 +79,21 @@ def plot_3dfunc(u): print(eq_stencil) # Create Operator -op = Operator([eq_stencil], name='DevitoOperator') -# Apply the operator for a number of timesteps -op.apply(time=nt, dt=dt, a=nu) -print("Devito Field norm is:", norm(u)) - -# Reset field -u.data[:, :, :, :] = 0 -u.data[:, :, :, int(nz/2)] = 1 -xdslop = XDSLOperator([eq_stencil], name='xDSLOperator') -# Apply the xdsl operator for a number of timesteps -xdslop.apply(time=nt, dt=dt, a=nu) - -if args.plot: - plot_3dfunc(u) - -print("XDSL Field norm is:", norm(u)) - -# import pdb;pdb.set_trace() +if args.devito: + op = Operator([eq_stencil], name='DevitoOperator') + # Apply the operator for a number of timesteps + op.apply(time=nt, dt=dt, a=nu) + print("Devito Field norm is:", norm(u)) + +if args.xdsl: + # Reset field + u.data[:, :, :, :] = 0 + u.data[:, :, :, int(nz/2)] = 1 + xdslop = XDSLOperator([eq_stencil], name='xDSLOperator') + # Apply the xdsl operator for a number of timesteps + xdslop.apply(time=nt, dt=dt, a=nu) + + if args.plot: + plot_3dfunc(u) + + print("XDSL Field norm is:", norm(u)) From 9a973cc2e2edb1202f2db69aee01f50701df9efd Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Fri, 4 Aug 2023 18:15:27 +0300 Subject: [PATCH 13/25] bench: Generalize benchmarking scripts --- fast/bench_utils.py | 8 ++++ fast/diffusion_2D_wBCs.py | 5 ++- fast/wave2d.py | 90 +++++++++++++++++---------------------- 3 files changed, 52 insertions(+), 51 deletions(-) create mode 100644 fast/bench_utils.py diff --git a/fast/bench_utils.py b/fast/bench_utils.py new file mode 100644 index 0000000000..fff49fefdb --- /dev/null +++ b/fast/bench_utils.py @@ -0,0 +1,8 @@ +from examples.seismic import plot_image + +__all__ = ['plot_2dfunc'] + + +def plot_2dfunc(u): + # Plot a 3D structured grid using pyvista + plot_image(u.data[0], cmap='seismic') diff --git a/fast/diffusion_2D_wBCs.py b/fast/diffusion_2D_wBCs.py index b5f850db85..869d6d2e7d 100644 --- a/fast/diffusion_2D_wBCs.py +++ b/fast/diffusion_2D_wBCs.py @@ -6,8 +6,8 @@ import numpy as np from devito import Grid, TimeFunction, Eq, solve, Operator, Constant, norm, XDSLOperator -from examples.seismic import plot_image from examples.cfd import init_hat +from fast.bench_utils import plot_2dfunc parser = argparse.ArgumentParser(description='Process arguments.') @@ -58,6 +58,9 @@ op.apply(time=nt, dt=dt, a=nu) print("Devito Field norm is:", norm(u)) + if args.plot: + plot_2dfunc(u) + # Reset data init_hat(field=u.data[0], dx=dx, dy=dy, value=1.) diff --git a/fast/wave2d.py b/fast/wave2d.py index 47d08a5d52..589802ba8c 100644 --- a/fast/wave2d.py +++ b/fast/wave2d.py @@ -6,8 +6,8 @@ from devito import (TimeFunction, Eq, Operator, solve, norm, XDSLOperator, configuration) from examples.seismic import RickerSource -from examples.seismic import Model, TimeAxis - +from examples.seismic import Model, TimeAxis, plot_image +from fast.bench_utils import plot_2dfunc from devito.tools import as_tuple import argparse @@ -27,26 +27,12 @@ parser.add_argument("-bls", "--blevels", default=1, type=int, nargs="+", help="Block levels") parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot2D") -parser.add_argument("-mode", "--mode", default='devito', type=str, help="Operator mode") +parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run") +parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run") args = parser.parse_args() -def plot_2dfunc(u): - # Plot a 3D structured grid using pyvista - - import matplotlib.pyplot as plt - import pyvista as pv - cmap = plt.colormaps["viridis"] - values = u.data[0, :, :, :] - vistagrid = pv.UniformGrid() - vistagrid.dimensions = np.array(values.shape) + 1 - vistagrid.spacing = (1, 1, 1) - vistagrid.origin = (0, 0, 0) # The bottom left corner of the data set - vistagrid.cell_data["values"] = values.flatten(order="F") - vistaslices = vistagrid.slice_orthogonal() - # vistagrid.plot(show_edges=True) - vistaslices.plot(cmap=cmap) - +mpiconf = configuration['mpi'] # Define a physical size # nx, ny, nz = args.shape @@ -122,47 +108,51 @@ def plot_2dfunc(u): if args.plot: plot_2dfunc(u) -#print("Init Devito linalg norm 0 :", np.linalg.norm(u.data[0])) -#print("Init Devito linalg norm 1 :", np.linalg.norm(u.data[1])) -#print("Init Devito linalg norm 2 :", np.linalg.norm(u.data[2])) +# print("Init Devito linalg norm 0 :", np.linalg.norm(u.data[0])) +# print("Init Devito linalg norm 1 :", np.linalg.norm(u.data[1])) +# print("Init Devito linalg norm 2 :", np.linalg.norm(u.data[2])) # print("Norm of initial data:", norm(u)) configuration['mpi'] = 0 u2.data[:] = u.data[:] -configuration['mpi'] = 'basic' +configuration['mpi'] = mpiconf -# Run more with no sources now (Not supported in xdsl) -op1 = Operator([stencil], name='DevitoOperator') -op1.apply(time=time_range.num-1, dt=model.critical_dt) - -configuration['mpi'] = 0 -ub.data[:] = u.data[:] -configuration['mpi'] = 'basic' +if args.devito: + # Run more with no sources now (Not supported in xdsl) + op1 = Operator([stencil], name='DevitoOperator') + op1.apply(time=time_range.num-1, dt=model.critical_dt) -#print("After Operator 1: Devito norm:", norm(u)) -#print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) -#print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) -#print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) + configuration['mpi'] = 0 + ub.data[:] = u.data[:] + configuration['mpi'] = mpiconf -# import pdb;pdb.set_trace() + if len(shape) == 2 and args.plot: + plot_2dfunc(u) + # print("After Operator 1: Devito norm:", norm(u)) + # print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) + # print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) + # print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) -# Reset initial data -configuration['mpi'] = 0 -u.data[:] = u2.data[:] -configuration['mpi'] = 'basic' -#v[:, ..., :] = 1 +if args.xdsl: + # Reset initial data + configuration['mpi'] = 0 + u.data[:] = u2.data[:] + configuration['mpi'] = mpiconf + # v[:, ..., :] = 1 + # print("Reinitialise data: Devito norm:", norm(u)) + # print("XDSL init linalg norm:", np.linalg.norm(u.data[0])) + # print("XDSL init linalg norm:", np.linalg.norm(u.data[1])) + # print("XDSL init linalg norm:", np.linalg.norm(u.data[2])) -#print("Reinitialise data: Devito norm:", norm(u)) -#print("Init XDSL linalg norm:", np.linalg.norm(u.data[0])) -#print("Init XDSL linalg norm:", np.linalg.norm(u.data[1])) -#print("Init XDSL linalg norm:", np.linalg.norm(u.data[2])) + # Run more with no sources now (Not supported in xdsl) + xdslop = Operator([stencil], name='xDSLOperator') + xdslop.apply(time=time_range.num-1, dt=model.critical_dt) -# Run more with no sources now (Not supported in xdsl) -xdslop = Operator([stencil], name='xDSLOperator') -xdslop.apply(time=time_range.num-1, dt=model.critical_dt) + if len(shape) == 2 and args.plot: + plot_2dfunc(u) -print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) -print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) -print("XDSL output norm 2:", np.linalg.norm(u.data[2]), "vs:", np.linalg.norm(ub.data[2])) + print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) + print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) + print("XDSL output norm 2:", np.linalg.norm(u.data[2]), "vs:", np.linalg.norm(ub.data[2])) From 4b34fc5b05f6b0894b99d57ee11aa65e22990db9 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Fri, 4 Aug 2023 18:24:03 +0300 Subject: [PATCH 14/25] bench: Generalize wave3d --- fast/bench_utils.py | 20 ++++++++++++++++++-- fast/wave3d.py | 25 +++++++------------------ 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/fast/bench_utils.py b/fast/bench_utils.py index fff49fefdb..1881a61806 100644 --- a/fast/bench_utils.py +++ b/fast/bench_utils.py @@ -1,8 +1,24 @@ +import numpy as np from examples.seismic import plot_image -__all__ = ['plot_2dfunc'] +__all__ = ['plot_2dfunc', 'plot_3dfunc'] def plot_2dfunc(u): - # Plot a 3D structured grid using pyvista + # Plot a 2D image using devito's machinery plot_image(u.data[0], cmap='seismic') + + +def plot_3dfunc(u): + # Plot a 3D structured grid using pyvista + import matplotlib.pyplot as plt + import pyvista as pv + cmap = plt.colormaps["viridis"] + values = u.data[0, :, :, :] + vistagrid = pv.ImageData() + vistagrid.dimensions = np.array(values.shape) + 1 + vistagrid.spacing = (1, 1, 1) + vistagrid.origin = (0, 0, 0) # The bottom left corner of the data set + vistagrid.cell_data["values"] = values.flatten(order="F") + vistaslices = vistagrid.slice_orthogonal() + vistaslices.plot(cmap=cmap) diff --git a/fast/wave3d.py b/fast/wave3d.py index d7db2ebeba..86b0d55afe 100644 --- a/fast/wave3d.py +++ b/fast/wave3d.py @@ -6,6 +6,7 @@ XDSLOperator, configuration) from examples.seismic import RickerSource from examples.seismic import Model, TimeAxis +from fast.bench_utils import plot_3dfunc from devito.tools import as_tuple @@ -26,24 +27,12 @@ parser.add_argument("-bls", "--blevels", default=2, type=int, nargs="+", help="Block levels") parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot3D") +parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run") +parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run") args = parser.parse_args() -def plot_3dfunc(u): - # Plot a 3D structured grid using pyvista - - import matplotlib.pyplot as plt - import pyvista as pv - cmap = plt.colormaps["viridis"] - values = u.data[0, :, :, :] - vistagrid = pv.ImageData() - vistagrid.dimensions = np.array(values.shape) + 1 - vistagrid.spacing = (1, 1, 1) - vistagrid.origin = (0, 0, 0) # The bottom left corner of the data set - vistagrid.cell_data["values"] = values.flatten(order="F") - vistaslices = vistagrid.slice_orthogonal() - # vistagrid.plot(show_edges=True) - vistaslices.plot(cmap=cmap) +mpiconf = configuration['mpi'] # Define a physical size @@ -121,7 +110,7 @@ def plot_3dfunc(u): configuration['mpi'] = 0 u2.data[:] = u.data[:] -configuration['mpi'] = 'basic' +configuration['mpi'] = mpiconf # Run more with no sources now (Not supported in xdsl) op1 = Operator([stencil], name='DevitoOperator') @@ -129,7 +118,7 @@ def plot_3dfunc(u): configuration['mpi'] = 0 ub.data[:] = u.data[:] -configuration['mpi'] = 'basic' +configuration['mpi'] = mpiconf if len(shape) == 3: if args.plot: @@ -143,7 +132,7 @@ def plot_3dfunc(u): # Reset initial data configuration['mpi'] = 0 u.data[:] = u2.data[:] -configuration['mpi'] = 'basic' +configuration['mpi'] = mpiconf # print("Reinitialise data for XDSL:", np.linalg.norm(u.data[:])) # print("Init XDSL linalg norm 0:", np.linalg.norm(u.data[0])) From 7d7e6394029391e63dd350ae666f92d18643c2c8 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Fri, 4 Aug 2023 19:50:26 +0300 Subject: [PATCH 15/25] wave: TryAdd example with no Operator --- fast/bench_utils.py | 1 + fast/wave2d.py | 1 + fast/wave2d_b.py | 168 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 170 insertions(+) create mode 100644 fast/wave2d_b.py diff --git a/fast/bench_utils.py b/fast/bench_utils.py index 1881a61806..fe1563bae5 100644 --- a/fast/bench_utils.py +++ b/fast/bench_utils.py @@ -7,6 +7,7 @@ def plot_2dfunc(u): # Plot a 2D image using devito's machinery plot_image(u.data[0], cmap='seismic') + plot_image(u.data[1], cmap='seismic') def plot_3dfunc(u): diff --git a/fast/wave2d.py b/fast/wave2d.py index 589802ba8c..2a7855a6e9 100644 --- a/fast/wave2d.py +++ b/fast/wave2d.py @@ -108,6 +108,7 @@ if args.plot: plot_2dfunc(u) +import pdb;pdb.set_trace() # print("Init Devito linalg norm 0 :", np.linalg.norm(u.data[0])) # print("Init Devito linalg norm 1 :", np.linalg.norm(u.data[1])) # print("Init Devito linalg norm 2 :", np.linalg.norm(u.data[2])) diff --git a/fast/wave2d_b.py b/fast/wave2d_b.py new file mode 100644 index 0000000000..aada12f24a --- /dev/null +++ b/fast/wave2d_b.py @@ -0,0 +1,168 @@ +# Based on the implementation of the Devito acoustic example implementation +# Not using Devito's source injection abstraction +import sys +import numpy as np + +from devito import (TimeFunction, Eq, Operator, solve, norm, + XDSLOperator, configuration, Grid) +from examples.seismic import RickerSource +from examples.seismic import Model, TimeAxis, plot_image +from fast.bench_utils import plot_2dfunc +from devito.tools import as_tuple + +import argparse +np.set_printoptions(threshold=np.inf) + + +parser = argparse.ArgumentParser(description='Process arguments.') + +parser.add_argument("-d", "--shape", default=(16, 16), type=int, nargs="+", + help="Number of grid points along each axis") +parser.add_argument("-so", "--space_order", default=4, + type=int, help="Space order of the simulation") +parser.add_argument("-to", "--time_order", default=2, + type=int, help="Time order of the simulation") +parser.add_argument("-nt", "--nt", default=20, + type=int, help="Simulation time in millisecond") +parser.add_argument("-bls", "--blevels", default=1, type=int, nargs="+", + help="Block levels") +parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot2D") +parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run") +parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run") +args = parser.parse_args() + + +mpiconf = configuration['mpi'] + +# Define a physical size +# nx, ny, nz = args.shape +nt = args.nt + +shape = (args.shape) # Number of grid point (nx, ny, nz) +spacing = as_tuple(10.0 for _ in range(len(shape))) # Grid spacing in m. The domain size is now 1km by 1km +origin = as_tuple(0.0 for _ in range(len(shape))) # What is the location of the top left corner. +domain_size = tuple((d-1) * s for d, s in zip(shape, spacing)) +grid = Grid(shape=shape, extent=(1000, 1000)) +# grid = Grid(shape=shape) +# This is necessary to define +# the absolute location of the source and receivers + +# Define a velocity profile. The velocity is in km/s +v = np.empty(shape, dtype=np.float32) +v[:, :] = 1 + +# With the velocity and model size defined, we can create the seismic model that +# encapsulates this properties. We also define the size of the absorbing layer as +# 10 grid points +so = args.space_order +to = args.time_order + +#model = Model(vp=v, origin=origin, shape=shape, spacing=spacing, +# space_order=so, nbl=0) + +# plot_velocity(model) + +t0 = 0. # Simulation starts a t=0 +tn = nt # Simulation last 1 second (1000 ms) +# dt = model.critical_dt # Time step from model grid spacing +# print("dt is:", dt) + +#time_range = TimeAxis(start=t0, stop=tn, step=dt) + +# The source is positioned at a $20m$ depth and at the middle of the +# $x$ axis ($x_{src}=500m$), +# with a peak wavelet frequency of $10Hz$. +f0 = 0.010 # Source peak frequency is 10Hz (0.010 kHz) +#src = RickerSource(name='src', grid=model.grid, f0=f0, +# npoint=1, time_range=time_range) + +# First, position source centrally in all dimensions, then set depth +#src.coordinates.data[0, :] = np.array(model.domain_size) * .5 + +# We can plot the time signature to see the wavelet +# src.show() + +# Define the wavefield with the size of the model and the time dimension +u = TimeFunction(name="u", grid=grid, time_order=to, space_order=so) +# Another one to clone data +u2 = TimeFunction(name="u", grid=grid, time_order=to, space_order=so) +ub = TimeFunction(name="ub", grid=grid, time_order=to, space_order=so) + +# We can now write the PDE +# pde = model.m * u.dt2 - u.laplace + model.damp * u.dt +# import pdb;pdb.set_trace() +pde = u.dt2 - u.laplace + +stencil = Eq(u.forward, solve(pde, u.forward)) +# stencil + +# Finally we define the source injection and receiver read function to generate +# the corresponding code +# print(time_range) + +#src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) +#op0 = Operator([stencil] + src_term, subs=model.spacing_map, name='SourceDevitoOperator') + +# Run with source and plot +#op0.apply(time=time_range.num-1, dt=model.critical_dt) + + + +# print("Init Devito linalg norm 0 :", np.linalg.norm(u.data[0])) +# print("Init Devito linalg norm 1 :", np.linalg.norm(u.data[1])) +# print("Init Devito linalg norm 2 :", np.linalg.norm(u.data[2])) +# print("Norm of initial data:", norm(u)) + +configuration['mpi'] = 0 +u2.data[:] = u.data[:] +configuration['mpi'] = mpiconf + +import pdb;pdb.set_trace() +u.data[:] = np.load("wave_dat2.npy", allow_pickle=True) + +if len(shape) == 2: + if args.plot: + plot_2dfunc(u) + +print("Init norm:", np.linalg.norm(u.data[:])) + +if args.devito: + # Run more with no sources now (Not supported in xdsl) + #op1 = Operator([stencil], name='DevitoOperator', subs=grid.spacing_map) + op1 = Operator([stencil], name='DevitoOperator') + op1.apply(time=10, dt=5.54) + + configuration['mpi'] = 0 + ub.data[:] = u.data[:] + configuration['mpi'] = mpiconf + + if len(shape) == 2 and args.plot: + plot_2dfunc(u) + + print("After Operator 1: Devito norm:", norm(u)) + # print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) + # print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) + # print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) + + +if args.xdsl: + # Reset initial data + configuration['mpi'] = 0 + u.data[:] = u2.data[:] + configuration['mpi'] = mpiconf + # v[:, ..., :] = 1 + # print("Reinitialise data: Devito norm:", norm(u)) + # print("XDSL init linalg norm:", np.linalg.norm(u.data[0])) + # print("XDSL init linalg norm:", np.linalg.norm(u.data[1])) + # print("XDSL init linalg norm:", np.linalg.norm(u.data[2])) + + # Run more with no sources now (Not supported in xdsl) + xdslop = Operator([stencil], name='xDSLOperator') + xdslop.apply(time=nt) + + if len(shape) == 2 and args.plot: + plot_2dfunc(u) + + print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) + print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) + print("XDSL output norm 2:", np.linalg.norm(u.data[2]), "vs:", np.linalg.norm(ub.data[2])) From 26abd2ed6a0912f780147d92af7943dfba14d3ef Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Fri, 4 Aug 2023 19:55:10 +0300 Subject: [PATCH 16/25] add datatest --- fast/wave_dat2.npy | Bin 0 -> 120128 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 fast/wave_dat2.npy diff --git a/fast/wave_dat2.npy b/fast/wave_dat2.npy new file mode 100644 index 0000000000000000000000000000000000000000..22e536f687c92d1691fa491c08b1d413dd02744c GIT binary patch literal 120128 zcmeI52Y43MwuWg^r1u&K36O*k5)w!tbg2?T6GT9oq9O_i2!hgk@4Y2~6d(-}Qc3x# z(o3icf>NZ3^dbn{@6E}7p2z|3(E|v{+~=8l#+kj>TC;a%@3p?~o1Zx?+IhC>kTJuc z4By!F_-Md~{cYSWZEAg7)uxiA&BuNFf7!o#pRRrT_xOnGjk|w7;3Hih(6jp&AL+ks zT`E~tadJxje|D8Dhgtsh|9dlo@qfk$7$abefH4Bb2pA(^jDRr$#t0ZAV2pq<0>%g! zBVdexF#^U27$abefH4Bb2pA(^jDRr$#t0ZA@TVhCD!9$;vO4~>k#X~-jX>Y;TNmlC zE0-{d0_z)r@^~%>S|BHXtMa67~4=lZEoQ^kjUQgF{>)w8~>7usu zxuw2GWF!5K;PZX*o}p7p29BJxJgWEmWzQAdzp`1;-^+*8+!}E4hNIo&Rt47tZMU`C zp<{;5ujpEs?zK>xpVc--edvpg=Nnb}v$kkWw}l?7CztHK{PA)hyB%eJ`lOYwTZdT{ zOe@E%TJvtvipRdFA9(VIDEp4JW><9AaaQLibS+5tHmi-h+MZS)`u3F#PuW@`o9S;P zpU)F}&n;eR;GEiCkzd(Z*+eJ$j+vd`W?;XftLysjt9WaLLr~MKg&i9H_CwIYbuEKp zb$IIBL)Y%;UIVrHN^Qrg4}GI$V~}k5$ma8jF8x@EzP@MSN8c4IRd&HtuPu8fJ`>kwX3!C|5 zyODhSZ&9ZG-d#%+TohM!M5*enFL{-!cgJeyrvZMW7h2l-+~~B%FS1>w4R!K`%=I7c zxbN-ZYqNb3;uPJvP?a4Q3pqA#e>}92j#!;@Ewk<&R2$mHsSkazF-o?uiR~Klktko? z zYp{I&y=c>JuQkKUFA8<`E%l2}5$|@H9jtd0&g4IC`()eZF-O-_?tH6KX4`K=@^|a* zRAAhSb)#o~<{UKdnp4N+k>Oo>c17wWUu1^~DCZ zJY{pdY~ur8zshI0douwt2!@eDc`2{Jo)JgU5Wee(SyutNz-(dNrStBdQ&IIHTH89m{lnN7tft z?@P6@R9pJoR$pvj3!B((EFXvDtBHJi$TxHH*P>2)t=TVFqj)cvixEAF_{SuL*m}kn z4$Ssm9f!oc&Y|n_UUjm|&>^f$mL4uy!up1fn0a@7hxi56a(NVS?HkvqT7l!!Bihz= zkL;+Uu+HDowZgjBPHnzX+c(swwE9+*jd0o8D4W=BAs?~wg-_;$_k^`8J?|&%-Ru{v z8@w0H#m$1H0^0|+bZA-EGxT7Kf>p|XGb=1OL$J$?b6M7x>d-7Ar|;XYBZhZ&9aN`8 z)Sg;vqsC3n9KH2m;pAhq&UfqDNZnhaHd)oSxcW?2-$2>0maRjwiEVt~i#g#vAqMsn z_V~2-m;dKaM14U{OCH7E&05YH#e2axCDmvjdZ&eTm0vOjhxwg3?K13O$Mv6tO^FDo zve~uQ&MdCW+)G3s+4ynvoFKR896P>@YJaJw8^^=#Uq+4AwPU&$pf-!u)<=ClR^M-A z!&kPh$mU(yJ|rLb;yq#QA_n$&@-A}W^i-cv>rr2jcalf3ceAFmj_|k3TsSt43hO>4 zW7YGvv)13Nw>o0)q@r#~!?Q;=tgs^5zKmz|o70NaI9=mFjrm6$+^xE|b1$o7fzD6r znuYFFP#aIR9i=|I)c3k7d#3Kgt|Mq z)>AD&{z)FiKFfZ=8pV6TT)1}eTEFjnOvJa&Z@H~2{U~z#;QP_Px3{aFt+Ty*P0v90 z!BO|!tK=wBbN-iQYcAHo`5;|$*1h~{lj%LX>PhM|LVdT$2DYloCUe4j!rH}tf^YIJ zVxc}EHtKk=m9&g^fenI3fq_sfQ|nO^kcX1{uy?axutxD-Fn)y#w{bH~YZ^7-pk+*8 zlgc$3XYhBAYI@JTf!$vBI}OLz*tPuo*eUs%Y*?rxTIY{+t(oqvP#aIR4N{+}>gy{T z`DF{6XJnhTi~WQ>9^b@3EXg&aVgp+N?*fB-{ojN?J^}m%csDpL*dTZm7zi~pbsDt+ z`6zi5dpG+9YZUJVb8)PJW6aH#qiQTUa@4)?jaBZAcQ@N$-?~ZM#x_Uee1Eo!KbQ1Q z{Es?(be^DV<#n&7+E}ZtpZeTZ-xIPCAX{5y6WfF2gFT+S3*W>*EW`x10^S7%2`)_R zX^B69OM-&{_XcYP8w8I6$Dr<})}y{4rzMYK?`FSXjpDsvE`BX{&i&ZcmNir<#`fb4khmG$gT+cWr;i?WhooCTC?wwbgx7GITUt89=tG;(+W3p^bm(4x0 z-BUg$$rrgW^$9T$3o(Irfk6^0F%vtS)*t`MUcdi>?}jghqk~I=zW_4_rv)1Xj{*as zhNf<#z96S1k7Dm;zhI5xycs;lz_zN&|a9Z#_FeWe%>TYU1>I-sO@+kIh_6yc1-V5g9 zS`ViL$M23O+%4TX@xup^iQO*!o_OqbLDNwkJ#^kf*SKe|Hn-K5IpIBVlnoErdMKM7 zvOQHkvdh;n`OG8V#4ug45Yt)3_E0gxpTKFsGs6AC$D@&W-QEoJQD}0|fuKEr?}lTA z*Mm!fzW_4_rv)1Xj{*as?xxnGz96S1k7Dm;zhI5xyhWmw&MwB%9j-RudhQT**Q7tQO)^g*=Q$QuCiH3whzn4aQP}EpE2@H3=b4bCdD*Eu@NIYBit{1JQ@jfD`<(( zccHm^)!!5Jv}icd4Wd;;AB83d9SGV3_-;5>I6QbGxC$_Ha9XfI@F*}4>TYU1>I-sO z@+kIh_6F7{-V5epe(X8Zrp|XwpH%qS)J;c%&KK+25Z$v;8`{z*v-+B3W4>$+m(3Hh zJzGAu%U2!wyd~csD28~&vOzJORc!Ao#Z`05M%#>@ z77Zu5L9}Yw0xe!&{W zd%;{3-4SCdSm{`zwT==xXHG11&sJ^TS6lY@O6rS^O|s=8n+s*zRzC2RMLrkFcSXfO zEE5$|pkm9U7@I2AQ;Kj_281=FTl*fX~71;qQF3?yQ%f4FUV=hqu9IIFIb~^FPMw% z-gfnO>&UG0p1Q_83$;0{w%gR_n)+UljYYDBO>Eba4}5i&&tvj!QVfF?OJT)yLa`C! ze#P2GF~=zOQF?|}dX|oQCbX&Oh0$=MtLIyXSLyo`Vj(7CBSvB+W;C4W2GOdak3y4! z4g~E1^8v>SuLqX|e*tE0o)7RSFc9i)YCY-;a$52z_HOnI)+pW!=7J-O&bh{Xlu?_5 zPnnNJDf7`#^RZp?;US;=Eo@c{#PXS9x~SMLDaMUxR#MGJ3q4DFJyR|{+tycjKFsec z?=J5!?=tT+?=~7vbc1Ns&_|)kK?lNm1>X(F3a;4i?;!D+z;!K1)HsJq#lsV~TB$)ni2nG@D1-V5epo~{khz0zvaMQx9& zPbT%+~|`8Xh7@5(2>KT-_j7UUS0w4W0rxd;3CCB zEqN4sH|rAz?*((wRQJep`E`s7q!JK1{6ujh78p2D=^8xl|#uLHhz~;c{!0N#4lC^57@IbvOG7dptQU2YWaB1#1-V1#?kIedej}D%tRpt-Z22E(K34DPQ;m z%PjZ=o;WZK;)&*Y0jo}pw}WMaX@hN}RYM;I<_-3Z_5i*cjul=HE(!hu%p9DS_XIqO z{e-%kT95jIgFK48oBe_{iuZ!KxUarfWn-FbO_0sqvJIvWwhzXCN4^UxhV5w_Z#T~e zTn?NL+zuQMTo0VjlYA815F8QO1Nd$@R(L(QB=`$3b8uR)LGUOr5bADfJq~hO@+kIh z_6yc1-V5g9q-=!B);8I^CEG*g1MbUSK1V)-kBUj7e3W@!;6C9%;X>hS;YL&QwQ#3! zsBo!h58zheSmE{HlHf1E%)x2F2En7iK&ZPps4vKA$)ni2*)Le5crTcXg0gi$HqXg6 zoG`xNh-aqowFSl3!XfWY<9w}oKF|iB5kM=DTF-!%08Ihf0yG9_4bU8*Jz!4YSmE{H zlHf1E%)x2F2En7iK&ZQ^^}ry>Y00D5yV);Tqj)cv3p5hgzA7Jw<;#?!XK0}q9!k%! zOL~TH(utnIJTGWz(A1!I-sO@+kIh_6yc1-V1D2k?ngadY>H9`_z`+XF@vE`(Q2+mOP5R8(XYVX1(Yw=|$t@^RRrsn~wFO=J`OojD{I4Gn!_!&CCh! z3A#b{6ZUxWE^=YC2k_l+tnhkpN$?k7=HRqogWyqMAk^K|dej%>wB%9j-Ru{vQM?z- z#VA8hZJv&Q|1Vm9H2=&A?+I%c`w4rz`FjKE6LcVG58%7uSmE{HlHf1E%)x2F2Dt|Y zLfuWRM}0v~OCE(y_6yc1-V5V90QvBK+_zgGh@2d4!aG=Hx~-Hi?E3vycHd$s=<`lnuCKVgq2?;;mQ z+f2sgq}HRpz_#%{=ToQPsTawG zsZWw?M%Bt-E8tyVkZ9Gw!_nlR13`NL-wnqKuLqX|e?c4CrvBcUx|>?h_}=-c!~V>R z)F;%8)bU^|;9X#l;KJwz(W=20p~*o9g7yHu8;%uTk9+VJVCLYo#_s_I_&u5b7}Do*0jwVW1WpT{5$+d09^D{XHS|#& z=s?gOaL@R?rROsAPc;I60;h%MAMO`E9*qRL6?B7W)i}`PpabC^e7Es?WKSK==XMe9 z7d{@11iBTpM9CUX>D17wp^rk7!#%VI#_#n#x1oQkA^3PS66jXY5~1%xbBBf#-5^>u zkA$BXBtP}h2{<&B-&K;!suy@-$Q=pkiC=} zXj9P(qv1AwFZ`tp`JXyLFN}s8T|I~Kd-{JmTrce=-#Qro-o#5A?*G>j-=Dln-^2XB z0UHgB5imx;7y)Ahj1e$Kz!(8z1dI_dM!*;WV+4#5Fh;-_0b>M=5imx;7y)Ahj1e$K zz!(8z1dI_dM!*;WV+4#5Fh;-_0b>M=5imx;7y)Ahj1e$Kz!(8z1dI_dM!*;WV+4#5 zFh;-_0b>M=5%>p3AZd3{-Xl8x!4M4BsYW2_m&iGXb)>pvoc;qNuzXQe{w$V3bIN-; z7q!)qRp%>o?H>rmaH$cbnT$-d8^IKjVk?m>)q3SflY4z;`KIvxYRYfQ|XAjTfAa* zRM&Y!U9;1@K5CPGN9g4|dDf9170NQ%c7e;lfD$`a=3iDbu}zupTD;*C<7j1-sN=ZK z>+0H0-K(lL+0-_>`n;Uz(w~=!mzpee-(R4lbM+sW`F&o-vct}TD=ICXVZG^xy1v_W zIO)8jt_A3xtJ>67Tl%EGk$SODoMnI0!drtHEzI9#bID>|2QD9z=cHw;Wj9xJdC#ld z&sBf)8(6eNg`qkw>%6P3wb8w?YEwdO3#iXF_01|9FE+~bW4&58t3?q<%9L7mEvwh8 z+oQ`4df$I#CeKObyB*FP(4fNaHrZpR1r*ZZq4N)Q?G4=kHfrXS-rRZP@5ZX z4Yet&^9#C`NB0h@&F^YUANux{jf=A7CYvuN+O%(V%hpYcN`HQ>RJCO%y)LdCW7T}p zeqYlUXDd|n^IvtUXR(S=&P4o!&Vgv zDmEvs()lAN1Gl#Dw9C@Q)9#VZOX%82-J7O1pQ^2=`UI$NKiR-mv}{JnwwruB->B2B zwUb{uFOKU!*5cZPao+RCRkQwJ(rUkv4S%w6opW&YTg4{Y?+yDY=~+qyVY-*3HVxIby86%;8w+F$n+dXgNsB3$-SEDtxmgTK9t;f{hye0fAe;-vf^oLJw2j98;pmLLL z4IDW*|4!Gw)4lv^GhJ0he~nTCx_hXc+PRpx)W=cx6WNVd#bZ?2;(Dt_a#H;T_*}&E>vWe}y@`0~9@>x&5$0~;Z?J3eu z@4nS%EPhx#PwAef3O*CR&R)*?4R8Nh0Xgjae|HRA)S!Zc-=5PUMS|Np70o+v-N)8f zs~`r#Q)Gz`DD%PxH3oP@A7I@EziQ3YqxcXvap=@Cj+xWoO zXY%=2zWXbNv@@cA?@4$sn2U8SbCupc-r2|h!=mLLdPe(~>0QCD=>BejjXwCmVM{if z(0i4Fo%*dfx9)=_vz>P?9pxN!&2D|0e12iK_KvH%tIqwZXMeG(c2UPpoj=gEHM+M( zZFZ?Gedvn~Y&DimY~ur8_}n4i>lFiYlAhla)+pW!=AuKEQspu}+Umc*kGoy)n_~kf zL{4>hXxSz7!(R8D>NvQEUB8*n#j}*N^Z9_C>)TW{Rju)GN453Y_Ezf}8d5E4Ps?gS zIy&oojjq+#z4~f%Mr}RRhrV@X16z+|bG&ThXqt7bh1sT zHu2t9uFb<2MQj}ZC^Blqz$o9Z7f1Q)h}8MFy7sB=Em51+YFkZxYN;E21T zX|A?Q)n}Xf-jR*{vW3lYvh5-t_`)ae2{8~0F_CwX3%|a9@&8@q)ZNs2)EDHNT8#+aPPMJv3uC1hVHJ%>^2-KlRYj+)jn}QZ(9_1R!0_{C+V8K?nSAMhuX%d z&sp_tCL6N+|0FB z@m0>wiFX`xEZ%!|riA4>>~vmP*JkP71-0>3+al_-L4AkI#@n)$TkpvX*H!B&WoSiyzC>fuj#me-kQgky!*gAam(05b<`1sep90t2D$re32)Ag3jd zV((_ZV2$FvU@n5~z1_ED%o?}(xk@m}O}~pBBdxTKHy6Y_6lX&d2K- z_wuXFZndS)B=wDzjf=9iS2oMZ_C@*FAzzPz?82_eH!&1bEW|`?#7L}Q^~BCIz%#=A zKK-xn)xC-i1RVi;)a zaK-rLS2H9uKf5Sl`4@!~OV;a`_(^Tw#5XP^BxcspQ|Das(!FMCi=j{*as?xxnGz96S1k7Dm;zhI5xybc1Ns&_|)kK?j2N0KOZJ6ue-pzi&8pV6TTy&}WU7|&{+9vP&6HSe$tTwGl zh%p7}@YH#LuG#3`hiVg}w!9~K)%S{Q9F?u*vbjUH!{p-!`C1^KZ^?I{@UD)EWwBx+ zHe%ebSf?mvxZfl_0~!f*D`<((ccHm^m35G4IMEHFRYM!cbZOp_min;+!<3r9p`oa zP}iF19&=JiZCSfAs4shb9ocFxo75*aR5NyvFR+z5^6jMeC^R|fK+qn*cf+T`>%kkr zUx1l|(}E3xM}dJ*cT?+8OOVr&N3nOaH?T(WUN9GHqJK6yJ=kWlT{_cLNr%18x%Rv6 ztx%fXvKb3 z&+vhs1$`HqJ9Lm}Q_%~f;eM5L^=SUlnWJq+Pm6{V-5^>u^igPX{G|OrdjQ{^%(05g zgG++H05eap~b!?D8a!6m_8fSH5Sf(?R4fq_tWQ|nP*kkgV!v3IjyutxD-Fc&oo&8lBr z$0?n2jeE2iueOiRCnUbDzCN;XNwylx=DV_uj}7v*Qa)G7_a()!RRHj%^Q{BlpYTo4Klz;&zKMZYh>6&U5p6Si zS~Q&K2GOdak3y4!4kYz_!0W*!!Cxd}<{BfgLGUOr5bADfJ?aZ`TJk9NZuSe-DBcU^ zqLz*%ojd7TUfrV&bN3)+K32-cW!dtOO>A4q2Xl|lee%uU#0tf7K=YAR^U+H45vf?u zwf3}&RqUR6hTJLhu|?0eLC+YiXC3nj&4>B@R<=y2S=3VBUMo)`|6Ww5--g)#< zXmVJu&>p~d!?D8a!6m_8fSI$FfDM92fq_tWQ|nP*uwIczv3IjJu}1M;Fc-h+oHgvY z?y;s_Rof5M=VSG?myP1GwJl}6x+@i zS8wZCcqX3hdp)DuE4E&l=YzeBJ&nDMJ&wJOJ&(N)JuMnebc1Ns*w4}Au%DwnfbUM` zShc6ZCBa{SnS;}^pMyuS7gKk$pHp9u(~?KAcQYrfQM?z-MQF-?zD#ZAs%=5_8LPfd zvTe0504cNe7OySo)r3&lo^ z+K3v7T8Wy8+KCzpZ8J3$8cuYBXw|5>(Bz;4L3;q-4aW+v z2bTnY0cH+P3pNNIMQumjO|8eAkkgV!nd?E;DBcU^qMO=8sBIJVDWJaC*ehF0WplS| zN6W|esp>&ucu%nqQ*cT>NbO1uTS>9=41@J7+Qt*=9DR>Dn5gVA!>yDS0=L4(_%nj@f3=S*~ zOb%=gj1H_0%nq%Z8Bat9g7yHu8;%vs4=xG*0?Zsd5o{1ViuZ)Ii&~HRf}EB-ioKit zf;Ecwg1N}AKAY6nLN;P$D@Ha~%66etc;egg4W>yf#FQFOY>@`=MDx6WO@mRJ@piCl zFl_X+XgJXgqE$2F?dU+z9>90QvBK-YCBa{SnS;}U4T493flzl->rr2DkVmn1vtO`A z@m??&v(y(GiL!;wN-22zcKNC#pJ4q1gtwP?0&gFcM(}p?e8Ann;lSm<>A>xn`6xIa zxF0hg1t$a_g$@Mm0em+cYcj7VE(!hu%$&6gY>+*kyoeb7hvY#v|xkaQD7j{-PC&27v!|$QS9C97pzgd7tF<1vK20y zaKT;V18x|O*v!|$9mCf?Na1T2rZK+OJRfKQ&;+0jNUdi;JAj4&EdiPWv;}5813D12 z2h0f^E4&_D68r_2IXErYAb1oQ2z58L9`yw|EqN4sH~R%^6z>IdkwZ3vWqXYD3|Uk3 z43AUk8LlX%gVHn9ke=a28r3tH=LPKy8XB}TXlhdHebC&Xy+MP676(lZIuNu6@ZE5% z?D23(@E6o4;Iv?a;89>8oKx#jUy#$1N3nOaUtojxg1Ip3ea5HgeJ09xap`?#OYbv3 z9q4_`^MSSujTu@qG-qhf(4eK(!=X__tHzw5$w3E#_5i*cj+I;(E(!hu%sd&V6&3{^ z1qMRhO|3_LK~75^MPF>NM)6)?+pLG{DLq_>^l+Ba!bRcLC;Je{i;q~B>;4i?;!D+z;xknq?QtMG) zkkev|y_@}lHH!Cwx#*P6zc(=(xOYJZ#y z{)_j7wTu0PJ)XRaT$uWV`iwf>{5>c7C^R|fK+qn*cf+y5>%k?#Ux1mX{+^S*)ZN&k zz96S1kHVMnJ?DSUhxrP7Jb4$nF!hP~duOl}G@Rgk=I@=+eD6Vl2?E!6#-|KsxGx0B5 z@QiT3@bPFQ(5;{)Lf?g+77Zr{S~c`hXmZej7{6!v{O2OIMKluVR?rfm??Q8jRuDZc z8cxp9s&NlZj`4f1sYmbmpQBqrON71)%^f;Ow5e#D(bJ;gL^sGi+8Dow{QT$OX%^Ak zp@T%5ie5NYx@NS^=xL4L3xC>(z1VB$Akn6x7e>R4t{$Da@q79&cHI9z>u9*q)$^@` z@$XIif1~x{@1v{dTL->BG5$T47eC(rtp~mdGXA}y|8}%qd^^6=`X|2=d+|?bcrZr5 z7y)Ahj1e$Kz!(8z1dI_dM!*;WV+4#5Fh;-_0b>M=5imx;7y)Ahj1e$Kz!(8z1dI_d zM!*;WV+4#5Fh;-_0b>M=5imx;7y)Ahj1e$Kz!(8z1dI_dM!*;WV+5W)0-HDFp7yPd zr{6WMzU&bgd9HN+F*;thU&G%YN1*0_GBdIdEu25x!f)<+9XWJfSJ(d7%J{FBH3AQJ z+7|SzQ+nR5canIxnVck99A*+Ps`- z(yy1gXMUP>`m|r6>3tm+miJjyB5CX8rFE?rm0t1Fey`O!8tD9puG#3`47E9{w)NE~ z{f^j+ee#NDk_soCT)gnpJwr;i%@MI|K*+i>kHV|^)Gd#O>FscY}*UT?KI zq_+0z^J1e)f7WwuJvjUB!2N~$d30Jds!YyON3xDxeld2Q<>bJYE5`r0)cWeWH+*mE zn5^@uy2ibuYBNS{`>W4A_026C>2Ks-%oF<$jxMrge#D~XIeki5*jRemu4`L%`2)9= zQ~lz~g@3=@Z$>Ns@=J6?>wKoJt`KZky zwWZHS_01w1gJjD}Hq*}tKHn!j5!<(D+uL&%pL~0U#ju{0y{81vuv&G%*Y|P0;T2Aw zzaOxEZce)+4?kJm{pj+5EINMFd9<$mtb1AV53jIBZKtZwaP>VU8}DLEHlJ^lY2VuP z1(tK~cI#7gY5(6A&wlH=MftHyy_eU@W9@HK)o)zKGB!VLk6$%7A!o&jj}NWM?c!DO zPR@7iymb`O`E*_TM)%IC%|x}OkC*yZkPU2k%jN^wo*^G;ertN5V&|&n@TrxN881fTx7pphqca^ z>e@Zs8>u#D)OLvaj8xyBWy4yw3d`pEvRzR=o=>!C*UH!aljd%S`&IA7kR@)*CzY=8 zsHx9K%MVyj`NYD%-Q9||KC7~<=~Zr;{jD4EK@+ZzbLh9vJ~(KFOK`J-GlM;J1n7K# zu1(UtyK3{f+S2D;^~J`=vgIzD*!Gr>c=;MCpJ{hQ|F3BK1bn)pk zea6=`+U` zB%Md-+FQCeT5Srd?Jo84QQu**fh}y7kZpWqldom+`M*V-_Ih7v7(chm&YH!RKia&c zYo@EE$DcI$WGK_ET>FWM{&pR2*cPvIdQFS37FG)Fn>!?AtF`0A57vh!_(z7c8r#z8 zcyUjsyJsgPA3x~4jjnyGdlS@Vo!YKdANrP%4Q#z9o7kQpANcZ?&u;R4Q8A>wk^Os5 zw6Nr7#kx1TxMbU&++{l3UH93Zdrr9lowEcS>+Wo4+1(}ZK#MF6KRuij($zZ7v0R@? zYwa!^cWl&qYL$l-9UVX2{O#K4H)pSl)G<=$x$Pa3@0C>>+ICT&GwPcp8`%0@HnF`! zKJfLee10h3OBBQ3i$3l2n%O*Q?!@c8iZyz4YsoFEvSm&V%)Mgw`0#Q)584G}J=M={ zaQLXe*sCoZWKCS&-ZIdA#z>iAOU zk-Fxnd$ghLE%o_IeX)TpY-0OU`M?)GZ_9Tx#V}2=%vMZkXN3RSlQ0*TYh_+Kv3JEX zWlNM?@!eZ{%7ssB7I5^d1$Oayg93w3%yu|4!7cRHoL){jTdZ2wc3fs>?~+HHyYAQ! zzW>6y@aX)B&dcuQtlA^2XVtzsmg@XBUHd`z-d7vi(ucm-z!o;+WE&q5@`X=p`JSv8 zt|%7f))FFf;Ecwg1Pv5g5`>j9FLbPeyLkP=VPIEMY1IYF1Jr`Sm8M$w5!_< zr=ZyUVSb|*INL9Ca+%p{!un4N{v1|weS>PVZq2VYsc2%g2Jvy#-qOMOOkHcFdo>IF z95zjDE22girrp`V{ESEWho_rtR6&gn8OPj8ns9t9#YWj|Q;Q0A(x!7dl8g;p;S zG3wTs$VXN?A~Ra=h#al+W4acgdneWAY_@|{1Jq}{`eLK0Z1t8+C)viwQTfUzpZF#Q z-V@d?_7nE_*Z(j7&!2$$f}EB-ioKh)oArXfQRd?F>~{h;4ZQ8J_uJ^u{`cOl(&_5# zu$Ip6xm-4lbD1AGGGb2d_-cVam|TCi9~>E0Wm(j;@j0SP_H~XP7GFR5YaN_V*0l|~ zSG3aL$QZR9qds%g_q=R0m#z0@bAW8`mygZzg-_;$7>I?K*yG8&$c58#eL}5AeL$fqqfoIn$1KmQxn%A7+;!|*(%NXaxh^@CgU4J-V+^vw?`Y4Yx z{i8dt{xbUW<~gGC?C$So@q3->1$0>G{4-r^pnGlA#zJk+t51~rj*ty(VH4X8<>N#7 zVorEZh=Kitn23!SsZXdG|NO7)^}7om1qMRhO&vykK~75^#lFf~&icW7!CbU$A03wa zeP0)kGhew>$#^GX!MN3~YwvitEwylXlytXg^Pzi}gX7)1 z=_s#rt_{(>scLgmZGF{eqx#N}4F}nJM>gM-ZG12%yeIf324W%aA{VAUAy(>mu$8og zcYzIpM}bpNcTk){4rwq3ovtVTChRzDDVpEZfZSh z0rF4sDE4ml3)U0f3+7_*zFV&E)EeZrf(?Qpfm={3Q=d^2kkgV!v3Ik^vqte= zFc;f1Ws5mjE+o45#ToANHO4o0m#(q>oA3=k9vT<7v3kb%obL7G+imU@-&#jjoyY4M z_oCEhw%X2ApRd&Sv}|CDIeAmI9pz&+`-yy#caaMd3o#KJF%l~nB)Blo0Dl6f^~b-u zSO34@yWvOS_281=FTl*fSiuItqrgC@nW@vLFUUv9qu9IIFIb~^FPMw4{$B2OB_Fx( z@M*JQi2lK@`G*?E-?+a!KHt0c33Wc|o8XXVT7ta}cb(_fwd=aKNNxJ6?Lqa)r@nEr zaZqRTmG{7NVZqX$65KBD4(n3+fgyJ zQY<$W6Bs18a57e}nBlZ|2Do3I35^806|_XJvc3zg8u}90QvBK-YCBa{S znS;H8_kl5iflzl->rr2j(~?KAce7uxM)6)S7a1hiu=F59*UB`Mjt)e!5~A~gErs(lOrcy(; znzl~fYf98%tMde1%lOM^(?qq2QQIHYr<3}!pKO+`9kThUY!8r+R`OL?b^HnWwiMpg zMX}sgO!*ZX{0Xt%QOv`|{hrk`ppo#_God9y--YH59VFV+S4l67wi!Jw8cuYB=+V$e zp~*o9g7yHu8;%uT4=xG*!dbXDI4#&9coY~2bvLyh^#wUCxe|Lfdje||?*(&ly<@D& z?dKz=ojLcL+WEzqy6M=Y^P0MLSNAe$PLkA?wd=b2UXzW}vNc0Cm&^87^3htpqU3X> zd^b=G-zyfd`X-9)sA8RHedq3_}u(Lth3MK6qo`=5SK`2GUTKRR=? z&FE>-aH1PTtA;)bO%6H`vQB*eMCW;QZJX{{sts-HtIr(utxqm2 zTbX5ZiEM+d6p$}`;(MxMn5S3Vri|I)+@H# zim`)YZKarZDfTUT23I}Hy_EU5rf00LXZ6!FultAR!~8z-Zt{-uuJX?E?(zgq0G_3}AbzQ-#DFU2xYF%4C0 zKPtxSiq&5+hbZcs;lz_zUxX4mJoL1qQ;NOsz+KK~75^#oo<+!5YPT z!CW-dHSV$J&sW<>^{J!2*w`dn*_Td;@z;L-P&vb3`NHQN`JST~t|*oW#neWzk#j^U zR&o$>k;zZ&=hO9!1NE$jU!ncnJTK&WpgXq(Z~qTxh0h*k}K6!|W> zCE5e{Za7wWJ-8(J3ovtVTJl}+DDqu$S86@#3vycWDBcD33)U#knTwmccTYL-F10

0M_i-r^3Ahj0yDC$9UAZQQZyWv>j_281=FTl*fX{qJFqrgC@ zyQ%f4FUV=hqu9IIFIb~^FPMvADfQrS^?6Tym&-=DY;BWGY>$+W_VSfQK98r=gB?@q z!ICNU;7G-aHkEqtjGiHro`q+!*RvhdGk*4p)`RAG0XqOg080Q<09ycKKxYo-0QP`} z6Ww5HyaXKx+5`A*FbsG-xFq-sFmrHPutBg7Fc9i)YCSVv!XD4w&3?fe#d|@U2N&!U zf6p-~;am0XCmR;ROE$=6OWD37A7C>P@_9K0F9Ev&FCn(j*BCD`&j;8V7#mm{m>bv| z7#vs}m>hT_7#&z0m>pU*@I*8@=s?gOz<0y3!t22$!C!!xgVTZyGACdltX=&$2ycf*H_v&M*zF?Bzi2<_trfh>*o|CUhDR?57Cb19`SSK;=M*I5&p7^Dn z$@?|I6V3Ai77Zp1HVsA%Rt;tic8#_fyd4cExKAcxpZB3l8!q_HOnI)+pW!<|0^qi^#?x+1f6fVEKvi5to9uPmymhe_|P| zn23!S&3OAL#SUliy5a5S`GBkO7iR-^1BU~b1E&MG1IGi`1LtGrqtL3sN1@3<2ZHth zz8mfcUJot_{sPP#oR<9rJPHhix|>>$`huL6Jc@hl7pzgd7t94VZphYd*({ohkAefk zcZg!R_5>d_RIyG=#Yg@4n(|TRd4b!6p~d!?D8a!6m_8fSH5S zf(?R4fq_tWQ|nP*kkgV!v3Ijyuts5vxj;*lD@E_~rF=Tex0hl-dz4*zpMlc*WRl*e zWE#`^nCAm68JaS*WoXRMnxQ#EdzM-chc*q38d|mFInm!RIuNu6@ZE5%@Op4b@E2g_ zU@PEV;8C1YcT?+8Uy#$1N3nNf16#Zo%tdNFoTKz`#ZvTey;AgWlN96oX;cqqo)@&2 zXfRXjMbT)Y)kL$2b`uRJbHaPVdV(eg9SGV3_-^VGcs;lz_zN&|a9XfI@F*}4>TYU1 z>I-sO@+fSuUtkm4X1(YZ`TS14yDElA#j-aY=ta%*fi@Y9GFoLc%c=F$XqwSBGbg+! ztX;|biS}DGIp{#p9>90QvBK-YC4oVLnS;}U4T493flzl->rr3OmpqES8=I_A_+T#j z$Y*}}E+swn57JZrC_VM|bfTv=&kI_9H2rA%(fFIcH(*Y9Pten{pRmW1ccG6$lY&_|)k zK?j2N0KOZJ6G^w3^SrQjv1YT!^DPXyFq(gAMzqafD`+^uAknI!k4o0$NC$%U0KOZJ6-i4kP%oW`r z{0aIfG&$%%&>p~d!?D8a!6nfK%$z={e-D7Vn_7?h0$;}O0X+M!k6f7ggqo2$9&81? z3k(vh6|5c&C%QqjYUrcT?0Is!-%VTk zz$F>KNA~Aq{@>gMtA{^<(}J6W`-P83BZ0OVJuMnebc1NsIOiHV5VQxhr4N1Kl8oQ$ z`)@}5*=^uY;I!Zw;eO%c(MX{AM|Xp^89gl;PIQBubB%lGKxhl!ZTz0)v&Z=F+=2Ur zk4GbcZUrq7`Yv?lXq(Z~a-bXJ8u!rT&=&21@q4d-XVjnF2RNQ@`A27N z{Cg8GYSf=?7rihVZglm0>%jLX#=pn%?D2Y8clg$U?@#z9$oTh)Ue>5T+YjHL@J$fk zX>l0;p4YR->Sf;H+d1RkyL*`<{^!2M=5imx;7y)Ah hj1e$Kz!(8z1dI_dM!*;WV+4#5Fh;-_ftND^{|AL-*M Date: Sat, 5 Aug 2023 12:05:00 +0300 Subject: [PATCH 17/25] setup: Add necessary data --- fast/setup_wave2d.py | 116 +++++++++++++++++++++++++++++++++++++++++++ fast/wave2d_b.py | 46 +++-------------- 2 files changed, 124 insertions(+), 38 deletions(-) create mode 100644 fast/setup_wave2d.py diff --git a/fast/setup_wave2d.py b/fast/setup_wave2d.py new file mode 100644 index 0000000000..74f5ae2b67 --- /dev/null +++ b/fast/setup_wave2d.py @@ -0,0 +1,116 @@ +# Script to save initial data for the Acoustic wave execution benchmark +# Based on the implementation of the Devito acoustic example implementation +# Not using Devito's source injection abstraction +import sys +import numpy as np + +from devito import (TimeFunction, Eq, Operator, solve, norm, + XDSLOperator, configuration) +from examples.seismic import RickerSource +from examples.seismic import Model, TimeAxis, plot_image +from fast.bench_utils import plot_2dfunc +from devito.tools import as_tuple + +import argparse +np.set_printoptions(threshold=np.inf) + + +parser = argparse.ArgumentParser(description='Process arguments.') + +parser.add_argument("-d", "--shape", default=(16, 16), type=int, nargs="+", + help="Number of grid points along each axis") +parser.add_argument("-so", "--space_order", default=4, + type=int, help="Space order of the simulation") +parser.add_argument("-to", "--time_order", default=2, + type=int, help="Time order of the simulation") +parser.add_argument("-nt", "--nt", default=20, + type=int, help="Simulation time in millisecond") +parser.add_argument("-bls", "--blevels", default=1, type=int, nargs="+", + help="Block levels") +parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot2D") +parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run") +parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run") +args = parser.parse_args() + + +mpiconf = configuration['mpi'] + +# Define a physical size +# nx, ny, nz = args.shape +nt = args.nt + +shape = (args.shape) # Number of grid point (nx, ny, nz) +spacing = as_tuple(10.0 for _ in range(len(shape))) # Grid spacing in m. The domain size is now 1km by 1km +origin = as_tuple(0.0 for _ in range(len(shape))) # What is the location of the top left corner. +# This is necessary to define +# the absolute location of the source and receivers + +# Define a velocity profile. The velocity is in km/s +v = np.empty(shape, dtype=np.float32) +v[:, :] = 1 + +# With the velocity and model size defined, we can create the seismic model that +# encapsulates this properties. We also define the size of the absorbing layer as +# 10 grid points +so = args.space_order +to = args.time_order + +model = Model(vp=v, origin=origin, shape=shape, spacing=spacing, + space_order=so, nbl=0) + +# plot_velocity(model) + +t0 = 0. # Simulation starts a t=0 +tn = nt # Simulation last 1 second (1000 ms) +dt = model.critical_dt # Time step from model grid spacing +print("dt is:", dt) + +time_range = TimeAxis(start=t0, stop=tn, step=dt) + +# The source is positioned at a $20m$ depth and at the middle of the +# $x$ axis ($x_{src}=500m$), +# with a peak wavelet frequency of $10Hz$. +f0 = 0.010 # Source peak frequency is 10Hz (0.010 kHz) +src = RickerSource(name='src', grid=model.grid, f0=f0, + npoint=1, time_range=time_range) + +# First, position source centrally in all dimensions, then set depth +src.coordinates.data[0, :] = np.array(model.domain_size) * .5 + +# We can plot the time signature to see the wavelet +# src.show() + +# Define the wavefield with the size of the model and the time dimension +u = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) +# Another one to clone data +u2 = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) +ub = TimeFunction(name="ub", grid=model.grid, time_order=to, space_order=so) + +# We can now write the PDE +# pde = model.m * u.dt2 - u.laplace + model.damp * u.dt +# import pdb;pdb.set_trace() +pde = u.dt2 - u.laplace + +stencil = Eq(u.forward, solve(pde, u.forward)) +# stencil + +# Finally we define the source injection and receiver read function to generate +# the corresponding code +# print(time_range) + +print("Init norm:", np.linalg.norm(u.data[:])) +src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) +op0 = Operator([stencil] + src_term, subs=model.spacing_map, name='SourceDevitoOperator') + +# Run with source and plot +op0.apply(time=time_range.num-1, dt=model.critical_dt) + +if len(shape) == 2: + if args.plot: + plot_2dfunc(u) + +import pdb;pdb.set_trace() +# Save Data here +shape_str = '_'.join(str(item) for item in shape) +np.save("critical_dt%s.npy" % shape_str, model.critical_dt, allow_pickle=True) +np.save("wave_dat%s.npy" % shape_str, u.data[:], allow_pickle=True) diff --git a/fast/wave2d_b.py b/fast/wave2d_b.py index aada12f24a..f0fc4c13b5 100644 --- a/fast/wave2d_b.py +++ b/fast/wave2d_b.py @@ -57,30 +57,8 @@ so = args.space_order to = args.time_order -#model = Model(vp=v, origin=origin, shape=shape, spacing=spacing, -# space_order=so, nbl=0) - -# plot_velocity(model) - t0 = 0. # Simulation starts a t=0 tn = nt # Simulation last 1 second (1000 ms) -# dt = model.critical_dt # Time step from model grid spacing -# print("dt is:", dt) - -#time_range = TimeAxis(start=t0, stop=tn, step=dt) - -# The source is positioned at a $20m$ depth and at the middle of the -# $x$ axis ($x_{src}=500m$), -# with a peak wavelet frequency of $10Hz$. -f0 = 0.010 # Source peak frequency is 10Hz (0.010 kHz) -#src = RickerSource(name='src', grid=model.grid, f0=f0, -# npoint=1, time_range=time_range) - -# First, position source centrally in all dimensions, then set depth -#src.coordinates.data[0, :] = np.array(model.domain_size) * .5 - -# We can plot the time signature to see the wavelet -# src.show() # Define the wavefield with the size of the model and the time dimension u = TimeFunction(name="u", grid=grid, time_order=to, space_order=so) @@ -94,19 +72,6 @@ pde = u.dt2 - u.laplace stencil = Eq(u.forward, solve(pde, u.forward)) -# stencil - -# Finally we define the source injection and receiver read function to generate -# the corresponding code -# print(time_range) - -#src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) -#op0 = Operator([stencil] + src_term, subs=model.spacing_map, name='SourceDevitoOperator') - -# Run with source and plot -#op0.apply(time=time_range.num-1, dt=model.critical_dt) - - # print("Init Devito linalg norm 0 :", np.linalg.norm(u.data[0])) # print("Init Devito linalg norm 1 :", np.linalg.norm(u.data[1])) @@ -117,8 +82,13 @@ u2.data[:] = u.data[:] configuration['mpi'] = mpiconf -import pdb;pdb.set_trace() -u.data[:] = np.load("wave_dat2.npy", allow_pickle=True) +# import pdb;pdb.set_trace() +shape_str = '_'.join(str(item) for item in shape) +u.data[:] = np.load("wave_dat%s.npy" % shape_str, allow_pickle=True) +dt = np.load("critical_dt%s.npy" % shape_str, allow_pickle=True) + +# np.save("critical_dt%s.npy" % shape_str, model.critical_dt, allow_pickle=True) +# np.save("wave_dat%s.npy" % shape_str, u.data[:], allow_pickle=True) if len(shape) == 2: if args.plot: @@ -130,7 +100,7 @@ # Run more with no sources now (Not supported in xdsl) #op1 = Operator([stencil], name='DevitoOperator', subs=grid.spacing_map) op1 = Operator([stencil], name='DevitoOperator') - op1.apply(time=10, dt=5.54) + op1.apply(time=nt, dt=dt) configuration['mpi'] = 0 ub.data[:] = u.data[:] From 8db89d76c8a0dec8a505ba5cf77ade28a3e1e533 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Sat, 5 Aug 2023 12:20:36 +0300 Subject: [PATCH 18/25] bench: Load dt to XDSL --- fast/diffusion_3D_wBCs.py | 2 +- fast/wave2d_b.py | 9 ++++----- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/fast/diffusion_3D_wBCs.py b/fast/diffusion_3D_wBCs.py index 9cb39d5179..9220e25247 100644 --- a/fast/diffusion_3D_wBCs.py +++ b/fast/diffusion_3D_wBCs.py @@ -76,7 +76,7 @@ def plot_3dfunc(u): x, y, z = grid.dimensions t = grid.stepping_dim -print(eq_stencil) +# print(eq_stencil) # Create Operator if args.devito: diff --git a/fast/wave2d_b.py b/fast/wave2d_b.py index f0fc4c13b5..ffc570578a 100644 --- a/fast/wave2d_b.py +++ b/fast/wave2d_b.py @@ -90,15 +90,14 @@ # np.save("critical_dt%s.npy" % shape_str, model.critical_dt, allow_pickle=True) # np.save("wave_dat%s.npy" % shape_str, u.data[:], allow_pickle=True) -if len(shape) == 2: - if args.plot: - plot_2dfunc(u) +if len(shape) == 2 and args.plot: + plot_2dfunc(u) print("Init norm:", np.linalg.norm(u.data[:])) if args.devito: # Run more with no sources now (Not supported in xdsl) - #op1 = Operator([stencil], name='DevitoOperator', subs=grid.spacing_map) + # op1 = Operator([stencil], name='DevitoOperator', subs=grid.spacing_map) op1 = Operator([stencil], name='DevitoOperator') op1.apply(time=nt, dt=dt) @@ -128,7 +127,7 @@ # Run more with no sources now (Not supported in xdsl) xdslop = Operator([stencil], name='xDSLOperator') - xdslop.apply(time=nt) + xdslop.apply(time=nt, dt=dt) if len(shape) == 2 and args.plot: plot_2dfunc(u) From bcd7a8d123dc6a9ba8d40f509ca61b3ad0ac28e3 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Sat, 5 Aug 2023 12:55:03 +0300 Subject: [PATCH 19/25] setup: Save extent --- fast/setup_wave2d.py | 4 ++-- fast/wave2d_b.py | 20 ++++++++------------ 2 files changed, 10 insertions(+), 14 deletions(-) diff --git a/fast/setup_wave2d.py b/fast/setup_wave2d.py index 74f5ae2b67..dafde35517 100644 --- a/fast/setup_wave2d.py +++ b/fast/setup_wave2d.py @@ -58,6 +58,7 @@ model = Model(vp=v, origin=origin, shape=shape, spacing=spacing, space_order=so, nbl=0) +import pdb;pdb.set_trace() # plot_velocity(model) t0 = 0. # Simulation starts a t=0 @@ -88,7 +89,6 @@ # We can now write the PDE # pde = model.m * u.dt2 - u.laplace + model.damp * u.dt -# import pdb;pdb.set_trace() pde = u.dt2 - u.laplace stencil = Eq(u.forward, solve(pde, u.forward)) @@ -109,8 +109,8 @@ if args.plot: plot_2dfunc(u) -import pdb;pdb.set_trace() # Save Data here shape_str = '_'.join(str(item) for item in shape) np.save("critical_dt%s.npy" % shape_str, model.critical_dt, allow_pickle=True) np.save("wave_dat%s.npy" % shape_str, u.data[:], allow_pickle=True) +np.save("grid_extent%s.npy" % shape_str, model.grid.extent, allow_pickle=True) diff --git a/fast/wave2d_b.py b/fast/wave2d_b.py index ffc570578a..96caab9c78 100644 --- a/fast/wave2d_b.py +++ b/fast/wave2d_b.py @@ -39,17 +39,13 @@ nt = args.nt shape = (args.shape) # Number of grid point (nx, ny, nz) +shape_str = '_'.join(str(item) for item in shape) spacing = as_tuple(10.0 for _ in range(len(shape))) # Grid spacing in m. The domain size is now 1km by 1km origin = as_tuple(0.0 for _ in range(len(shape))) # What is the location of the top left corner. domain_size = tuple((d-1) * s for d, s in zip(shape, spacing)) -grid = Grid(shape=shape, extent=(1000, 1000)) -# grid = Grid(shape=shape) -# This is necessary to define -# the absolute location of the source and receivers +extent = np.load("grid_extent%s.npy" % shape_str, allow_pickle=True) -# Define a velocity profile. The velocity is in km/s -v = np.empty(shape, dtype=np.float32) -v[:, :] = 1 +grid = Grid(shape=shape, extent=as_tuple(extent)) # With the velocity and model size defined, we can create the seismic model that # encapsulates this properties. We also define the size of the absorbing layer as @@ -82,8 +78,6 @@ u2.data[:] = u.data[:] configuration['mpi'] = mpiconf -# import pdb;pdb.set_trace() -shape_str = '_'.join(str(item) for item in shape) u.data[:] = np.load("wave_dat%s.npy" % shape_str, allow_pickle=True) dt = np.load("critical_dt%s.npy" % shape_str, allow_pickle=True) @@ -132,6 +126,8 @@ if len(shape) == 2 and args.plot: plot_2dfunc(u) - print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) - print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) - print("XDSL output norm 2:", np.linalg.norm(u.data[2]), "vs:", np.linalg.norm(ub.data[2])) + print("After Operator 1: Devito norm:", norm(u)) + + # print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) + # print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) + # print("XDSL output norm 2:", np.linalg.norm(u.data[2]), "vs:", np.linalg.norm(ub.data[2])) From 8c8858ac68aed62205584e538c29f032efad9d05 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Sat, 5 Aug 2023 13:12:29 +0300 Subject: [PATCH 20/25] bench: Add so to saved data --- fast/setup_wave2d.py | 7 +++---- fast/wave2d_b.py | 16 ++++++---------- 2 files changed, 9 insertions(+), 14 deletions(-) diff --git a/fast/setup_wave2d.py b/fast/setup_wave2d.py index dafde35517..5492d742fc 100644 --- a/fast/setup_wave2d.py +++ b/fast/setup_wave2d.py @@ -58,7 +58,6 @@ model = Model(vp=v, origin=origin, shape=shape, spacing=spacing, space_order=so, nbl=0) -import pdb;pdb.set_trace() # plot_velocity(model) t0 = 0. # Simulation starts a t=0 @@ -111,6 +110,6 @@ # Save Data here shape_str = '_'.join(str(item) for item in shape) -np.save("critical_dt%s.npy" % shape_str, model.critical_dt, allow_pickle=True) -np.save("wave_dat%s.npy" % shape_str, u.data[:], allow_pickle=True) -np.save("grid_extent%s.npy" % shape_str, model.grid.extent, allow_pickle=True) +np.save("so%s_critical_dt%s.npy" % (so, shape_str), model.critical_dt, allow_pickle=True) +np.save("so%s_wave_dat%s.npy" % (so, shape_str), u.data[:], allow_pickle=True) +np.save("so%s_grid_extent%s.npy" % (so, shape_str), model.grid.extent, allow_pickle=True) diff --git a/fast/wave2d_b.py b/fast/wave2d_b.py index 96caab9c78..bfdaa995ec 100644 --- a/fast/wave2d_b.py +++ b/fast/wave2d_b.py @@ -37,13 +37,14 @@ # Define a physical size # nx, ny, nz = args.shape nt = args.nt +so = args.so shape = (args.shape) # Number of grid point (nx, ny, nz) shape_str = '_'.join(str(item) for item in shape) spacing = as_tuple(10.0 for _ in range(len(shape))) # Grid spacing in m. The domain size is now 1km by 1km origin = as_tuple(0.0 for _ in range(len(shape))) # What is the location of the top left corner. domain_size = tuple((d-1) * s for d, s in zip(shape, spacing)) -extent = np.load("grid_extent%s.npy" % shape_str, allow_pickle=True) +extent = np.load("so%s_grid_extent%s.npy" % (so, shape_str), allow_pickle=True) grid = Grid(shape=shape, extent=as_tuple(extent)) @@ -78,8 +79,8 @@ u2.data[:] = u.data[:] configuration['mpi'] = mpiconf -u.data[:] = np.load("wave_dat%s.npy" % shape_str, allow_pickle=True) -dt = np.load("critical_dt%s.npy" % shape_str, allow_pickle=True) +u.data[:] = np.load("so%s_wave_dat%s.npy" % (so, shape_str), allow_pickle=True) +dt = np.load("so%s_critical_dt%s.npy" % (so, shape_str), allow_pickle=True) # np.save("critical_dt%s.npy" % shape_str, model.critical_dt, allow_pickle=True) # np.save("wave_dat%s.npy" % shape_str, u.data[:], allow_pickle=True) @@ -102,18 +103,13 @@ if len(shape) == 2 and args.plot: plot_2dfunc(u) - print("After Operator 1: Devito norm:", norm(u)) + print("Devito norm:", norm(u)) # print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) # print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) # print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) if args.xdsl: - # Reset initial data - configuration['mpi'] = 0 - u.data[:] = u2.data[:] - configuration['mpi'] = mpiconf - # v[:, ..., :] = 1 # print("Reinitialise data: Devito norm:", norm(u)) # print("XDSL init linalg norm:", np.linalg.norm(u.data[0])) # print("XDSL init linalg norm:", np.linalg.norm(u.data[1])) @@ -126,7 +122,7 @@ if len(shape) == 2 and args.plot: plot_2dfunc(u) - print("After Operator 1: Devito norm:", norm(u)) + print("XDSL norm:", norm(u)) # print("XDSL output norm 0:", np.linalg.norm(u.data[0]), "vs:", np.linalg.norm(ub.data[0])) # print("XDSL output norm 1:", np.linalg.norm(u.data[1]), "vs:", np.linalg.norm(ub.data[1])) From 72fb27d7a4006c26bbe0280ca441da7ea582ffb3 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Sat, 5 Aug 2023 13:46:53 +0300 Subject: [PATCH 21/25] bench: Add wave3d setup --- fast/setup_wave3d.py | 115 ++++++++++++++++++++++++++++++++++++++ fast/wave2d_b.py | 2 +- fast/wave3d_b.py | 129 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 245 insertions(+), 1 deletion(-) create mode 100644 fast/setup_wave3d.py create mode 100644 fast/wave3d_b.py diff --git a/fast/setup_wave3d.py b/fast/setup_wave3d.py new file mode 100644 index 0000000000..795c96b284 --- /dev/null +++ b/fast/setup_wave3d.py @@ -0,0 +1,115 @@ +# Script to save initial data for the Acoustic wave execution benchmark +# Based on the implementation of the Devito acoustic example implementation +# Not using Devito's source injection abstraction +import sys +import numpy as np + +from devito import (TimeFunction, Eq, Operator, solve, norm, + configuration) +from examples.seismic import RickerSource +from examples.seismic import Model, TimeAxis +from fast.bench_utils import plot_3dfunc +from devito.tools import as_tuple + +import argparse +np.set_printoptions(threshold=np.inf) + + +parser = argparse.ArgumentParser(description='Process arguments.') + +parser.add_argument("-d", "--shape", default=(16, 16, 16), type=int, nargs="+", + help="Number of grid points along each axis") +parser.add_argument("-so", "--space_order", default=4, + type=int, help="Space order of the simulation") +parser.add_argument("-to", "--time_order", default=2, + type=int, help="Time order of the simulation") +parser.add_argument("-nt", "--nt", default=20, + type=int, help="Simulation time in millisecond") +parser.add_argument("-bls", "--blevels", default=1, type=int, nargs="+", + help="Block levels") +parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot2D") +parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run") +parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run") +args = parser.parse_args() + + +mpiconf = configuration['mpi'] + +# Define a physical size +# nx, ny, nz = args.shape +nt = args.nt + +shape = (args.shape) # Number of grid point (nx, ny, nz) +spacing = as_tuple(10.0 for _ in range(len(shape))) # Grid spacing in m. The domain size is now 1km by 1km +origin = as_tuple(0.0 for _ in range(len(shape))) # What is the location of the top left corner. +# This is necessary to define +# the absolute location of the source and receivers + +# Define a velocity profile. The velocity is in km/s +v = np.empty(shape, dtype=np.float32) +v[:, :, :] = 1 + +# With the velocity and model size defined, we can create the seismic model that +# encapsulates this properties. We also define the size of the absorbing layer as +# 10 grid points +so = args.space_order +to = args.time_order + +model = Model(vp=v, origin=origin, shape=shape, spacing=spacing, + space_order=so, nbl=0) + +# plot_velocity(model) + +t0 = 0. # Simulation starts a t=0 +tn = nt # Simulation last 1 second (1000 ms) +dt = model.critical_dt # Time step from model grid spacing +print("dt is:", dt) + +time_range = TimeAxis(start=t0, stop=tn, step=dt) + +# The source is positioned at a $20m$ depth and at the middle of the +# $x$ axis ($x_{src}=500m$), +# with a peak wavelet frequency of $10Hz$. +f0 = 0.010 # Source peak frequency is 10Hz (0.010 kHz) +src = RickerSource(name='src', grid=model.grid, f0=f0, + npoint=1, time_range=time_range) + +# First, position source centrally in all dimensions, then set depth +src.coordinates.data[0, :] = np.array(model.domain_size) * .5 + +# We can plot the time signature to see the wavelet +# src.show() + +# Define the wavefield with the size of the model and the time dimension +u = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) +# Another one to clone data +u2 = TimeFunction(name="u", grid=model.grid, time_order=to, space_order=so) +ub = TimeFunction(name="ub", grid=model.grid, time_order=to, space_order=so) + +# We can now write the PDE +# pde = model.m * u.dt2 - u.laplace + model.damp * u.dt +pde = u.dt2 - u.laplace + +stencil = Eq(u.forward, solve(pde, u.forward)) +# stencil + +# Finally we define the source injection and receiver read function to generate +# the corresponding code +# print(time_range) + +print("Init norm:", np.linalg.norm(u.data[:])) +src_term = src.inject(field=u.forward, expr=src * dt**2 / model.m) +op0 = Operator([stencil] + src_term, subs=model.spacing_map, name='SourceDevitoOperator') + +# Run with source and plot +op0.apply(time=time_range.num-1, dt=model.critical_dt) + +if len(shape) == 3: + if args.plot: + plot_3dfunc(u) + +# Save Data here +shape_str = '_'.join(str(item) for item in shape) +np.save("so%s_critical_dt%s.npy" % (so, shape_str), model.critical_dt, allow_pickle=True) +np.save("so%s_wave_dat%s.npy" % (so, shape_str), u.data[:], allow_pickle=True) +np.save("so%s_grid_extent%s.npy" % (so, shape_str), model.grid.extent, allow_pickle=True) diff --git a/fast/wave2d_b.py b/fast/wave2d_b.py index bfdaa995ec..411d61ef14 100644 --- a/fast/wave2d_b.py +++ b/fast/wave2d_b.py @@ -37,7 +37,7 @@ # Define a physical size # nx, ny, nz = args.shape nt = args.nt -so = args.so +so = args.space_order shape = (args.shape) # Number of grid point (nx, ny, nz) shape_str = '_'.join(str(item) for item in shape) diff --git a/fast/wave3d_b.py b/fast/wave3d_b.py new file mode 100644 index 0000000000..9e0aa32190 --- /dev/null +++ b/fast/wave3d_b.py @@ -0,0 +1,129 @@ +# Based on the implementation of the Devito acoustic example implementation +# Not using Devito's source injection abstraction +import sys +import numpy as np + +from devito import (TimeFunction, Eq, Operator, solve, norm, + XDSLOperator, configuration, Grid) +from examples.seismic import RickerSource +from examples.seismic import Model, TimeAxis, plot_image +from fast.bench_utils import plot_3dfunc +from devito.tools import as_tuple + +import argparse +np.set_printoptions(threshold=np.inf) + + +parser = argparse.ArgumentParser(description='Process arguments.') + +parser.add_argument("-d", "--shape", default=(16, 16, 16), type=int, nargs="+", + help="Number of grid points along each axis") +parser.add_argument("-so", "--space_order", default=4, + type=int, help="Space order of the simulation") +parser.add_argument("-to", "--time_order", default=2, + type=int, help="Time order of the simulation") +parser.add_argument("-nt", "--nt", default=20, + type=int, help="Simulation time in millisecond") +parser.add_argument("-bls", "--blevels", default=1, type=int, nargs="+", + help="Block levels") +parser.add_argument("-plot", "--plot", default=False, type=bool, help="Plot2D") +parser.add_argument("-devito", "--devito", default=False, type=bool, help="Devito run") +parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run") +args = parser.parse_args() + + +mpiconf = configuration['mpi'] + +# Define a physical size +# nx, ny, nz = args.shape +nt = args.nt +so = args.space_order + +shape = (args.shape) # Number of grid point (nx, ny, nz) +shape_str = '_'.join(str(item) for item in shape) +spacing = as_tuple(10.0 for _ in range(len(shape))) # Grid spacing in m. The domain size is now 1km by 1km +origin = as_tuple(0.0 for _ in range(len(shape))) # What is the location of the top left corner. +domain_size = tuple((d-1) * s for d, s in zip(shape, spacing)) +extent = np.load("so%s_grid_extent%s.npy" % (so, shape_str), allow_pickle=True) + +grid = Grid(shape=shape, extent=as_tuple(extent)) + +# With the velocity and model size defined, we can create the seismic model that +# encapsulates this properties. We also define the size of the absorbing layer as +# 10 grid points +so = args.space_order +to = args.time_order + +t0 = 0. # Simulation starts a t=0 +tn = nt # Simulation last 1 second (1000 ms) + +# Define the wavefield with the size of the model and the time dimension +u = TimeFunction(name="u", grid=grid, time_order=to, space_order=so) +# Another one to clone data +u2 = TimeFunction(name="u", grid=grid, time_order=to, space_order=so) +ub = TimeFunction(name="ub", grid=grid, time_order=to, space_order=so) + +# We can now write the PDE +# pde = model.m * u.dt2 - u.laplace + model.damp * u.dt +# import pdb;pdb.set_trace() +pde = u.dt2 - u.laplace + +stencil = Eq(u.forward, solve(pde, u.forward)) + +# print("Init Devito linalg norm 0 :", np.linalg.norm(u.data[0])) +# print("Init Devito linalg norm 1 :", np.linalg.norm(u.data[1])) +# print("Init Devito linalg norm 2 :", np.linalg.norm(u.data[2])) +# print("Norm of initial data:", norm(u)) + +configuration['mpi'] = 0 +u2.data[:] = u.data[:] +configuration['mpi'] = mpiconf + +u.data[:] = np.load("so%s_wave_dat%s.npy" % (so, shape_str), allow_pickle=True) +dt = np.load("so%s_critical_dt%s.npy" % (so, shape_str), allow_pickle=True) + +# np.save("critical_dt%s.npy" % shape_str, model.critical_dt, allow_pickle=True) +# np.save("wave_dat%s.npy" % shape_str, u.data[:], allow_pickle=True) + +if len(shape) == 3 and args.plot: + plot_3dfunc(u) + +print("Init norm:", np.linalg.norm(u.data[:])) +# print("Init linalg norm:", np.linalg.norm(u.data[0])) +# print("Init linalg norm:", np.linalg.norm(u.data[1])) +# print("Init linalg norm:", np.linalg.norm(u.data[2])) + + +if args.devito: + # Run more with no sources now (Not supported in xdsl) + # op1 = Operator([stencil], name='DevitoOperator', subs=grid.spacing_map) + op1 = Operator([stencil], name='DevitoOperator') + op1.apply(time=nt, dt=dt) + + configuration['mpi'] = 0 + ub.data[:] = u.data[:] + configuration['mpi'] = mpiconf + + if len(shape) == 3 and args.plot: + plot_3dfunc(u) + + print("Devito norm:", norm(u)) + # print("Devito linalg norm 0:", np.linalg.norm(u.data[0])) + # print("Devito linalg norm 1:", np.linalg.norm(u.data[1])) + # print("Devito linalg norm 2:", np.linalg.norm(u.data[2])) + + +if args.xdsl: + + # Run more with no sources now (Not supported in xdsl) + xdslop = XDSLOperator([stencil], name='xDSLOperator') + xdslop.apply(time=nt, dt=dt) + + if len(shape) == 3 and args.plot: + plot_3dfunc(u) + + print("XDSL norm:", norm(u)) + + # print("XDSL output norm 0:", np.linalg.norm(u.data[0])) + # print("XDSL output norm 1:", np.linalg.norm(u.data[1])) + # print("XDSL output norm 2:", np.linalg.norm(u.data[2])) From 722e99810cd78fe631ea881f412e59071f191ba0 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Sat, 5 Aug 2023 20:48:58 +0300 Subject: [PATCH 22/25] bench: compress saved data --- fast/setup_wave3d.py | 4 +- fast/temp1 | 201 ---------------------------------------- fast/temp2 | 212 ------------------------------------------- fast/wave3d_b.py | 3 +- 4 files changed, 4 insertions(+), 416 deletions(-) delete mode 100644 fast/temp1 delete mode 100644 fast/temp2 diff --git a/fast/setup_wave3d.py b/fast/setup_wave3d.py index 795c96b284..b5c27b229d 100644 --- a/fast/setup_wave3d.py +++ b/fast/setup_wave3d.py @@ -112,4 +112,6 @@ shape_str = '_'.join(str(item) for item in shape) np.save("so%s_critical_dt%s.npy" % (so, shape_str), model.critical_dt, allow_pickle=True) np.save("so%s_wave_dat%s.npy" % (so, shape_str), u.data[:], allow_pickle=True) -np.save("so%s_grid_extent%s.npy" % (so, shape_str), model.grid.extent, allow_pickle=True) + +np.savez_compressed("so%s_grid_extent%s" % (so, shape_str), model.grid.extent, + allow_pickle=True) diff --git a/fast/temp1 b/fast/temp1 deleted file mode 100644 index 46f207a419..0000000000 --- a/fast/temp1 +++ /dev/null @@ -1,201 +0,0 @@ -module { - func.func @apply_kernel(%arg0: memref<260x260xf32>, %arg1: memref<260x260xf32>) -> memref<260x260xf32> attributes {param_names = ["u_vec_0", "u_vec_1"]} { - %c28_i64 = arith.constant 28 : i64 - %c12_i64 = arith.constant 12 : i64 - %c24_i64 = arith.constant 24 : i64 - %c8_i64 = arith.constant 8 : i64 - %c20_i64 = arith.constant 20 : i64 - %c16_i64 = arith.constant 16 : i64 - %c257 = arith.constant 257 : index - %cst = arith.constant 1.000000e-01 : f32 - %cst_0 = arith.constant -2.000000e+00 : f32 - %c-2_i64 = arith.constant -2 : i64 - %cst_1 = arith.constant 0.00392156886 : f32 - %c-1_i64 = arith.constant -1 : i64 - %cst_2 = arith.constant 3.075740e-05 : f32 - %cst_3 = arith.constant 0.00999999977 : f32 - %c-1 = arith.constant -1 : index - %c64 = arith.constant 64 : index - %c738197504_i32 = arith.constant 738197504 : i32 - %c4_i64 = arith.constant 4 : i64 - %c-1_i32 = arith.constant -1 : i32 - %c4_i32 = arith.constant 4 : i32 - %c1_i32 = arith.constant 1 : i32 - %c0_i32 = arith.constant 0 : i32 - %c1275069450_i32 = arith.constant 1275069450 : i32 - %c66_i32 = arith.constant 66 : i32 - %c1_i64 = arith.constant 1 : i64 - %c1140850688_i32 = arith.constant 1140850688 : i32 - %c8_i32 = arith.constant 8 : i32 - %c0 = arith.constant 0 : index - %c1 = arith.constant 1 : index - %0 = llvm.alloca %c8_i32 x i32 {alignment = 32 : i64} : (i32) -> !llvm.ptr - %1 = llvm.alloca %c1_i64 x i32 {alignment = 32 : i64} : (i64) -> !llvm.ptr - %2 = call @MPI_Comm_rank(%c1140850688_i32, %1) : (i32, !llvm.ptr) -> i32 - %3 = llvm.load %1 : !llvm.ptr - %alloc = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr = memref.extract_aligned_pointer_as_index %alloc : memref<66xf32> -> index - %4 = arith.index_cast %intptr : index to i64 - %5 = llvm.inttoptr %4 : i64 to !llvm.ptr - %alloc_4 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_5 = memref.extract_aligned_pointer_as_index %alloc_4 : memref<66xf32> -> index - %6 = arith.index_cast %intptr_5 : index to i64 - %7 = llvm.inttoptr %6 : i64 to !llvm.ptr - %alloc_6 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_7 = memref.extract_aligned_pointer_as_index %alloc_6 : memref<66xf32> -> index - %8 = arith.index_cast %intptr_7 : index to i64 - %9 = llvm.inttoptr %8 : i64 to !llvm.ptr - %alloc_8 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_9 = memref.extract_aligned_pointer_as_index %alloc_8 : memref<66xf32> -> index - %10 = arith.index_cast %intptr_9 : index to i64 - %11 = llvm.inttoptr %10 : i64 to !llvm.ptr - %alloc_10 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_11 = memref.extract_aligned_pointer_as_index %alloc_10 : memref<66xf32> -> index - %12 = arith.index_cast %intptr_11 : index to i64 - %13 = llvm.inttoptr %12 : i64 to !llvm.ptr - %alloc_12 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_13 = memref.extract_aligned_pointer_as_index %alloc_12 : memref<66xf32> -> index - %14 = arith.index_cast %intptr_13 : index to i64 - %15 = llvm.inttoptr %14 : i64 to !llvm.ptr - %alloc_14 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_15 = memref.extract_aligned_pointer_as_index %alloc_14 : memref<66xf32> -> index - %16 = arith.index_cast %intptr_15 : index to i64 - %17 = llvm.inttoptr %16 : i64 to !llvm.ptr - %alloc_16 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_17 = memref.extract_aligned_pointer_as_index %alloc_16 : memref<66xf32> -> index - %18 = arith.index_cast %intptr_17 : index to i64 - %19 = llvm.inttoptr %18 : i64 to !llvm.ptr - %20 = arith.remui %3, %c4_i32 : i32 - %21 = arith.divui %3, %c4_i32 : i32 - %22 = arith.remui %21, %c4_i32 : i32 - %23 = arith.addi %22, %c-1_i32 : i32 - %24 = arith.cmpi sge, %23, %c0_i32 : i32 - %25 = arith.muli %23, %c4_i32 : i32 - %26 = arith.addi %20, %25 : i32 - %27 = llvm.ptrtoint %0 : !llvm.ptr to i64 - %28 = llvm.inttoptr %27 : i64 to !llvm.ptr - %29 = arith.addi %27, %c16_i64 : i64 - %30 = llvm.inttoptr %29 : i64 to !llvm.ptr - %31 = arith.addi %22, %c1_i32 : i32 - %32 = arith.cmpi slt, %31, %c4_i32 : i32 - %33 = arith.muli %31, %c4_i32 : i32 - %34 = arith.addi %20, %33 : i32 - %35 = arith.addi %27, %c4_i64 : i64 - %36 = llvm.inttoptr %35 : i64 to !llvm.ptr - %37 = arith.addi %27, %c20_i64 : i64 - %38 = llvm.inttoptr %37 : i64 to !llvm.ptr - %39 = arith.addi %20, %c-1_i32 : i32 - %40 = arith.cmpi sge, %39, %c0_i32 : i32 - %41 = arith.muli %22, %c4_i32 : i32 - %42 = arith.addi %39, %41 : i32 - %43 = arith.addi %27, %c8_i64 : i64 - %44 = llvm.inttoptr %43 : i64 to !llvm.ptr - %45 = arith.addi %27, %c24_i64 : i64 - %46 = llvm.inttoptr %45 : i64 to !llvm.ptr - %47 = arith.addi %20, %c1_i32 : i32 - %48 = arith.cmpi slt, %47, %c4_i32 : i32 - %49 = arith.addi %47, %41 : i32 - %50 = arith.addi %27, %c12_i64 : i64 - %51 = llvm.inttoptr %50 : i64 to !llvm.ptr - %52 = arith.addi %27, %c28_i64 : i64 - %53 = llvm.inttoptr %52 : i64 to !llvm.ptr - %54 = llvm.inttoptr %c1_i64 : i64 to !llvm.ptr - %55 = math.fpowi %cst_2, %c-1_i64 : f32, i64 - %56 = math.fpowi %cst_1, %c-2_i64 : f32, i64 - %57 = arith.mulf %56, %cst_0 : f32 - %58:2 = scf.for %arg2 = %c0 to %c257 step %c1 iter_args(%arg3 = %arg0, %arg4 = %arg1) -> (memref<260x260xf32>, memref<260x260xf32>) { - %subview = memref.subview %arg4[2, 2] [64, 64] [1, 1] : memref<260x260xf32> to memref<64x64xf32, strided<[260, 1], offset: 522>> - %subview_18 = memref.subview %arg3[2, 2] [66, 66] [1, 1] : memref<260x260xf32> to memref<66x66xf32, strided<[260, 1], offset: 522>> - scf.if %24 { - %subview_19 = memref.subview %subview_18[-1, 0] [66, 1] [1, 1] : memref<66x66xf32, strided<[260, 1], offset: 522>> to memref<66xf32, strided<[260], offset: 262>> - memref.copy %subview_19, %alloc : memref<66xf32, strided<[260], offset: 262>> to memref<66xf32> - %60 = func.call @MPI_Isend(%5, %c66_i32, %c1275069450_i32, %26, %c0_i32, %c1140850688_i32, %28) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - %61 = func.call @MPI_Irecv(%7, %c66_i32, %c1275069450_i32, %26, %c0_i32, %c1140850688_i32, %30) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - } else { - llvm.store %c738197504_i32, %28 : !llvm.ptr - llvm.store %c738197504_i32, %30 : !llvm.ptr - } - scf.if %32 { - %subview_19 = memref.subview %subview_18[-1, 63] [66, 1] [1, 1] : memref<66x66xf32, strided<[260, 1], offset: 522>> to memref<66xf32, strided<[260], offset: 325>> - memref.copy %subview_19, %alloc_6 : memref<66xf32, strided<[260], offset: 325>> to memref<66xf32> - %60 = func.call @MPI_Isend(%9, %c66_i32, %c1275069450_i32, %34, %c0_i32, %c1140850688_i32, %36) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - %61 = func.call @MPI_Irecv(%11, %c66_i32, %c1275069450_i32, %34, %c0_i32, %c1140850688_i32, %38) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - } else { - llvm.store %c738197504_i32, %36 : !llvm.ptr - llvm.store %c738197504_i32, %38 : !llvm.ptr - } - scf.if %40 { - %subview_19 = memref.subview %subview_18[0, -1] [1, 66] [1, 1] : memref<66x66xf32, strided<[260, 1], offset: 522>> to memref<66xf32, strided<[1], offset: 521>> - memref.copy %subview_19, %alloc_10 : memref<66xf32, strided<[1], offset: 521>> to memref<66xf32> - %60 = func.call @MPI_Isend(%13, %c66_i32, %c1275069450_i32, %42, %c0_i32, %c1140850688_i32, %44) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - %61 = func.call @MPI_Irecv(%15, %c66_i32, %c1275069450_i32, %42, %c0_i32, %c1140850688_i32, %46) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - } else { - llvm.store %c738197504_i32, %44 : !llvm.ptr - llvm.store %c738197504_i32, %46 : !llvm.ptr - } - scf.if %48 { - %subview_19 = memref.subview %subview_18[63, -1] [1, 66] [1, 1] : memref<66x66xf32, strided<[260, 1], offset: 522>> to memref<66xf32, strided<[1], offset: 16901>> - memref.copy %subview_19, %alloc_14 : memref<66xf32, strided<[1], offset: 16901>> to memref<66xf32> - %60 = func.call @MPI_Isend(%17, %c66_i32, %c1275069450_i32, %49, %c0_i32, %c1140850688_i32, %51) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - %61 = func.call @MPI_Irecv(%19, %c66_i32, %c1275069450_i32, %49, %c0_i32, %c1140850688_i32, %53) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - } else { - llvm.store %c738197504_i32, %51 : !llvm.ptr - llvm.store %c738197504_i32, %53 : !llvm.ptr - } - %59 = func.call @MPI_Waitall(%c8_i32, %0, %54) : (i32, !llvm.ptr, !llvm.ptr) -> i32 - scf.if %24 { - %subview_19 = memref.subview %subview_18[-1, -1] [66, 1] [1, 1] : memref<66x66xf32, strided<[260, 1], offset: 522>> to memref<66xf32, strided<[260], offset: 261>> - memref.copy %subview_19, %alloc_4 : memref<66xf32, strided<[260], offset: 261>> to memref<66xf32> - } - scf.if %32 { - %subview_19 = memref.subview %subview_18[-1, 64] [66, 1] [1, 1] : memref<66x66xf32, strided<[260, 1], offset: 522>> to memref<66xf32, strided<[260], offset: 326>> - memref.copy %subview_19, %alloc_8 : memref<66xf32, strided<[260], offset: 326>> to memref<66xf32> - } - scf.if %40 { - %subview_19 = memref.subview %subview_18[-1, -1] [1, 66] [1, 1] : memref<66x66xf32, strided<[260, 1], offset: 522>> to memref<66xf32, strided<[1], offset: 261>> - memref.copy %subview_19, %alloc_12 : memref<66xf32, strided<[1], offset: 261>> to memref<66xf32> - } - scf.if %48 { - %subview_19 = memref.subview %subview_18[64, -1] [1, 66] [1, 1] : memref<66x66xf32, strided<[260, 1], offset: 522>> to memref<66xf32, strided<[1], offset: 17161>> - memref.copy %subview_19, %alloc_16 : memref<66xf32, strided<[1], offset: 17161>> to memref<66xf32> - } - scf.parallel (%arg5) = (%c0) to (%c64) step (%c1) { - %60 = arith.addi %arg5, %c-1 : index - %61 = arith.addi %arg5, %c1 : index - scf.for %arg6 = %c0 to %c64 step %c1 { - %62 = memref.load %subview_18[%arg5, %arg6] : memref<66x66xf32, strided<[260, 1], offset: 522>> - %63 = memref.load %subview_18[%60, %arg6] : memref<66x66xf32, strided<[260, 1], offset: 522>> - %64 = memref.load %subview_18[%61, %arg6] : memref<66x66xf32, strided<[260, 1], offset: 522>> - %65 = arith.addi %arg6, %c-1 : index - %66 = memref.load %subview_18[%arg5, %65] : memref<66x66xf32, strided<[260, 1], offset: 522>> - %67 = arith.addi %arg6, %c1 : index - %68 = memref.load %subview_18[%arg5, %67] : memref<66x66xf32, strided<[260, 1], offset: 522>> - %69 = arith.mulf %55, %62 : f32 - %70 = arith.mulf %56, %63 : f32 - %71 = arith.mulf %56, %64 : f32 - %72 = arith.mulf %57, %62 : f32 - %73 = arith.addf %70, %71 : f32 - %74 = arith.addf %73, %72 : f32 - %75 = arith.mulf %56, %66 : f32 - %76 = arith.mulf %56, %68 : f32 - %77 = arith.addf %75, %76 : f32 - %78 = arith.addf %77, %72 : f32 - %79 = arith.addf %74, %78 : f32 - %80 = arith.mulf %79, %cst : f32 - %81 = arith.addf %69, %cst_3 : f32 - %82 = arith.addf %81, %80 : f32 - %83 = arith.mulf %82, %cst_2 : f32 - memref.store %83, %subview[%arg5, %arg6] : memref<64x64xf32, strided<[260, 1], offset: 522>> - } - scf.yield - } - scf.yield %arg4, %arg3 : memref<260x260xf32>, memref<260x260xf32> - } - return %58#0 : memref<260x260xf32> - } - func.func private @MPI_Comm_rank(i32, !llvm.ptr) -> i32 - func.func private @MPI_Isend(!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - func.func private @MPI_Irecv(!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - func.func private @MPI_Waitall(i32, !llvm.ptr, !llvm.ptr) -> i32 -} - diff --git a/fast/temp2 b/fast/temp2 deleted file mode 100644 index 40783a347a..0000000000 --- a/fast/temp2 +++ /dev/null @@ -1,212 +0,0 @@ -#map = affine_map<()[s0] -> (s0 + 2)> -module { - func.func @apply_kernel(%arg0: memref<260x260xf32>, %arg1: memref<260x260xf32>) -> memref<260x260xf32> attributes {param_names = ["u_vec_0", "u_vec_1"]} { - %c28_i64 = arith.constant 28 : i64 - %c12_i64 = arith.constant 12 : i64 - %c24_i64 = arith.constant 24 : i64 - %c8_i64 = arith.constant 8 : i64 - %c20_i64 = arith.constant 20 : i64 - %c16_i64 = arith.constant 16 : i64 - %c257 = arith.constant 257 : index - %cst = arith.constant 1.000000e-01 : f32 - %cst_0 = arith.constant -2.000000e+00 : f32 - %c-2_i64 = arith.constant -2 : i64 - %cst_1 = arith.constant 0.00392156886 : f32 - %c-1_i64 = arith.constant -1 : i64 - %cst_2 = arith.constant 3.075740e-05 : f32 - %cst_3 = arith.constant 0.00999999977 : f32 - %c-1 = arith.constant -1 : index - %c64 = arith.constant 64 : index - %c738197504_i32 = arith.constant 738197504 : i32 - %c4_i64 = arith.constant 4 : i64 - %c-1_i32 = arith.constant -1 : i32 - %c4_i32 = arith.constant 4 : i32 - %c1_i32 = arith.constant 1 : i32 - %c0_i32 = arith.constant 0 : i32 - %c1275069450_i32 = arith.constant 1275069450 : i32 - %c66_i32 = arith.constant 66 : i32 - %c1_i64 = arith.constant 1 : i64 - %c1140850688_i32 = arith.constant 1140850688 : i32 - %c8_i32 = arith.constant 8 : i32 - %c0 = arith.constant 0 : index - %c1 = arith.constant 1 : index - %0 = llvm.alloca %c8_i32 x i32 {alignment = 32 : i64} : (i32) -> !llvm.ptr - %1 = llvm.alloca %c1_i64 x i32 {alignment = 32 : i64} : (i64) -> !llvm.ptr - %2 = call @MPI_Comm_rank(%c1140850688_i32, %1) : (i32, !llvm.ptr) -> i32 - %3 = llvm.load %1 : !llvm.ptr - %alloc = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr = memref.extract_aligned_pointer_as_index %alloc : memref<66xf32> -> index - %4 = arith.index_cast %intptr : index to i64 - %5 = llvm.inttoptr %4 : i64 to !llvm.ptr - %alloc_4 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_5 = memref.extract_aligned_pointer_as_index %alloc_4 : memref<66xf32> -> index - %6 = arith.index_cast %intptr_5 : index to i64 - %7 = llvm.inttoptr %6 : i64 to !llvm.ptr - %alloc_6 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_7 = memref.extract_aligned_pointer_as_index %alloc_6 : memref<66xf32> -> index - %8 = arith.index_cast %intptr_7 : index to i64 - %9 = llvm.inttoptr %8 : i64 to !llvm.ptr - %alloc_8 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_9 = memref.extract_aligned_pointer_as_index %alloc_8 : memref<66xf32> -> index - %10 = arith.index_cast %intptr_9 : index to i64 - %11 = llvm.inttoptr %10 : i64 to !llvm.ptr - %alloc_10 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_11 = memref.extract_aligned_pointer_as_index %alloc_10 : memref<66xf32> -> index - %12 = arith.index_cast %intptr_11 : index to i64 - %13 = llvm.inttoptr %12 : i64 to !llvm.ptr - %alloc_12 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_13 = memref.extract_aligned_pointer_as_index %alloc_12 : memref<66xf32> -> index - %14 = arith.index_cast %intptr_13 : index to i64 - %15 = llvm.inttoptr %14 : i64 to !llvm.ptr - %alloc_14 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_15 = memref.extract_aligned_pointer_as_index %alloc_14 : memref<66xf32> -> index - %16 = arith.index_cast %intptr_15 : index to i64 - %17 = llvm.inttoptr %16 : i64 to !llvm.ptr - %alloc_16 = memref.alloc() {alignment = 64 : i64} : memref<66xf32> - %intptr_17 = memref.extract_aligned_pointer_as_index %alloc_16 : memref<66xf32> -> index - %18 = arith.index_cast %intptr_17 : index to i64 - %19 = llvm.inttoptr %18 : i64 to !llvm.ptr - %20 = arith.remui %3, %c4_i32 : i32 - %21 = arith.divui %3, %c4_i32 : i32 - %22 = arith.remui %21, %c4_i32 : i32 - %23 = arith.addi %22, %c-1_i32 : i32 - %24 = arith.cmpi sge, %23, %c0_i32 : i32 - %25 = arith.muli %23, %c4_i32 : i32 - %26 = arith.addi %20, %25 : i32 - %27 = llvm.ptrtoint %0 : !llvm.ptr to i64 - %28 = llvm.inttoptr %27 : i64 to !llvm.ptr - %29 = arith.addi %27, %c16_i64 : i64 - %30 = llvm.inttoptr %29 : i64 to !llvm.ptr - %31 = arith.addi %22, %c1_i32 : i32 - %32 = arith.cmpi slt, %31, %c4_i32 : i32 - %33 = arith.muli %31, %c4_i32 : i32 - %34 = arith.addi %20, %33 : i32 - %35 = arith.addi %27, %c4_i64 : i64 - %36 = llvm.inttoptr %35 : i64 to !llvm.ptr - %37 = arith.addi %27, %c20_i64 : i64 - %38 = llvm.inttoptr %37 : i64 to !llvm.ptr - %39 = arith.addi %20, %c-1_i32 : i32 - %40 = arith.cmpi sge, %39, %c0_i32 : i32 - %41 = arith.muli %22, %c4_i32 : i32 - %42 = arith.addi %39, %41 : i32 - %43 = arith.addi %27, %c8_i64 : i64 - %44 = llvm.inttoptr %43 : i64 to !llvm.ptr - %45 = arith.addi %27, %c24_i64 : i64 - %46 = llvm.inttoptr %45 : i64 to !llvm.ptr - %47 = arith.addi %20, %c1_i32 : i32 - %48 = arith.cmpi slt, %47, %c4_i32 : i32 - %49 = arith.addi %47, %41 : i32 - %50 = arith.addi %27, %c12_i64 : i64 - %51 = llvm.inttoptr %50 : i64 to !llvm.ptr - %52 = arith.addi %27, %c28_i64 : i64 - %53 = llvm.inttoptr %52 : i64 to !llvm.ptr - %54 = llvm.inttoptr %c1_i64 : i64 to !llvm.ptr - %55 = math.fpowi %cst_2, %c-1_i64 : f32, i64 - %56 = math.fpowi %cst_1, %c-2_i64 : f32, i64 - %57 = arith.mulf %56, %cst_0 : f32 - %58:2 = scf.for %arg2 = %c0 to %c257 step %c1 iter_args(%arg3 = %arg0, %arg4 = %arg1) -> (memref<260x260xf32>, memref<260x260xf32>) { - scf.if %24 { - %subview = memref.subview %arg3[1, 2] [66, 1] [1, 1] : memref<260x260xf32> to memref<66xf32, strided<[260], offset: 262>> - memref.copy %subview, %alloc : memref<66xf32, strided<[260], offset: 262>> to memref<66xf32> - %60 = func.call @MPI_Isend(%5, %c66_i32, %c1275069450_i32, %26, %c0_i32, %c1140850688_i32, %28) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - %61 = func.call @MPI_Irecv(%7, %c66_i32, %c1275069450_i32, %26, %c0_i32, %c1140850688_i32, %30) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - } else { - llvm.store %c738197504_i32, %28 : !llvm.ptr - llvm.store %c738197504_i32, %30 : !llvm.ptr - } - scf.if %32 { - %subview = memref.subview %arg3[1, 65] [66, 1] [1, 1] : memref<260x260xf32> to memref<66xf32, strided<[260], offset: 325>> - memref.copy %subview, %alloc_6 : memref<66xf32, strided<[260], offset: 325>> to memref<66xf32> - %60 = func.call @MPI_Isend(%9, %c66_i32, %c1275069450_i32, %34, %c0_i32, %c1140850688_i32, %36) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - %61 = func.call @MPI_Irecv(%11, %c66_i32, %c1275069450_i32, %34, %c0_i32, %c1140850688_i32, %38) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - } else { - llvm.store %c738197504_i32, %36 : !llvm.ptr - llvm.store %c738197504_i32, %38 : !llvm.ptr - } - scf.if %40 { - %subview = memref.subview %arg3[2, 1] [1, 66] [1, 1] : memref<260x260xf32> to memref<66xf32, strided<[1], offset: 521>> - memref.copy %subview, %alloc_10 : memref<66xf32, strided<[1], offset: 521>> to memref<66xf32> - %60 = func.call @MPI_Isend(%13, %c66_i32, %c1275069450_i32, %42, %c0_i32, %c1140850688_i32, %44) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - %61 = func.call @MPI_Irecv(%15, %c66_i32, %c1275069450_i32, %42, %c0_i32, %c1140850688_i32, %46) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - } else { - llvm.store %c738197504_i32, %44 : !llvm.ptr - llvm.store %c738197504_i32, %46 : !llvm.ptr - } - scf.if %48 { - %subview = memref.subview %arg3[65, 1] [1, 66] [1, 1] : memref<260x260xf32> to memref<66xf32, strided<[1], offset: 16901>> - memref.copy %subview, %alloc_14 : memref<66xf32, strided<[1], offset: 16901>> to memref<66xf32> - %60 = func.call @MPI_Isend(%17, %c66_i32, %c1275069450_i32, %49, %c0_i32, %c1140850688_i32, %51) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - %61 = func.call @MPI_Irecv(%19, %c66_i32, %c1275069450_i32, %49, %c0_i32, %c1140850688_i32, %53) : (!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - } else { - llvm.store %c738197504_i32, %51 : !llvm.ptr - llvm.store %c738197504_i32, %53 : !llvm.ptr - } - %59 = func.call @MPI_Waitall(%c8_i32, %0, %54) : (i32, !llvm.ptr, !llvm.ptr) -> i32 - scf.if %24 { - %subview = memref.subview %arg3[1, 1] [66, 1] [1, 1] : memref<260x260xf32> to memref<66xf32, strided<[260], offset: 261>> - memref.copy %subview, %alloc_4 : memref<66xf32, strided<[260], offset: 261>> to memref<66xf32> - } - scf.if %32 { - %subview = memref.subview %arg3[1, 66] [66, 1] [1, 1] : memref<260x260xf32> to memref<66xf32, strided<[260], offset: 326>> - memref.copy %subview, %alloc_8 : memref<66xf32, strided<[260], offset: 326>> to memref<66xf32> - } - scf.if %40 { - %subview = memref.subview %arg3[1, 1] [1, 66] [1, 1] : memref<260x260xf32> to memref<66xf32, strided<[1], offset: 261>> - memref.copy %subview, %alloc_12 : memref<66xf32, strided<[1], offset: 261>> to memref<66xf32> - } - scf.if %48 { - %subview = memref.subview %arg3[66, 1] [1, 66] [1, 1] : memref<260x260xf32> to memref<66xf32, strided<[1], offset: 17161>> - memref.copy %subview, %alloc_16 : memref<66xf32, strided<[1], offset: 17161>> to memref<66xf32> - } - scf.parallel (%arg5) = (%c0) to (%c64) step (%c1) { - %60 = arith.addi %arg5, %c-1 : index - %61 = arith.addi %arg5, %c1 : index - scf.for %arg6 = %c0 to %c64 step %c1 { - %62 = affine.apply #map()[%arg5] - %63 = affine.apply #map()[%arg6] - %64 = memref.load %arg3[%62, %63] : memref<260x260xf32> - %65 = affine.apply #map()[%60] - %66 = affine.apply #map()[%arg6] - %67 = memref.load %arg3[%65, %66] : memref<260x260xf32> - %68 = affine.apply #map()[%61] - %69 = affine.apply #map()[%arg6] - %70 = memref.load %arg3[%68, %69] : memref<260x260xf32> - %71 = arith.addi %arg6, %c-1 : index - %72 = affine.apply #map()[%arg5] - %73 = affine.apply #map()[%71] - %74 = memref.load %arg3[%72, %73] : memref<260x260xf32> - %75 = arith.addi %arg6, %c1 : index - %76 = affine.apply #map()[%arg5] - %77 = affine.apply #map()[%75] - %78 = memref.load %arg3[%76, %77] : memref<260x260xf32> - %79 = arith.mulf %55, %64 : f32 - %80 = arith.mulf %56, %67 : f32 - %81 = arith.mulf %56, %70 : f32 - %82 = arith.mulf %57, %64 : f32 - %83 = arith.addf %80, %81 : f32 - %84 = arith.addf %83, %82 : f32 - %85 = arith.mulf %56, %74 : f32 - %86 = arith.mulf %56, %78 : f32 - %87 = arith.addf %85, %86 : f32 - %88 = arith.addf %87, %82 : f32 - %89 = arith.addf %84, %88 : f32 - %90 = arith.mulf %89, %cst : f32 - %91 = arith.addf %79, %cst_3 : f32 - %92 = arith.addf %91, %90 : f32 - %93 = arith.mulf %92, %cst_2 : f32 - %94 = affine.apply #map()[%arg5] - %95 = affine.apply #map()[%arg6] - memref.store %93, %arg4[%94, %95] : memref<260x260xf32> - } - scf.yield - } - scf.yield %arg4, %arg3 : memref<260x260xf32>, memref<260x260xf32> - } - return %58#0 : memref<260x260xf32> - } - func.func private @MPI_Comm_rank(i32, !llvm.ptr) -> i32 - func.func private @MPI_Isend(!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - func.func private @MPI_Irecv(!llvm.ptr, i32, i32, i32, i32, i32, !llvm.ptr) -> i32 - func.func private @MPI_Waitall(i32, !llvm.ptr, !llvm.ptr) -> i32 -} - diff --git a/fast/wave3d_b.py b/fast/wave3d_b.py index 9e0aa32190..a2a412e8ef 100644 --- a/fast/wave3d_b.py +++ b/fast/wave3d_b.py @@ -44,8 +44,7 @@ spacing = as_tuple(10.0 for _ in range(len(shape))) # Grid spacing in m. The domain size is now 1km by 1km origin = as_tuple(0.0 for _ in range(len(shape))) # What is the location of the top left corner. domain_size = tuple((d-1) * s for d, s in zip(shape, spacing)) -extent = np.load("so%s_grid_extent%s.npy" % (so, shape_str), allow_pickle=True) - +extent = np.load("so%s_grid_extent%s.npz" % (so, shape_str), allow_pickle=True)['arr_0'] grid = Grid(shape=shape, extent=as_tuple(extent)) # With the velocity and model size defined, we can create the seismic model that From aa73e333a9769c2f059cf2b3360e08951e501cb8 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Sat, 5 Aug 2023 21:02:16 +0300 Subject: [PATCH 23/25] bench: compress properly u.data[:] --- fast/setup_wave3d.py | 2 +- fast/wave3d_b.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/fast/setup_wave3d.py b/fast/setup_wave3d.py index b5c27b229d..797d8257a8 100644 --- a/fast/setup_wave3d.py +++ b/fast/setup_wave3d.py @@ -111,7 +111,7 @@ # Save Data here shape_str = '_'.join(str(item) for item in shape) np.save("so%s_critical_dt%s.npy" % (so, shape_str), model.critical_dt, allow_pickle=True) -np.save("so%s_wave_dat%s.npy" % (so, shape_str), u.data[:], allow_pickle=True) +np.savez_compressed("so%s_wave_dat%s" % (so, shape_str), u.data[:], allow_pickle=True) np.savez_compressed("so%s_grid_extent%s" % (so, shape_str), model.grid.extent, allow_pickle=True) diff --git a/fast/wave3d_b.py b/fast/wave3d_b.py index a2a412e8ef..3e3081f597 100644 --- a/fast/wave3d_b.py +++ b/fast/wave3d_b.py @@ -44,7 +44,7 @@ spacing = as_tuple(10.0 for _ in range(len(shape))) # Grid spacing in m. The domain size is now 1km by 1km origin = as_tuple(0.0 for _ in range(len(shape))) # What is the location of the top left corner. domain_size = tuple((d-1) * s for d, s in zip(shape, spacing)) -extent = np.load("so%s_grid_extent%s.npz" % (so, shape_str), allow_pickle=True)['arr_0'] +extent = np.load("so%s_grid_extent%s.npz" % (so, shape_str))['arr_0'] grid = Grid(shape=shape, extent=as_tuple(extent)) # With the velocity and model size defined, we can create the seismic model that @@ -78,7 +78,7 @@ u2.data[:] = u.data[:] configuration['mpi'] = mpiconf -u.data[:] = np.load("so%s_wave_dat%s.npy" % (so, shape_str), allow_pickle=True) +u.data[:] = np.load("so%s_wave_dat%s.npz" % (so, shape_str), allow_pickle=True)['arr_0'] dt = np.load("so%s_critical_dt%s.npy" % (so, shape_str), allow_pickle=True) # np.save("critical_dt%s.npy" % shape_str, model.critical_dt, allow_pickle=True) From ad55988736aea65abc362d73ed8fb3ccf328ff20 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Sun, 6 Aug 2023 14:01:19 +0300 Subject: [PATCH 24/25] bench: More cleanup and tiling merge --- fast/bench_utils.py | 5 +++-- fast/diffusion_2D_wBCs.py | 6 +++--- fast/diffusion_3D_wBCs.py | 40 ++++++++------------------------------- 3 files changed, 14 insertions(+), 37 deletions(-) diff --git a/fast/bench_utils.py b/fast/bench_utils.py index fe1563bae5..a01a959759 100644 --- a/fast/bench_utils.py +++ b/fast/bench_utils.py @@ -1,3 +1,5 @@ +import matplotlib.pyplot as plt +import pyvista as pv import numpy as np from examples.seismic import plot_image @@ -12,8 +14,7 @@ def plot_2dfunc(u): def plot_3dfunc(u): # Plot a 3D structured grid using pyvista - import matplotlib.pyplot as plt - import pyvista as pv + cmap = plt.colormaps["viridis"] values = u.data[0, :, :, :] vistagrid = pv.ImageData() diff --git a/fast/diffusion_2D_wBCs.py b/fast/diffusion_2D_wBCs.py index 869d6d2e7d..917454487d 100644 --- a/fast/diffusion_2D_wBCs.py +++ b/fast/diffusion_2D_wBCs.py @@ -43,9 +43,6 @@ grid = Grid(shape=(nx, ny), extent=(2., 2.)) u = TimeFunction(name='u', grid=grid, space_order=so) -# Reset our data field and ICs -init_hat(field=u.data[0], dx=dx, dy=dy, value=1.) - a = Constant(name='a') # Create an equation with second-order derivatives # eq = Eq(u.dt, a * u.laplace, subdomain=grid.interior) @@ -53,6 +50,9 @@ stencil = solve(eq, u.forward) eq_stencil = Eq(u.forward, stencil) +# Reset our data field and ICs +init_hat(field=u.data[0], dx=dx, dy=dy, value=1.) + if args.devito: op = Operator([eq_stencil], name='DevitoOperator') op.apply(time=nt, dt=dt, a=nu) diff --git a/fast/diffusion_3D_wBCs.py b/fast/diffusion_3D_wBCs.py index 9220e25247..9d0715f340 100644 --- a/fast/diffusion_3D_wBCs.py +++ b/fast/diffusion_3D_wBCs.py @@ -5,7 +5,9 @@ import argparse import numpy as np -from devito import Grid, TimeFunction, Eq, solve, Operator, Constant, norm, XDSLOperator +from devito import (Grid, TimeFunction, Eq, solve, Operator, Constant, + norm, XDSLOperator) +from fast.bench_utils import plot_3dfunc parser = argparse.ArgumentParser(description='Process arguments.') @@ -24,25 +26,6 @@ parser.add_argument("-xdsl", "--xdsl", default=False, type=bool, help="xDSL run") args = parser.parse_args() - -def plot_3dfunc(u): - # Plot a 3D structured grid using pyvista - - import matplotlib.pyplot as plt - import pyvista as pv - - cmap = plt.cm.get_cmap("viridis") - values = u.data[0, :, :, :] - vistagrid = pv.ImageData() - vistagrid.dimensions = np.array(values.shape) + 1 - vistagrid.spacing = (1, 1, 1) - vistagrid.origin = (0, 0, 0) # The bottom left corner of the data set - vistagrid.cell_data["values"] = values.flatten(order="F") - vistaslices = vistagrid.slice_orthogonal() - # vistagrid.plot(show_edges=True) - vistaslices.plot(cmap=cmap) - - # Some variable declarations nx, ny, nz = args.shape nt = args.nt @@ -61,29 +44,24 @@ def plot_3dfunc(u): grid = Grid(shape=(nx, ny, nz), extent=(2., 2., 2.)) u = TimeFunction(name='u', grid=grid, space_order=so) -# init_hat(field=u.data[0], dx=dx, dy=dy, value=2.) -u.data[:, :, :, :] = 0 -u.data[:, :, :, int(nz/2)] = 1 a = Constant(name='a') # Create an equation with second-order derivatives eq = Eq(u.dt, a * u.laplace) - stencil = solve(eq, u.forward) eq_stencil = Eq(u.forward, stencil) -# Create boundary condition expressions -x, y, z = grid.dimensions -t = grid.stepping_dim - -# print(eq_stencil) # Create Operator if args.devito: + u.data[:, :, :, :] = 0 + u.data[:, :, :, int(nz/2)] = 1 op = Operator([eq_stencil], name='DevitoOperator') # Apply the operator for a number of timesteps op.apply(time=nt, dt=dt, a=nu) print("Devito Field norm is:", norm(u)) + if args.plot: + plot_3dfunc(u) if args.xdsl: # Reset field @@ -92,8 +70,6 @@ def plot_3dfunc(u): xdslop = XDSLOperator([eq_stencil], name='xDSLOperator') # Apply the xdsl operator for a number of timesteps xdslop.apply(time=nt, dt=dt, a=nu) - + print("XDSL Field norm is:", norm(u)) if args.plot: plot_3dfunc(u) - - print("XDSL Field norm is:", norm(u)) From 6cfe5694a1d356e9bf56800ecffb388ae46e65a3 Mon Sep 17 00:00:00 2001 From: George Bisbas Date: Sun, 6 Aug 2023 14:15:33 +0300 Subject: [PATCH 25/25] bench: Hide pyvista req --- fast/bench_utils.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/fast/bench_utils.py b/fast/bench_utils.py index a01a959759..2efc9fa478 100644 --- a/fast/bench_utils.py +++ b/fast/bench_utils.py @@ -1,5 +1,4 @@ import matplotlib.pyplot as plt -import pyvista as pv import numpy as np from examples.seismic import plot_image @@ -14,7 +13,7 @@ def plot_2dfunc(u): def plot_3dfunc(u): # Plot a 3D structured grid using pyvista - + import pyvista as pv cmap = plt.colormaps["viridis"] values = u.data[0, :, :, :] vistagrid = pv.ImageData()