Permalink
Browse files

Merge pull request #1002 from lahnerml/bounce_back_fix

potential bugfix for bounce back in CPU lb
  • Loading branch information...
2 parents 8d50f0d + b724456 commit b36e54dec8c990bd4142825538ebb8656b916cf7 @rempferg rempferg committed on GitHub Jan 31, 2017
@@ -558,12 +558,13 @@ void lb_bounce_back() {
population_shift = 0;
for (l = 0; l < 3; l++) {
population_shift -=
- lbpar.agrid * lbpar.agrid * lbpar.agrid * lbpar.agrid *
- lbpar.agrid * lbpar.rho[0] * 2 * lbmodel.c[i][l] *
+ lbpar.agrid * lbpar.agrid * lbpar.agrid *
+ lbpar.rho[0] * 2 * lbmodel.c[i][l] *
lbmodel.w[i] *
lb_boundaries[lbfields[k].boundary - 1].velocity[l] /
lbmodel.c_sound_sq;
}
+
if (x - lbmodel.c[i][0] > 0 &&
x - lbmodel.c[i][0] < lblattice.grid[0] + 1 &&
y - lbmodel.c[i][1] > 0 &&
View
@@ -1323,13 +1323,13 @@ int lb_lbnode_get_pi_neq(int* ind, double* p_pi) {
index = get_linear_index(ind_shifted[0],ind_shifted[1],ind_shifted[2],lblattice.halo_grid);
mpi_recv_fluid(node,index,&rho,j,pi);
- // unit conversion // TODO: Check Unit Conversion!
- p_pi[0] = pi[0]/lbpar.tau/lbpar.tau/lbpar.agrid/lbpar.agrid/lbpar.agrid;
- p_pi[1] = pi[1]/lbpar.tau/lbpar.tau/lbpar.agrid/lbpar.agrid/lbpar.agrid;
- p_pi[2] = pi[2]/lbpar.tau/lbpar.tau/lbpar.agrid/lbpar.agrid/lbpar.agrid;
- p_pi[3] = pi[3]/lbpar.tau/lbpar.tau/lbpar.agrid/lbpar.agrid/lbpar.agrid;
- p_pi[4] = pi[4]/lbpar.tau/lbpar.tau/lbpar.agrid/lbpar.agrid/lbpar.agrid;
- p_pi[5] = pi[5]/lbpar.tau/lbpar.tau/lbpar.agrid/lbpar.agrid/lbpar.agrid;
+ // unit conversion
+ p_pi[0] = pi[0]/lbpar.tau/lbpar.tau/lbpar.agrid;
+ p_pi[1] = pi[1]/lbpar.tau/lbpar.tau/lbpar.agrid;
+ p_pi[2] = pi[2]/lbpar.tau/lbpar.tau/lbpar.agrid;
+ p_pi[3] = pi[3]/lbpar.tau/lbpar.tau/lbpar.agrid;
+ p_pi[4] = pi[4]/lbpar.tau/lbpar.tau/lbpar.agrid;
+ p_pi[5] = pi[5]/lbpar.tau/lbpar.tau/lbpar.agrid;
#endif // LB
}
return 0;
@@ -1,17 +1,17 @@
# Copyright (C) 2011,2012,2013,2014,2015,2016 The ESPResSo project
-#
+#
# This file is part of ESPResSo.
-#
+#
# ESPResSo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
-#
+#
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
-#
+#
source "tests_common.tcl"
@@ -33,16 +33,16 @@ setmd box_l $l $l $l
setmd time_step 0.01
thermostat lb 0
-set agrid 1.00
+set agrid 0.5
set visc 7.
set rho 2.
set tau 0.04
setmd skin [ expr 0.4*$agrid ]
-lbfluid cpu agrid $agrid visc $visc dens $rho friction 1 tau $tau
+lbfluid cpu agrid $agrid visc $visc dens $rho friction 1 tau $tau
-# Wall positions are set to $agrid and $l-$agrid, leaving one layer of boundary nodes
+# Wall positions are set to $agrid and $l-$agrid, leaving one layer of boundary nodes
# on each side of the box
#
# We start with couette flow setting the velocity on boundary 2 to $v_couette
@@ -61,7 +61,7 @@ lset v_boundary $couette_flow_direction $v_couette
lbboundary wall normal [ lindex $normal1 0 ] [ lindex $normal1 1 ] [ lindex $normal1 2 ] dist $agrid
lbboundary wall normal [ lindex $normal2 0 ] [ lindex $normal2 1 ] [ lindex $normal2 2 ] dist [ expr -$l +$agrid ] \
- velocity [ lindex $v_boundary 0 ] [ lindex $v_boundary 1 ] [ lindex $v_boundary 2 ]
+ velocity [ lindex $v_boundary 0 ] [ lindex $v_boundary 1 ] [ lindex $v_boundary 2 ]
set dist [ expr $l - 2*$agrid ]
@@ -76,7 +76,7 @@ for { set i 2 } { $i < int(floor($l/$agrid))-2 } { incr i } {
set u_theory [ expr $v_couette*($pos-$agrid)/$dist ]
set p_xy_theory [ expr -$v_couette/$dist*$rho*$visc ]
- ## We go though the box in normal direction
+ ## We go through the box in normal direction
set m [ expr $i * [ lindex $normal1 0 ] ]
set n [ expr $i * [ lindex $normal1 1 ] ]
set o [ expr $i * [ lindex $normal1 2 ] ]
@@ -102,8 +102,8 @@ set couette_u_accuracy [ expr $accuracy_u / $meanabs_u ]
set couette_p_accuracy [ expr $accuracy_p / $meanabs_p ]
puts "Couette flow result:"
-puts "flow accuary $couette_u_accuracy"
-puts "pressure accuary $couette_p_accuracy"
+puts "flow accuracy $couette_u_accuracy"
+puts "pressure accuracy $couette_p_accuracy"
puts "----------"
# Now we add a force density in normal direction, and compress the flow.
@@ -151,13 +151,13 @@ for { set i 2 } { $i < int(floor($l/$agrid))-2 } { incr i } {
}
set hydrostatic_p_accuracy [ expr $accuracy_p / $meanabs_p ]
puts "Hydrostatic test result:"
-puts "pressure accuary $hydrostatic_p_accuracy"
+puts "pressure accuracy $hydrostatic_p_accuracy"
puts "-------------"
# Now we add a force density in the direction of the Couette flow
# We'd like to switch of the Couette flow, but changing the velocity of the BCs is not implemented
# yet. So we just make a much larger force density.
-#
+#
# We expect a flow profile u=f/eta/2*(x-agrid)*($l-$agrid-$x)
# and the pressure tensor has to be the derivative of that: p_xy = f*(x-$l/2)
@@ -204,8 +204,8 @@ set poisseuille_u_accuracy [ expr $accuracy_u / $meanabs_u ]
set poisseuille_p_accuracy [ expr $accuracy_p / $meanabs_p ]
puts "Poisseuille flow result:"
-puts "flow accuary $poisseuille_u_accuracy"
-puts "pressure accuary $poisseuille_p_accuracy"
+puts "flow accuracy $poisseuille_u_accuracy"
+puts "pressure accuracy $poisseuille_p_accuracy"
puts "----------"
if { $couette_u_accuracy > 1e-5 } {
@@ -214,13 +214,13 @@ if { $couette_u_accuracy > 1e-5 } {
if { $couette_p_accuracy > 1e-5 } {
error_exit "Couette pressure accuracy not achieved"
}
-if { $hydrostatic_p_accuracy > 1e-3 } {
+if { $hydrostatic_p_accuracy > 5e-3 } {
error_exit "hydrostatic pressure accuracy not achieved"
}
if { $poisseuille_u_accuracy > 3e-2 } {
error_exit "Poisseuille flow accuracy not achieved"
}
-if { $poisseuille_p_accuracy > 1e-2 } {
+if { $poisseuille_p_accuracy > 1.5e-2 } {
error_exit "Poisseuille pressure accuracy not achieved"
}
@@ -1,17 +1,17 @@
# Copyright (C) 2011,2012,2013,2014,2015,2016 The ESPResSo project
-#
+#
# This file is part of ESPResSo.
-#
+#
# ESPResSo is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
-#
+#
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
-#
+#
source "tests_common.tcl"
@@ -32,7 +32,7 @@ set l 12
setmd box_l $l $l $l
setmd time_step 0.01
-set agrid 1.0
+set agrid 0.5
set visc 7.
set rho 2.
set tau 0.04
@@ -47,7 +47,7 @@ for {set i 0} {$i < 10} {incr i} {
part $i pos $i 0 0
}
-# Wall positions are set to $agrid and $l-$agrid, leaving one layer of boundary nodes
+# Wall positions are set to $agrid and $l-$agrid, leaving one layer of boundary nodes
# on each side of the box
#
# We start with couette flow setting the velocity on boundary 2 to $v_couette
@@ -66,7 +66,7 @@ lset v_boundary $couette_flow_direction $v_couette
lbboundary wall normal [ lindex $normal1 0 ] [ lindex $normal1 1 ] [ lindex $normal1 2 ] dist $agrid
lbboundary wall normal [ lindex $normal2 0 ] [ lindex $normal2 1 ] [ lindex $normal2 2 ] dist [ expr -$l +$agrid ] \
- velocity [ lindex $v_boundary 0 ] [ lindex $v_boundary 1 ] [ lindex $v_boundary 2 ]
+ velocity [ lindex $v_boundary 0 ] [ lindex $v_boundary 1 ] [ lindex $v_boundary 2 ]
set dist [ expr $l - 2*$agrid ]
@@ -82,7 +82,7 @@ for { set i 2 } { $i < int(floor($l/$agrid))-2 } { incr i } {
set u_theory [ expr $v_couette*($pos-$agrid)/$dist ]
set p_xy_theory [ expr -$v_couette/$dist*$rho*$visc ]
- ## We go though the box in normal direction
+ ## We go through the box in normal direction
set m [ expr $i * [ lindex $normal1 0 ] ]
set n [ expr $i * [ lindex $normal1 1 ] ]
set o [ expr $i * [ lindex $normal1 2 ] ]
@@ -121,7 +121,7 @@ set fx [ expr $f_hydrostatic*[ lindex $normal1 0 ] ]
set fy [ expr $f_hydrostatic*[ lindex $normal1 1 ] ]
set fz [ expr $f_hydrostatic*[ lindex $normal1 2 ] ]
-lbfluid agrid $agrid visc $visc dens $rho friction 1 tau $tau ext_force $fx $fy $fz
+lbfluid cpu agrid $agrid visc $visc dens $rho friction 1 tau $tau ext_force $fx $fy $fz
integrate 1000
# Equation of state: p = rho*c_2**2
@@ -163,7 +163,7 @@ puts "-------------"
# Now we add a force density in the direction of the Couette flow
# We'd like to switch of the Couette flow, but changing the velocity of the BCs is not implemented
# yet. So we just make a much larger force density.
-#
+#
# We expect a flow profile u=f/eta/2*(x-agrid)*($l-$agrid-$x)
# and the pressure tensor has to be the derivative of that: p_xy = f*(x-$l/2)
@@ -172,7 +172,7 @@ set f_body 0.1
set f_body_vec [ list 0 0 0 ]
lset f_body_vec $couette_flow_direction $f_body
-lbfluid agrid $agrid visc $visc dens $rho friction 1 tau $tau ext_force \
+lbfluid cpu agrid $agrid visc $visc dens $rho friction 1 tau $tau ext_force \
[ lindex $f_body_vec 0 ] [ lindex $f_body_vec 1 ] [ lindex $f_body_vec 2 ]
integrate 2000
@@ -220,13 +220,13 @@ if { $couette_u_accuracy > 1e-5 } {
if { $couette_p_accuracy > 1e-5 } {
error_exit "Couette pressure accuracy not achieved"
}
-if { $hydrostatic_p_accuracy > 1e-3 } {
+if { $hydrostatic_p_accuracy > 5e-3 } {
error_exit "hydrostatic pressure accuracy not achieved"
}
if { $poisseuille_u_accuracy > 3e-2 } {
error_exit "Poisseuille flow accuracy not achieved"
}
-if { $poisseuille_p_accuracy > 1e-2 } {
+if { $poisseuille_p_accuracy > 1.5e-2 } {
error_exit "Poisseuille pressure accuracy not achieved"
}
@@ -41,7 +41,7 @@ cellsystem domain_decomposition -no_verlet_list
set kinematic_visc 10.
set dens 1.
-lbfluid cpu visc [ expr $kinematic_visc / $dens ] dens $dens friction 1. agrid 1.0 tau .01
+lbfluid cpu visc [ expr $kinematic_visc / $dens ] dens $dens friction 1. agrid 0.5 tau .01
set v1 1.0
lbboundary wall normal -1. 0. 0. dist [ expr -(+0.5+$w) ] velocity 0.00 $v1 0.

0 comments on commit b36e54d

Please sign in to comment.