Commit 8c946751 authored by Daniel Brown's avatar Daniel Brown
Browse files

Fixing rom q settings to deal with 3 map separations

parent c3b7e949
......@@ -51,7 +51,7 @@ maxintop 400000 # maximum number of DCUHRE integration calls
mapintmethod 3 # integration method for mirror maps: 1 Rieman,
# 2 Cuba serial, 3 Cuba parallel
cuba_numprocs 1 # maximum number for parallel Cuba processes
maxintcuba 1000000 # maximum number of Cuba integration calls
maxintcuba 10000000 # maximum number of Cuba integration calls
mapinterpmethod 2 # interploation: 1 neareast, 2 linear, 2 cubic
mapinterpsize 3 # size of interpolation array
abserr 1E-3 # absolute error for Cuba intergations
......
......@@ -994,65 +994,66 @@ int set_k_mirror(int mirror_index) {
mirror->knm_order[2] = DEFAULT_KNM_ORDER_THIRD;
}
complex_t qLx_11 = complex_0, qLx_12 = complex_0, qLx_21 = complex_0, qLx_22 = complex_0;
complex_t qLy_11 = complex_0, qLy_12 = complex_0, qLy_21 = complex_0, qLy_22 = complex_0;
// compute what q_L is, the beam parameter of the unit expansion inserted into the
// coupling coefficient inner product.
if(mirror->knm_change_q == 1){
qLx_11 = knm_q.qxt1_11;
qLy_11 = knm_q.qyt1_11;
qLx_12 = knm_q.qxt1_12;
qLy_12 = knm_q.qyt1_12;
qLx_21 = knm_q.qxt1_21;
qLy_21 = knm_q.qyt1_21;
qLx_22 = knm_q.qxt1_22;
qLy_22 = knm_q.qyt1_22;
}else if(mirror->knm_change_q == 2){
qLx_11 = knm_q.qxt2_11;
qLy_11 = knm_q.qyt2_11;
qLx_12 = knm_q.qxt2_12;
qLy_12 = knm_q.qyt2_12;
qLx_21 = knm_q.qxt2_21;
qLy_21 = knm_q.qyt2_21;
qLx_22 = knm_q.qxt2_22;
qLy_22 = knm_q.qyt2_22;
} else
bug_error("value of mirror->knm_change_q was %i, which is not handled.", mirror->knm_change_q);
// the inner product calculations uses the mirrors own knm_q structure
// for storing the incoming q(x/y)t1_(direction) and outgoing q(x/y)t2_(direction)
// beam parameters. These need to be set before the solver are called below
// here we need to set the initial q values used for the first coupling coefficient
// matrix this always goes from incoming q_L to outgoing q_2
if (CALC_MR_KNM(mirror,11)){
mirror->knm_q.qxt1_11 = qLx_11;
mirror->knm_q.qyt1_11 = qLy_11;
mirror->knm_q.qxt2_11 = knm_q.qxt2_11;
mirror->knm_q.qyt2_11 = knm_q.qyt2_11;
mirror->knm_q.qxt1_11 = knm_q.qxt1_11;
mirror->knm_q.qyt1_11 = knm_q.qxt1_11;
// if first solver should handle q change
if(mirror->knm_change_q == 1) {
mirror->knm_q.qxt2_11 = knm_q.qxt2_11;
mirror->knm_q.qyt2_11 = knm_q.qyt2_11;
} else {
mirror->knm_q.qxt2_11 = knm_q.qxt1_11;
mirror->knm_q.qyt2_11 = knm_q.qyt1_11;
}
}
if (CALC_MR_KNM(mirror,12)){
mirror->knm_q.qxt1_12 = qLx_12;
mirror->knm_q.qyt1_12 = qLy_12;
mirror->knm_q.qxt2_12 = knm_q.qxt2_12;
mirror->knm_q.qyt2_12 = knm_q.qyt2_12;
mirror->knm_q.qxt1_12 = knm_q.qxt1_12;
mirror->knm_q.qyt1_12 = knm_q.qxt1_12;
// if first solver should handle q change
if(mirror->knm_change_q == 1) {
mirror->knm_q.qxt2_12 = knm_q.qxt2_12;
mirror->knm_q.qyt2_12 = knm_q.qyt2_12;
} else {
mirror->knm_q.qxt2_12 = knm_q.qxt1_12;
mirror->knm_q.qyt2_12 = knm_q.qyt1_12;
}
}
if (CALC_MR_KNM(mirror,21)){
mirror->knm_q.qxt1_21 = qLx_21;
mirror->knm_q.qyt1_21 = qLy_21;
mirror->knm_q.qxt2_21 = knm_q.qxt2_21;
mirror->knm_q.qyt2_21 = knm_q.qyt2_21;
mirror->knm_q.qxt1_21 = knm_q.qxt1_21;
mirror->knm_q.qyt1_21 = knm_q.qxt1_21;
// if first solver should handle q change
if(mirror->knm_change_q == 1) {
mirror->knm_q.qxt2_21 = knm_q.qxt2_21;
mirror->knm_q.qyt2_21 = knm_q.qyt2_21;
} else {
mirror->knm_q.qxt2_21 = knm_q.qxt1_21;
mirror->knm_q.qyt2_21 = knm_q.qyt1_21;
}
}
if (CALC_MR_KNM(mirror,22)){
mirror->knm_q.qxt1_22 = qLx_22;
mirror->knm_q.qyt1_22 = qLy_22;
mirror->knm_q.qxt2_22 = knm_q.qxt2_22;
mirror->knm_q.qyt2_22 = knm_q.qyt2_22;
mirror->knm_q.qxt1_22 = knm_q.qxt1_22;
mirror->knm_q.qyt1_22 = knm_q.qxt1_22;
// if first solver should handle q change
if(mirror->knm_change_q == 1) {
mirror->knm_q.qxt2_22 = knm_q.qxt2_22;
mirror->knm_q.qyt2_22 = knm_q.qyt2_22;
} else {
mirror->knm_q.qxt2_22 = knm_q.qxt1_22;
mirror->knm_q.qyt2_22 = knm_q.qyt1_22;
}
}
// set tmp matrix to identity, this is the matrix we will apply each knm too
......@@ -1073,6 +1074,11 @@ int set_k_mirror(int mirror_index) {
case 2:
message(" %i Bayer-Helms\n", i+1);
break;
case 4:
message(" %i ROMHOM\n", i+1);
break;
default:
bug_error("not handled");
}
}
}
......@@ -1096,6 +1102,9 @@ int set_k_mirror(int mirror_index) {
case 4:
message(" - ROMHOM\n");
break;
default:
bug_error("not handled");
}
message(" 11: q_in = x:%s y:%s\n",complex_form(mirror->knm_q.qxt1_11),complex_form(mirror->knm_q.qyt1_11));
......@@ -1193,35 +1202,48 @@ int set_k_mirror(int mirror_index) {
bug_error("Could not handle knm_order digit %i",mirror->knm_order[i]);
}
if (i==0){
// Once we have computed the first matrix using q_L -> q_2 we now have to compute
// q'_1 -> q_L, seeing as we only have 2 separations currently
if (i == 0 && mirror->knm_change_q == 1){
// If the q changes at the start then we the first solver
// has handled any mode mismatch so we're not only using q2
if (CALC_MR_KNM(mirror,11)){
mirror->knm_q.qxt1_11 = knm_q.qxt1_11;
mirror->knm_q.qyt1_11 = knm_q.qyt1_11;
mirror->knm_q.qxt2_11 = qLx_11;
mirror->knm_q.qyt2_11 = qLy_11;
mirror->knm_q.qxt1_11 = knm_q.qxt2_11;
mirror->knm_q.qyt1_11 = knm_q.qyt2_11;
}
if (CALC_MR_KNM(mirror,12)){
mirror->knm_q.qxt1_12 = knm_q.qxt1_12;
mirror->knm_q.qyt1_12 = knm_q.qyt1_12;
mirror->knm_q.qxt2_12 = qLx_12;
mirror->knm_q.qyt2_12 = qLy_12;
mirror->knm_q.qxt1_12 = knm_q.qxt2_12;
mirror->knm_q.qyt1_12 = knm_q.qyt2_12;
}
if (CALC_MR_KNM(mirror,21)){
mirror->knm_q.qxt1_21 = knm_q.qxt1_21;
mirror->knm_q.qyt1_21 = knm_q.qyt1_21;
mirror->knm_q.qxt2_21 = qLx_21;
mirror->knm_q.qyt2_21 = qLy_21;
mirror->knm_q.qxt1_21 = knm_q.qxt2_21;
mirror->knm_q.qyt1_21 = knm_q.qyt2_21;
}
if (CALC_MR_KNM(mirror,22)){
mirror->knm_q.qxt1_22 = knm_q.qxt1_22;
mirror->knm_q.qyt1_22 = knm_q.qyt1_22;
mirror->knm_q.qxt2_22 = qLx_22;
mirror->knm_q.qyt2_22 = qLy_22;
mirror->knm_q.qxt1_22 = knm_q.qxt2_22;
mirror->knm_q.qyt1_22 = knm_q.qyt2_22;
}
} else if (i == NUM_KNM_TYPES-1 && mirror->knm_change_q == 2){
// If the final solver should handle mode mismatch...
if (CALC_MR_KNM(mirror,11)){
mirror->knm_q.qxt2_11 = knm_q.qxt2_11;
mirror->knm_q.qyt2_11 = knm_q.qyt2_11;
}
if (CALC_MR_KNM(mirror,12)){
mirror->knm_q.qxt2_12 = knm_q.qxt2_12;
mirror->knm_q.qyt2_12 = knm_q.qyt2_12;
}
if (CALC_MR_KNM(mirror,21)){
mirror->knm_q.qxt2_21 = knm_q.qxt2_21;
mirror->knm_q.qyt2_21 = knm_q.qyt2_21;
}
if (CALC_MR_KNM(mirror,22)){
mirror->knm_q.qxt2_22 = knm_q.qxt2_22;
mirror->knm_q.qyt2_22 = knm_q.qyt2_22;
}
}
}
......@@ -1287,8 +1309,6 @@ int set_k_mirror(int mirror_index) {
for (m = 0; m < inter.num_fields; m++) {
get_tem_modes_from_field_index(&n2, &m2, m);
//printf("%i%i->%i%i %s\n",n1,m1,n2,m2,complex_form15(mirror->knm_no_rgouy.k11[n][m]));
// 1) the old map code had wrong qx1/qx2 for k12/k21 at this place and I assume
// this was copied over. Fixed this now.
// 2) the old code uses the `turned' q's to compute thre reverse gouy, so
......
......@@ -885,9 +885,6 @@ void fill_unn_cache(unn_cache_t *c, u_n_accel_t *acc1, u_n_accel_t *acc2){
assert(acc1 != NULL);
assert(acc2 != NULL);
assert(ceq(acc1->negK_2q, acc2->negK_2q));
assert(acc1->sqrt2_wz == acc2->sqrt2_wz);
c->acc1 = acc1;
c->acc2 = acc2;
......@@ -986,7 +983,7 @@ complex_t do_romhom_real_int(romhom_weights_t *rom, double** u_xy, unn_cache_t *
// integrand, the conjugates do no matter
for(p=0; p<cx->num_nodes; p++) {
for(q=0; q<cy->num_nodes; q++) {
u_xy[p][q] = cx->values[idx][p] * cy->values[idy][q];
u_xy[q][p] = cx->values[idx][q] * cy->values[idy][p];
}
}
......
......@@ -869,7 +869,7 @@ void compute_mirror_knm_romhom(mirror_t *mirror, double nr1, double nr2, int mis
return;
romhom_weights_t *rom = mirror->map_rom;
mirror_knm_t *knm = &(mirror->knm_romhom);
mirror_knm_t *knm = &(mirror->knm_romhom);
if (!IS_MIRROR_KNM_ALLOCD(knm))
bug_error("mirror romhom knm has not been allocated");
......@@ -887,30 +887,47 @@ void compute_mirror_knm_romhom(mirror_t *mirror, double nr1, double nr2, int mis
mirror_knm_q_t *kq = &(mirror->knm_q);
if(!ceq(kq->qxt1_11, kq->qxt2_11) || !ceq(kq->qyt1_11, kq->qyt2_11)) {
warn("Mismatch on side 1 reflection at %s, can't compute coupling using ROM\n", mirror->name);
warn("1->1: qx = %s, qx' = %s\n", complex_form15(kq->qxt1_11), complex_form15(kq->qxt2_11));
warn("1->1: qy = %s, qy' = %s\n", complex_form15(kq->qyt1_11), complex_form15(kq->qyt2_11));
}
if(!ceq(kq->qxt1_22, kq->qxt2_22) || !ceq(kq->qyt1_22, kq->qyt2_22)) {
warn("Mismatch on side 2 reflection at %s, can't compute coupling using ROM\n", mirror->name);
warn("2->2: qx = %s, qx' = %s\n", complex_form15(kq->qxt1_22), complex_form15(kq->qxt2_22));
warn("2->2: qy = %s, qy' = %s\n", complex_form15(kq->qyt1_22), complex_form15(kq->qyt2_22));
}
if(!ceq(kq->qxt1_12, kq->qxt2_12) || !ceq(kq->qxt1_21, kq->qxt2_21)
|| !ceq(kq->qyt1_12, kq->qyt2_12) || !ceq(kq->qyt1_21, kq->qyt2_21)) {
warn("Mismatch on transmission at %s, can't compute coupling using ROM\n", mirror->name);
warn("1->2: qx = %s, qx' = %s\n", complex_form15(kq->qxt1_12), complex_form15(kq->qxt2_12));
warn("1->2: qy = %s, qy' = %s\n", complex_form15(kq->qyt1_12), complex_form15(kq->qyt2_12));
warn("2->1: qx = %s, qx' = %s\n", complex_form15(kq->qxt1_21), complex_form15(kq->qxt2_21));
warn("2->1: qy = %s, qy' = %s\n", complex_form15(kq->qyt1_21), complex_form15(kq->qyt2_21));
}
// K11
if (CALC_MR_KNM(mirror,11)) {
assert(ceq(kq->qxt1_11, kq->qyt2_11));
u_nm_accel_get(acc_11_nr1_1, max_n, max_m, kq->qxt1_11, kq->qyt1_11, nr1);
u_nm_accel_get(acc_11_nr1_2, max_n, max_m, kq->qxt2_11, kq->qyt2_11, nr1);
}
// K22
if (CALC_MR_KNM(mirror,22)) {
assert(ceq(kq->qxt1_22, kq->qyt2_22));
u_nm_accel_get(acc_22_nr2_1, max_n, max_m, kq->qxt1_22, kq->qyt1_22, nr2);
u_nm_accel_get(acc_22_nr2_2, max_n, max_m, kq->qxt2_22, kq->qyt2_22, nr2);
}
// K12
if (CALC_MR_KNM(mirror,12)) {
assert(ceq(kq->qxt1_12, kq->qyt2_12));
u_nm_accel_get(acc_12_nr1_1, max_n, max_m, kq->qxt1_12, kq->qyt1_12, nr1);
u_nm_accel_get(acc_12_nr2_2, max_n, max_m, kq->qxt2_12, kq->qyt2_12, nr2);
}
// K21
if (CALC_MR_KNM(mirror,21)) {
assert(ceq(kq->qxt1_21, kq->qyt2_21));
u_nm_accel_get(acc_21_nr2_1, max_n, max_m, kq->qxt1_21, kq->qyt1_21, nr2);
u_nm_accel_get(acc_21_nr1_2, max_n, max_m, kq->qxt2_21, kq->qyt2_21, nr1);
}
......
......@@ -6952,8 +6952,8 @@ void read_conf(const char *command_string) {
for(j=0; j<NUM_KNM_TYPES; j++){
int i = config_setting[j] - '0';
if(!(i==1 || i==2))
gerror("knm_order can only contain digits 1 or 2 not '%c'",i);
if(!(i==1 || i==2 || i==4))
gerror("knm_order can only contain digits 1, 2 or 4 not '%c'",i);
mirror->knm_order[j] = i;
}
......@@ -7732,16 +7732,16 @@ void read_rom_map(romhom_weights_t *rom){
for(n=0; n<xQ1; n++){
for(m=0; m<yQ1; m++){
rom->w_Q1Q2[n][m] = z_pl_z(wQ1[n][m], wQ2[n][m]);
rom->w_Q1Q3[n][m] = z_pl_z(wQ1[n][m], wQ3[n][m]);
rom->w_Q1Q4[n][m] = z_pl_z(wQ1[n][m], wQ4[n][m]);
rom->w_Q1Q2[m][n] = z_pl_z(wQ1[m][n], wQ2[m][n]);
rom->w_Q1Q3[m][n] = z_pl_z(wQ1[m][n], wQ3[m][n]);
rom->w_Q1Q4[m][n] = z_pl_z(wQ1[m][n], wQ4[m][n]);
rom->w_Q2Q3[n][m] = z_pl_z(wQ2[n][m], wQ3[n][m]);
rom->w_Q2Q4[n][m] = z_pl_z(wQ2[n][m], wQ4[n][m]);
rom->w_Q2Q3[m][n] = z_pl_z(wQ2[m][n], wQ3[m][n]);
rom->w_Q2Q4[m][n] = z_pl_z(wQ2[m][n], wQ4[m][n]);
rom->w_Q3Q4[n][m] = z_pl_z(wQ3[n][m], wQ4[n][m]);
rom->w_Q3Q4[m][n] = z_pl_z(wQ3[m][n], wQ4[m][n]);
rom->w_Q1Q2Q3Q4[n][m] = z_pl_z(rom->w_Q1Q2[n][m], rom->w_Q3Q4[n][m]);
rom->w_Q1Q2Q3Q4[m][n] = z_pl_z(rom->w_Q1Q2[m][n], rom->w_Q3Q4[m][n]);
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment