From 9de2a09796876ee55093cb32f052ea638ec660c9 Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Wed, 28 Aug 2024 16:03:34 +0100 Subject: [PATCH 1/6] reduce quadrature order for compressibly flow pressure terms from 6 to 4 --- utils/compressible_flow/forms.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/utils/compressible_flow/forms.py b/utils/compressible_flow/forms.py index 4c2c95a8..eb3db34e 100644 --- a/utils/compressible_flow/forms.py +++ b/utils/compressible_flow/forms.py @@ -47,7 +47,7 @@ def u_mass(u, w): def u_tendency(w, n, u, theta, rho, gas, Up, mu=None, f=None, F=None, - pi_degree=6): + pi_degree=4): """ Written in a dimension agnostic way """ @@ -86,12 +86,13 @@ def form_mass(u, rho, theta, du, drho, dtheta): def get_form_function(n, Up, c_pen, gas, mu, f=None, F=None, - viscosity=None, diffusivity=None): + viscosity=None, diffusivity=None, + pi_degree=4): def form_function(u, rho, theta, du, drho, dtheta, t): eqn = theta_tendency(dtheta, u, theta, n, Up, c_pen) eqn += rho_tendency(drho, rho, u, n) eqn += u_tendency(du, n, u, theta, rho, - gas, Up, mu, f, F) + gas, Up, mu, f, F, pi_degree=pi_degree) if viscosity: eqn += form_viscosity(u, du, viscosity) if diffusivity: From 6a7dd1a2e3655d77b1152dc6d84951c54e6a95ec Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Wed, 28 Aug 2024 16:07:43 +0100 Subject: [PATCH 2/6] add reference state specifically used for auxblockpc --- asQ/preconditioners/auxiliary_blockpc.py | 40 +++++++++++++++++------- 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/asQ/preconditioners/auxiliary_blockpc.py b/asQ/preconditioners/auxiliary_blockpc.py index 4444b970..48c8329a 100644 --- a/asQ/preconditioners/auxiliary_blockpc.py +++ b/asQ/preconditioners/auxiliary_blockpc.py @@ -14,8 +14,11 @@ def _setup(self, pc, v, u): self.prefix = pc.getOptionsPrefix() + self._prefix self.options = PETSc.Options(self.prefix) - self.uref = appctx.get('uref') - assert self.uref.function_space() == u.function_space() + uref = appctx.get('uref') + self.uref = appctx.get('aux_uref', uref) + self.u = fd.Function(u.function_space()) + + self.frozen = self.options.getBool("frozen", default=False) self.bcs = appctx['bcs'] self.tref = appctx['tref'] @@ -26,6 +29,12 @@ def _setup(self, pc, v, u): self.form_mass = appctx.get('aux_form_mass', form_mass) self.form_function = appctx.get('aux_form_function', form_function) + def update(self, pc): + if not self.frozen: + self.update_state(pc) + super().update(pc) + return + class AuxiliaryRealBlockPC(AuxiliaryBlockPCBase): """ @@ -57,8 +66,12 @@ class AuxiliaryRealBlockPC(AuxiliaryBlockPCBase): 'aux_theta': Alternative implicit theta parameter. """ + def update_state(self, pc): + self.u.assign(self.uref) + def form(self, pc, v, u): self._setup(pc, v, u) + self.update_state(pc) dt = self.appctx['dt'] theta = self.appctx['theta'] @@ -66,7 +79,7 @@ def form(self, pc, v, u): dt = self.options.getReal('dt', default=dt) theta = self.options.getReal('theta', default=theta) - us = fd.split(self.uref) + us = fd.split(self.u) vs = fd.split(v) dt1 = fd.Constant(1/dt) @@ -75,7 +88,7 @@ def form(self, pc, v, u): M = self.form_mass(*fd.split(u), *vs) F = self.form_function(*us, *vs, self.tref) - K = fd.derivative(F, self.uref) + K = fd.derivative(F, self.u) a = dt1*M + thet*K @@ -118,10 +131,18 @@ class AuxiliaryComplexBlockPC(AuxiliaryBlockPCBase): 'aux_d2i': Imaginary part of an alternative complex coefficient on the stiffness matrix. """ + def update_state(self, pc): + if self.u.function_space() == self.uref.function_space(): + self.u.assign(self.uref) + else: + self.cpx.set_real(self.u, self.uref) + self.cpx.set_imag(self.u, self.uref) + def form(self, pc, v, u): self._setup(pc, v, u) - cpx = self.appctx['cpx'] + self.cpx = cpx + self.update_state(pc) d1 = self.appctx['d1'] d2 = self.appctx['d2'] @@ -138,15 +159,10 @@ def form(self, pc, v, u): # complex and real valued function spaces W = v.function_space() - V = W.sub(0) - - bcs = tuple((cb - for bc in self.bcs - for cb in cpx.DirichletBC(W, V, bc, 0*bc.function_arg))) M = cpx.BilinearForm(W, d1, self.form_mass) - K = cpx.derivative(d2, partial(self.form_function, t=self.tref), self.uref) + K = cpx.derivative(d2, partial(self.form_function, t=self.tref), self.u) a = M + K - return (a, bcs) + return (a, self.bcs) From 6b5901f1f39d948c69d8fa50d2171617b258a8dd Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Wed, 28 Aug 2024 17:03:21 +0100 Subject: [PATCH 3/6] serial non-hydrostatic gravity wave vertical slice script --- .../vertical_slice/gravity_wave_nh_serial.py | 318 ++++++++++++++++++ 1 file changed, 318 insertions(+) create mode 100644 case_studies/vertical_slice/gravity_wave_nh_serial.py diff --git a/case_studies/vertical_slice/gravity_wave_nh_serial.py b/case_studies/vertical_slice/gravity_wave_nh_serial.py new file mode 100644 index 00000000..2da80bcb --- /dev/null +++ b/case_studies/vertical_slice/gravity_wave_nh_serial.py @@ -0,0 +1,318 @@ +import firedrake as fd +from firedrake.petsc import PETSc +from firedrake.output import VTKFile +from pyop2.mpi import MPI +from utils.diagnostics import convective_cfl_calculator +from utils.serial import SerialMiniApp +from utils import compressible_flow as euler +import numpy as np + + +def initial_conditions(mesh, W, Vv, gas, H, perturbation=True, hydrostatic=False): + if hydrostatic: + return NotImplementedError + + x, z = fd.SpatialCoordinate(mesh) + V2 = W.subfunctions[1] + up = fd.as_vector([fd.Constant(0.0), fd.Constant(1.0)]) # up direction + + Un = fd.Function(W) + un, rhon, thetan = Un.subfunctions + + Tsurf = fd.Constant(300.) + thetab = Tsurf*fd.exp(gas.N**2*z/gas.g) + + thetan.interpolate(thetab) + + Pi = fd.Function(V2) + + euler.hydrostatic_rho(Vv, V2, mesh, thetan, rhon, + pi_boundary=fd.Constant(1.0), + gas=gas, Up=up, top=True, Pi=Pi) + + un.project(fd.as_vector([20.0, 0.0])) + + Uback = Un.copy(deepcopy=True) + + if perturbation: + a = fd.Constant(5e3) + dtheta = fd.Constant(1e-2) + + theta_pert = dtheta*fd.sin(np.pi*z/H)/(1 + x**2/a**2) + thetan.interpolate(thetan + theta_pert) + + return Un, Uback + + +import argparse +parser = argparse.ArgumentParser(description='Nonhydrostatic gravity wave.') +parser.add_argument('--nlayers', type=int, default=5, help='Number of layers.') +parser.add_argument('--ncolumns', type=int, default=150, help='Number of columns.') +parser.add_argument('--nt', type=int, default=10, help='Number of timesteps to solve.') +parser.add_argument('--dt', type=float, default=12, help='Timestep in seconds.') +parser.add_argument('--atol', type=float, default=1e-4, help='Average absolute tolerance for each timestep') +parser.add_argument('--rtol', type=float, default=1e-8, help='Relative absolute tolerance for each timestep') +parser.add_argument('--degree', type=int, default=1, help='Degree of finite element space (the DG space). Default 1.') +parser.add_argument('--filename', type=str, default='gravity_wave_nh', help='Name of vtk file.') +parser.add_argument('--write_file', action='store_true', help='Write vtk file at end of each timestep.') +parser.add_argument('--output_freq', type=int, default=10, help='How often to write to file.') +parser.add_argument('--show_args', action='store_true', help='Output all the arguments.') + +args = parser.parse_known_args() +args = args[0] + +if args.show_args: + PETSc.Sys.Print(args) + +PETSc.Sys.Print("Setting up problem") + +comm = fd.COMM_WORLD + +# parameters + +output_freq = args.output_freq +nt = args.nt +dt = args.dt + +# set up the mesh +L = 300e3 +H = 10e3 + +distribution_parameters = { + "partition": True, + "overlap_type": (fd.DistributedMeshOverlapType.VERTEX, 2) +} + +base_mesh = fd.PeriodicIntervalMesh( + args.ncolumns, float(L), + distribution_parameters=distribution_parameters, + comm=comm) +base_mesh.coordinates.dat.data[:] -= float(L)/2 + +mesh = fd.ExtrudedMesh(base_mesh, + layers=args.nlayers, + layer_height=float(H/args.nlayers)) +n = fd.FacetNormal(mesh) + +# function spaces + +W, Vv = euler.function_space(mesh, horizontal_degree=args.degree, + vertical_degree=args.degree, + vertical_velocity_space=True) +V1, V2, Vt = W.subfunctions # velocity, density, temperature + +PETSc.Sys.Print(f"DoFs: {W.dim()}") +PETSc.Sys.Print(f"DoFs/core: {W.dim()/comm.size}") +PETSc.Sys.Print("") + +PETSc.Sys.Print("Calculating initial condiions") + +gas = euler.StandardAtmosphere(N=0.01) + +Un, Uback = initial_conditions(mesh, W, Vv, gas, H) +uback, rho_back, theta_back = Uback.subfunctions + +U0 = Uback.copy(deepcopy=True) # background state at rest +U0.subfunctions[0].assign(0) + +form_mass = euler.get_form_mass() + +up = fd.as_vector([fd.Constant(0.0), fd.Constant(1.0)]) # up direction + +form_function = euler.get_form_function( + n, up, c_pen=fd.Constant(2.0**(-7./2)), gas=gas, mu=None) + +bcs = [fd.DirichletBC(W.sub(0), 0., "bottom"), + fd.DirichletBC(W.sub(0), 0., "top")] +for bc in bcs: + bc.apply(Un) + +# Parameters for the newton iterations +patch_parameters = { + 'pc_type': 'python', + 'pc_python_type': 'firedrake.AssembledPC', + 'assembled': { + 'pc_type': 'python', + 'pc_python_type': 'firedrake.ASMStarPC', + 'pc_star': { + 'construct_dim': 0, + 'sub_sub_pc_type': 'lu', + 'sub_sub_pc_factor_mat_solver_type': 'petsc', + }, + }, +} + +aux_parameters = { + 'pc_type': 'python', + 'pc_python_type': 'asQ.AuxiliaryRealBlockPC', + 'aux': { + 'frozen': None, + 'pc_type': 'lu', + 'pc_factor_mat_solver_type': 'mumps', + } +} + +composite_parameters = { + 'pc_type': 'composite', + 'pc_composite_type': 'multiplicative', + 'pc_composite_pcs': 'ksp,ksp', + 'sub_0_ksp': aux_parameters, + 'sub_0_ksp_ksp': { + 'type': 'gmres', + 'max_it': 2, + 'convergence_test': 'skip', + 'converged_maxits': None, + }, + 'sub_1_ksp': patch_parameters, + 'sub_1_ksp_ksp': { + 'type': 'gmres', + 'max_it': 2, + 'convergence_test': 'skip', + 'converged_maxits': None, + }, +} + +solver_parameters = { + 'snes': { + 'monitor': None, + 'converged_reason': None, + 'atol': args.atol, + "rtol": args.rtol, + 'ksp_ew': None, + 'ksp_ew_version': 2, + 'ksp_ew_rtol0': 1e-2, + 'ksp_ew_threshold': 1e-5, + 'lag_preconditioner': -2, + 'lag_preconditioner_persists': None, + 'lag_jacobian': -2, + 'lag_jacobian_persists': None, + }, + 'ksp_type': 'fgmres', + # 'ksp_pc_side': 'right', + # 'ksp_norm_type': 'unpreconditioned', + 'ksp': { + 'monitor': None, + 'converged_rate': None, + 'atol': 1.0*args.atol, + 'max_it': 50, + "min_it": 1, + # 'view': None, + }, +} + +solver_parameters.update(composite_parameters) + +appctx = {'aux_uref': U0} + +theta = 0.5 + +miniapp = SerialMiniApp(dt=dt, theta=theta, w_initial=Un, + form_mass=form_mass, + form_function=form_function, + solver_parameters=solver_parameters, + bcs=bcs, appctx=appctx) + +PETSc.Sys.Print("Solving problem") + +uout = fd.Function(V1, name='velocity') +thetaout = fd.Function(Vt, name='temperature') +rhoout = fd.Function(V2, name='density') + +ofile = VTKFile(f'output/{args.filename}.pvd', comm=comm) + + +def assign_out_functions(): + uout.assign(miniapp.w0.subfunctions[0]) + rhoout.assign(miniapp.w0.subfunctions[1]) + thetaout.assign(miniapp.w0.subfunctions[2]) + + rhoout.assign(rhoout - rho_back) + thetaout.assign(thetaout - theta_back) + + +def write_to_file(time): + ofile.write(uout, rhoout, thetaout, time=time) + + +PETSc.Sys.Print('### === --- Timestepping loop --- === ###') +PETSc.Sys.Print('') +linear_its = 0 +nonlinear_its = 0 +solver_time = [] + +cfl_calc = convective_cfl_calculator(mesh) +cfl_series = [] + + +def max_cfl(u, dt): + with cfl_calc(u, dt).dat.vec_ro as v: + return v.max()[1] + + +assign_out_functions() +cfl = max_cfl(uout, dt) +cfl_series.append(cfl) +PETSc.Sys.Print(f'Maximum initial CFL = {round(cfl, 4)}') +PETSc.Sys.Print('') + + +def preproc(app, it, time): + PETSc.Sys.Print('') + PETSc.Sys.Print(f'### === --- Calculating time-step {it} --- === ###') + PETSc.Sys.Print('') + comm.Barrier() + stime = MPI.Wtime() + solver_time.append(stime) + + +def postproc(app, it, time): + global linear_its + global nonlinear_its + + etime = MPI.Wtime() + stime = solver_time[-1] + duration = etime - stime + solver_time[-1] = duration + PETSc.Sys.Print('') + PETSc.Sys.Print(f'Timestep solution time: {duration}\n') + PETSc.Sys.Print('') + + linear_its += miniapp.nlsolver.snes.getLinearSolveIterations() + nonlinear_its += miniapp.nlsolver.snes.getIterationNumber() + + if (it % output_freq) == 0: + assign_out_functions() + if args.write_file: + write_to_file(time=time) + + cfl = max_cfl(uout, dt) + cfl_series.append(cfl) + PETSc.Sys.Print(f'Time = {time}') + PETSc.Sys.Print(f'Maximum CFL = {round(cfl, 4)}') + + +# solve for each window +miniapp.solve(nt=nt, + preproc=preproc, + postproc=postproc) + +PETSc.Sys.Print('') +PETSc.Sys.Print('### === --- Iteration counts --- === ###') +PETSc.Sys.Print('') + +PETSc.Sys.Print(f'linear iterations: {linear_its} | iterations per timestep: {linear_its/nt}') +PETSc.Sys.Print(f'nonlinear iterations: {nonlinear_its} | iterations per timestep: {nonlinear_its/nt}') +PETSc.Sys.Print('') + +PETSc.Sys.Print(f'Total DoFs: {W.dim()}') +PETSc.Sys.Print(f'Number of MPI ranks: {mesh.comm.size} ') +PETSc.Sys.Print(f'DoFs/rank: {W.dim()/mesh.comm.size}') +PETSc.Sys.Print('') + +if len(solver_time) > 1: + # solver_time = solver_time[1:] + solver_time[0] = solver_time[1] + +PETSc.Sys.Print(f'Total solution time: {sum(solver_time)}') +PETSc.Sys.Print(f'Average timestep solution time: {sum(solver_time)/len(solver_time)}') +PETSc.Sys.Print('') From 82b43b75932b8c51f053b9c1ef95090cacf06657 Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Thu, 29 Aug 2024 10:08:22 +0100 Subject: [PATCH 4/6] paradiag non-hydrostatic gravity wave vertical slice script --- .../vertical_slice/gravity_wave_nh.py | 358 ++++++++++++++++++ .../vertical_slice/gravity_wave_nh_serial.py | 9 +- 2 files changed, 364 insertions(+), 3 deletions(-) create mode 100644 case_studies/vertical_slice/gravity_wave_nh.py diff --git a/case_studies/vertical_slice/gravity_wave_nh.py b/case_studies/vertical_slice/gravity_wave_nh.py new file mode 100644 index 00000000..2dbb569b --- /dev/null +++ b/case_studies/vertical_slice/gravity_wave_nh.py @@ -0,0 +1,358 @@ +import firedrake as fd +from firedrake.petsc import PETSc +from firedrake.output import VTKFile +from pyop2.mpi import MPI +from utils.diagnostics import convective_cfl_calculator +from utils import compressible_flow as euler +from math import sqrt +import numpy as np +import asQ + + +def initial_conditions(mesh, W, Vv, gas, H, perturbation=True, hydrostatic=False): + if hydrostatic: + return NotImplementedError + + x, z = fd.SpatialCoordinate(mesh) + V2 = W.subfunctions[1] + up = fd.as_vector([fd.Constant(0.0), fd.Constant(1.0)]) # up direction + + Un = fd.Function(W) + un, rhon, thetan = Un.subfunctions + + Tsurf = fd.Constant(300.) + thetab = Tsurf*fd.exp(gas.N**2*z/gas.g) + + thetan.interpolate(thetab) + + Pi = fd.Function(V2) + + euler.hydrostatic_rho(Vv, V2, mesh, thetan, rhon, + pi_boundary=fd.Constant(1.0), + gas=gas, Up=up, top=True, Pi=Pi) + + un.project(fd.as_vector([20.0, 0.0])) + + Uback = Un.copy(deepcopy=True) + + if perturbation: + a = fd.Constant(5e3) + dtheta = fd.Constant(1e-2) + + theta_pert = dtheta*fd.sin(np.pi*z/H)/(1 + x**2/a**2) + thetan.interpolate(thetan + theta_pert) + + return Un, Uback + + +import argparse +parser = argparse.ArgumentParser(description='Nonhydrostatic gravity wave.') +parser.add_argument('--nlayers', type=int, default=5, help='Number of layers.') +parser.add_argument('--ncolumns', type=int, default=150, help='Number of columns.') +parser.add_argument('--dt', type=float, default=12, help='Timestep in seconds.') +parser.add_argument('--atol', type=float, default=1e-4, help='Average absolute tolerance for each timestep') +parser.add_argument('--nwindows', type=int, default=1, help='Number of windows to solve.') +parser.add_argument('--nslices', type=int, default=2, help='Number of slices in the all-at-once system.') +parser.add_argument('--slice_length', type=int, default=2, help='Number of timesteps in each slice of the all-at-once system.') +parser.add_argument('--alpha', type=float, default=1e-4, help='Circulant parameter') +parser.add_argument('--degree', type=int, default=1, help='Degree of finite element space (the DG space).') +parser.add_argument('--filename', type=str, default='gravity_wave_nh', help='Name of vtk file.') +parser.add_argument('--write_file', action='store_true', help='Write vtk file at end of each window.') +parser.add_argument('--metrics_dir', type=str, default='output', help='Directory to save paradiag metrics and vtk to.') +parser.add_argument('--show_args', action='store_true', help='Output all the arguments.') + +args = parser.parse_known_args() +args = args[0] + +if args.show_args: + PETSc.Sys.Print(args) + +PETSc.Sys.Print("Setting up problem") + +# set up the ensemble communicator for space-time parallelism +time_partition = tuple((args.slice_length for _ in range(args.nslices))) +window_length = sum(time_partition) + +global_comm = fd.COMM_WORLD +ensemble = asQ.create_ensemble(time_partition, comm=global_comm) + +dt = args.dt + +# set up the mesh +L = 300e3 +H = 10e3 + +distribution_parameters = { + "partition": True, + "overlap_type": (fd.DistributedMeshOverlapType.VERTEX, 2) +} + +base_mesh = fd.PeriodicIntervalMesh( + args.ncolumns, float(L), + distribution_parameters=distribution_parameters, + comm=ensemble.comm) +base_mesh.coordinates.dat.data[:] -= float(L)/2 + +mesh = fd.ExtrudedMesh(base_mesh, + layers=args.nlayers, + layer_height=float(H/args.nlayers)) +n = fd.FacetNormal(mesh) + +# function spaces + +W, Vv = euler.function_space(mesh, horizontal_degree=args.degree, + vertical_degree=args.degree, + vertical_velocity_space=True) +V1, V2, Vt = W.subfunctions # velocity, density, temperature + +PETSc.Sys.Print(f"DoFs/timestep: {W.dim()}") +PETSc.Sys.Print(f"DoFs/core: {W.dim()*window_length/global_comm.size}") +PETSc.Sys.Print(f"Block DoFs/core: {2*W.dim()/ensemble.comm.size}") +PETSc.Sys.Print("") + +PETSc.Sys.Print("Calculating initial condiions") + +gas = euler.StandardAtmosphere(N=0.01) + +Un, Uback = initial_conditions(mesh, W, Vv, gas, H) +uback, rho_back, theta_back = Uback.subfunctions + +U0 = Uback.copy(deepcopy=True) # background state at rest +U0.subfunctions[0].assign(0) + +# finite element forms + +form_mass = euler.get_form_mass() + +up = fd.as_vector([fd.Constant(0.0), fd.Constant(1.0)]) # up direction + +form_function = euler.get_form_function( + n, up, c_pen=fd.Constant(2.0**(-7./2)), gas=gas, mu=None) + +zv = fd.as_vector([fd.Constant(0.), fd.Constant(0.)]) +bcs = [fd.DirichletBC(W.sub(0), zv, "bottom"), + fd.DirichletBC(W.sub(0), zv, "top")] +for bc in bcs: + bc.apply(Un) + +# Parameters for the blocks +patch_parameters = { + 'pc_type': 'python', + 'pc_python_type': 'firedrake.AssembledPC', + 'assembled': { + 'pc_type': 'python', + 'pc_python_type': 'firedrake.ASMStarPC', + 'pc_star': { + 'construct_dim': 0, + 'sub_sub_pc_type': 'lu', + 'sub_sub_pc_factor_mat_solver_type': 'petsc', + }, + }, +} + +aux_parameters = { + 'pc_type': 'python', + 'pc_python_type': 'asQ.AuxiliaryComplexBlockPC', + 'aux': { + 'frozen': None, + 'pc_type': 'lu', + 'pc_factor_mat_solver_type': 'mumps', + } +} + +composite_parameters = { + 'pc_type': 'composite', + 'pc_composite_type': 'multiplicative', + 'pc_composite_pcs': 'ksp,ksp', + 'sub_0_ksp': aux_parameters, + 'sub_0_ksp_ksp': { + 'type': 'gmres', + 'max_it': 2, + 'convergence_test': 'skip', + 'converged_maxits': None, + }, + 'sub_1_ksp': patch_parameters, + 'sub_1_ksp_ksp': { + 'type': 'gmres', + 'max_it': 2, + 'convergence_test': 'skip', + 'converged_maxits': None, + }, +} + +block_parameters = { + 'mat_type': 'matfree', + 'ksp_type': 'fgmres', + # 'ksp_pc_side': 'right', + # 'ksp_norm_type': 'unpreconditioned', + 'ksp': { + 'rtol': args.alpha, + 'max_it': 8, + 'converged_maxits': None, + 'convergence_test': 'skip', + }, + 'snes_lag_preconditioner': -2, + 'snes_lag_preconditioner_persists': None, + 'snes_lag_jacobian': -2, + 'snes_lag_jacobian_persists': None, +} + +block_parameters.update(composite_parameters) + +atol = sqrt(window_length)*args.atol +solver_parameters_diag = { + "snes": { + "monitor": None, + "converged_reason": None, + "rtol": 1e-8, + "atol": atol, + "ksp_ew": None, + "ksp_ew_version": 1, + "ksp_ew_rtol0": 1e-1, + "ksp_ew_threshold": 1e-1, + 'lag_preconditioner': -2, + 'lag_preconditioner_persists': None, + # 'lag_jacobian': -2, + # 'lag_jacobian_persists': None, + }, + "ksp_type": "fgmres", + "mat_type": "matfree", + "ksp": { + "monitor": None, + "converged_rate": None, + "atol": atol, + }, + "pc_type": "python", + "pc_python_type": "asQ.CirculantPC", + "circulant_alpha": args.alpha, + "circulant_state": "window", + "circulant_block": block_parameters +} + +appctx = {'block_appctx': {'aux_uref': U0}} + +theta = 0.5 + +pdg = asQ.Paradiag(ensemble=ensemble, + time_partition=time_partition, + form_function=form_function, + form_mass=form_mass, + ics=Un, dt=dt, theta=theta, + bcs=bcs, appctx=appctx, + solver_parameters=solver_parameters_diag) + +aaofunc = pdg.aaofunc +is_last_slice = pdg.layout.is_local(-1) + +PETSc.Sys.Print("Solving problem") + +# only last slice does diagnostics/output +if is_last_slice: + uout = fd.Function(V1, name='velocity') + if args.write_file: + thetaout = fd.Function(Vt, name='temperature') + rhoout = fd.Function(V2, name='density') + + ofile = VTKFile(f'output/{args.filename}.pvd', + comm=ensemble.comm) + + def assign_out_functions(): + uout.assign(aaofunc[-1].subfunctions[0]) + if args.write_file: + rhoout.assign(aaofunc[-1].subfunctions[1]) + thetaout.assign(aaofunc[-1].subfunctions[2]) + + rhoout.assign(rhoout - rho_back) + thetaout.assign(thetaout - theta_back) + + def write_to_file(): + if args.write_file: + ofile.write(uout, rhoout, thetaout) + + cfl_calc = convective_cfl_calculator(mesh) + cfl_series = [] + + def max_cfl(u, dt): + with cfl_calc(u, dt).dat.vec_ro as v: + return v.max()[1] + +solver_time = [] + + +def window_preproc(pdg, wndw, rhs): + PETSc.Sys.Print('') + PETSc.Sys.Print(f'### === --- Calculating time-window {wndw} --- === ###') + PETSc.Sys.Print('') + ensemble.global_comm.Barrier() + stime = MPI.Wtime() + solver_time.append(stime) + + +def window_postproc(pdg, wndw, rhs): + etime = MPI.Wtime() + stime = solver_time[-1] + duration = etime - stime + solver_time[-1] = duration + PETSc.Sys.Print('', comm=global_comm) + PETSc.Sys.Print(f'Window solution time: {duration}', comm=global_comm) + PETSc.Sys.Print('', comm=global_comm) + + # postprocess this timeslice + if is_last_slice: + assign_out_functions() + if args.write_file: + write_to_file() + PETSc.Sys.Print('', comm=ensemble.comm) + + cfl = max_cfl(uout, dt) + cfl_series.append(cfl) + PETSc.Sys.Print(f'Maximum CFL = {cfl}', comm=ensemble.comm) + + +# solve for each window +pdg.solve(nwindows=args.nwindows, + preproc=window_preproc, + postproc=window_postproc) + +PETSc.Sys.Print('') +PETSc.Sys.Print('### === --- Iteration counts --- === ###') +PETSc.Sys.Print('') + +from asQ import write_paradiag_metrics +write_paradiag_metrics(pdg, directory=f'{args.metrics_dir}') + +nw = pdg.total_windows +nt = pdg.total_timesteps +PETSc.Sys.Print(f'windows: {nw}') +PETSc.Sys.Print(f'timesteps: {nt}') +PETSc.Sys.Print('') + +lits = pdg.linear_iterations +nlits = pdg.nonlinear_iterations +blits = pdg.block_iterations.data() + +PETSc.Sys.Print(f'linear iterations: {lits} | iterations per window: {lits/nw}') +PETSc.Sys.Print(f'nonlinear iterations: {nlits} | iterations per window: {nlits/nw}') +PETSc.Sys.Print(f'block linear iterations: {blits} | iterations per block solve: {blits/lits}') +PETSc.Sys.Print('') + +ensemble.global_comm.Barrier() +if is_last_slice: + PETSc.Sys.Print(f'Maximum CFL = {max(cfl_series)}', comm=ensemble.comm) + PETSc.Sys.Print(f'Minimum CFL = {min(cfl_series)}', comm=ensemble.comm) + PETSc.Sys.Print('', comm=ensemble.comm) +ensemble.global_comm.Barrier() + +PETSc.Sys.Print(f'DoFs per timestep: {W.dim()}', comm=global_comm) +PETSc.Sys.Print(f'Number of MPI ranks per timestep: {mesh.comm.size} ', comm=global_comm) +PETSc.Sys.Print(f'DoFs/rank: {W.dim()/mesh.comm.size}', comm=global_comm) +PETSc.Sys.Print(f'Block DoFs/rank: {2*W.dim()/mesh.comm.size}', comm=global_comm) +PETSc.Sys.Print('') + +if len(solver_time) > 1: + solver_time[0] = solver_time[1] + +PETSc.Sys.Print(f'Total solution time: {sum(solver_time)}', comm=global_comm) +PETSc.Sys.Print(f'Average window solution time: {sum(solver_time)/len(solver_time)}', comm=global_comm) +PETSc.Sys.Print(f'Average timestep solution time: {sum(solver_time)/(window_length*len(solver_time))}', comm=global_comm) +PETSc.Sys.Print('', comm=global_comm) diff --git a/case_studies/vertical_slice/gravity_wave_nh_serial.py b/case_studies/vertical_slice/gravity_wave_nh_serial.py index 2da80bcb..0a94bca7 100644 --- a/case_studies/vertical_slice/gravity_wave_nh_serial.py +++ b/case_studies/vertical_slice/gravity_wave_nh_serial.py @@ -3,8 +3,8 @@ from firedrake.output import VTKFile from pyop2.mpi import MPI from utils.diagnostics import convective_cfl_calculator -from utils.serial import SerialMiniApp from utils import compressible_flow as euler +from utils.serial import SerialMiniApp import numpy as np @@ -115,6 +115,8 @@ def initial_conditions(mesh, W, Vv, gas, H, perturbation=True, hydrostatic=False U0 = Uback.copy(deepcopy=True) # background state at rest U0.subfunctions[0].assign(0) +# finite element forms + form_mass = euler.get_form_mass() up = fd.as_vector([fd.Constant(0.0), fd.Constant(1.0)]) # up direction @@ -122,8 +124,9 @@ def initial_conditions(mesh, W, Vv, gas, H, perturbation=True, hydrostatic=False form_function = euler.get_form_function( n, up, c_pen=fd.Constant(2.0**(-7./2)), gas=gas, mu=None) -bcs = [fd.DirichletBC(W.sub(0), 0., "bottom"), - fd.DirichletBC(W.sub(0), 0., "top")] +zv = fd.as_vector([fd.Constant(0.), fd.Constant(0.)]) +bcs = [fd.DirichletBC(W.sub(0), zv, "bottom"), + fd.DirichletBC(W.sub(0), zv, "top")] for bc in bcs: bc.apply(Un) From 9a6c5f6d898ce0b1c92375859e4068819de95c9b Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Tue, 3 Sep 2024 10:19:23 +0100 Subject: [PATCH 5/6] update gw_nh params --- .../vertical_slice/gravity_wave_nh.py | 23 +++++++------------ .../vertical_slice/gravity_wave_nh_serial.py | 11 ++++----- 2 files changed, 12 insertions(+), 22 deletions(-) diff --git a/case_studies/vertical_slice/gravity_wave_nh.py b/case_studies/vertical_slice/gravity_wave_nh.py index 2dbb569b..6aaef78a 100644 --- a/case_studies/vertical_slice/gravity_wave_nh.py +++ b/case_studies/vertical_slice/gravity_wave_nh.py @@ -145,7 +145,7 @@ def initial_conditions(mesh, W, Vv, gas, H, perturbation=True, hydrostatic=False 'pc_star': { 'construct_dim': 0, 'sub_sub_pc_type': 'lu', - 'sub_sub_pc_factor_mat_solver_type': 'petsc', + 'sub_sub_pc_factor_mat_solver_type': 'mumps', }, }, } @@ -183,49 +183,41 @@ def initial_conditions(mesh, W, Vv, gas, H, perturbation=True, hydrostatic=False block_parameters = { 'mat_type': 'matfree', 'ksp_type': 'fgmres', - # 'ksp_pc_side': 'right', - # 'ksp_norm_type': 'unpreconditioned', 'ksp': { 'rtol': args.alpha, - 'max_it': 8, + 'max_it': 50, 'converged_maxits': None, 'convergence_test': 'skip', }, - 'snes_lag_preconditioner': -2, - 'snes_lag_preconditioner_persists': None, - 'snes_lag_jacobian': -2, - 'snes_lag_jacobian_persists': None, } block_parameters.update(composite_parameters) -atol = sqrt(window_length)*args.atol +patol = sqrt(window_length)*args.atol solver_parameters_diag = { "snes": { "monitor": None, "converged_reason": None, "rtol": 1e-8, - "atol": atol, + "atol": patol, "ksp_ew": None, "ksp_ew_version": 1, "ksp_ew_rtol0": 1e-1, "ksp_ew_threshold": 1e-1, 'lag_preconditioner': -2, 'lag_preconditioner_persists': None, - # 'lag_jacobian': -2, - # 'lag_jacobian_persists': None, }, - "ksp_type": "fgmres", "mat_type": "matfree", + "ksp_type": "fgmres", "ksp": { "monitor": None, "converged_rate": None, - "atol": atol, + "atol": patol, }, "pc_type": "python", "pc_python_type": "asQ.CirculantPC", + "circulant_state": "reference", "circulant_alpha": args.alpha, - "circulant_state": "window", "circulant_block": block_parameters } @@ -238,6 +230,7 @@ def initial_conditions(mesh, W, Vv, gas, H, perturbation=True, hydrostatic=False form_function=form_function, form_mass=form_mass, ics=Un, dt=dt, theta=theta, + reference_state=Uback, bcs=bcs, appctx=appctx, solver_parameters=solver_parameters_diag) diff --git a/case_studies/vertical_slice/gravity_wave_nh_serial.py b/case_studies/vertical_slice/gravity_wave_nh_serial.py index 0a94bca7..56a7f151 100644 --- a/case_studies/vertical_slice/gravity_wave_nh_serial.py +++ b/case_studies/vertical_slice/gravity_wave_nh_serial.py @@ -140,7 +140,7 @@ def initial_conditions(mesh, W, Vv, gas, H, perturbation=True, hydrostatic=False 'pc_star': { 'construct_dim': 0, 'sub_sub_pc_type': 'lu', - 'sub_sub_pc_factor_mat_solver_type': 'petsc', + 'sub_sub_pc_factor_mat_solver_type': 'mumps', }, }, } @@ -182,21 +182,18 @@ def initial_conditions(mesh, W, Vv, gas, H, perturbation=True, hydrostatic=False 'atol': args.atol, "rtol": args.rtol, 'ksp_ew': None, - 'ksp_ew_version': 2, + 'ksp_ew_version': 1, 'ksp_ew_rtol0': 1e-2, 'ksp_ew_threshold': 1e-5, 'lag_preconditioner': -2, 'lag_preconditioner_persists': None, - 'lag_jacobian': -2, - 'lag_jacobian_persists': None, }, + 'mat_type': 'matfree', 'ksp_type': 'fgmres', - # 'ksp_pc_side': 'right', - # 'ksp_norm_type': 'unpreconditioned', 'ksp': { 'monitor': None, 'converged_rate': None, - 'atol': 1.0*args.atol, + 'atol': args.atol, 'max_it': 50, "min_it": 1, # 'view': None, From 85f2b0a06c2ca56aa726fddf1b112ca06c852940 Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Tue, 3 Sep 2024 10:19:57 +0100 Subject: [PATCH 6/6] rcm ordering makes hydrostatic balance calculation unstable for some parameter combinations --- utils/compressible_flow/compressible_flow.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/compressible_flow/compressible_flow.py b/utils/compressible_flow/compressible_flow.py index aff10d25..b7fb358c 100644 --- a/utils/compressible_flow/compressible_flow.py +++ b/utils/compressible_flow/compressible_flow.py @@ -98,7 +98,7 @@ def hydrostatic_rho(Vv, V2, mesh, thetan, rhon, pi_boundary, 'mat_type': 'aij', 'ksp_type': 'preonly', 'pc_type': 'lu', - "pc_factor_mat_ordering_type": "rcm", + # "pc_factor_mat_ordering_type": "rcm", "pc_factor_mat_solver_type": "mumps", }