Using CFL-limited timestep for non VectorField; Compound bases ; Parallel with $mpiexec

[mohammed ali]

Hello Everyone, 

I am new to Dedalus and I have a couple of questions that hopefully someone can help me with.

I am trying to solve an NCC IVP PDE, so I started with a simple PDE to get familiar with the libraries, and my questions are :

1. Can we use CFL to integrate with different time steps when the Field is not VectorField because add_velocity(u)  requires that u be a vector?

Error: 

ValueError: Velocity must be a vector

2. Can we use Compound bases  in Dedalus  v3 (because according to the documentations it can be used in  Dedalus  v2  )?

3. Can we integrate in parallel using $mpiexec -n  like in the example Rayleigh-Benard convection (2D IVP) but for 1D problem?  By the way, I ran the Rayleigh-Benard convection 2D example and it worked,  but it does not work for my 1D PDE.

Error: 

ValueError: Mesh ([4]) must have lower dimension than distributor (1)

Thanks,

Mohammed

This is the Code:

"""

 Moahmmed Ali Mohammed

Date 07/19/2023

Dedalus v3 script simulating the 1D nonuniform transport PDE.

This script demonstrates solving a 1D initial value problem and produces

a space-time plot of the solution.

We use a Chyb basis to solve the IVP:

    dt(u) + c(x)*dx(u) = 0   , where c(x)=1/(x^2+1)

with ICs

    u((0,x)=1/((x+3)^2+1)

To run:

    $ python3 NCC_Example_1.py

"""

import numpy as np

import matplotlib.pyplot as plt

import dedalus.public as d3

import logging

import pathlib

import subprocess

import os

import h5py # HDF5 file manipulation

import scipy.io

from mpi4py import MPI

import time

import logging

logger = logging.getLogger(__name__)

# Parameters

L = 20

Nx = 801

dealias = 3/2

stop_sim_time = 10

timestepper = d3.RK222

max_timestep = 1e-1

dtype = np.float64

# Bases

x_coord = d3.Coordinate('x')

dist = d3.Distributor(x_coord, dtype=dtype)

xbasis = d3.Chebyshev(x_coord, size=Nx, bounds=(-L, L), dealias=dealias)

# Fields

u = dist.Field(name='u', bases=xbasis)

tau1 = dist.Field(name = 'tau1')

# tau polynomials

tau_basis = xbasis.derivative_basis(1)

p1 = dist.Field(bases = tau_basis)

p1['c'][-1] = 1

# NCC

x = dist.local_grid(xbasis)

c = dist.Field(bases=xbasis)

c['g'] = 1/(x**2+1)

# Substitutions

dx = lambda A: d3.Differentiate(A, x_coord)

# Problem

problem = d3.IVP([u,tau1], namespace=locals())

problem.add_equation("dt(u)+ tau1*p1 =-c*dx(u)")

# Add BC

problem.add_equation("u(x = -L) = 1/((-L+3)**2+1)")

# Initial conditions

u['g'] = 1/((x+3)**2+1)

# Solver

solver = problem.build_solver(timestepper)

solver.stop_sim_time = stop_sim_time

#storage

u_list = [np.copy(u['g'])]

t_list = [solver.sim_time]

# analysis

# Output solution fields

full_solution = solver.evaluator.add_file_handler('NCC_Example_2_Chy_solution2', iter =1 , max_writes =800000000)

full_solution.add_task(u, layout = 'g', name='u')

# Analysis

snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.1, max_writes=800000000)

snapshots.add_task(u, layout = 'g', name='u')

# CFL

CFL = d3.CFL(solver, initial_dt=max_timestep, cadence=10, safety=0.2, threshold=0.1,

             max_change=1.5, min_change=0.5, max_dt=max_timestep)

CFL.add_velocity(u)

# Solving/iterating

while solver.proceed:

    time_step = CFL.compute_timestep()

    solver.step(time_step)

    u.change_scales(1)

    u_list.append(np.copy(u['g']))

    t_list.append(solver.sim_time)

    if solver.iteration % 100 == 0:

        logger.info('Iteration = %i, time = %e, dt = %e' %(solver.iteration, solver.sim_time, time_step))

[Daniel Lecoanet]

Hi Mohammed,

1. I would define u to be a VectorField so it can be used with CFL. Since your problem is 1D, vector fields are equivalent to scalar fields.

2. We are working on implementing this feature in v3, but it is not possible to use yet.

3. 1D problems can only be run in serial.

Hope that helps,

Daniel

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/6b4c31d0-061b-4520-9a7b-2147ecbecdeen%40googlegroups.com.

[mohammed ali]

Hi Daniel,

Thanks a lot for the prompt response.

I followed your suggestion about changing u to be a vector field, and you were right it does not affect the behavior of the system. However, I am still running  into problems with 

CFL.add_velocity(u) 

***************************************************************************************************************

This is the error I am getting :

Traceback (most recent call last):
  File "/home/mohammed/Class_Dedalus_v1/NCC_Example/NCC_CFL_Example.py", line 98, in <module>
    CFL.add_velocity(u)
  File "/home/mohammed/mambaforge/envs/dedalus3/lib/python3.11/site-packages/dedalus/extras/flow_tools.py", line 232, in add_velocity
    cfl_operator = operators.AdvectiveCFL(velocity, coords[0])
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mohammed/mambaforge/envs/dedalus3/lib/python3.11/site-packages/dedalus/tools/dispatch.py", line 35, in __call__
    raise NotImplementedError("No subclasses of {} found for the supplied arguments: {}, {}".format(cls, args, kw))
NotImplementedError: No subclasses of <class 'dedalus.core.operators.AdvectiveCFL'> found for the supplied arguments: (<Field 140227416922960>, <dedalus.core.coords.Coordinate object at 0x7f8946a78690>), {}

*********************************************************************************************************************

And this is the full code

import numpy as np
import matplotlib.pyplot as plt
import dedalus.public as d3
import logging
import pathlib
import subprocess
import os
import h5py # HDF5 file manipulation
import scipy.io
from mpi4py import MPI
import time
import logging
logger = logging.getLogger(__name__)


# Parameters
L = 20
Nx = 801
dealias = 3/2
stop_sim_time = 10
timestepper = d3.RK222

max_timestep = 1e-3

dtype = np.float64

# Bases
x_coord = d3.Coordinate('x')
dist = d3.Distributor(x_coord, dtype=dtype)
xbasis = d3.Chebyshev(x_coord, size=Nx, bounds=(-L, L), dealias=dealias)


# Fields

u = dist.VectorField(x_coord,name='u', bases=xbasis)
tau1 = dist.VectorField(x_coord,name = 'tau1')

# tau polynomials
tau_basis = xbasis.derivative_basis(1)
p1 = dist.Field(bases = tau_basis)

p1['c'][-1] = 2

# NCC
x = dist.local_grid(xbasis)
c = dist.Field(bases=xbasis)
c['g'] = 1/(x**2+1)

# Substitutions
dx = lambda A: d3.Differentiate(A, x_coord)

# Problem
problem = d3.IVP([u,tau1], namespace=locals())
problem.add_equation("dt(u)+ tau1*p1 =-c*dx(u)")

# Add BC
problem.add_equation("u(x = -L) = 1/((-L+3)**2+1)")

# Initial conditions
u['g'] = 1/((x+3)**2+1)

# Solver
solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_sim_time

#storage
u_list = [np.copy(u['g'])]
t_list = [solver.sim_time]

# Analysis
snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.1, max_writes=800000000)
snapshots.add_task(u, layout = 'g', name='u')

# CFL

CFL = d3.CFL(solver, initial_dt=max_timestep, cadence=10, safety=2, threshold=0.1, max_change=1.5, min_change=0.5, max_dt=max_timestep)

CFL.add_velocity(u)

# Solving/iterating
while solver.proceed:
    time_step = CFL.compute_timestep()
    solver.step(time_step)
    u.change_scales(1)
    u_list.append(np.copy(u['g']))
    t_list.append(solver.sim_time)
    if solver.iteration % 100 == 0:
        logger.info('Iteration = %i, time = %e, dt = %e' %(solver.iteration, solver.sim_time, time_step))

***********************************************************************************************************************

Thanks,

Mohammed