Skip to content

Commit

Permalink
fixed a bug on the version that imposes water height at inflow
Browse files Browse the repository at this point in the history 
git-svn-id: svn://svn.sissa.it/openship@567 958231e9-8e66-0410-a4a0-f664109d2741
  • Loading branch information
@andreamola
andreamola committed Dec 28, 2015
1 parent dd4c4d0 commit a90b19a
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 56 deletions.
1 change: 1 addition & 0 deletions include/restart_nonlinear_problem_diff.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,7 @@ public NewtonArgument
if ( (comp_dom.flags[i] & water) &&
(comp_dom.flags[i] & near_boat)==0 &&
(comp_dom.flags[i] & transom_on_water)==0 &&
(comp_dom.flags[i] & near_inflow)==0 &&
(comp_dom.vector_constraints.is_constrained(3*i+2))==0 )
{
//cout<<"W "<<i<<" ("<<3*i+2<<")"<<endl;
Expand Down
116 changes: 62 additions & 54 deletions source/free_surface.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4411,15 +4411,17 @@ std::cout<<"Preparing interpolated solution for restart"<<std::endl;
loc_y_smooth_res[i] += loc_stiffness_matrix[i][j]*(y_displs[j]-(1-blend_factor)*comp_dom.old_map_points(3*local_dof_indices[j]+1));
}
if ( !constraints.is_constrained(local_dof_indices[i]) &&
!(comp_dom.flags[local_dof_indices[i]] & transom_on_water) )
!(comp_dom.flags[local_dof_indices[i]] & transom_on_water) &&
!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
{
unsigned int ii = local_dof_indices[i];
eta_res[ii] += loc_eta_res[i].val();
phi_res[ii] += loc_phi_res[i].val();

}
if ( !(constraints.is_constrained(local_dof_indices[i])) &&
!(comp_dom.flags[local_dof_indices[i]] & edge) )
!(comp_dom.flags[local_dof_indices[i]] & edge) &&
!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
{
unsigned int ii = local_dof_indices[i];
x_smoothing_res[ii] += loc_x_smooth_res[i].val();
Expand Down Expand Up @@ -4619,10 +4621,10 @@ std::cout<<"Preparing interpolated solution for restart"<<std::endl;
RestartNonlinearProblemDiff rest_nonlin_prob_diff(*this,comp_dom,t,y,yp,jacobian_dot_matrix);
std::map<unsigned int,unsigned int> &map_diff = rest_nonlin_prob_diff.indices_map;

/*

// these lines test the correctness of the jacobian for the
// restart (reduced) nonlinear problem

/*
Vector<double> restart_prob_solution(rest_nonlin_prob_diff.n_dofs());
Vector<double> restart_prob_residual(rest_nonlin_prob_diff.n_dofs());
for (std::map<unsigned int, unsigned int>::iterator it = map_diff.begin(); it != map_diff.end(); ++it)
Expand All @@ -4646,7 +4648,7 @@ std::cout<<"Preparing interpolated solution for restart"<<std::endl;
rest_nonlin_prob_diff.residual(restart_prob_residual,restart_prob_solution);
cout<<"----------Test---------"<<endl;
for (unsigned int i=0; i<rest_nonlin_prob_diff.n_dofs(); ++i)
if (fabs(restart_prob_residual(i)-delta_res(i)) > 1e-10)
if (fabs(restart_prob_residual(i)-delta_res(i)) > 1e-15)
cout<<i<<" "<<delta_res(i)<<" vs "<<restart_prob_residual(i)<<" err "<<restart_prob_residual(i)-delta_res(i)<<" "<<sys_comp(i)<<endl;
delta_res*=-1;
delta_res.add(restart_prob_residual);
Expand Down Expand Up @@ -4681,9 +4683,9 @@ std::cout<<"Preparing interpolated solution for restart"<<std::endl;
//output_results(filename1, t, y, yp);



/*
// these lines test the jacobian of the DAE system
/*

Vector<double> delta_y(this->n_dofs());
Vector<double> delta_res(this->n_dofs());

Expand Down Expand Up @@ -4969,11 +4971,11 @@ int FreeSurface<dim>::jacobian(const double t,
bem.solve_system(bem_phi, bem_dphi_dn, bem_bc);

/*
cout<<"44!!!!!!"<<endl;
for (SparsityPattern::iterator col=jacobian_sparsity_pattern.begin(1484); col!=jacobian_sparsity_pattern.end(1484); ++col)
cout<<"JACOBIANO!!!!!!"<<endl;
for (SparsityPattern::iterator col=jacobian_sparsity_pattern.begin(80); col!=jacobian_sparsity_pattern.end(80); ++col)
{
unsigned int j = col->column();
cout<<j<<"("<<jacobian_matrix(1484,j)<<") ";
cout<<j<<"("<<jacobian_matrix(80,j)<<") ";
}
cout<<endl;
//*/
Expand Down Expand Up @@ -5304,7 +5306,7 @@ int FreeSurface<dim>::residual_and_jacobian(const double t,
if (comp_dom.flags[i] & near_inflow)
{
working_map_points(3*i+2) = comp_dom.old_map_points(3*i+2);
jacobian_matrix.add(3*i+2,3*i+2,1.0);
jacobian_matrix.set(3*i+2,3*i+2,1.0);
}
//cout<<"& "<<3*i<<" "<<3*i+1<<endl;
}
Expand Down Expand Up @@ -7035,26 +7037,29 @@ int FreeSurface<dim>::residual_and_jacobian(const double t,
loc_y_smooth_res[i] += loc_stiffness_matrix[i][j]*(y_displs[j]-(1-blend_factor)*comp_dom.old_map_points(3*local_dof_indices[j]+1));
}
if ( !constraints.is_constrained(local_dof_indices[i]) &&
!(comp_dom.flags[local_dof_indices[i]] & transom_on_water))
!(comp_dom.flags[local_dof_indices[i]] & transom_on_water) )
{
unsigned int ii = local_dof_indices[i];
eta_res[ii] += loc_eta_res[i].val();
if (!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
eta_res[ii] += loc_eta_res[i].val();
phi_res[ii] += loc_phi_res[i].val();
//cout<<"* "<<cell<<" "<<local_dof_indices[i]<<" "<<loc_phi_res[i]<<endl;
for (unsigned int j=0;j<dofs_per_cell;++j)
{
unsigned int jj = local_dof_indices[j];
for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(3*ii+2,
3*jj+k,
loc_eta_res[i].fastAccessDx(3*j+k));
jacobian_matrix.add(3*ii+2,
jj+comp_dom.vector_dh.n_dofs(),
loc_eta_res[i].fastAccessDx(3*dofs_per_cell+j));
jacobian_matrix.add(3*ii+2,
jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(4*dofs_per_cell+j));

if (!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
{
for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(3*ii+2,
3*jj+k,
loc_eta_res[i].fastAccessDx(3*j+k));
jacobian_matrix.add(3*ii+2,
jj+comp_dom.vector_dh.n_dofs(),
loc_eta_res[i].fastAccessDx(3*dofs_per_cell+j));
jacobian_matrix.add(3*ii+2,
jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(4*dofs_per_cell+j));
}
for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(ii+comp_dom.vector_dh.n_dofs(),
3*jj+k,
Expand All @@ -7065,17 +7070,19 @@ int FreeSurface<dim>::residual_and_jacobian(const double t,
jacobian_matrix.add(ii+comp_dom.vector_dh.n_dofs(),
jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_phi_res[i].fastAccessDx(4*dofs_per_cell+j));

for (unsigned int k=0; k<dim; ++k)
jacobian_dot_matrix.add(3*ii+2,
3*jj+k,
loc_eta_res[i].fastAccessDx(5*dofs_per_cell+3*j+k));
jacobian_dot_matrix.add(3*ii+2,
jj+comp_dom.vector_dh.n_dofs(),
loc_eta_res[i].fastAccessDx(8*dofs_per_cell+j));
jacobian_dot_matrix.add(3*ii+2,
jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(9*dofs_per_cell+j));
if (!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
{
for (unsigned int k=0; k<dim; ++k)
jacobian_dot_matrix.add(3*ii+2,
3*jj+k,
loc_eta_res[i].fastAccessDx(5*dofs_per_cell+3*j+k));
jacobian_dot_matrix.add(3*ii+2,
jj+comp_dom.vector_dh.n_dofs(),
loc_eta_res[i].fastAccessDx(8*dofs_per_cell+j));
jacobian_dot_matrix.add(3*ii+2,
jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(9*dofs_per_cell+j));
}
for (unsigned int k=0; k<dim; ++k)
jacobian_dot_matrix.add(ii+comp_dom.vector_dh.n_dofs(),
3*jj+k,
Expand All @@ -7087,23 +7094,25 @@ int FreeSurface<dim>::residual_and_jacobian(const double t,
jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_phi_res[i].fastAccessDx(9*dofs_per_cell+j));
}

for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(3*ii+2,
k+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(k+10*dofs_per_cell));
for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(3*ii+2,
k+3+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(k+3+10*dofs_per_cell));
for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(3*ii+2,
k+6+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(k+6+10*dofs_per_cell));
for (unsigned int k=0; k<4; ++k)
jacobian_matrix.add(3*ii+2,
k+9+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(k+9+10*dofs_per_cell));
if (!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
{
for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(3*ii+2,
k+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(k+10*dofs_per_cell));
for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(3*ii+2,
k+3+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(k+3+10*dofs_per_cell));
for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(3*ii+2,
k+6+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(k+6+10*dofs_per_cell));
for (unsigned int k=0; k<4; ++k)
jacobian_matrix.add(3*ii+2,
k+9+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
loc_eta_res[i].fastAccessDx(k+9+10*dofs_per_cell));
}
for (unsigned int k=0; k<dim; ++k)
jacobian_matrix.add(ii+comp_dom.vector_dh.n_dofs(),
k+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
Expand Down Expand Up @@ -8112,8 +8121,7 @@ for (unsigned int i=0; i<this->n_dofs(); ++i)
{
unsigned int j = col->column();
//cout<<" "<<i<<" "<<j<<" "<<jacobian_matrix(i,j)+alpha*jacobian_dot_matrix(i,j)<<endl;
//if ( (j > (int)i-preconditioner_band) &&
// (j > (int)i+preconditioner_band) )

jacobian_preconditioner_matrix.set(i,j,jacobian_matrix(i,j)+alpha*jacobian_dot_matrix(i,j));
}

Expand Down
4 changes: 2 additions & 2 deletions waveBem.prm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,8 +156,8 @@ subsection Boundary Conditions ID Numbers
set Wall Surface 1 ID = 3
set Wall Surface 2 ID = 4
set Wall Surface 3 ID = 0
set Inflow Surface 1 ID = 0
set Inflow Surface 2 ID = 0
set Inflow Surface 1 ID = 9
set Inflow Surface 2 ID = 10
set Inflow Surface 3 ID = 0
set Pressure Surface ID = 0
set Free Surface Edge On Boat ID = 20
Expand Down

0 comments on commit a90b19a

Please sign in to comment.