On Thu, Sep 21, 2017 at 11:28 AM, Maximilian Hartig < [email protected]> wrote:
> On 20. Sep 2017, at 22:51, Matthew Knepley <[email protected]> wrote: > > On Wed, Sep 20, 2017 at 3:46 PM, Maximilian Hartig <imilian.hartig@gmail. > com> wrote: > >> >> On 20. Sep 2017, at 19:05, Matthew Knepley <[email protected]> wrote: >> >> On Wed, Sep 20, 2017 at 12:57 PM, Maximilian Hartig < >> [email protected]> wrote: >> >>> On 20. Sep 2017, at 18:17, Matthew Knepley <[email protected]> wrote: >>> >>> On Wed, Sep 20, 2017 at 11:46 AM, Maximilian Hartig < >>> [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. > Okay, I see. Lets leave that for the time being. > > >> 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. > Right, but if the stress is local, then you can just compute it in the kernel (this is what I do). Even if the relationship is implicit, you can compute that solve local to the (Gauss) point. > 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. > Okay, I got this. Can you tell me how to run it and why you think stuff is wrong? Thanks, Matt > 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/ >>> >>> >>> >> >> >> -- >> 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/ >> >> >> > > > -- > 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/ > > > -- 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/
