> On 20. Sep 2017, at 22:51, Matthew Knepley <[email protected]> wrote:
>
> On Wed, Sep 20, 2017 at 3:46 PM, Maximilian Hartig <[email protected]
> <mailto:[email protected]>> wrote:
>
>> On 20. Sep 2017, at 19:05, Matthew Knepley <[email protected]
>> <mailto:[email protected]>> wrote:
>>
>> On Wed, Sep 20, 2017 at 12:57 PM, Maximilian Hartig
>> <[email protected] <mailto:[email protected]>> wrote:
>>> On 20. Sep 2017, at 18:17, Matthew Knepley <[email protected]
>>> <mailto:[email protected]>> wrote:
>>>
>>> On Wed, Sep 20, 2017 at 11:46 AM, Maximilian Hartig
>>> <[email protected] <mailto:[email protected]>> wrote:
>>> Hello,
>>>
>>> I’m trying to implement plasticity using petscFE but I am quite stuck since
>>> a while. Here’s what I’m trying to do:
>>>
>>> I have a TS which solves the following equation:
>>> gradient(stress) +Forces = density*acceleration
>>> where at the moment stress is a linear function of the strain and hence the
>>> gradient of the displacement. This works fine. Now I want to compare the
>>> stress to a reference value and if it lies above this yield stress, I have
>>> to reevaluate the stress at the respective location. Then I need to update
>>> the plastic strain / yield stress at this location.
>>> I tried doing that first by solving three fields at the same time:
>>> displacements, stresses and yield stress. This failed.
>>> Then, I tried solving only for displacement increments, storing the
>>> displacements, stresses and yield stress from the past time step in an
>>> auxiliary field. The auxiliary fields are updated after each time step with
>>> a second SNES, using the displacement increments from the current,
>>> converged time step. This also failed.
>>> In both cases the code had problems converging and when it did, I ended up
>>> with negative plastic strain. This is not possible and I don’t know how it
>>> happens because I explicitly only increment the plastic strain when the
>>> increment is positive.
>>>
>>> I am sure there is an easy solution to how I can update the internal
>>> variables and determine the correct stress for the residual but I just
>>> cannot figure it out. I’d be thankful for any hints.
>>>
>>> It looks like there are two problems above:
>>>
>>> 1) Convergence
>>>
>>> For any convergence question, we at minimum need to see the output of
>>>
>>> -snes_view -snes_converged_reason -snes_monitor
>>> -ksp_monitor_true_residual -snes_linesearch_monitor
>>>
>>> However, this does not seem to be the main issue.
>>>
>>> 2) Negative plastic strain
>>
>> This is what I’m mainly concerned with.
>>>
>>> If the system really converged (I cannot tell without other information),
>>> then the system formulation is wrong. Of course, its
>>> really easy to check by just plugging your solution into the residual
>>> function too. I do not understand your explanation above
>>> completely however. Do you solve for the plastic strain or the increment?
>>
>> I am trying to find a formulation that works and I think there is a core
>> concept I am just not “getting”.
>> I want to solve for the displacements.
>> This works fine in an elastic case. When plasticity is involved, I need to
>> determine the actual stress for my residual evaluation and I have not found
>> a way to do that.
>> All formulations for stress I found in literature use strain increments so I
>> tried to just solve for increments each timestep and then add them together
>> in tspoststep. But I still need to somehow evaluate the stress for my
>> displacement increment residuals. So currently, I have auxiliary fields with
>> the stress and the plastic strain.
>>
>> First question: Don't you get stress by just applying a local operator,
>> rather than a solve?
> That depends on the type of plasticity.
>
> What type of plasticity is not local?
I did not express myself correctly. The plasticity is local, but for nonlinear
hardening I would have to iteratively solve to get the correct stress and
plastic strain from the displacement increments.
>
> For a linear hardening formulation it is correct that I could just apply a
> local operator. I’d be happy with that for now. But I’d still need to save
> stress state and plastic strain to determine whether or not I’m dealing with
> a plasticity. I don’t know how to do that inside the residual evaluation.
>
> I do not know what you mean by this, meaning why you can't just save these as
> auxiliary fields. Also, it would seem to be enough to have the old
> displacement and the plastic strain.
Yes, I can update the auxiliary fields but only after I solved for the
displacements in my understanding. I still need the stress though as I have to
determine whether I am in the plastic or the elastic domain.
>
> Plus DMProjectField seems to have problems evaluating the gradient when
> boundary conditions are imposed.
>
> There are several examples where we do exactly this. Can you show me what you
> mean by this?
Yes, of course. I sent an example a week ago but maybe there was a problem with
the attached files. I’ll copy and paste the code and the gmsh file as text
below.
Thanks,
Max
#include <petscdm.h>
#include <petscdmlabel.h>
#include <petscds.h>
#include <petscdmplex.h>
#include <petscksp.h>
#include <petscsnes.h>
#include <petscts.h>
#include <math.h>
#include <petscsys.h>
/* define the pointwise functions we'd like to project */
void projectstress(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt
uOff[],
const PetscInt uOff_x[], const PetscScalar u[], const
PetscScalar u_t[],
const PetscScalar u_x[], const PetscInt aOff[], const
PetscInt aOff_x[],
const PetscScalar a[], const PetscScalar a_t[], const
PetscScalar a_x[],
PetscReal t, const PetscScalar x[], PetscInt
numConstants,
const PetscScalar constants[], PetscScalar v[]){
const PetscReal mu =76.923076923, lbda=115.384615385;
PetscInt Ncomp = dim;
PetscInt comp,d;
PetscReal sigma[dim*dim];
for(comp=0;comp<Ncomp;comp++){
for(d=0;d<dim;d++){
sigma[comp*dim+d]=mu*(u_x[comp*dim+d]+u_x[d*dim+comp]);
}
for(d=0;d<dim;d++){
sigma[comp*dim+comp]+=lbda*u_x[d*dim+d];
}
}
for(d=0;d<dim;d++){
v[d] = sigma[d*dim+d];
}
v[3] = sigma[0*dim+1];
v[4] = sigma[0*dim+2];
v[5] = sigma[1*dim+2];
}
void projectdisplacement(PetscInt dim, PetscInt Nf, PetscInt NfAux, const
PetscInt uOff[],
const PetscInt uOff_x[], const PetscScalar u[], const
PetscScalar u_t[],
const PetscScalar u_x[], const PetscInt aOff[], const
PetscInt aOff_x[],
const PetscScalar a[], const PetscScalar a_t[], const
PetscScalar a_x[],
PetscReal t, const PetscScalar x[], PetscInt
numConstants,
const PetscScalar constants[], PetscScalar v[]){
PetscInt d;
for(d=0;d<dim;d++){
v[d] =u[d];
}
}
static PetscErrorCode initial_displacement_vector(PetscInt dim, PetscReal time,
const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
u[0]=0.0;
u[1]=0.1*pow(x[2],2);
u[2]=0.1*x[2];
return 0;
}
static PetscErrorCode expected_stress_vector(PetscInt dim, PetscReal time,
const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
const PetscReal mu =76.923076923, lbda=115.384615385;
PetscReal strain[dim*dim];
PetscReal gradu[dim*dim];
PetscReal stress[dim*dim];
PetscInt i,j;
/* gradient of the displacement field: */
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
gradu[i*dim+j]=0.0;
}
}
gradu[1*dim+2]=0.2*x[2];
gradu[2*dim+2]=0.1;
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
strain[i*dim+j]=0.5*(gradu[i*dim+j]+gradu[j*dim+i]);
}
}
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
stress[i*dim+j]=2.0*mu*strain[i*dim+j];
}
for(j=0;j<dim;j++){
stress[i*dim+i]+=lbda*strain[j*dim+j];
}
}
for(i=0;i<dim;i++){
u[i] = stress[i*dim+i];
}
u[3] = stress[0*dim+1];
u[4] = stress[0*dim+2];
u[5] = stress[1*dim+2];
return 0;
}
static PetscErrorCode zero_stress(PetscInt dim, PetscReal time, const PetscReal
x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
const PetscInt Ncomp = 2*dim;
PetscInt comp;
for (comp=0;comp<Ncomp;comp++) u[comp]=0.0;
return 0;
}
static PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const PetscReal
x[], PetscInt Nf, PetscScalar *u, void *ctx)
{
PetscInt comp;
for (comp=0;comp<dim;comp++) u[comp]=0.0;
return 0;
}
static PetscErrorCode SetupDiscretization(DM dm){
PetscInt dim = 3;
PetscFE fe_displacement, fe_stress;
PetscDS prob;
PetscErrorCode ierr;
PetscBool simplex = PETSC_TRUE;
PetscQuadrature q;
PetscInt order;
/* get the dimension of the problem from the DM */
ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr);
/* Creating the FE for the displacement */
ierr = PetscFECreateDefault(dm, dim, dim,
simplex,"disp_",PETSC_DEFAULT,&fe_displacement);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) fe_displacement,
"displacement");CHKERRQ(ierr);
ierr = PetscFEGetQuadrature(fe_displacement,&q);CHKERRQ(ierr);
ierr = PetscQuadratureGetOrder(q,&order);CHKERRQ(ierr);
/* Creating the FE for the stress */
ierr =
PetscFECreateDefault(dm,dim,2*dim,simplex,"stress_",PETSC_DEFAULT,&fe_stress);CHKERRQ(ierr);
ierr = PetscFESetQuadrature(fe_stress,q);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) fe_stress,
"cauchy_stress");CHKERRQ(ierr);
/* Discretization and boundary conditons: */
ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe_displacement);
CHKERRQ(ierr);
ierr = PetscDSSetDiscretization(prob, 1, (PetscObject)
fe_stress);CHKERRQ(ierr);
ierr = DMSetDS(dm, prob); CHKERRQ(ierr);
/* Define the boundaries */
const PetscInt Ncomp = dim;
const PetscInt Nfid = 1;
PetscInt fid[Nfid]; /* fixed faces [numer of fixed faces] */
PetscInt restrictall[3] = {0, 1, 2}; /* restricting all movement */
fid[0] = 3; /* fixed face */
ierr = DMAddBoundary(dm, PETSC_TRUE, "fixed", "Face
Sets",0,Ncomp,restrictall,(void (*)()) zero_vector,
Nfid,fid,NULL);CHKERRQ(ierr);
ierr = PetscFEDestroy(&fe_displacement); CHKERRQ(ierr);
ierr = PetscFEDestroy(&fe_stress); CHKERRQ(ierr);
return(0);
}
int main(int argc, char *argv[])
{
DM dm, distributeddm; /* problem definition */
Vec u,expected_solution,projected_solution;
PetscViewer viewer;
int dim; /* dimension of the anlysis */
PetscErrorCode ierr;
PetscMPIInt rank, numProcs;
PetscPartitioner part;
ierr = PetscInitialize(&argc,&argv,NULL,NULL);
ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD,"testcube.msh",
PETSC_TRUE,&(dm));
ierr = DMGetDimension(dm,&(dim));
/* distribute the mesh */
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &numProcs);
DMPlexGetPartitioner(dm, &part);
PetscPartitionerSetType(part, PETSCPARTITIONERPARMETIS);
ierr = DMPlexDistribute(dm,0,NULL,&distributeddm);
if (distributeddm) {
ierr=DMDestroy(&(dm));
dm = distributeddm;
}
ierr = DMView(dm,PETSC_VIEWER_STDOUT_WORLD);
ierr = DMSetMatType(dm, MATAIJ);CHKERRQ(ierr);
ierr = SetupDiscretization(dm);CHKERRQ(ierr);
ierr = DMPlexCreateClosureIndex(dm,NULL);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(dm, &(u)); CHKERRQ(ierr);
ierr = VecDuplicate(u,&expected_solution); CHKERRQ(ierr);
ierr = VecDuplicate(u,&projected_solution); CHKERRQ(ierr);
/* intitialize the fields: */
PetscErrorCode (*initial[2])(PetscInt dim, PetscReal time, const PetscReal
x[], PetscInt Nf, PetscScalar u[], void* ctx) =
{initial_displacement_vector,zero_stress};
ierr = DMProjectFunction(dm,0.0,initial,NULL, INSERT_ALL_VALUES, u);
CHKERRQ(ierr);
PetscErrorCode (*expected_sol[2])(PetscInt dim, PetscReal time, const
PetscReal x[], PetscInt Nf, PetscScalar u[], void* ctx) =
{initial_displacement_vector,expected_stress_vector};
ierr = DMProjectFunction(dm,0.0,expected_sol,NULL, INSERT_ALL_VALUES,
expected_solution); CHKERRQ(ierr);
ierr =
PetscViewerVTKOpen(PETSC_COMM_WORLD,"expected_solution.vtu",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) expected_solution,
"expected_fields"); CHKERRQ(ierr);
ierr = VecView(expected_solution,viewer);CHKERRQ(ierr);
ierr= PetscViewerDestroy(&viewer); CHKERRQ(ierr);
/* project the fields: */
void (*projection[2])(PetscInt dim, PetscInt Nf, PetscInt NfAux, const
PetscInt uOff[],
const PetscInt uOff_x[], const PetscScalar u[], const
PetscScalar u_t[],
const PetscScalar u_x[], const PetscInt aOff[], const
PetscInt aOff_x[],
const PetscScalar a[], const PetscScalar a_t[], const
PetscScalar a_x[],
PetscReal t, const PetscScalar x[], PetscInt
numConstants,
const PetscScalar constants[], PetscScalar v[]) =
{projectdisplacement, projectstress};
ierr = DMProjectField(dm, 0.0, u,
projection,INSERT_ALL_VALUES,projected_solution); CHKERRQ(ierr);
ierr =
PetscViewerVTKOpen(PETSC_COMM_WORLD,"projected_solution.vtu",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) projected_solution,
"projected_fields"); CHKERRQ(ierr);
ierr = VecView(projected_solution,viewer);CHKERRQ(ierr);
ierr= PetscViewerDestroy(&viewer); CHKERRQ(ierr);
VecDestroy(&u);
VecDestroy(&expected_solution);
VecDestroy(&projected_solution);
DMDestroy(&dm);
PetscFinalize();
return 0;
}
testcube.msh:
$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
6
2 1 "back"
2 2 "front"
2 3 "bottom"
2 4 "right"
2 5 "top"
2 6 "left"
$EndPhysicalNames
$Nodes
10
1 0 -0.05 0
2 0.1 -0.05 0
3 0.1 0.05 0
4 0 0.05 0
5 0 -0.05 0.1
6 0.1 -0.05 0.1
7 0.1 0.05 0.1
8 0 0.05 0.1
9 0.05 0 0
10 0.05 0 0.1
$EndNodes
$Elements
28
1 2 2 1 6 1 2 9
2 2 2 1 6 1 9 4
3 2 2 1 6 2 3 9
4 2 2 1 6 3 4 9
5 2 2 3 15 5 1 6
6 2 2 3 15 1 2 6
7 2 2 4 19 6 2 7
8 2 2 4 19 2 3 7
9 2 2 5 23 7 3 8
10 2 2 5 23 3 4 8
11 2 2 6 27 8 4 1
12 2 2 6 27 8 1 5
13 2 2 2 28 5 6 10
14 2 2 2 28 5 10 8
15 2 2 2 28 6 7 10
16 2 2 2 28 7 8 10
17 4 2 29 1 1 2 9 6
18 4 2 29 1 6 5 10 9
19 4 2 29 1 1 5 6 9
20 4 2 29 1 1 9 4 8
21 4 2 29 1 10 5 8 9
22 4 2 29 1 9 5 8 1
23 4 2 29 1 2 3 9 7
24 4 2 29 1 7 6 10 9
25 4 2 29 1 2 6 7 9
26 4 2 29 1 3 4 9 8
27 4 2 29 1 8 7 10 9
28 4 2 29 1 3 7 8 9
$EndElements
>
> Thanks,
>
> Matt
>
> Thanks,
> Max
>>
>> Thanks,
>>
>> Matt
>>
>> I evaluate the current trial stress by adding a stress increment assuming
>> elastic behaviour. If the trial stress lies beyond the yield stress I
>> calculate the corrected stress to evaluate my residual for the
>> displacements. But now I somehow need to update my plastic strain and the
>> stress in the auxiliary fields. So in tspoststep I created another SNES to
>> now calculate the stress and plastic strain while the displacement is the
>> auxiliary field.
>>
>> I’m sure there’s an elegant solution on how to update internal variables but
>> I have not found it.
>>
>> Thanks,
>> Max
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>> Thanks,
>>> Max
>>>
>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>> http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/~mk51/>
>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their experiments
>> is infinitely more interesting than any results to which their experiments
>> lead.
>> -- Norbert Wiener
>>
>> http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/~mk51/>
>
>
>
> --
> What most experimenters take for granted before they begin their experiments
> is infinitely more interesting than any results to which their experiments
> lead.
> -- Norbert Wiener
>
> http://www.caam.rice.edu/~mk51/ <http://www.caam.rice.edu/~mk51/>
