mulimoen/SummationByParts
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

move advance out of Fields

Browse Source
This commit is contained in:
Magnus Ulimoen
2019-12-14 13:59:20 +01:00
parent 946caf5e1b
commit 066c7619f7
2 changed files with 106 additions and 110 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/lib.rs
View File

@@ -61,7 +61,7 @@ impl Universe {
/// Using artifical dissipation with the upwind operator
pub fn advance_upwind(&mut self, dt: f32) {
Field::advance_upwind::<operators::Upwind4>(
maxwell::advance_upwind::<operators::Upwind4>(
&self.sys.0,
&mut self.sys.1,
dt,
@@ -72,7 +72,7 @@ impl Universe {
}
pub fn advance(&mut self, dt: f32) {
Field::advance::<operators::Upwind4>(
maxwell::advance::<operators::Upwind4>(
&self.sys.0,
&mut self.sys.1,
dt,

212
src/maxwell.rs
View File

@@ -68,124 +68,120 @@ impl Field {
ey.into_shape((ny, nx)).unwrap(),
)
}
}
pub(crate) fn advance_upwind<UO>(
&self,
fut: &mut Self,
dt: f32,
grid: &Grid<UO>,
work_buffers: Option<&mut WorkBuffers>,
) where
UO: UpwindOperator,
{
assert_eq!(self.0.shape(), fut.0.shape());
pub(crate) fn advance_upwind<UO>(
prev: &Field,
fut: &mut Field,
dt: f32,
grid: &Grid<UO>,
work_buffers: Option<&mut WorkBuffers>,
) where
UO: UpwindOperator,
{
assert_eq!(prev.0.shape(), fut.0.shape());
let mut wb: WorkBuffers;
let (y, k, tmp) = if let Some(x) = work_buffers {
(&mut x.y, &mut x.buf, &mut x.tmp)
} else {
wb = WorkBuffers::new(self.nx(), self.ny());
(&mut wb.y, &mut wb.buf, &mut wb.tmp)
let mut wb: WorkBuffers;
let (y, k, tmp) = if let Some(x) = work_buffers {
(&mut x.y, &mut x.buf, &mut x.tmp)
} else {
wb = WorkBuffers::new(prev.nx(), prev.ny());
(&mut wb.y, &mut wb.buf, &mut wb.tmp)
};
let boundaries = BoundaryTerms {
north: Boundary::This,
south: Boundary::This,
west: Boundary::This,
east: Boundary::This,
};
for i in 0..4 {
// y = y0 + c*kn
y.assign(&prev);
match i {
0 => {}
1 | 2 => {
y.scaled_add(1.0 / 2.0 * dt, &k[i - 1]);
}
3 => {
y.scaled_add(dt, &k[i - 1]);
}
_ => {
unreachable!();
}
};
let boundaries = BoundaryTerms {
north: Boundary::This,
south: Boundary::This,
west: Boundary::This,
east: Boundary::This,
};
for i in 0..4 {
// y = y0 + c*kn
y.assign(&self);
match i {
0 => {}
1 | 2 => {
y.scaled_add(1.0 / 2.0 * dt, &k[i - 1]);
}
3 => {
y.scaled_add(dt, &k[i - 1]);
}
_ => {
unreachable!();
}
};
RHS_upwind(&mut k[i], &y, grid, &boundaries, tmp);
}
Zip::from(&mut fut.0)
.and(&self.0)
.and(&*k[0])
.and(&*k[1])
.and(&*k[2])
.and(&*k[3])
.apply(|y1, &y0, &k1, &k2, &k3, &k4| {
*y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
});
RHS_upwind(&mut k[i], &y, grid, &boundaries, tmp);
}
/// Solving (Au)_x + (Bu)_y
/// with:
/// A B
/// [ 0, 0, 0] [ 0, 1, 0]
/// [ 0, 0, -1] [ 1, 0, 0]
/// [ 0, -1, 0] [ 0, 0, 0]
pub(crate) fn advance<SBP>(
&self,
fut: &mut Self,
dt: f32,
grid: &Grid<SBP>,
work_buffers: Option<&mut WorkBuffers>,
) where
SBP: SbpOperator,
{
assert_eq!(self.0.shape(), fut.0.shape());
Zip::from(&mut fut.0)
.and(&prev.0)
.and(&*k[0])
.and(&*k[1])
.and(&*k[2])
.and(&*k[3])
.apply(|y1, &y0, &k1, &k2, &k3, &k4| *y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4));
}
let mut wb: WorkBuffers;
let (y, k, tmp) = if let Some(x) = work_buffers {
(&mut x.y, &mut x.buf, &mut x.tmp)
} else {
wb = WorkBuffers::new(self.nx(), self.ny());
(&mut wb.y, &mut wb.buf, &mut wb.tmp)
/// Solving (Au)_x + (Bu)_y
/// with:
/// A B
/// [ 0, 0, 0] [ 0, 1, 0]
/// [ 0, 0, -1] [ 1, 0, 0]
/// [ 0, -1, 0] [ 0, 0, 0]
pub(crate) fn advance<SBP>(
prev: &Field,
fut: &mut Field,
dt: f32,
grid: &Grid<SBP>,
work_buffers: Option<&mut WorkBuffers>,
) where
SBP: SbpOperator,
{
assert_eq!(prev.0.shape(), fut.0.shape());
let mut wb: WorkBuffers;
let (y, k, tmp) = if let Some(x) = work_buffers {
(&mut x.y, &mut x.buf, &mut x.tmp)
} else {
wb = WorkBuffers::new(prev.nx(), prev.ny());
(&mut wb.y, &mut wb.buf, &mut wb.tmp)
};
let boundaries = BoundaryTerms {
north: Boundary::This,
south: Boundary::This,
west: Boundary::This,
east: Boundary::This,
};
for i in 0..4 {
// y = y0 + c*kn
y.assign(&prev);
match i {
0 => {}
1 | 2 => {
y.scaled_add(1.0 / 2.0 * dt, &k[i - 1]);
}
3 => {
y.scaled_add(dt, &k[i - 1]);
}
_ => {
unreachable!();
}
};
let boundaries = BoundaryTerms {
north: Boundary::This,
south: Boundary::This,
west: Boundary::This,
east: Boundary::This,
};
for i in 0..4 {
// y = y0 + c*kn
y.assign(&self);
match i {
0 => {}
1 | 2 => {
y.scaled_add(1.0 / 2.0 * dt, &k[i - 1]);
}
3 => {
y.scaled_add(dt, &k[i - 1]);
}
_ => {
unreachable!();
}
};
RHS(&mut k[i], &y, grid, &boundaries, tmp);
}
Zip::from(&mut fut.0)
.and(&self.0)
.and(&*k[0])
.and(&*k[1])
.and(&*k[2])
.and(&*k[3])
.apply(|y1, &y0, &k1, &k2, &k3, &k4| {
*y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
});
RHS(&mut k[i], &y, grid, &boundaries, tmp);
}
Zip::from(&mut fut.0)
.and(&prev.0)
.and(&*k[0])
.and(&*k[1])
.and(&*k[2])
.and(&*k[3])
.apply(|y1, &y0, &k1, &k2, &k3, &k4| *y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4));
}
#[allow(non_snake_case)]