From 998c80873d667a91feb9fdfda04e5209318e857e Mon Sep 17 00:00:00 2001 From: Srinath Kailasa Date: Fri, 17 Nov 2023 17:50:30 +0000 Subject: [PATCH] Some kind of bug with empty boxes --- fmm/src/field_translation.rs | 4 +- fmm/src/fmm.rs | 128 ++++++++++++++++++----------------- 2 files changed, 67 insertions(+), 65 deletions(-) diff --git a/fmm/src/field_translation.rs b/fmm/src/field_translation.rs index 8e70d53e..f824d71f 100644 --- a/fmm/src/field_translation.rs +++ b/fmm/src/field_translation.rs @@ -206,7 +206,6 @@ where ); // Now compute the multipole expansions, with each of chunk_size boxes at a time. - // let chunk_size = 8; // This has to be a divisor of nleaves, 8 is the minimum possible answer for adaptive case let chunk_size = find_chunk_size(nleaves, P2M_MAX_CHUNK_SIZE); check_potentials @@ -254,8 +253,7 @@ where max_chunk_size = P2M_MAX_CHUNK_SIZE; } let chunk_size = find_chunk_size(nsources, max_chunk_size); - // println!("LEVEL {:?} CHUNK SIZE {:?}", level, chunk_size); - + multipoles .par_chunks_exact(nsiblings * ncoeffs*chunk_size) .zip(self.level_multipoles[(level - 1) as usize].par_chunks_exact(chunk_size)) diff --git a/fmm/src/fmm.rs b/fmm/src/fmm.rs index 56dc77cc..7328e19d 100644 --- a/fmm/src/fmm.rs +++ b/fmm/src/fmm.rs @@ -1294,6 +1294,8 @@ mod test { let alpha_outer = 2.95; let adaptive = false; let ncrit = 150; + + // TODO: There is a bug for when boxes are empty ... let depth = 5; let kernel = Laplace3dKernel::default(); @@ -1326,68 +1328,70 @@ mod test { println!("linear p2m {:?}", s.elapsed()); - let kernel = Laplace3dKernel::default(); - - let tree = SingleNodeTree::new( - points.data(), - adaptive, - Some(ncrit), - Some(depth), - &global_idxs[..], - ); - - let m2l_data_fft = - FftFieldTranslationKiFmm::new(kernel.clone(), order, *tree.get_domain(), alpha_inner); - let fmm = KiFmm::new(order, alpha_inner, alpha_outer, kernel, tree, m2l_data_fft); - - // Form charge dict, matching charges with their associated global indices - let charge_dict = build_charge_dict(&global_idxs[..], &charges[..]); - let old_datatree = FmmData::new(fmm, &charge_dict); - - let &idx = datatree.fmm.tree().key_to_index.get(&ROOT).unwrap(); - let old_leaf = old_datatree.fmm.tree().get_all_leaves().unwrap()[idx]; - let old_key = old_datatree.fmm.tree().get_all_keys().unwrap()[idx]; - // let old_points = old_datatree.points.get(&old_leaf).unwrap(); - // let old_points = old_points.iter().map(|p| p.coordinate).flat_map(|[x, y, z]| vec![x, y, z]).collect_vec(); - - let new_leaf = datatree.fmm.tree().get_all_leaves().unwrap()[idx]; - let new_key = datatree.fmm.tree().get_all_keys().unwrap()[idx]; - println!("old {:?} new {:?} keys", old_key, new_key); - - let (l, r) = datatree.charge_index_pointer[idx]; - // let new_points = &datatree.fmm.tree().get_all_coordinates().unwrap()[l*3..r*3]; - - let s = Instant::now(); - old_datatree.p2m(); - for level in (1..=depth).rev() { - old_datatree.m2m(level) - } - println!("old p2m {:?}", s.elapsed()); - - // Check potentials - let midx = datatree.fmm.tree().key_to_index.get(&new_key).unwrap(); - // let (l, r) = datatree.expansion_index_pointer[*midx]; - let ncoeffs = datatree.fmm.m2l.ncoeffs(datatree.fmm.order); - let new_multipole = &datatree.multipoles[midx * ncoeffs..(midx + 1) * ncoeffs]; - let old_multipole = old_datatree - .multipoles - .get(&old_key) - .unwrap() - .deref() - .lock() - .unwrap(); - - // println!("HERE {:?} {:?}", old_key, old_multipole.data()); - // println!("HERE {:?} {:?}", new_key, new_multipole); - let abs_error: f64 = old_multipole - .data() - .iter() - .zip(new_multipole.iter()) - .map(|(a, b)| (a - b).abs()) - .sum(); - - let rel_error = abs_error / (old_multipole.data().iter().sum::()); - println!("rel error {:?}", rel_error); + // let kernel = Laplace3dKernel::default(); + + // let tree = SingleNodeTree::new( + // points.data(), + // adaptive, + // Some(ncrit), + // Some(depth), + // &global_idxs[..], + // ); + + // let m2l_data_fft = + // FftFieldTranslationKiFmm::new(kernel.clone(), order, *tree.get_domain(), alpha_inner); + // let fmm = KiFmm::new(order, alpha_inner, alpha_outer, kernel, tree, m2l_data_fft); + + // // Form charge dict, matching charges with their associated global indices + // let charge_dict = build_charge_dict(&global_idxs[..], &charges[..]); + // let s = Instant::now(); + // let old_datatree = FmmData::new(fmm, &charge_dict); + // println!("old data tree setup {:?}", s.elapsed()); + + // let &idx = datatree.fmm.tree().key_to_index.get(&ROOT).unwrap(); + // let old_leaf = old_datatree.fmm.tree().get_all_leaves().unwrap()[idx]; + // let old_key = old_datatree.fmm.tree().get_all_keys().unwrap()[idx]; + // // let old_points = old_datatree.points.get(&old_leaf).unwrap(); + // // let old_points = old_points.iter().map(|p| p.coordinate).flat_map(|[x, y, z]| vec![x, y, z]).collect_vec(); + + // let new_leaf = datatree.fmm.tree().get_all_leaves().unwrap()[idx]; + // let new_key = datatree.fmm.tree().get_all_keys().unwrap()[idx]; + // println!("old {:?} new {:?} keys", old_key, new_key); + + // let (l, r) = datatree.charge_index_pointer[idx]; + // // let new_points = &datatree.fmm.tree().get_all_coordinates().unwrap()[l*3..r*3]; + + // let s = Instant::now(); + // old_datatree.p2m(); + // for level in (1..=depth).rev() { + // old_datatree.m2m(level) + // } + // println!("old p2m {:?}", s.elapsed()); + + // // Check potentials + // let midx = datatree.fmm.tree().key_to_index.get(&new_key).unwrap(); + // // let (l, r) = datatree.expansion_index_pointer[*midx]; + // let ncoeffs = datatree.fmm.m2l.ncoeffs(datatree.fmm.order); + // let new_multipole = &datatree.multipoles[midx * ncoeffs..(midx + 1) * ncoeffs]; + // let old_multipole = old_datatree + // .multipoles + // .get(&old_key) + // .unwrap() + // .deref() + // .lock() + // .unwrap(); + + // // println!("HERE {:?} {:?}", old_key, old_multipole.data()); + // // println!("HERE {:?} {:?}", new_key, new_multipole); + // let abs_error: f64 = old_multipole + // .data() + // .iter() + // .zip(new_multipole.iter()) + // .map(|(a, b)| (a - b).abs()) + // .sum(); + + // let rel_error = abs_error / (old_multipole.data().iter().sum::()); + // println!("rel error {:?}", rel_error); assert!(false) }