GeneralFunction problem

[nimrod gavriel]

Hi all, 

I am trying to implement a GeneralFunction as a forcing for my model, but I have a problem when running in parallel. When the code receives my 'y_field', which is just a field with the y values, it allways gets the same lower bottom slice of the grid.

e.g., running on 2 cores, if my grid is -1 to 1 in x and y and the resolution is 128x128, 
when I print the y values inside my storms_fun(x, y, t), I always get y values from -1 to zero with a resolution 128x64. Why don't I ever get the other half? 

This also happened when I tried to get the y values inside the function with dist.local_grid.

Thanks!

My code is:

days_in_simulation = 15
# Resulution should be 2^n for best performance
# : 32, 64, 128, 256, 512, 1024
res = 256
Nx, Ny = res, res
dealias = 3 / 2
timestepper = d3.RKSMR
initial_timestep = 0.001
dtype = np.float64
Logger_cadence, results_cadence, max_dt_days = (
150,
0.03,
0.005,
) # result cadence in days
nrows, ncols = 1, 1 # For plot
CFL_cadence = 40
coords = d3.CartesianCoordinates("x", "y")
dist = d3.Distributor(coords, dtype=dtype)
xbasis = d3.RealFourier(
coords["x"], size=Nx, bounds=(-L_domain, L_domain), dealias=dealias
)
ybasis = d3.RealFourier(
coords["y"], size=Ny, bounds=(-L_domain, L_domain), dealias=dealias
)
# Fields
xi = dist.Field(name="xi", bases=(xbasis, ybasis))
psi = dist.Field(name="psi", bases=(xbasis, ybasis))
tau_psi = dist.Field(name="tau_psi")
u = dist.VectorField(coords, name="u", bases=(xbasis, ybasis))
f = dist.Field(name="f", bases=(xbasis, ybasis))
x_field = dist.Field(name="x_field", bases=(xbasis, ybasis))
y_field = dist.Field(name="y_field", bases=(xbasis, ybasis))
# Substitutions
x, y = dist.local_grids(xbasis, ybasis)
# domain = d3.Domain(dist, (xbasis, ybasis))
domain = Domain(dist, (xbasis, ybasis))
t = dist.Field()
ex, ey = coords.unit_vector_fields(dist)

x_field["g"] = x
y_field["g"] = y

# Define forcing
def storm_function(*args):
result = storms_fun(
x=args[0]["g"], y=args[1]["g"], t=args[2]["g"].flatten()
)
return result


def S(*args, domain=domain, F=storm_function):
return d3.GeneralFunction(
dist=dist,
domain=domain,
tensorsig=(),
dtype=np.float64,
layout="g",
func=F,
args=args,
)


storms = S(x_field, y_field, t)

u = d3.skew(d3.grad(psi))
xi = -d3.div(d3.skew(u))

# Problem
problem = d3.IVP([psi, tau_psi], time=t, namespace=locals())
problem.add_equation(
"dt(lap(psi)- Ld_fact*psi) + tau_psi - Ek*lap(lap(psi)) = -u@grad(f/Ro + lap(psi)- Ld_fact*psi) + S(x_field, y_field, t)"
)
problem.add_equation("integ(psi) = 0")

# Solver
solver = problem.build_solver(timestepper)
solver.stop_sim_time = stop_sim_time
file_handler_mode = "overwrite"

# Analysis
snapshots = solver.evaluator.add_file_handler(
"snapshots", sim_dt=results_cadence, max_writes=48
)
snapshots.add_task(xi / t_scale, name="vorticity", layout="g")
snapshots.add_task(storms, name="storms", layout="g")

# CFL
CFL = d3.CFL(
solver,
initial_dt=initial_timestep,
cadence=CFL_cadence,
safety=3,
max_change=1.2,
min_change=0.89,
max_dt=max_dt,
threshold=0.1,
)
CFL.add_velocity(u)
# Flow properties
flow = d3.GlobalFlowProperty(solver, cadence=Logger_cadence)
flow.add_property(((u @ ey * U_scale) ** 2), name="v")
flow.add_property(u @ u, name="KE")

# Main loop
try:
logger.info("Starting main loop")
while solver.proceed:
timestep = CFL.compute_timestep()
# storms["g"] = sim_case.storms_fun(x_dealias, y_dealias, solver.sim_time)
solver.step(timestep)
if (solver.iteration - 1) % Logger_cadence == 0:
max_w = np.sqrt(flow.max("v"))
tot_KE = flow.volume_integral("KE")
logger.info(
"Iteration=%i, Time=%e, dt=%e, max(v)=%f, Tot_KE =%f"
% (
solver.iteration,
solver.sim_time * t_scale_days,
timestep * t_scale_days,
max_w,
tot_KE,
)
)
if np.isnan(max_w):
break

except:
logger.error("Exception raised, triggering end of main loop.")
raise
finally:
solver.log_stats()

[Daniel Lecoanet]

Hi Nimrod,

I tried running your code in parallel (taking L_domain=1). If I print the values of y_field:

print(np.min(y_field['g']))
print(np.max(y_field['g']))

and run on 2 cores, I get:

-1.0

-0.0078125

0.0

0.9921875

So everything seems to be working correctly to me. By the way, I’m not sure why you are making the x_field and y_field variables, I think you should be able to use the x and y numpy arrays directly in your function.

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/709689c6-928f-4844-994d-06fa5652ac8fn%40googlegroups.com.

[nimrod gavriel]

Thanks, Daniel!

I tried again with a simpler version (attached), and I now understand my problem better.

I do indeed get the whole domain; the problem is that the time ("t") variable only passes on the first core (related to the lower half in the two cores test). 

See the attached log file I got as a result of this run.
There, a value for "t" only appears where the ymin, ymax = -1, 0.
As I need "t" for evaluating my source term, I only got it for the lower half in the previous attempts. 

Do you have an idea of how to solve this issue? 

Thanks again! 

Nimrod

[Daniel Lecoanet]

Hi Nimrod,

Sorry for the slow response. The easiest way to do this is to use the allgather_data() method, e.g., time = args[2].allgather_data(). See attached code showing this.

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/70cdcc3b-7284-4d4c-9560-092a05d6970cn%40googlegroups.com.
<try_general_function.py><model_log.txt>

[nimrod gavriel]

Thanks, Danial!

I tried this but for some reason the parallelization makes a mass and I get certain cores activated in the wrong times and in a mirrored fashion when I run on multiple cores.

My solution ended up being to add the source in the "while solver.proceed" loop with the load_from_global_grid_data() function, although it probably has a bigger computational cost in parallel.

Thanks, 

Nimrod