Or: qb.dat.data[:] = np.nan_to_num(qb.dat.data_ro) On Fri, 17 Jun 2016 at 08:47 Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
I think the more direct approach is to edit the numpy array qb.dat.data.
E.g. something like qb.dat.data[np.isnan(qb.dat.data)] = 0
On 17 June 2016 at 03:52, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
I want to define an initial condition for a field that I am pretty sure will have nan's. It does in python so I'm pretty sure that it will in firedrake as well.
In python I know how to remove that nan's by doing the following.
qb[np.isnan(qb)] = 0.
Is there a way of doing this in firedrake?
Cheers, Francis
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
Hi FD's Trying to get my feet wet. Simple query. Trying to make a DG conservation law/system code, starting with linear advection, backward Euler in steps. How do I make the code snippet below forward Euler or Crank-Nicolson? I used all the old time step terms in L but that does not work. (Now commented out -thet factor terms) Also, why does the nonlinear solver not work? I think it is because the printing can't be done in DG0 so there needs to assign between the DG solution of the prime unknown. Would be nice to have a standard DG example for U_t + div F = 0 forward Euler or RK3. Thanks for the reply as then I can keep playing in the plane! Onno P.S. Code below tend = 4.211 timestep = tend/40 V = FunctionSpace(tmesh, "DG", 0) W = VectorFunctionSpace(tmesh, "CG", 1) # Advection velocity velocity = Expression(("1.0", "0.0")) u0 = project(velocity, W) #Inflow BC inflow = Expression("(x[0] < 0.0) ? 1.0 : 0") inflow2 = Expression("(x[0] < 0.5) ? 1.0 : 0") inflow3 = Expression("(x[0] < 1.0) ? 1.0 : 0") #Inflow Variable D0 = Function(V) D0.interpolate(inflow) m1 = Function(V) m1.interpolate(inflow2) m2 = Function(V) m2.interpolate(inflow3) #Initial Condition ic = Expression("(x[0] < .5) ? 1.0 : 0") ic2 = Expression("(x[0] < .5) ? 1.0 : 0") ic3 = Expression("(x[0] < .5) ? 1.0 : 0") #Set up element normal n = FacetNormal(tmesh) # un = 0.5*(dot(u0, n) + abs(dot(u0, n))) # Previous timestep for BE timestep D_ = Function(V) m1_ = Function(V) m2_ = Function(V) D_.interpolate(ic) m1_.interpolate(ic2) m2_.interpolate(ic3) D = TrialFunction(V) phi = TestFunction(V) # # Setting up Bilinear form # thet = 0; a1 = ((D/timestep)*phi - (1-thet)*D*dot(u0, grad(phi)))*dx # (1-thet) interior term at new time step a2 = (1-thet)*dot(jump(phi), un('+')*D('+') - un('-')*D('-'))*dS # Internal flux a3 = (1-thet)*dot(phi, un*D)*(ds(112)) # outflow at top wall a = a1 + a2 + a3 # L1 = (D_/timestep)*phi*dx - D0*phi*dot(u0, n)*(ds(114)) # previous timestep & given inflow at bottom wall # L2 = thet*(D_*dot(u0, grad(phi)))*dx # thet times interior term at old time stepL # L3 = -thet*dot(jump(phi), un('+')*D_('+') - un('-')*D_('-'))*dS # Internal flux # L4 = -thet*dot(phi, un*D_)*(ds(112)) # outflow at top wall L = L1 F = a-L1 out = Function(V) outfile = File("./Results/swe2d.pvd") # Write IC to file # outfile << D_ D_.assign(out) outfile.write(out) # outfile.wrote(project(out, W, name="D")) t = 0.0 while (t <= tend): #Solve problem solve(a == L, out) #solve(F == 0, out) #Update previous timestep D_.assign(out) #Update time t+=timestep # outfile << out outfile.write(out) # outfile.write(project(out, W, name="D"))
Hello Onno, I have a 1D, one-layer SW model without rotation in Fenics that uses Crank-Nicholson. As Colin pointed out already, the QG code uses RK3, but if my Fenics code might help at all please let me know an I will send it along. Safe travels, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Onno Bokhove [O.Bokhove@leeds.ac.uk] Sent: Friday, June 17, 2016 7:56 AM To: firedrake Subject: [firedrake] feet wet Hi FD's Trying to get my feet wet. Simple query. Trying to make a DG conservation law/system code, starting with linear advection, backward Euler in steps. How do I make the code snippet below forward Euler or Crank-Nicolson? I used all the old time step terms in L but that does not work. (Now commented out -thet factor terms) Also, why does the nonlinear solver not work? I think it is because the printing can't be done in DG0 so there needs to assign between the DG solution of the prime unknown. Would be nice to have a standard DG example for U_t + div F = 0 forward Euler or RK3. Thanks for the reply as then I can keep playing in the plane! Onno P.S. Code below tend = 4.211 timestep = tend/40 V = FunctionSpace(tmesh, "DG", 0) W = VectorFunctionSpace(tmesh, "CG", 1) # Advection velocity velocity = Expression(("1.0", "0.0")) u0 = project(velocity, W) #Inflow BC inflow = Expression("(x[0] < 0.0) ? 1.0 : 0") inflow2 = Expression("(x[0] < 0.5) ? 1.0 : 0") inflow3 = Expression("(x[0] < 1.0) ? 1.0 : 0") #Inflow Variable D0 = Function(V) D0.interpolate(inflow) m1 = Function(V) m1.interpolate(inflow2) m2 = Function(V) m2.interpolate(inflow3) #Initial Condition ic = Expression("(x[0] < .5) ? 1.0 : 0") ic2 = Expression("(x[0] < .5) ? 1.0 : 0") ic3 = Expression("(x[0] < .5) ? 1.0 : 0") #Set up element normal n = FacetNormal(tmesh) # un = 0.5*(dot(u0, n) + abs(dot(u0, n))) # Previous timestep for BE timestep D_ = Function(V) m1_ = Function(V) m2_ = Function(V) D_.interpolate(ic) m1_.interpolate(ic2) m2_.interpolate(ic3) D = TrialFunction(V) phi = TestFunction(V) # # Setting up Bilinear form # thet = 0; a1 = ((D/timestep)*phi - (1-thet)*D*dot(u0, grad(phi)))*dx # (1-thet) interior term at new time step a2 = (1-thet)*dot(jump(phi), un('+')*D('+') - un('-')*D('-'))*dS # Internal flux a3 = (1-thet)*dot(phi, un*D)*(ds(112)) # outflow at top wall a = a1 + a2 + a3 # L1 = (D_/timestep)*phi*dx - D0*phi*dot(u0, n)*(ds(114)) # previous timestep & given inflow at bottom wall # L2 = thet*(D_*dot(u0, grad(phi)))*dx # thet times interior term at old time stepL # L3 = -thet*dot(jump(phi), un('+')*D_('+') - un('-')*D_('-'))*dS # Internal flux # L4 = -thet*dot(phi, un*D_)*(ds(112)) # outflow at top wall L = L1 F = a-L1 out = Function(V) outfile = File("./Results/swe2d.pvd") # Write IC to file # outfile << D_ D_.assign(out) outfile.write(out) # outfile.wrote(project(out, W, name="D")) t = 0.0 while (t <= tend): #Solve problem solve(a == L, out) #solve(F == 0, out) #Update previous timestep D_.assign(out) #Update time t+=timestep # outfile << out outfile.write(out) # outfile.write(project(out, W, name="D"))
participants (3)
-
David Ham
-
Francis Poulin
-
Onno Bokhove