Skip to content

Commit c69fafc

Browse files
committed
1. add Exx::set_Ds_delta()
1 parent c476389 commit c69fafc

File tree

4 files changed

+114
-16
lines changed

4 files changed

+114
-16
lines changed

include/RI/global/Map_Operator.h

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ namespace RI
1212

1313
namespace Map_Operator
1414
{
15+
// m1+m2
1516
template<typename Tkey, typename Tvalue>
1617
std::map<Tkey,Tvalue> operator+(const std::map<Tkey,Tvalue> &m1, const std::map<Tkey,Tvalue> &m2)
1718
{
@@ -51,6 +52,51 @@ namespace Map_Operator
5152
m[item.first] = item.second;
5253
return m;
5354
}
55+
56+
// m1-m2
57+
template<typename Tkey, typename Tvalue>
58+
std::map<Tkey,Tvalue> operator-(const std::map<Tkey,Tvalue> &m1, const std::map<Tkey,Tvalue> &m2)
59+
{
60+
std::map<Tkey,Tvalue> m;
61+
auto ptr1 = m1.begin();
62+
auto ptr2 = m2.begin();
63+
while(ptr1!=m1.end() && ptr2!=m2.end())
64+
{
65+
if(ptr1->first == ptr2->first)
66+
{
67+
m.emplace_hint(m.end(), ptr1->first, ptr1->second - ptr2->second);
68+
++ptr1;
69+
++ptr2;
70+
}
71+
else if(ptr1->first < ptr2->first)
72+
{
73+
m.emplace_hint(m.end(), ptr1->first, ptr1->second);
74+
++ptr1;
75+
}
76+
else
77+
{
78+
m.emplace_hint(m.end(), ptr2->first, -ptr2->second);
79+
++ptr2;
80+
}
81+
}
82+
m.insert(ptr1, m1.end());
83+
while(ptr2!=m2.end())
84+
{
85+
m.emplace_hint(m.end(), ptr2->first, -ptr2->second);
86+
++ptr2;
87+
}
88+
return m;
89+
}
90+
91+
// -m_in
92+
template<typename Tkey, typename Tvalue>
93+
std::map<Tkey,Tvalue> operator-(const std::map<Tkey,Tvalue> &m_in)
94+
{
95+
std::map<Tkey,Tvalue> m;
96+
for(const auto &item : m_in)
97+
m.emplace_hint(m.end(), item.first, -item.second);
98+
return m;
99+
}
54100
}
55101

56102
}

include/RI/physics/Exx.h

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,10 @@ class Exx
4646
const std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds,
4747
const Tdata_real &threshold_D,
4848
const std::string &save_name_suffix="");
49+
void set_Ds_delta(
50+
const std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds,
51+
const Tdata_real &threshold_D,
52+
const std::string &save_name_suffix="");
4953
void set_dCs(
5054
const std::array<std::map<TA, std::map<TAC, Tensor<Tdata>>>,Npos> &dCs,
5155
const Tdata_real &threshold_dC,
@@ -80,6 +84,7 @@ class Exx
8084
bool C=false;
8185
bool V=false;
8286
bool D=false;
87+
bool D_delta=false;
8388
bool dC=false;
8489
bool dV=false;
8590
};

include/RI/physics/Exx.hpp

Lines changed: 61 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include "Exx.h"
99
#include "../ri/Cell_Nearest.h"
1010
#include "../ri/Label.h"
11+
#include "../global/Map_Operator.h"
1112

1213
#include <cassert>
1314

@@ -64,6 +65,41 @@ void Exx<TA,Tcell,Ndim,Tdata>::set_Ds(
6465
this->lri.set_tensors_map2( Ds, Label::ab::a2b1, {{"threshold_filter", threshold_D}}, "Ds_a2b1_"+save_name_suffix );
6566
this->lri.set_tensors_map2( Ds, Label::ab::a2b2, {{"threshold_filter", threshold_D}}, "Ds_a2b2_"+save_name_suffix );
6667
this->flag_finish.D = true;
68+
this->flag_finish.D_delta = false;
69+
70+
//if()
71+
this->post_2D.saves["Ds_"+save_name_suffix] = this->post_2D.set_tensors_map2(Ds);
72+
}
73+
74+
template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
75+
void Exx<TA,Tcell,Ndim,Tdata>::set_Ds_delta(
76+
const std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds,
77+
const Tdata_real &threshold_D,
78+
const std::string &save_name_suffix)
79+
{
80+
using namespace Map_Operator;
81+
82+
assert(flag_finish.D);
83+
for(const Label::ab label : {Label::ab::a1b1, Label::ab::a1b2, Label::ab::a2b1, Label::ab::a2b2})
84+
{
85+
this->lri.set_tensors_map2(
86+
Ds,
87+
label,
88+
{{"flag_filter", false}},
89+
"Ds_tmp" );
90+
this->lri.set_tensors_map2(
91+
this->lri.data_pool["Ds_tmp"].Ds_ab - this->lri.data_pool["Ds_"+Label_Tools::get_name(label)+"_"+save_name_suffix].Ds_ab,
92+
label,
93+
{{"flag_period", false}, {"flag_comm", false}, {"flag_filter", true}, {"threshold_filter", threshold_D}},
94+
"Ds_delta_"+Label_Tools::get_name(label)+"_"+save_name_suffix);
95+
this->lri.set_tensors_map2(
96+
this->lri.data_pool["Ds_delta_"+Label_Tools::get_name(label)+"_"+save_name_suffix].Ds_ab + this->lri.data_pool["Ds_"+Label_Tools::get_name(label)+"_"+save_name_suffix].Ds_ab,
97+
label,
98+
{{"flag_period", false}, {"flag_comm", false}, {"flag_filter", false}},
99+
"Ds_"+Label_Tools::get_name(label)+"_"+save_name_suffix);
100+
}
101+
this->flag_finish.D_delta = true;
102+
this->flag_finish.D = true;
67103

68104
//if()
69105
this->post_2D.saves["Ds_"+save_name_suffix] = this->post_2D.set_tensors_map2(Ds);
@@ -102,18 +138,32 @@ template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
102138
void Exx<TA,Tcell,Ndim,Tdata>::cal_Hs(
103139
const std::array<std::string,3> &save_names_suffix) // "Cs","Vs","Ds"
104140
{
141+
using namespace Map_Operator;
142+
105143
assert(this->flag_finish.stru);
106-
assert(this->flag_finish.C);
107-
assert(this->flag_finish.V);
108-
assert(this->flag_finish.D);
109144

145+
assert(this->flag_finish.C);
110146
this->lri.data_ab_name[Label::ab::a ] = "Cs_a_" +save_names_suffix[0];
111147
this->lri.data_ab_name[Label::ab::b ] = "Cs_b_" +save_names_suffix[0];
148+
149+
assert(this->flag_finish.V);
112150
this->lri.data_ab_name[Label::ab::a0b0] = "Vs_" +save_names_suffix[1];
113-
this->lri.data_ab_name[Label::ab::a1b1] = "Ds_a1b1_"+save_names_suffix[2];
114-
this->lri.data_ab_name[Label::ab::a1b2] = "Ds_a1b2_"+save_names_suffix[2];
115-
this->lri.data_ab_name[Label::ab::a2b1] = "Ds_a2b1_"+save_names_suffix[2];
116-
this->lri.data_ab_name[Label::ab::a2b2] = "Ds_a2b2_"+save_names_suffix[2];
151+
152+
if(!this->flag_finish.D_delta)
153+
{
154+
assert(this->flag_finish.D);
155+
this->lri.data_ab_name[Label::ab::a1b1] = "Ds_a1b1_"+save_names_suffix[2];
156+
this->lri.data_ab_name[Label::ab::a1b2] = "Ds_a1b2_"+save_names_suffix[2];
157+
this->lri.data_ab_name[Label::ab::a2b1] = "Ds_a2b1_"+save_names_suffix[2];
158+
this->lri.data_ab_name[Label::ab::a2b2] = "Ds_a2b2_"+save_names_suffix[2];
159+
}
160+
else
161+
{
162+
this->lri.data_ab_name[Label::ab::a1b1] = "Ds_delta_a1b1_"+save_names_suffix[2];
163+
this->lri.data_ab_name[Label::ab::a1b2] = "Ds_delta_a1b2_"+save_names_suffix[2];
164+
this->lri.data_ab_name[Label::ab::a2b1] = "Ds_delta_a2b1_"+save_names_suffix[2];
165+
this->lri.data_ab_name[Label::ab::a2b2] = "Ds_delta_a2b2_"+save_names_suffix[2];
166+
}
117167

118168
std::vector<std::map<TA, std::map<TAC, Tensor<Tdata>>>> Hs_vec(1);
119169
this->lri.coefficients = {nullptr};
@@ -123,20 +173,15 @@ void Exx<TA,Tcell,Ndim,Tdata>::cal_Hs(
123173
Label::ab_ab::a0b0_a2b1,
124174
Label::ab_ab::a0b0_a2b2},
125175
Hs_vec);
126-
this->Hs = std::move(Hs_vec[0]);
176+
if(this->flag_finish.D_delta)
177+
this->Hs = this->Hs + Hs_vec[0];
178+
else
179+
this->Hs = std::move(Hs_vec[0]);
127180

128181
//if()
129182
this->energy = this->post_2D.cal_energy(
130183
this->post_2D.saves["Ds_"+save_names_suffix[2]],
131184
this->post_2D.set_tensors_map2(this->Hs) );
132-
133-
this->lri.data_ab_name[Label::ab::a ] = "Cs_a_" +save_names_suffix[0];
134-
this->lri.data_ab_name[Label::ab::b ] = "Cs_b_" +save_names_suffix[0];
135-
this->lri.data_ab_name[Label::ab::a0b0] = "Vs_" +save_names_suffix[1];
136-
this->lri.data_ab_name[Label::ab::a1b1] = "Ds_a1b1_"+save_names_suffix[2];
137-
this->lri.data_ab_name[Label::ab::a1b2] = "Ds_a1b2_"+save_names_suffix[2];
138-
this->lri.data_ab_name[Label::ab::a2b1] = "Ds_a2b1_"+save_names_suffix[2];
139-
this->lri.data_ab_name[Label::ab::a2b2] = "Ds_a2b2_"+save_names_suffix[2];
140185
}
141186

142187
template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>

unittests/physics/Exx-test.hpp

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,8 @@ namespace Exx_Test
2727
exx.cal_Hs();
2828
exx.cal_force();
2929
exx.cal_stress();
30+
exx.set_Ds_delta({}, 1E-4);
31+
exx.cal_Hs();
3032

3133
MPI_Finalize();
3234
}

0 commit comments

Comments
 (0)